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

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

お願いします！！

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

No.15 ベストアンサー 回答者： lycaon 回答日時： 2009/04/05 14:42

【正誤表】

回答No3.
　(誤)(6)D50に =(2*C50+2*D49)/44　を記入、緑に塗る。
　(正)(6)D50に =(2*C50+2*D49)/4　を記入、緑に塗る。

回答No5.
　(誤)B≦1/2 が必須条件。
　(正)B≦1/4 のときのみ意味ある近似解を得る。

回答No9.
　(誤)1回計算前後の max(i,j){T<τ+1>－T<τ-Δτ>}<0.005℃で打切り
　(正)...　T<τ>－T<τ-Δτ>　...

回答No11.
　(誤)現在圧倒的に不足
　(正)現在絶対的に不足

回答No13　添付図.
　(誤)格子の右端が119,120　
　(正)格子の右端は109,110

回答No14　添付図.
　(誤)右下回帰式　T-20＝(140.05-20)*...　
　(正)T-20＝(160.05-20)*...

No.14 回答者： lycaon 回答日時： 2009/04/05 01:27

Dim T(122),TT(122),境界(4),p As Long

'--
Private Sub Command1_Click()
Fil="c:\..(フルパス)..\デスクトップ\1.csv"
Open Fil For Output As #1
ρh=7900:λh=15:Ch=500:αh=λh/ρh/Ch
ρc=2100:λc=1:Cc=900:αc=λc/ρc/Cc
ρa=1.293:λa=0.024:μa=0.0000171
D=0.3 '直径(m)
U=15 '(m/s)
Re=ρa*U*D/μa:Pr=0.72
h強=(λa/D)*0.036*Re^0.8*Pr^0.333
Δτ=0.02 '(s)
Δx=0.0005 '(m)0.5mm
q=3000/(3.1415/4*D^2*0.005)
Bh=αh*Δτ/Δx^2 '0.30
Bc=αc*Δτ/Δx^2 '0.04
ω=q*Δτ/ρh/Ch
Ts=20 '遠方
Flag=1
境界(1)="断熱":境界(2)="自然対流":境界(3)="強制対流":境界(4)=Ts & "℃"
10 Print #1, 境界(Flag) & Time '☆
For i=0 To 110:T(i)=Ts:TT(i)=Ts:Next
For p=0 To 270000 '★90分
TT(0)=2*Bh*T(1)+(1-2*Bh)*T(0)+ω
For i=1 To 10
TT(i)=Bh*(T(i-1)+T(i+1))+(1-2*Bh)*T(i)+ω
Next i
For i=11 To 109
TT(i)=Bc*(T(i-1)+T(i+1))+(1-2*Bc)*T(i)
Next i
Select Case Flag
Case 1
TT(110)=2*Bc*T(109)+(1-2*Bc)*T(110)
Case 2
h自=1.27*((T(110)-Ts)/0.3)^0.25
B1=Δx*h自/λc
TT(110)=2*Bc*T(109)+(1-2*Bc)*T(110)+2*B1*Bc*(Ts-T(110))
Case 3
B1=Δx*h強/λc
TT(110)=2*Bc*T(109)+(1-2*Bc)*T(110)+2*B1*Bc*(Ts-T(110))
Case 4
TT(110)=Ts
End Select
If p Mod 15000=0 Then '5分
Print #1, p & ",,";
For i=0 To 110
Print #1, TT(i) & ",";
Next
Print #1, ""
End If
For i=0 To 110:T(i)=TT(i):Next
DoEvents
Next p '★
20 Print #1, ""
If Flag<4 Then Flag=Flag+1:GoTo 10 '☆
Print #1, "終了" & Time
MsgBox ("終了")
End
End Sub
'--
'D/Tに1.csv作成。編集後1.xls名で保管。
'CPU333MHzで24分。
'コード短縮のためSelect Case判定27万×4回。冗長でも10～20行を4回書き並べる方が速い。

参考URL：http://okwave.jp/qa4821596.html

No.13 回答者： lycaon 回答日時： 2009/04/05 01:18

系内発熱＆対流熱伝達

【例題7】30cmΦ*5mmのニクロムNI両面を銀板AGで挟み3kW発熱体HEとし底面側面断熱。通電直後出力安定.
[1]HE上面も断熱したら通電10分後のHEは何℃か.
[2]HE上面に30cmΦ*50mmのコンクリート円柱COを載せ密着し側面断熱。初期温度は全て20℃。上面が
　(a)断熱 (b)自然対流 (c)強制対流(15m/s) (d)定温
のとき軸上温度推移を調べよ

