熱伝導方程式および初期条件について

質問

温度実験の物理モデルをたててシミュレーションを行おうと試みたのですが上手くいかず分からなかったので投稿させていただきました．少しでも何かお分かりになる方がいましたら，ぜひよろしくお願いします．

実験内容：
温度が室温に保たれたアルミ棒(10cm)の右端をt=0において-5℃に引き下げ，以後その温度で維持した場合の，左端の温度変化について熱伝導方程式で物理モデルを表したいと思います．なお，左端には外乱は一切加えないものとし，一次元で考えます．

時間をt(s)，温度をθ，温度棒の一次元座標をxとすると，

初期条件　θ(t，0)= 0　
境界条件θ=(0，x)=-5
と考えたのですが

自分の解法では
∂θ/∂t=a^2*(∂^2)θ/∂(x^2)について
θ(t，x)= f(x)g(t)とし
{1/(a^2)*g(t)}*{∂g(t)/∂t}=λ
{1/f(x)}*{(∂^2)f(x)/(∂x^2)}=λ
と変数分離を行ったのですが
λの正・零・負の判定を行おうとしたところ初期条件が足らず出来ませんでした．
この方法ではすでに間違っているのかもしれませんが，
どなたか，初期条件・境界条件の指摘やこの実験の物理モデルがたてられる方がいらっしゃいましたら，厚かましいようですがよろしくお願いします．

回答No.2

inara1 · ベストアンサー

棒の長さが与えられているので、有限長の棒の加熱（棒の他端は断熱）の問題だと思います。

元の微分方程式は
　　　∂θ/∂t = a*(∂^2)θ/∂(x^2) --- (1)
だと思います（a は熱伝達率で a > 0 ）。θ= f(x)*g(t) とおくと
　　∂θ/∂t = f(x)*dg(t)/dt
　　(∂^2)θ/∂(x^2) = g(t)*d^2f(x)/dx^2
なので、これらを式(1)に代入すると
　　　f(x)*dg(t)/dt = a*g(t)*d^2f(x)/dx^2
となります。この両辺を a*f(x)*g(t) で割ると
　　　[ 1/{ a*g(t) } ]*dg(t)/dt = { 1/f(x) }*d^2f(x)/dx^2
この左辺は t だけの関数で、右辺は x だけの関数なので、等式が成り立つには右辺も左辺も定数でなくてはなりません。その定数の値をλとすれば
　　　[ 1/{ a*g(t) } ]*dg(t)/dt = { 1/f(x) }*d^2f(x)/dx^2 = λ
　　　→ dg(t)/dt - a*λ*g(t) = 0 --- (2)
　　　　d^2f(x)/dx^2 - λ*f(x) = 0 --- (3)
という2つの微分方程式が得られます。この段階ではλの符号はまだ決まりません。

式(2)は簡単に解けて
　　　g(t) = C1*exp( a*λ*t ) --- (4)
となります（C1 は定数）。もしλ= 0 なら、g(t) = C1 となるので、θ= f(x)*g(t) = C1*f(x) となって棒の温度 は時間変化しないことになってしまいます。したがってλ= 0 は間違った仮定です。

λが正か負の判定は、式(4)で t →∞ としたときを考えます。λ> 0 の場合、g(t) は無限に大きくなってしまう（ a > 0 なので a*λ> 0 ）ので、これも間違った仮定です。したがって、λが 0 でも正でもなければλ< 0 でなければなりません。

λ< 0 のとき、式(4)で t →∞ としたとき g(t) → 0 になりますが、0 は絶対零度ではありません。θは規格化温度で、Ts を端部の温度、T0 を初期温度としたときの θ= ( T - Ts )/( T0 - Ts ) ですから、θ= 0 というのは T = Ts の意味です。「冷却している端部以外の部分で棒が断熱されていれば」、時間の経過とともに、棒全体の温度は端部の温度 Ts になっていくので辻褄があっています。

λが負なので、λ= -p^2 と置き換えれば式(4)は
　　　g(t) = C1*exp( - a*p^2*t ) --- (5)
となり、また式(3)は
　　　d^2f(x)/dx^2 + p^2*f(x) = 0 --- (6)
となります。p^2 > 0 ですから式(6)の解は
　　　f(x) = C2*sin( p*x) + C3*cos( p*x ) --- (7)
となります（C2 と C3 は定数）。実際
　　　df(x)/dt = p*C2*cos( p*x ) - p*C3*sin( p*x )
　　　d^2f(x)/dx^2 = -p^2*C2*sin( p*x ) - p^2*C3*cos( p*x ) = -p^2*{ C2*sin( p*x) + C3*cos( p*x ) } = -p^2*f(x)
　　　→ d^2f(x)/dx^2 + p^2*f(x) = 0
ですから、式(7) は式(6) の解になっています。

