断面積の違う熱伝導方程式

長さL、断面積S(x)の棒での熱伝導方程式の解が知りたいです。
偏微分の形で構いません。
どの教科書を見ても面積が一定で分かりません。

お願いします！！

もし境界条件が必要なら、
T(0,x)=a
T(t,L)=a
T(t,0)=bt+a
でといていただけるとめちゃめちゃありがたいです。

No.5 回答者： lycaon 回答日時： 2009/03/19 18:22

【例題3】例題2の円錐台。

側面断熱、初期は全体が20℃。
上面20℃のままτ＝0 に下面加熱開始、一定昇温速度で下面を100℃まで上げ、以後は下面100℃、上面20℃に保つ。
加熱開始後の温度分布推移を求めよ。

(4)式でΘとζの項を消し a=λ/(Cvρ)=温度伝導率 と置き

∂T/∂τ＝a{∂^2T/∂r^2+(1/r)(∂T/∂r)+∂^2T/∂x^2 }・・・(7)

差分方程式は見つからず自分で導出。長時間経過後に例題2の分布と一致、まず正しいでしょう。

T<ｎ+1>(i,j)＝(1-4B)T(i,j)+B{T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1)}+B{T(i+1,j)-T(i-1,j)}/2i}・・・(8)

Δx=Δr、B=aΔτ/Δx^2。
iは半径、jは軸方向。
右辺は時刻τ=nΔτの各部温度、左辺はτ=(n+1)Δτの点(i,j)の温度。

(8)式はオイラー陽解法で時間発展的に解けます。
B≦1/2 が必須条件。

円錐台高0.1m、底面0.05mR、上面0.01mR
Δx=Δr=0.001ｍ。格子数x方向100、r方向50
a=9.2*10^-8 m2/s（アクリル）

最初全体20℃。時刻0から0.02秒/sで下面加熱、100℃到達後横這い。
中央温度推移（10分刻み）の太い青線まで左端下底の間隔均等。
太い青線は70分経過後下底が100℃に揃った瞬間。太い赤線は前回のもの。
青線と赤線の間は定常状態に移行中。

昇温中は左端上昇に引っ張られ下に凸の曲線。
昇温が止まると、フーリエの法則
　Q＝A・λ (T1-T2) /L
で、断熱ゆえ任意断面でQ一定。
Aが半分になれば(T1-T2)/L は2倍、上下底の半径比大ほど大きく上に凸へ変化。
(終)

No.4 回答者： lycaon 回答日時： 2009/03/19 18:14

【例題2】台形を、上底と下底の各中点を結ぶ中心軸回りに回転させて得る円錐台。

側面断熱、下面を100℃、上面を20℃に保持したときの中心軸を通る平面内の定常温度分布を求めよ。

円筒座標の熱伝導方程式

Cvρ∂T/∂τ＝λ{∂^2T/∂r^2 + (1/r)(∂T/∂r) + (1/r^2)(∂^2T/∂Θ^2 +∂^2T/∂x^2 }+ζ　・・・(4)

Cv、ρ、λ、T、τ、ζの意味は(1)式に同じ。
物性値が一様で系内に発熱がない、軸対称の円錐台の定常状態では、τ、Θ、ζの項が消え、次のラプラス方程式になります。

　　∂^2T/∂r^2 + (1/r)(∂T/∂r)+∂^2T/∂x^2 ＝0　・・・(5)

この差分形は
http://oacis.lib.kaiyodai.ac.jp/dspace/bitstream …
の巻末（62頁以降）で算出されています。

　　T(i,j)＝{T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1)+{T(i+1,j)-T(i-1,j)}/2i}/4　・・・(6)

「第5項のiは半径方向の格子番号であり，相対的な半径を意味する．このiが差分公式中に表れることが円筒座標を用いるときの特徴となる」との記載があります。

Excelのセルに式を埋め込む方法でも可能ですが、階段状の輪郭内に式をコピペするのは厄介ゆえ、今回はＶＢ6 を用いて解きました。
円錐台高さ100mm、底面 50mmR（100℃）、上面 10mmR（20℃）
CVSファイルに結果を落とし Excelでグラフ化しました。
プログラムはこのサイトの800字制限に全然納まらず、割愛します。

No.3 回答者： lycaon 回答日時： 2009/03/19 18:05

【例題1の解法手順】

(1)新しい Excelを開く。
(2)ツール－オプション－計算方法－手動－反復計算レ－最大反復回数12000－変化の最大幅0.001－OK
(3)A1～A51セルに100を記入、これらのセルを赤で塗る。

(4)B1に =(2*A1+2*B2)/4 を記入、緑に塗る。
(5)C2に =(B2+D2+2*C3)/4 を記入、緑に塗る。
(6)D2に =(2*C2+2*D3)/4　を記入、緑に塗る。
(7)C2とD2の2つのセルをマウスで選択し Ctrl+C でコピーし、E3,G4・・・と、2つづつ階段状に、AM20まで Ctrl+V でペーストする。