(NI)ρ:7900kg/m^3,λ:15W/m℃, Cp:500J/kg℃
(CO)ρ:2100,λ:1, Cp:900
(Air)ρ:1.293,λ:0.024,μ:1.71*10^-5kg/ms,Pr数=0.72
全て一定

解
[1]HE完全断熱→T(℃)-20 ∝τ(s)。NI面積0.0707m2、体積353cc、重量2.79kg。10分後3000*10*60=500*2.79*(T-20)→T=1310℃([2]ではCOに奪熱され(a)断熱でも90分後250℃以下)

[2]熱流//軸,1次元.電極銀でNI均一発熱.AG温度勾配≒0→差分格子でAG省略.

1次元非定常
　∂T/∂τ＝a*∂^2T/∂x^2 + q/ρ/C・・・(12)
差分式
　T<n+1>(i)=B*{T(i-1)+T(i+1)}+(1-2*B)*T(i)+ω・・・(13)
　B= a*Δτ/Δx^2 （≦1/2 )、ω=q*Δτ/ρ/C
q[J/sm^3]:発熱量,ρ:密度,λ:熱伝導率,C:比熱,a[m^2/s]=λ/ρ/C:熱拡散率

CO,AG区間はq=0.HE区間:q=3000/353E-6＝8.5*10^6(J/sm^3)
Δx=0.0005(m)としHEを10分割,COを100分割。B≦1/2よりΔτ=0.02(s)と設定.
対流熱伝達の式はURLによる.(下図とコード中に記載)

次回答図(a)→(d)順に上面放熱増
(b)90分後上面138℃の対流損熱 hAΔT=[1.27*(118/0.3)^0.25]*0.0707*118=47W
CO輻射損熱 0.0707*4.88*0.87*(4.11^4-2.93^4}=63kcal/h=73W >対流だがここでは輻射無視

(d)定常到達後HE表面 3000=1*0.0707*(T-20)/0.05→2142℃
HE表面T予測式:T=20+(T∞-20)*(1-exp(-ｋτ))、T∞=2142℃最適ｋなし
T∞=160.05(物理意味不明),ｋ=0.0008 が適合,図示

VB6だがVBAで動く
1次元は解析解あり.定常HE温度:放物線、非定常CO温度:xの正弦*τの指数関数の無限級数

参考URL：http://chemeng.on.coocan.jp/ce/ce4.html

No.12 回答者： lycaon 回答日時： 2009/03/28 03:44

Dim T(50,50) As Single,Ta(50,50) As Single,TT(50,50) As Single

'----
Sub test4()
'Application.ScreenUpdating=False
t1=Time
'**初期
With Sheets("Sheet1")
For j=1 To 50: For i=1 To 50
T(i,j)=.Cells(i,j).Value
Next: Next
Tc=-1
T(25,25)=Tc 'τ=0で変更
End With
Δx=0.01
Δτ=60
a=0.0000001
B=a*Δτ/Δx^2
C=1-4*B
Sheets("Sheet4").Activate
With ActiveSheet
.Cells.Clear
.Cells(3,52).Value="●25行": .Cells(4,52).Value="分"
.Cells(3,104).Value="●25列": .Cells(4,104).Value="分"
.Cells(28,52).Value="●24行"
.Cells(28,104).Value="●24列"
回=0
For p=1 To 6000 '★
For j=2 To 49: For i=2 To 49
Ta(i,j)=(1-4*B)*T(i,j)+B*(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1))
Next: Next j
Ta(25,25)=Tc
For j=2 To 49: For i=2 To 49: T(i,j)=Ta(i,j): Next: Next

'** 出力
v=20
If p=1 Or p Mod v=0 Then '☆
If p=1 Then k=1 Else k=53*(p/v)
.Cells(k,1).Value="■" & p*Δτ & "秒後"
For j=1 To 50: For i=1 To 50: .Cells(k+i,j).Value=T(i,j): Next: Next
.Cells(回+5,52)=p*Δτ/60
For j=1 To 50: .Cells(回+5,52+j)=T(25,j): .Cells(回+30,52+j)=T(24,j): Next
For i=1 To 50: .Cells(回+5,104+i)=T(i,25): .Cells(回+30,104+i)=T(i,24): Next

