極座標系の3次元熱伝導方程式の解法

質問

上記の式dT/dt=(a/r^2)(d/dr)(r^2(dT/dr))の解法をご存知の方がいらっしゃいましたら、教えてください。“d”はラウンド・ディーのつもりです。“a”は熱伝導率などを含んだ定数と扱ってよろしいです。

回答No.6

inara1 · ベストアンサー

温度 Ts の球の表面が、遠方温度 T∞ の流体で冷却または加熱されている場合、熱伝達係数を h (W/m^2/K）とすれば、球の表面からの伝熱量 q [W/m^2] は、ニュートンの冷却則により
　　　q = h*( Ts - T∞ )
ですが、これは球の表面での温度勾配による伝熱量に等しい（フーリエの法則）ので
　　　q = -k*∂Ts/∂r
が、全ての時間 t に対して、球の表面（ r = D/2 ）で成り立ちます。また、時間 t = 0 では、0≦r≦D/2 の範囲で
　　　T(r,0) = T0
という初期条件を仮定します（最初は球内部の温度は一定）。

この境界条件と初期条件から球内部の温度分布 T( r , t ) を求めることができますが、手元の書籍 [3] にその解が出ているで、それを引用します。
　　　T(r,t) = 2*T0*∑(n=1～∞) [ { sin (νn) - νn*cos (νn) }/{ νn - sin (νn)*cos (νn) } ]*exp( -4*νn^2*a*t/D^2 )*sin( 2*νn*r/D )/( 2*νn*r/D )
a は球の温度伝達率 = 熱伝導率λ/密度ρ/比熱c です。またνn は以下の式の解です（0 ≦ ν1 < ν2 <ν3 <・・・）。
　　　ν*cos(ν) = ( 1 - h*D/2 )*sin(ν)
h*D/2 がどのような値であっても（0～無限大）、0≦ν1≦π、ν1≦ν2≦2*π、ν2≦ν3≦3*π の範囲にあります。「実用上はν1～ν4までで間に合う」とこの書籍には書かれています（無限級数計算しなくてもn = 1～4 でいいということ）。

[3] 谷下市松「伝熱工学」p.65、裳華房（1986）

お礼

詳しいご回答有難うございました。先日、問題が解決しました(答えが出ました)と書きましたが、私も谷下市松著「伝熱工学」をインターネットで見つけ、半径Rの3次元非定常伝熱計算が扱われていることが載っていたので購入しました。そしてP65ページのTを求める無限級数の式と、超越関数の方程式の解νを用いて計算したために解決できました。一点だけ気になる仮定を自分成りに立てたのですが、それは熱伝達率が外に熱を逃がす場合と、外から熱を貰う場合とで同じになるという仮定です。熱伝達率(係数)hは境界条件を代表するものだと思いますので、熱抵抗という考えからは、どちらの方向でも同じになりそうで、講演会の質問の際に、むかし東北大学の先生に聞いたときは、同じだと答えてくれましたが、輻射熱を考えると光の進行方向が逆になるような気がして、エネルギーが移動する方向側にある物体の吸収率とか透過率が向きによって違うような気がします。そこだけが気になっています。流体力学については詳しくないですので、ご回答の内容を勉強してみようと思います。上のご回答内容についてはほとんど理解できます。有難うございます。

回答No.7

inara1

加熱する場合と冷却する場合で熱伝達係数が違うのではという点ですが、加熱する場合と冷却する場合で流体の温度が違いますが、それらの温度での物性値の違いが無視できるのなら同じはずです。固体や流体を介した伝熱では、接触している部分だけでの熱のやりとりなので、物性値が同じなら方向性はないはずです。

しかし輻射の場合はちょっと複雑です。輻射には放射率という量がありますが、放射率は、輻射エネルギーの透過率のことなので、放射率が1より小さければ、吸収されずに表面で反射される熱があるということです。その熱が最終的にどこに行くのかは物体の形状によるので、輻射熱を吸収する場合と放出する場合では対象性が失われると思います。輻射に関しては勉強不足で充分にお答えできません。

補足ですが、輻射によって伝わる熱エネルギーは2つの物体の絶対温度の4乗の差に比例しますが、温度差が小さければ、温度差に比例するとみなせるので、熱伝達（温度差に比例）と同じ扱いになります。つまり、輻射による熱移動を、ある熱伝達係数を持った伝熱として表わすことができます。

お礼

ご回答有難うございました。輻射の問題は量子力学的な分野になると思いますので、理論が極めて複雑になると思います。2物体の絶対温度の4乗の差になるのはステファン＝ボルツマンの法則だったと思います。伝熱学の本では、この輻射熱による移動も熱伝達係数として一緒くたにして扱うので、方向性があるというのが正しいような気はします。

回答No.5

inara1

手元にある伝熱の書籍 [1] を見てみたところ、自然対流による球の冷却では、熱伝達率は一定ではないようです。

自然対流による球のNu数は次式で表わされます[2]。
　　Nu = 2 + 0.589*Ra^(1/4)/[ 1 + ( 0.469/Pr )^(9/16) ]^(4/9)
これは Pr>0.7、Ra < 10^11 で有効な経験式です。Ra はレイリー数、Pr はプラントル数です。Ra は
　　Ra = g*β*( Ts - T∞ )*D^3/( ν*α)
g：重力加速度(m/s^2)、β：流体の体積膨張率、Ts：球の表面温度(℃)、T∞：遠方流体温度(℃)、D：球の直径(m)、ν：流体の動粘性係数(m^2/s)、α：流体の温度伝達率(m^2/s)
で表わされます。冷却媒体が空気なら、βは 1/Tf [1/K] になります（Tf = ( Ts+T∞)/2 )。冷却媒体が温度300Kでの空気なら ν=1.568×10^(-5) m^2/s、α=2.213×10^(-5) m^2/s、Pr = 0.708 です。

