浸出水の予測シミュレーション

一日の侵入水量I(m3/日　面積×日降雨量)から、地表に流出する地下水量Q(m3/日)を予測しています。地下にタンクがあると仮定して、Qはタンクの水位Hに比例して、抵抗Rに反比例すると仮定するとQ=H/R 水位Hの変化はdH/dt=I-Qとして、t=0の時Q=Q0の条件で、Q=I+(Q0-I0)exp(-t/R)・・・式1となります。
式1を使って、当日降雨量Ijと前日降雨量I'j-1'、前日排水量Q'j-1'によりj日の排出量Qjを求めたいのですが、教えていただけませんか?
参考書では、Qj=Q'j-1'exp(-1/R)+{Ij-R[Ij-I'j-1'][1-exp(-1/R)]-I'j-1'exp(-1/R)}となっているのですが、なぜ途中で積分を使っているのかもわかりません。どうかよろしくお願いします。

No.1 回答者： inara 回答日時： 2007/03/10 20:48

式(1)がそもそもおかしいと思われます。

・流路抵抗（Ｒ）の単位が合っていません
・I は時間変化するはずなので、その具体的な形とか I の微分方程式がないと計算（積分）ができません。

【単位をあわせて途中まで考えた経過】
タンクから排出される水の流量を Q [m^3/day]、タンクに入る水の流量（進入水量）を I [m^3/day] 、タンクの水位（排水口からの）を H [m]、タンク内の水面の断面積を A [m^2] とします。このとき、タンクに流入している水量（流量）は I - Q です。これは水位の変化率 dH/dt に水面の断面積 A をかけたものに等しいので、

A*dH/dt = I - Q --- [1]

が成り立ちます（ここですでに sargassosea さんの質問中の「dH/dt=I-Qとして」と違っています）。

一方、Q は水位 h に比例しまので、Q = H/R --- [2]