'** 定常確認
For j=2 To 49: For i=2 To 49
TT(i,j)=(1-4*B)*T(i,j)+B*(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1))
Next: Next
TT(25,25)=Tc: 差=0
For i=2 To 49: For j=2 To 49
Buf=Abs(TT(i,j)-T(i,j))
If Buf > 差 Then 差=Buf
Next: Next
回=回+1
If 差 < 0.1 Then Exit For
End If '☆
Next p '★
.Cells(1,52).Value="計算=" & p-1 & "回 " & t1 & "～" & Time
End With
MsgBox ("終了")
End Sub

No.11 回答者： lycaon 回答日時： 2009/03/28 03:38

【例題6】例題4の温度分布を初期条件とし、時刻τ＝0 で格子点(25,25)の温度を-1℃に変え、そのまま保つとき、面内温度分布の時間推移を調べよ。

Sheet1 に例題4（冷点なし）の Sub test2 の結果のT(i,j)を置き、VB EditorのSheet4に次回答のコードを転写、実行して下さい。

コード中で、時刻τ=0に(25,25)格子点の値を-1(℃)にして、Δτ=60秒で(9)式を計算し、20分経過毎に面内温度分布をSheet4 のA～AX列に書き出します。
同時にSheet1の25・24行と25・24列の線状温度推移をSheet4のAZ～EX列に書き出します。
最長τ=6000分を設定。1ループ前後のmax｜T<n+1>-T<n>差｜<0.1℃で終了させます。
実際にはτ=120分経過後に条件を満たし計算終了。結果は下図。

等温等高線図は、20分より120分の方が陥没の太さが拡大しており、120分で例題5（定常）の結果とほぼ同じ。
正方形板中央十字線上の縦横温度分布は、開始後 20分まで一気に進行。60分以降はゆっくり。

今回は前問の関数は呼び出さず。冷点でも構わず(9)式を計算し、一周計算毎に-1℃に戻しました。（途中を問わない収束計算向き。関数呼出は途中を問う非定常向きで逆になりましたが、最初の回答から何度もデバッグし、今はやり直す気も出ず。）

以上の知識があれば、端面が一定速度で昇温する、太さが変わる棒の非定常熱伝導問題は解けると思います。
----

非定常熱伝導を解析的および差分近似で解く文献(pdf)がネットで見つかりました。
解析解法は簡単な形状でも相当な数学素養を要することを改めて確認しました。
数値解法（差分法など）を多くの日本人が身につけられ、現在圧倒的に不足すると言われるプログラマや、インド人に負けない数学・物理の達人たちが日本に輩出することを願って、回答を終わります。

参考URL：http://www.kz.tsukuba.ac.jp/~abe/ohp-heat/kougi0 …

No.10 回答者： lycaon 回答日時： 2009/03/28 03:31

【ワークシート(W/S)とVBA】

Excel W/S は2次元。
1次元定常と非定常、2次元定常はセルに式を埋め込み可能。
3次元(xyz)定常はW/Sに垂直方向のセル必要。z方向20分割ならW/S20枚。
第nシートの点(i,j)の4隣と、n-1、n+1シートの点(i,j)の計6セルの平均をnシートの点(i,j)へ。作成煩雑。
2次元非定常も3軸(xyτ)でW/Sでは困難。
3次元非定常は4軸で作成不可能。

MS Office装備のVBAを使えばこれらが可能に。
例題4・5で動くExcel VBAを例示。

【非定常】
熱伝導方程式(1)の2次元差分式は定常では(3)。非定常は