(8)B51に =(2*A51+2*B50)/4 を記入、緑に塗る。
(5)C50に =(B50+D50+2*C49)/4 を記入、緑に塗る。
(6)D50に =(2*C50+2*D49)/44　を記入、緑に塗る。
(7)C50とD50の2つのセルをマウスで選択して Ctrl+C でコピーし、E49,G48・・・と、2つづつ階段状に、AM32まで Ctrl+V でペーストする。

(8)AN20に数字20を記入、青く塗り直してから、AN21～AN32までコピーする。

(9)これで上底13cm（定常ではmでもmmでも同じ結果。非定常は単位で違う結果)、20℃、下底51cm、100℃、側面断熱条件の高さ40cmの台形の枠が完成。枠内の全セルに、「隣接する上下左右のセルの値の合計/4」の式を書き込む。

(10)ツール－オプション－計算方法－手動－再計算実行。
(11)結果の面内温度分布を等高線グラフ（5℃刻み）で、上底と下底の中点を結ぶ線上の温度分布を線グラフで作図。（前回示しました。）

No.2 回答者： lycaon 回答日時： 2009/03/19 17:56

【例題1】上底と下底以外は断熱で、下底を100℃、上底を20℃に保持した２次元台形板の板内定常温度分布を求めよ。

直交座標の熱伝導方程式
Cv・ρ・∂T/∂τ＝λ{∂^2T/∂x^2 + ∂^2T/∂y^2 + ∂^2T/∂z^2 }+ζ　・・・(1)

　Cv(J/kg・℃)：定積比熱（固体では≒Cp、定圧比熱）
　ρ(kg/m^3)：密度
　λ(J/m・s・℃)：熱伝導率
　Ｔ(℃)：温度
　τ(s)：時間
　ζ(J/m^3・s)：発熱量
物性値一様で面内発熱がない定常２次元では、τ、ｚ、ζの項が消え、物性値も消え去って次のラプラス方程式になります。

∂^2T/∂x^2 + ∂^2T/∂y^2 ＝0　・・・(2)

差分法では微小有限長さΔx、Δy を用い、∂T/∂x等を｛T(i+1,j)-T(i,j)}/Δx 等と近似します。
∂^2T/∂x^2 は、[{T(i+1,j)-T(i,j)}/Δx-{T(i,j)-T(i-1,j)}/Δx]/Δx＝{T(i+1,j)+T(i-1,j)-2T(i,j)}/Δx^2
同様に ∂^2T/∂y^2 ＝{T(i,j+1)+T(i,j-1)-2T(i,j)}/Δy^2
ΔxやΔyは計算者が決め、普通はΔx＝Δy にとります。(2)式に入れ整理し

T(i,j)＝{T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1)}/4　・・・(3)

Excelをイメージすれば i行j列セルの温度T(i,j)は、そのセルに隣接する上下左右の4つのセルの単純平均になります。
境界条件を与え(3)式を収束するまで計算。
断熱境界条件は、そのセルの外側に架空のセルを想定、境界の一つ内側のセルと等温にします。境界セルの内外が同温なら境界に直角な方向に熱は流れません。
結果を下に添付。計算法の説明は次回答に。

No.1 回答者： lycaon 回答日時： 2009/03/19 17:49

こんにちは。

lycaonです。
化学工学便覧(3版)の記述では、
「初期および境界条件が簡単な場合は、非定常熱伝導方程式を解くことができる。しかしその解は複雑で実用上不便であるので、これらを線図化したものを用いるのが便利」とし、矩形平面板・太さ一様の円柱・球のHottel線図等が掲載されています。

後出する等温線は、果たして定式化できるでしょうか？
定常２次元ですら、矩形でない板の解析解は（流体や電気も同様ですが）あったとしても相当複雑でしょう。

しかし太さが変わる棒の非定常温度分布を求めるのは、偏微分方程式を差分法で解くのに、まさに絶好の練習問題です。老人ボケ防止にもなりそうで、久しぶりに取り組んでみました。

折角ですんで、プログラムは組めないが EXCELは使える多くの方々のため、EXCELだけで解く解法も例題1で掲示します。

【例題1】上底と下底以外は断熱で、下底を100℃、上底を20℃に保持した２次元台形板の板内定常温度分布を求めよ。

【例題2】台形を、上底と下底の各中点を結ぶ中心軸回りに回転させて得る円錐台。
側面断熱、下面を100℃、上面を20℃に保持したときの中心軸を通る平面内の、定常温度分布を求めよ。

【例題3】例題2の円錐台。側面断熱、初期は全体が20℃。
上面は20℃に保ったまま時刻τ＝0 に下面の加熱を開始し、一定の昇温速度で下面を100℃まで上げる。それ以後は下面 100℃、上面 20℃に保つ。
加熱開始直後からの温度分布の時間推移を求めよ。

以下、順に取り上げます。
基礎事項を取り上げたい上に添付図が多いので、回答を分けます。悪しからず。

この回答へのお礼

ほんとにありがとうございます！！
めちゃめちゃうれしいです！
参考にさせてもらいます！！

お礼日時：2009/03/19 21:59