数値流体力学入門  

 

 理想流体・ポテンシャル流れ

 

理想流体(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とする。入口、出口ではψを直線的に変化させる。

BASICプログラム例(一部)

10 DIM P(NJ,NI) 格子点のψ値
------------
3000 FOR J=SY TO EY
3010 FOR I=SX TO EX
3020 PP=(P(J,I+1)+P(J,,I-1)+P(J+1,I)+P(J-1,I))/4
3030 DD=PP-P(J,I): P(J,I)=PP: IF ABS(DD).>EPSI THEN K=1
3040 NEXT I: NEXT J
-----------------

(注:実際はこの程度の問題ではプログラムを組まずとも表計算ソフトで解ける。)

 

 差分法によるナビエ・ストークス式の数値解法

 

粘性を考慮したナビエ・ストークスの式も上記、差分法・有限要素法を基礎に解かれる。現在ではパソコン上でも動く汎用の流体数値計算プログラム(コード、多くは有限要素法)があり、ユーザーがプログラムを書くことは無くなっている。

差分法による粘性流体の数値計算法のひとつの方法は概略以下のようである。(生井・井上「粘性流体の力学」p.53)簡単のため非圧縮性2次元流れにおいて、ナビエ・ストークス式の代わりに渦度ζ(Zeta)による渦度輸送方程式を用いる。

流れ関数ψの定義を代入すると、

となる。一方、渦度と流れ関数の関係は次式。

この二つの式を差分式にすると、次式となる。

さらに、上の式のnn+1に置き換えて整理すると、

である。

以上の3つの差分式により、
ある時刻t =t0+nΔtにおける各格子点の流れ関数の値ψi,j(n) が与えられたとすると、第一式より渦度ζi,j(n)が得られる。
それらと第二式から次の時間の渦度ζi,j(n+1) が計算できる。
最後に第3式から次の時間の流れ関数ψi,j(n+1)を計算するのであるが、この式は最初のポテンシャル流れの差分式同様、全格子点にわたる連立方程式となるので、繰り返し計算法が必要である。

 

 有限要素法によるポテンシャル流れの数値解法

 

有限要素法(finite element method, FEM)は、流れの領域を差分法よりも粗い大きさの異なる要素(多くは三角形)に分割することで値を離散化し、基礎式における残差を要素の節点で最小にすることで解を求める方法である。有限要素法の理論的根拠を与えているのが重み付き残差法(weighted resisual method)とガレルキン法(Galerkin’s method)である。FEMは差分法に比べて繰り返し計算を必要としない、計算点数が少ない、複雑な形状に適用できるなどの特徴によりいろいろの分野で数値計算に用いられている。汎用のコードも多い。

先ず領域注の三角形要素一つを考えて、求める流れ関数ψ(x,y)を既知の関数Ni(x,y)と未知の節点量ψ123により離散化する。(各頂点の座標x1,y1,x2,y2,x3,y3 は既知とする。添え字1.2.3は各頂点をあらわす。)すなわち要素内のψが、一次式:


で表わされるとして、連立方程式:(i =1,2,3) を解いて、αiを求め、ψiでまとめると、次式となる。

       

       

       

=2×(三角要素面積)

これで要素内の任意の点Pにおけるψは節点値ψiの重み付け平均として求められる。

Niは形状関数あるいは内挿関数といい、連続かつ微分可能なものである。要素形状、節点数により種々ある。要素内の点の位置は直交座標(x, y)を使う代わりに(N1,N2,N3)で与えることもでき、これを面積座標とも呼ぶ。上式はψが面積座標による荷重平均で与えられることを示している。

このように節点で離散化したψiを基礎式に代入すると、これらが解でない場合は残差Rがでる。

重みつき残差法は基礎式の厳密解を求めることはあきらめ、この残差が最小となるような節点値の近似解を求めようとするものである。

誤差の取り方はいろいろ考えられるが(例えば二乗和)誤差に重みつけをおこなって加算することにする。その重みつけの分布を示す重み関数をW(x, y)とすると、重みつけ残差の総和を0とする条件、すなわちψが解ならば、

が成立する。そこで未知節点値ψiと同数の互いに独立な関数Wiを選び、それらのもとに上式が成立すれば、

である。これはn個のψiに関する同次n元連連立方程式である。Wiの選び方によりこの式は種々の形になる。

ガレルキン法はWiに要素の形状関数Niそのものを用いる方法である。すなわち、

部分積分して、

これがガレルキン法の基礎式となる。この式中の各項は要素について次のようになる。

ただし、

これを基礎式に代入して、N1について書くと、左辺は次式である。

[]内はx, yの関数ではないのでこれは簡単に積分できて、

となる。(他の形状関数の場合は、被積分関数がxyの関数であるので、積分はもっと複雑になり、数値積分公式を用いる必要がある。)よって1要素についての基礎式はN1,N2,N2について全て書いて、

である。左辺をnと書く。は節点座標のみから決まり、要素に固有な行列である。よって基礎式は各要素の面積分をすべての要素について加えあわせることで計算できることになった。

境界要素の取り扱い

基礎式の右辺 は自然境界S2に対してのみ存在する。たとえば2-3辺上のqnの値が一定の場合は、s-n系とx-y系の関係が、

したがってS2上の積分すなわち、

中のxys (2-3上ではn=0となる)を代入して積分すると次のようになる。

したがってこの要素の平衡方程式は

((3×3)(3×1)=(3×1))

となる。

流速u,vの計算

すなわち、 (右辺の係数を[B]マトリックスと呼ぶ)

以上で要素のマトリックス およびが得られたので、後は全体のマトリックスを組み立てて、連立一次方程式を解く。

【例題とプログラム例】円柱まわりの理想流体流れ。

 

(解析解は

 

【格子点】x=0-100, y=0-80の場に半円柱を置く。節点総数32、要素総数43。

節点番号 座標x,y
10 0
2 0 17
3 0 37
4 0 55
5 0 80
6 30 0
7 30 24
8 30 42
9 30 59
10 30 80
11 60 0
12 56 16
13 53 32
14 51 47
15 50 62
16 50 80
17 66 22
18 76 33
19 74 41
20 72 52
21 71 66
22 70 80
23 90 39
24 88 48
25 86 56
26 85 67
27 85 80
28 100 40
29 100 50
30 100 59
31 100 69
32 100 80

要素番号 頂点の節点番号(反時計まわり)
1 : 1 6 2
2: 2 6 7
3: 2 7 8
4: 2 8 3
5: 3 8 4
6: 8 9 4
7: 4 9 10
8: 4 10 5
9: 6 11 12
10: 6 12 13
11: 6 13 7
12: 7 13 8
13: 13 4 8
14: 8 14 15
15: 8 15 9
16: 9 15 10
17: 15 16 10
18: 11 17 12
19: 12 17 13
20: 17 18 13
21: 18 19 13
22: 13 19 20
23: 13 20 14
24: 14 20 15
25: 20 21 15
26: 15 21 22
27: 15 22 16
28: 18 23 24
29: 18 24 19
30: 19 24 20
31: 24 25 20
32: 25 26 20
33: 20 26 21
34: 21 26 22
35: 26 27 22
36: 23 28 24
37: 28 29 24
38: 29 30 24
39: 24 30 25
40: 25 30 26
41: 30 31 26
42: 26 31 32
43: 26 32 27

【境界条件】
節点 1 6 11 17 18 23 28; ψ=0 (基本境界条件S1)

節点 5 10 16 22 27 32; ψ=2 (S1)

節点 28 29 30 31 32の線分; - (∂ψ/∂x)=0 (自然境界条件S2)

【自然境界条件の取り扱い方】

各要素の部分マトリックスは、

[42]:

[41]:

[38]:

[37]:

【BASICプログラム】

[変数リスト]

NJ:節点総数(32), NE:要素総数(43), MTJ(I,J):要素番号Iの三角形要素の3頂点の節点番号(43,3), LANK(I):全体の剛性マトリックス[K]等の未知節点番号の並び, KAIK(I):[KAA] の未知節点番号の並び

EK(I,J):要素の影響係数マトリックス[Ke], SK(I,J):全体の影響係数マトリックス[K], COM(I,J):[B]マトリックス

F(I):ψの値 [SK()]{F()}={PB()}, P(I):座標 P(IJ):x座標、P(IJ+NJ):y座標

[マトリックスの構成]

⇒LANK

⇒KAIK

[プログラムの構成]

“INPUT”入力
“CLASFY”全体のマトリックスがどんな順で並んでいるかの指標をつくる
“CALSYK”全体の影響係数マトリックス[K]を計算して[KAA]をつくる
“CALELK”
“SYSK”
“BCON”
“DINV”

[INPUT]
100’
110 READ THICK, NJ, NE 厚さ、節点総数、要素総数の入力
111 DATA 1,32,43
120 FOR JE=1 TO NE: READ KE: FOR K=1 TO 3: READ MTJ(JE,K):NEXT K:NEXT JE
要素番号 IEの3頂点の節点番号入力
130 FOR IJ=1 TO NJ: READ KJ, P(IJ),P(IJ+NJ):NEXT I
x,y座標の入力
140 READ NL
141 DATA 5
150 FOR I=1 TO NL: READ ML(I):NEXT I: FOR I=1 TO NL: READ QN(I):NEXT I
自然境界条件 節点数、番号、値
160 READ NB
161 DATA 16
170 FOR I=1 TO NB: READ MB(I):NEXT I: FOR I=1 TO NL: READ FK(I):NEXT I
値指定節点 数、番号、値ψ

[CLASFY]
200’
210 FOR I=1 TO NJ: LANK(I)=I:NEXT
220 IJ=0: FOR I=1 TO NJ: FLAG=0: FOR K=1 TO NB: IF LANK(I)=MB(K) THEN FLAG=1
230 NEXT K: IF FLAG=0 THEN IJ=IJ+1: KAIK(IJ)=LANK(I)
240 NEXT I
250 IX=NJ-NB:FOR K=IX+1 TO NJ: KAIK(K)=MB(K-IXX):NENXT K
[CALSK]
900’
910 FOR I=I TO NJ
920 FOR J=1 TO NL
921 IF ML(J)=I THEN PB(I)=QN(J)
922 NEXT J: NEXT I
940 FOR JJ=1 TO NB
945 IF MB(JJ)=1 THEN F(I)=FK(JJ)
950 NEXT JJ
970 FOR JE=1 TO NE
980 GOSUB 600 要素[Ke]を作って全体[K]に加え込む
985 NEXT JE
990 GOSUB 800全体のマトリックスから既知節点を抜く
995 RETURN

[CALEK] AND [SYSK]
600 ‘
610 IX=MTJ(JE,1):JX=MTJ(JE,2):KX=MTJ(JE,3):IY=IX+NJ:JY=JX+NJ:KY=KX+NJ
620 S=(P(JX)*P(KY)+P(KX)*P(IY)+P(IX)*P(JY)-P(KX)*P(JY)-P(IX)*P(KY)-P(JX)*P(IY))/2
630 GOSUB 400
  400’
  410 COM(1,1)=(P(JY)-P(KY))/2/S:COMM(2,1)=(P(KX)-P(JX))/2/S
  420 COM(1,2)=(P(KY)-P(IY))/2/S:COMM(2,2)=(P(IX)-P(KX))/2/S
  430 COM(1,3)=(P(IY)-P(JY))/2/S:COMM(2,3)=(P(JX)-P(IX))/2/S
  440 RETURN
640 FOR I=1 TO 3: FOR K=1 TO 3: T=0: FOR J=1 TO 2: T=T+COM(J,I)*COM(J,K):NEXT J: EK(I,K)=T*S*TH:NEXT K: NEXT I
[Ke]の計算([要素マトリックの例] 要素21 S=93.0
  [ 1.40 -1.32 -0.00] [F(18)]
  [-1.32 1.42 -0.10] [F(19)]
  [-0.00 -0.10 0.10] [F(13)])

670 FOR I=1 TO 3
680 IF I=1 THEN IT=IX
690 IF I=2 THEN IT=JX
700 IF I=3 THEN IT=KX
710 FOR J=1 TO 3
720 IF J=1 THEN JT=IX
730 IF J=2 THEN JT=JX
740 IF J=3 THEN JT=KX
750 SK(IT,JT)=SK(IT,JT)+EK(I,J)
全体マトリックッス[SK]への加え込み
760 NEXT J: NEXT I
770 RETURN

[BCON]
800 ‘
810 FOR I=1 TO NJ: FOR K=1 TO NJ
815 IF LANK(I)<>KAIK(K) THEN 836
820 F1(K)=F(I):P1(K)=PB(I)
825 FOR J=1 TO NJ: FOR L=1 TO NJ
830 ID LANK(J)<>KAIK(L) THEN 835
831 B(K,L)=SK(I,J)
835 NEXT L: NEXT J
836 NEXT K: NEXT I
840 FOR I=1 TO NJ:F(I)=F1(I):PB(I)=P1(I):NEXT
845 IA=NJ-NB
850 FOR I=1 TO IA
855 T=0 : IB=IA+1: FOR J=IB TO NJ: T=T+B(I,J)*F(J):NEXT J
860 PB(I)=PB(I)-T: FOR J=1 TO IA:SK(I,J)=B(I,J):NEXT J
865 NEXT I
870 RETURN

[BCON]でおこなったこと

⇒(並び替え)
⇒(公式より) ([KAA]がSK(),右辺がPB())

最終的に解くべき式は次のようになる。[KAA]{ψ}={P}


演習レポート 12

【差分式】一次元非定常の伝熱の式は、時間 t と位置 x に関する偏微分方程式である。

(θ:温度[K],λs:固体の熱伝導度[J/(m・s・K)]) (αsssCsp:固体の熱拡散率[m2/s])

整数p, n により時間をt =pΔt, 位置を x =nΔx で区切り,Tnp [K]を数値解における温度(節点値)として基礎式を差分化せよ。