T<n+1>(i,j)＝(1-4*B)*T(i,j) + B*{T(i+1,j)+T(i-1,j)+(T(i,j+1)+T(i,j-1)}・・・(9)

Δx=Δy、B=aΔτ/(Δx)^2、a：熱拡散率（温度伝導率）
右辺は時刻τ=nΔτの各部温度、左辺はτ=(n+1)Δτの点(i,j)の温度。

もしB=1/4 に決めれば(9)式は

T<n+1>(i,j)＝{T(i+1,j)+T(i-1,j)+(T(i,j+1)+T(i,j-1)}/4・・・(10)

定常(3)式と同形だが左辺は現在、右辺は1ステップ過去の温度。

積み重ねたW/Sのイメージが判り易い。
nシートの或る温度T(i,j)は、n-1シートの点(i,j)回りの4～5点の値から(9)または(10)式で算出。3次元定常と違い nとn+1シートは使わず。

非定常熱伝導は初期条件と境界条件が必要。
1枚目W/Sに時刻0の全格子点の値を与えれば、Δτ時間後の全部のT(i,j)は(9)式で算出でき、収束計算なしで2枚目のW/Sが完成。
同様に3枚目、4枚目、と時間発展的に作成。

Δx≠Δy なら(9)式は

T<n+1>(i,j)＝(1-2*r-2*s)*T(i,j) + r*{T(i+1,j)+T(i-1,j)} + s*{(T(i,j+1)+T(i,j-1)}・・・(11)

ただし r=aΔτ/(Δx^2)、s=aΔτ/(Δy^2)
(8)(9)式で B>1/4、(11)で r+s>1/2 なら右辺T(i,j)の係数が負。
----
回答No5「B≦1/2 が必須条件」は1次元なら正だが 2次元なので誤り。
【訂正】
B≦1/4 のときのみ意味ある近似解を得る。

No.9 回答者： lycaon 回答日時： 2009/03/28 03:23

【例題5】例題4板のほぼ中央の格子点（25,25）が－１℃に保たれるときの定常温度分布を求めよ。

他の条件は例題4に同じとする。

Excel VBAの例続き。

結果下図。四周＆冷点で系外と熱の出入あり、特に冷点で出熱大。
λを与え冷点を閉曲線で囲めば線上温度勾配からQ(J/m2・s)も計算可能の筈。

配列宣言は各Subの外。
収束計算ではi,jを各々1から格子数まで増す。(3)式右辺の T(i+1,j)とT(i,j+1)は1ループ前の値、T(i-1,j)とT(i,j-1)は直前に更新した値。(ガウス・ザイデル法)
冷点（=境界値）も計算すれば、冷点通過以前のTは正しく通過後は不正。
2つの関数LL1、LL2 で点(25,25)を避け計算。
全セルで if i<>25 AND j<>25 ～ と判定するより速い。

1回計算前後の max(i,j){T<τ+1>－T<τ-Δτ>}<0.005℃で打切り→600回で終了。(下のコード）
コードを少し変え単純に10万回計算したら、中央陥没やや拡がる。添付図は後者。

Dim T(50,50),TT(50,50)
'---
Sub test3()
t1=Time
'**初期
For i=1 To 50:T(i,1)=200:T(i,50)=30:Next
For j=2 To 49:T(1,j)=120:T(50,j)=80:Next
For j=2 To 49:For i=2 To 49: T(i,j)=120:Next:Next
T(25,25)=-1
123 '**計算(20回毎)
For k=1 To 20
Call LL1(2,2,24,49)
Call LL1(25,2,25,24)
Call LL1(25,26,25,49)
Call LL1(26,2,49,49)
Next k
'**収束
For j=2 To 49:For i=2 To 49:TT(i,j)=T(i,j):Next:Next
TT(25,25)=-1
Call LL2(2,2,24,49)
Call LL2(25,2,25,24)
Call LL2(25,26,25,49)
Call LL2(26,2,49,49)
s=0
　(この間例題4と同じ)
End Sub

'---
Sub LL1(i1,j1,i2,j2)
For j=j1 To j2:For i=i1 To i2
T(i,j)=(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1))/4
Next i:Next j
End Sub

'---
Sub LL2(i1,j1,i2,j2)
For j=j1 To j2:For i=i1 To i2
TT(i,j)=(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1))/4
Next i:Next j
End Sub

No.8 回答者： lycaon 回答日時： 2009/03/22 11:43

Sub test2()

Dim T(50,50),TT(50,50)
t1=Time
For i=1 To 50:T(i,1)=200:T(i,50)=30:Next
For j=2 To 49:T(1,j)=120:T(50,j)=80:Next
For j=2 To 49:For i=2 To 49: T(i,j)=120:Next:Next '■仮温
'** 計算
123 For k=1 To 20
For j=2 To 49:For i=2 To 49
T(i,j)=(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1))/4
Next:Next
Next k
'** 収束確認
For j=2 To 49:For i=2 To 49:TT(i,j)=T(i,j):Next:Next
For j=2 To 49:For i=2 To 49
T(i,j)=(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1))/4
Next:Next
s=0
For j=2 To 49:For i=2 To 49
b=Abs(TT(i,j)-T(i,j))
If b>s Then s=b
Next:Next
If s>0.01 Then '■
r=r+1
DoEvents
GoTo 123
End If
'** 結果出力
With ActiveSheet
For j=1 To 50:For i=1 To 50
.Cells(i,j).Value=T(i,j)
Next:Next
t2=Time
.Cells(52,1).Value="計算=" & r*20 & "回 " & t1 & "～" & t2
End With
MsgBox ("終了")
End Sub