これで g(t) も f(x) も求められたので、結局
　　　θ= f(x)*g(t) = C1*exp( - a*p^2*t ) *{ C2*sin( p*x ) + C3*cos( p*x ) }
となります。C1*C2 = A1、C1*C3 = A2 と書き換えれば
　　　θ= f(x)*g(t) = exp( - a*p^2*t ) *{ A1*sin( p* x) + A2*cos( p*x ) } --- (8)
となります。A1 と A2 は初期条件と境界条件から求められますが、実は、p の値は1つでなく無限にあります。それぞれの p の値に対応して A1 と A2 の値が決まるので、A1 と A2 の組み合わせも無限にあります。そのようなときの解は
　　　θ =∑[ n = 1, 2, ・・, ∞ ] exp( - a*p[n]^2*t ) *{ A1[n]*sin( p[n]*x ) + A2[n]*cos( p[n]*x ) } --- (9)
となります。p[n]、A1[n]、A2[n] のいうのは n の関数という意味です。ある n に対する p[n]、A1[n]、A2[n] を式(8)に代入したものが解ならば、それらを足し合わせたものも解になるからです。

これ以降は、初期条件や境界条件を使って p[n]、A1[n]、A2[n] を求めることになりますが、このまま説明していくと問題を全部解いてしまうことになります。ご質問は λ の符号の決め方ですが、ここまでの説明で理解されたと思います。

お礼

回答していただき，ありがとうございました．
大変分かり易く理解できました

係数A1 A2 pについてなのですが
θ(x,t)= exp( - a*p^2*t ) *{ A1*sin( p* x) + A2*cos( p*x ) }

θ(0,t)=-5より
θ(0,t)=exp( - a*p^2*t ) *{ A1*sin( p* 0) + A2*cos( p*0 ) }
=exp( - a*p^2*t ) *{A2*1} =-5
A2=(-5)/exp( - a*p^2*t )
よってθ= exp( - a*p^2*t ) * A1*sin( p* x) + (-5)*cos( p*x )

θ(x,0)=0より
θ(x,0)=exp( - a*p^2*0 ) * A1*sin( p* x) + (-5)*cos( p*x )
=1*A1+(-5)*cos( p*x ) = 0
A1= 5*cos( p*x )/sin( p* x)
よって
θ= exp( - a*p^2*t ) * 5*cos( p*x ) + (-5)*cos( p*x )
　 =5*cos(p*x)*{exp( - a*p^2*t )-1}

左端までの距離をLとすると
dθ(L,t)/dx =0より
dθ(L,t)/dx =-5*p*sin(p*L)*{exp( - a*p^2*t )-1}
これを満たすのは
・-5*p=0 → p=0
・{exp( - a*p^2*t )-1}=0 → p=0
・sin(p*L)=0 → Lは定数より p =n*pi (n=1.2..3.4........)
よって p=n*pi (n=0.1.2.3.....)

これを代入して
θ =∑[ n =0, 1, 2, ・・, ∞ ] 5*cos(n*pi*x)*(exp( - a*(n*pi)^2*t ) -1)

と計算したのですが，実際にt.nに数値を入れてみても上手くいかず分かりませんでした．

この係数決定に間違いなどございましたら，ご教授お願いします．

回答No.1

eatern27

>初期条件　θ(t，0) = 0
>境界条件θ=(0，x) =-5
初期条件の方がθ(0, x)、境界条件の方がθ(t, 0)ですね。
左端の方に境界条件を課していませんが、そうすると無限に長い場合を考えている事になってしまいますので、左端の影響を考えるのであれば何かの境界条件は課す事は必要です。
多分、左端に熱流が流れない(θのx微分が左端で0)という境界条件で良いかと。

λの符号の判定って何をやっているのか分かりませんが、境界条件が増えれば判定できる話なのでしょうか。とりあえず考えてみて下さい。

※正確にはこのままだとθ(0,0)の値がいくらか分からないので問題と言えば問題ですが、１点だけなのであまり気にならないはず。

お礼

遅くなってすみませんでした．
回答のおかげで先に進むことが出来ました．
また，詰まってしまいましたが前進することが出来，感謝しております