自然対流による球の熱伝達係数 h (W/m^2/K) は
　　　h = k*Nu/D
となりますが、Nu が球の表面温度 Ts の関数なので h は一定ではありません。ただし、Ts がある温度（元の球の温度）から室温まで変化してもNuがあまり変わらないのであれば h を定数とみなすことができます（実際の球の大きさや温度で計算してみてください）。h を定数とみなせる場合の計算法は後で紹介します。

[1] F.P.Incropera and D.P.DeWitt, Fundamentals of Heat and Mass Transfer ( 5th Ed.), John Wiley & Sons, Inc.(2001)
[2] S.W.Churchill, "Free Convection Around Immersed Bodies", Heat Exchange Design Handbook, Section 2.5.7, Hemisphere Publishing, New York, 1983.

回答No.4

inara1

熱伝達率が一定の場合の解を考えてみます（ちょっとお時間を下さい）。
輻射による熱伝達を考慮すると、微分方程式を解析的に解くことができなくなるので、とりあえず輻射を無視して計算してみます。

お礼

よろしくお願いします。出来れば私も後追いで確認計算してみたいです。

回答No.3

inara1

＞球形の高温の液滴が室温に持ってこられたときの冷却問題です。初期条件は液滴の初期温度と室温、境界条件は熱伝達率です。

状況は分かりました。液滴の初期温度が非常に高温だと、液滴周囲の対流が層流でなく乱流となって、液滴と周囲（室温）との間の熱伝達率は一定とみなせないかもしれません。熱伝達率が一定とみなしていいのなら解を求めることはできます。平日は仕事なので少し時間を下さい。

お礼

計算時の熱伝達率は文献値を用いました。鉄板を炉内に入れるというもので、840Wm^2Kでした。乱流とか層流の影響は分かりませんが、実験値だろうと思います。値が比較的高いので、放射伝熱(輻射熱)も考慮されているのではないかとは感じました(定量的な判断はできませんが)。お時間を取ってもらってよろしいのでしょうか。恐縮しております。有難うございます。

回答No.2

inara1

ANo.1です。
変数分離法は熱伝導問題でよく使われます。C1 > 0 だと時間の経過とともに温度が発散してしまうため、C1 < 0 という条件をつけています。

考えている領域が有限の場合、その境界で温度一定とか、熱流束一定（断熱も含む）とか、熱伝達率一定などの条件がつきますが、その境界条件を満たす C1 の値は１つでなく複数になります。その場合、解は
　　　T(r,t) = ∑( n = 0～∞ ) ( An/r )*exp( a*Bn*t )*[ sin { r*√(-Bn) } + cos { r*√(-Bn) } ]
というように、無限級数であらわされます。An と Bn は n の関数ですが、その形は境界条件から決定できます。

お礼

再度ご回答頂き有難うございました。文献でも無限級数の式を見つけました。球形の高温の液滴が室温に持ってこられたときの冷却問題です。初期条件は液滴の初期温度と室温、境界条件は熱伝達率です。この条件で解かれた文献の式は非常に複雑で、自力では導き出せそうにありません。しかし、問題は一応解決しました(答えが出ました)。inara1さんに教えていただいた式にある程度似ています。もし解法を詳しく説明されている本が有りましたら、教えていただけると助かります。ラプラシアンの角度成分を扱った球面調和関数などは解いてみたことがあるのですが、今回の式はほとんど扱っている文献が見当たりませんでした。

回答No.1

inara1

T(r,t) = R(r)*T(t) と置いて、元の微分方程式に代入すれば
　　　dT/dt = a*C1*T(t)
　　　d^2R/dr^2 = C1*R - (2/r)*dR/dr
なので
　　　T(t) = C2*exp( a*C1*t )
　　　R(r) = (C3/r)*[ sin { r*√(-C1) } + cos { r*√(-C1) } ]
C1 (< 0 )、C2、C3 は境界条件から求められるはずです。

お礼

ご回答有難うございました。元々熱伝導の問題だったので、昨日までは伝熱学と微分方程式の本ばかり探して見つからなかったのですが、今朝は量子力学の本を見ていたら水素原子の電子軌道に関して結構見つかりました。会社の就業時間になりましたので、仕事が終わってからじっくり考えて見たいと思います。とり急ぎお礼申し上げます。