偏微分方程式の数値解法


差分法による理想流体流れの数値解法

理想流体(ideal fluid)とは流体の粘性と圧縮性を無視した流れである。理想流体を仮定すれば、運動方程式(ナビエ・スト−クス式)と連続の式の組を簡略化できる。理想流体の渦なし流れ(irrotational flow)をポテンシャル流れ(potential flow)という。特に2次元のポテンシャル流れの解析は複素関数を用いた数学的に整った形で流れを記述することができるので、本格的な流体力学の教科書では多くの部分をこれに割いている。(たとえば航空機の翼の断面形状と揚力との関係が純数学的に導かれる。最も「流体力学」らしい部分である。)

いま2次元速度成分u, vと次の関係にある流れ関数(stream function)ψ(Psi)を導入する。

これを連続の式: に代入すると、となるので、連続の式は自動的に満足されている。さらに渦なしの条件: に代入しても、

となり、満足する。この形の式をラプラス(Laplace)の式といい、ポテンシャル流れの基礎式である。

流れ関数一定の曲線は流線(stream line)である。(流線とはその線上の接線がその点における速度の方向と一致する曲線をいう。)また、二つの流線間の流量は流れ関数の差に等しい。流れ関数はポテンシャル流れに限らず粘性流れでも存在する。

境界条件 ラプラスの式の境界条件は基本境界条件S1 (ψ)と自然境界条件 S2(ψの微分値)がある。

(1)S1:壁では境界上でψが一定値となるのでその値を指定する。
(2)S2: 領域の入り口および出口では法線速度が与えられる。ここでは
 を用いる。

 

差分法(finite difference mehtod, FDM)は解析領域を等間隔メッシュに分割し、基礎式を差分方程式に変換して解く方法である。各格子点での流れ関数ψ(x,y)の値をψi,j であらわすと、各方向の二次微分は、

のように差分化される。これをラプラスの式に代入して、

が得られる。ラプラスの式は2階の偏微分方程式であるが、差分方程式で書くと「ある点の関数値はその周囲の4点の値の平均値」という簡単な関係となる。この関係を全格子点について解く。すなわち基本的には各格子点でのψに関する連立方程式となる。

計算上はこのような大規模な連立方程式を直接・同時に解くことは困難なので、仮の初期値から出発して近似計算をおこない、収束した値を解とする。繰り返し計算はi,jの小さい値からおこない、n回目の繰り返し計算での値をψnとする。するとψi-1,jとψi,j-1はn回目の値であり、他はn-1回目の値である。

この計算法で繰り返し計算が進行する。この、繰り返し計算が先行した格子については直ちに新しい値を使って近似計算をおこなう方法をガウス・ザイデル法という。

例えば右のような格子点について解析をおこなう。境界条件は壁の境界では流れ関数一定で、下側でψ=0、上側でψ=1とする。入口、出口ではψを直線的に変化させる。

【実習】表計算によるポテンシャル流れ pd02.xls, 雛形 pd02_temp.xls (セルコピーのマクロつき)

 

偏微分方程式の数値解法

 

1 熱伝導の基礎式(1次元)

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

 (θ:温度[K],λs:固体の熱伝導度[J/(msK)])
同じく x+Δx 面を通過する熱流束は,
 
である.両者の差がこの体積の温度上昇Δθに必要な熱:

ρsCpsΔθΔxΔyΔz に等しいことから,

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

2 差分化

 いま,θを厳密解における温度(連続値)[Kまたは℃]、T [Kまたは℃]を数値解における温度(節点値)とする。整数p,n により時間をt=pΔt, 位置を x=nΔx で区切ると,Eq.(1) の各項が以下のように差分化される。(この差分化法は Explicit)

これより、よって差分化された基礎式が次式となる。
TP+1n    =Θx (T Pn+1 + T Pn-1) + (1-2Θx)T Pn  (4)
(次の時間)    (今の時間)
()
この式により下図のような手順で時間方向に次々と節点の値(温度)が求まる.なお,解が発散しないためにはΘx < 1/2の条件が必要である.

3 境界条件

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

@境界面温度を与える
 θx=0 = θ0(t)    (5)

A断熱境界条件:    ()x=0=0
T P+10 = 2ΘxT P1 + (1-2Θx)T P0  (6)
B境界面熱流束を与える:   -λs()x=0= q0(t)
T P+10 = 2ΘxT P1 + (1-2Θx)T P0 + 2A1Θxq(t)   (7)
(A1= Δx/λs)
C線形熱伝達: -λs()x=0= h(θ-θx=0) (h:伝熱係数[J/(m2sK)])
T P+10 = 2ΘxT P1 + (1-2Θx)T P0 +2B1Θx(T-T P0)   (8)
(B1=Δxh/λs)

【実習例題】) アクリル板の冷却

 図の形状(板幅 l =0.2m ×0.2 m,厚み 0.02 m) のアクリル板の冷却を考える。初期温度は一様で 373.2 K(=100)とする。冷却法を次の3種で計算する。
A: 十分にかく拌した273.2 K(=0)の水に漬ける
B: 0℃の静止空気中に放置
C: 層流(Re<1x105)の空気流れ中に置く(うちわで扇ぐ).u [m/s]は各自推定.


使用する物性値は以下の値を使用。温度依存性を無視する。
アクリルの物性値:(添え字s)  密度 ρs=1180 kg/m3, 熱容量 Cps=1.38×103 J/(Kkg),
熱伝導度 λs =0.15 J/(msK), 熱拡散率 αs=9.21×10-8 m2/s
空気の物性値:(添え字なし)  密度 ρ=1.293 kg/m3, 粘性係数 μ=1.71×10-5 kg/(ms),
熱伝導度 λ=2.44×10-2 J/(msK), Pr=0.72

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

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

【実習】表計算による計算 pd01.xls、雛形:pd01_temp.xls