非定常熱伝導の数値解法　−偏微分方程式−

　物質，熱などの移動現象を厳密に扱うにはナビエ・ストークスの方程式などの偏微分方程式を直接解く必要がある．特殊な例以外は偏微分方程式を解析的に解くことはできない。普通は求める値を離散化して，偏微分方程式を連立方程式に還元して解くことがおこなわれる．その手法には各種あって、有限差分法，有限要素法等がある．偏微分方程式の数値解法は流れのシミュレーション手法として発展し，気象予測等に応用されている．
　ここでは最も簡単な例として熱伝導の偏微分方程式を差分化して解いてみる．

実習例題　ｱｸﾘﾙ板の冷却

　図の形状(板幅 l =0.2m × 0.2m,厚み 0.02m) のｱｸﾘﾙ板の冷却を考える。初期温度は一様で 373.2K(=100℃)とする。冷却法を次の３種で計算する。

A:　十分にかく拌した273.2K(=0℃)の水に漬ける
B:　0℃の静止空気中に放置
C:　層流(Re <1x105)の空気流れ中に置く(うちわで扇ぐ)．u [m/s]は各自推定．

使用する物性値は以下の値を使用。温度依存性を無視する。
ｱｸﾘﾙ板：(添え字s) 密度 ρs=1180kg/m3, 熱容量 Cps=1.38x103J/(K･kg),
熱伝導度 λs =0.15J/(m･s･K), 熱拡散率 αs=9.21x10-8m2/s
空気:(添え字なし) 密度 ρ=1.293kg/m3, 粘性係数 μ=1.71x10-5kg/(m･s),
熱伝導度 λ=2.44x10-2J/(m･s･K), Pr =0.72

熱伝導の基礎式（１次元）

　x 方向の一次元の熱伝導を考える．Δx,Δy,Δzの立体要素につき，時間Δt の間に x 面を通過する熱流束は，

(θ:温度[K],λs:固体の熱伝導度[J/(m･s･K)])
同じく x+Δx 面を通過する熱流束は，

である．両者の差がこの体積の温度上昇Δθに必要な熱: ρsCpsΔθΔxΔyΔz に等しいことから，

すなわち、

(αs=λs/ρsCsp :固体の熱拡散率[m2/s])
これが１次元の熱伝導の基礎式で，時間 t と位置 x に関する偏微分方程式である。

偏微分方程式の差分法による数値解法


　いま，θを厳密解における温度（連続値）[Kまたは℃]，T [Kまたは℃]を数値解における温度（節点値）とする。整数p, n　により時間をt =pΔt, 位置を　x =nΔx で区切ると，Eq.(1) は，


(この差分化法は Explicit 型という)
よって差分化された基礎式が次式となる．

この式により図のような手順で時間方向に次々と節点の値（温度）が求まる．なお，解が発散しないためにはΘx < 1/2の条件が必要である．

境界条件

　数値計算で最も重要で困難な問題は境界条件の取り扱いである．熱伝導における代表的な境界条件と差分式は以下のようである。

例題の解

各冷却法における板表面の境界条件は：
A:　表面温度一定
B:　線形熱伝達，自然対流　（水平に置かれた正方形板について)
伝熱係数 h　[J/(m2･s･K)]=1.27(Δθ/l )1/4
(l :板長さ(=0.2m), Δθ:空気温度と表面温度の差(θx=0-θ∞))
C:　線形熱伝達，強制対流，層流
h =(λ/l )0.664Re1/2Pr1/3 (Re=ρul /μ)
(D:　線形熱伝達，強制対流，乱流 h=(λ/l )0.036Re0.8Pr1/3)

温度分布は板中心面に対象と考え，板中央は断熱境界条件となる。
区間分割は片側(1cm)についてのみn=10,Δx =0.001m(=1mm)とする。時間幅Δxは解が発散しないための条件
Θx=αsΔt ／(Δx )2 < 1/2
からΔt=5sとする。(なお h >100では計算不能）

プログラム

100 '= ﾋﾃｲｼﾞｮｳ ﾃﾞﾝﾈﾂ =
110 '"UHC"
120 '-------[ﾃｲｽｳ]
　　　　　　　n Δx l Δt　λs Cps ρs 最大時間
130 READ N, DX, L, DT, RS, CPS, , DS, TM
140 DATA 10, 0.001,0.2, 5, 0.15, 1.38E3, 1180, 3600
150 DIM TP(N), T1(N)
160 '-------[ｼｮｷﾁ]
170 T0=100
180 FOR I=0 TO N
190 TP(I)=T0
200 NEXT
205 AS=RS/CPS/DS
206 OX=AS*DT/DX/DX
209 P=0
210 '
220 FOR I=1 TO N-1
230 T1(I)=OX*(TP(I+1)+TP(I-1))+(1-2*OX)*TP(I)
240 NEXT
250 '-------[ｷｮｳｶｲ 0]
260 T1(0)=2*OX*TP(1)+(1-2*OX)*TP(0)
270 '-------[ｷｮｳｶｲ n]
271 TB=0
275 '--------------- (A) ｲｯﾃｲ ｵﾝﾄﾞ
276 ' T1(N)=TB
280 '--------------- (B) ｾﾝｹｲﾈﾂﾃﾞﾝﾀﾂ　(ｼｾﾞﾝﾀｲﾘｭｳ)
281 ' H=1.27*((TP(N)-TB)/L)^(1/4)
282 ' B1=DX*H/RS
283 ' T1(N)=2*OX*TP(N-1)+(1-2*OX)*TP(N)+2*B1*OX*(TB-TP(N))
285 '----------------(C,D) ｾﾝｹｲﾈﾂﾃﾞﾝﾀﾂ (ｷｮｳｾｲﾀｲﾘｭｳ)
286 H=100 (各自変更）
287 B1=DX*H/RS
288 T1(N)=2*OX*TP(N-1)+(1-2*OX)*TP(N)+2*B1*OX*(TB-TP(N))
290 '
300 FOR I=0 TO N
310 TP(I)=T1(I)
320 NEXT
330 '
340 P=P+1
341 TIME=P*DT
350 IF (TIME/600)<>INT(TIME/600) THEN 400
360 '--------[ｼｭﾂﾘｮｸ]
370 PRINT "RETｦｵｾ";:BEEP 3:PRINT TIME/60;" min";
380 FOR I=0 TO N:PRINT USING " T(#)=###.#";I;TP(I); :NEXT:PRINT "WAITﾁｭｳ[RET]ｦｵｽ":WAIT
400 '
410 IF TIME>=TM THEN END
420 GOTO 210

実習レポート様式　（A4方眼紙にて提出）

10分おきに内部の分布を描く。（図１(A),(B),(C)）
各場合の中心温度の径時変化を描く（図２）

　

inserted by FC2 system