No.7 回答者： lycaon 回答日時： 2009/03/22 11:30

【EXCEL VBA】

非定常の前に駆け足で VBAの初歩を。

(1)新しいExcelを開く。
(2)ツール－マクロ－Visual Basic Editor(下図)－左欄のSheet1(Sheet1) をクリック。
(3)右欄空白に、次のコードを記述。
----
Sub test1()
x=InputBox("x の y乗を求める。x＝", x)
y=InputBox(x & " のy乗を求める。y＝", y)
MsgBox (x & "^" & y & "=" & x ^ y)
End Sub
----
(4)右向き▲を押す。（大きなx,yを入れるとエラー。気にせず終了。）

マクロ実行無効のエラーが出たら
http://www.ne.jp/asahi/hishidama/home/tech/excel …

　　　　☆★☆★☆★☆★☆

【例題4】50cm角の2次元板を、左辺200℃、上辺120℃、下辺80℃、右辺30℃に保つときの板内温度分布を、Excel VBA で求めよ。

上の(1)～(4)手順で次回答のコードをVB Editor へコピーし実行。
難しい関数は不使用で文法解説省略。VBAは本もサイトも豊富です。

333MHz CPU（遅）のPCで収束計算540回、約1.5分。
(仮温の行がないと980回、約2.5分。収束条件：Max｜T(ｎ+1回目)-T(ｎ回目)｜<0.01℃を0.5℃に変えると40回で打ち切り収束不良。）

結果のグラフ（次回）は Excel仕様で上下逆。

・Excel VBAからcsvファイルに書き出し可能。
・Excel VBAはExcelと関係なく技術計算可能。
・VBでも、配列T(i,j)をcsvに書き込めば同じコードを使用可。
・VBの中からOLE経由でExcelセルに配列を書き込みグラフ化等も可能。

No.6 回答者： lycaon 回答日時： 2009/03/22 10:49

lycaon、African wild dogです。

続けます。

【境界条件について】
直交座標で境界がx軸やy軸に平行なら、断熱条件（∂T/∂n＝0、ｎは境界に垂直な方向）は境界格子の内外で等温、として解きますが、台形や円錐台、円や球のように境界が軸に平行でない場合、差分法では境界条件を正しく設定できません。

下図で点Ａが境界の内側なら、その温度は点Ｂ、Ｃ，Ｄ，Ｅ各点温度の平均。
点Ａがx、y軸に対して斜めの境界線（赤線）上の場合、例えば点ＢとＤの温度を各2倍して4で割ります。これは境界外に架空の格子Ｃ，Ｅを置き、ＣはＢ，ＥはＤと等温としたのと同じ。
ところが本来は境界に垂直な線上に点を選ぶべき。つまり点Ｐ、Ｑという、格子上にない点を使うべきで、当然誤差が出ます。

対策の一つとしてΔx、Δyを極力小さくします（格子を細かくする）。
下底100℃、上底20℃の台形で高さを100分割すれば、図のＣ点とＱ点（図では≒Ｆ点）に平均0.8℃程度の温度差があります。1000分割なら平均0.08℃の差になるので、まあ、Ｃ点とＱ点の温度は同じと見ていいでしょう。不安なら更に分割を増せばよい。
ところが分割数を2倍にすると、計算時間は2次元で概ね4倍（3次元で8倍）、10倍で100倍、100倍で1万倍（3次元なら100万倍!）必要です。

対策とし、境界付近だけ格子を細かくする方法等があります。
これは Excel のセル内に式を埋め込む方法では実現困難で、Excel VBA などのプログラムで可能となります。

（境界付近は熟練プログラマでなければ全体の格子数を増し、計算中は一晩寝て待つ方が、トータル時間を節約できるでしょう。
プログラムは一旦完成すれば次回からは楽。ところが最初はデバッグで一晩どころか何日も潰す可能性があります。）