(この定義だと R の単位は動粘性係数の逆数 [h/m^2] になります）。式[2]を式[1]に代入すると、

A*dH/dt = I - H/R --- [3]

となります。進入水量 I が時間 t の関数ならばこの解は

H(t) = e^{ -t/(R*A) }*[ C -1/A*∫{I(t)*e^{ -t/(R*A) } dt --- [4]

となります（I(t)の形が分からないので積分のまま）。もし仮に、I(t) = I0 なら、H(0) = R*Q0 と式[2]を使えば、式 [3] の解は

Q = H(t)/R = -I0 + ( Q0 + I0 )*e^{ -t/( R*A ) } --- [5]

となりますが、この式は sargassoseaさんの式1とも違っています。

物理量の式を立てるときは、＝ の左右の単位（次元）が合っていないと式として成り立ちません（質量＝時間は数字が合っていても無意味）。 sargassoseaさんの式1の左辺は m^3/day、右辺の I+(Q0-I0)もm^3/day ですが、expの中は本来無次元であるべきなので、R の単位は時間 [day] ということになります。しかし、R の定義は Q = H/R ですので、これだと R の単位は m/(m^3/day) = day/m^2 となって合いません。また式[1]が sargassoseaさんの「dH/dt=I-Q」と違いますが、単位をあわせるなら、式[1]が正しいはずです（だから expの中も違ってくる）。

この質問は問題文を省略して書かれていると思われますが、
浸入水量 I について、問題に何か書いてないですか？　
参考書に書かれているという「Qj=Q'j-1'exp(-1/R)+{Ij-R[Ij-I'j-1'][1-exp(-1/R)]-I'j-1'exp(-1/R)}」では、Ｒが無次元ですが、質問文にある Q = H/R というのは、本当は Q = B*H/R （ B は排水口の断面積[m^2]）になっていませんか？だとしたら Ｒ は無次元でいいです。

この回答への補足

すみません。おっしゃるとおり単位が間違っていました。
Q=H/R 　　　　　　A'dH/dt=AI-AQ
Q : mm/d 単位面積あたりの浸出水量
R : タンク時定数　d
A : 集水面積　

浸入水量Iは、日降雨量(mm/d)により変化します。日降雨量は気象台のデータです。回答中の4式中の
∫{I(t)*e^{ -t/(R*A) } dt について解説してくれませんか？
データが日降雨量であるから、積分するということでしょうか？
なにかヒントでもお願いします。よろしくお願いします。

補足日時：2007/03/10 21:43

この回答へのお礼

お礼を回答への補足に入れてしまいました。
inaraさますみませんでした。そして回答ありがとうございます。
内容についてかなり早とちりしたようです。積分の意味がわかりました。どうもありがとうごさいました。sargassosea

お礼日時：2007/03/10 22:18

No.2 回答者： noname#26313 回答日時： 2007/03/10 22:49

Q=H/Rを用い、dH/dt=I-Qから、R・(dQ/dt)=I-Q、

そして(dQ/dt)+(Q/R)=I/Rが導けます。

積分因子、exp(t/R)を乗じると、
(d/dt){Q・exp(t/R)}=(I/R)・exp(t/R) ・・・・・（2式）
となります。
Q・exp(t/R)-Q_0=∫{(I/R)・exp(t/R)}dtですが、Iがtの函数であるようなので、
（2式）をΔ{Q・exp(t/R)}=(I/R)・exp(t/R)(Δt)と変形して、漸化式にします。

(Δt)を一日、Δ{ }を日毎の差とし、右辺は、t_'j-1'～t_'j'の間、I_'j-1'≒I_'j'で一定であるとして積分すると、
Q_'j'・exp(t_'j'/R)-Q_'j-1'・exp(t_'j-1'/R)=(I_'j'/R)・R・{(exp(t_'j'/R)-exp(t_'j-1'/R)}
となります。両辺をexp(t_'j'/R)で割ると、
Q_'j'-Q_'j-1'・exp(-1/R)=(I_'j'/R)・R・{1-exp(-1/R)}
∴　Q_'j'=Q_'j-1'・exp(-1/R)+(I_'j'/R)・R・{1-exp(-1/R)}
=Q_'j-1'・exp(-1/R)+I_'j'・[1-exp(-1/R)]

ここまで進めましたが、これから先は、I_'j'・[1-exp(-1/R)]の部分を丁寧な近似法で進め、
Q_'j'=Q_'j-1'・exp(-1/R)+{I_'j'-R[I_'j'-I_'j-1'][1-exp(-1/R)]-I_'j'-1・exp(-1/R)}
に持っていくのでしょう。

No.3 ベストアンサー 回答者： noname#26313 回答日時： 2007/03/11 10:22

#2の回答者です。

(I/R)・exp(t/R)(Δt)の部分の積分を、部分積分でやってみました。

'_j'などの記号の「'」は微分と紛らわしいので、省略します。
t_j-1～t_jの間で部分積分すると
∫(I/R)・exp(t/R)dt=(I/R)・[I・R・exp(t/R)-∫I'・R・exp(t/R)dt]
ですので
∫[t_j-1→t_j](I/R)・exp(t/R)dtは
(I/R)・[{I_j・R・exp(t_j/R)-I_j-1・R・exp(t_j-1/R)}
-(I_j-I_j-1)・R・R・{exp(t_j/R)-exp(t_j-1/R)}]
（なお、I'に相当する所を、I_j-I_j-1とおき、これを定数と見て
∫R・exp(t/R)dtの積分を行っています）
exp(t_j/R)で割ったものは、
(I/R)・[{I_j・R-I_j-1・R・exp(-1/R)}-(I_j-I_j-1)・R・R・{1-exp(-1/R)}]
={I_j-I_j-1・exp(-1/R)}-(I_j-I_j-1)・R・{1-exp(-1/R)}
=I_j-R・(I_j-I_j-1)・{1-exp(-1/R)}-{I_j-1・exp(-1/R)}
となり、全体は、
Q_j=Q_j-1・exp(-1/R)+I_j-R[I_j-I_j-1][1-exp(-1/R)]-I_j-1・exp(-1/R)
となり、参考書に書かれているとされているものと一致します。

この回答へのお礼

ていねいな回答ありがとうございます。
降雨量Iがtの関数である意味をもう一度全体を検討します。
どうもありがとうございました。

お礼日時：2007/03/12 12:16