1
目次
第 1 章 差分解法・数値流束
基礎方程式 : : : : : : : : : : :
特性方程式 : : : : : : : : : : :
差分法 : : : : : : : : : : : : : :
差分 : : : : : : : : : : :
陽解法と陰解法 : : : : :
数値流束 : : : : : : : : : : : : :
精度 : : : : : : : : : : : : : : :
スキーム : : : : :
スキーム : : : : : :
法 : : : :
空間1次精度のスキーム
安定性 : : : : : : : : : : : : : :
数値流束に関するまとめ : : : :
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.3.1
1.3.2
1.5.1 FTCS
1.5.2 Lax
1.5.3 Lax-Wendro
1.5.4
1.7.1 FTCS : : : : : : : : : :
1.7.2 Lax-Wendro スキーム :
1.7.3 Lax スキーム : : : : : :
1.7.4 空間1次精度のスキーム
3
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
第 2 章 ポリト ロピック・ガスの流体力学数値解法
流束ヤコビアン A : : : : : : : : : : : : : :
基本量に対する流束ヤコビアン M :
対角化 : : : : : : : : : : : : : : : :
法 : : : :
法の場合の結果 : : : : : : : :
による空間2次精度化 : : : : : :
内挿 : : : : : : : : : : : :
法の空間2次精度化 : : : : : : : : : :
流束制限関数 : : : : : : : : : : : : : : : :
時間2次精度化 : : : : : : : : : : : : : : :
多次元問題 : : : : : : : : : : : : : : : : :
2.1
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
: : : : : : : : :
2.1.1
2.1.2
2.2 FDS (Flux Dierence Splitting)
2.2.1 FDS
2.3 MUSCL
2.3.1 MUSCL
2.4 FDS
2.5
2.6
2.7
2.8 FVS (Flux Vector Splitting)
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
3
4
6
6
7
7
8
8
10
10
11
13
15
16
16
17
17
19
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
19
20
20
21
23
24
24
26
27
28
28
31
2
第 3 章 流体力学の風上差分法
断熱気体の流体力学の差分解法
オイラー方程式 : : : : :
特性速度 : : : : : : : : :
法: : : : : : : : : : : : : :
法: : : : : : : : : : : : : :
3.1
3.1.1
3.1.2
3.2 FDS
3.3 FVS
33
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
33
33
34
37
41
付 録 A 衝撃波管の解析解
45
等温衝撃波 : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
等温の場合の
不変量 : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
等温の場合の衝撃波管の解析解 : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
A.1
A.2
A.3
45
45
46
Riemann
付 録 B 計算のこと
B.1 Makele : : : : : : : : : : : : : : : : : :
B.1.1 基本的な Makele を書く方法 : :
B.2 Gnuplot : : : : : : : : : : : : : : : : : :
B.2.1 データの作成法 : : : : : : : : : :
B.2.2 プロット法 : : : : : : : : : : : :
49
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
49
49
51
51
51
3
第 1 章 差分解法・数値流束
1.1
基礎方程式
まず最初の基礎方程式は、ある体積の中に含まれる流体の質量が単位時間に流れ込む質量流速によっ
て増減するという連続の式
R
Z
@
dV
0 vdS
@t
S
から得られる。右辺を体積積分に変換し、V として微小体積を考えれば、
=
V
@
@t
+ div(v) = 0
(1.1)
(1.2)
が得られる。
つぎに、ある体積の中に含まれる流体の運動量は、質量と同じように、単位時間に流れ込む運動量流
速によって増減するが、それに加えて運動量の場合は、この流体の体積に加わっている「力」によって
も増減する。流体の中でかならず考慮しなければならない圧力による力は、圧力勾配に比例するので、
@
R
V
vdV
@t
=0
Z
S
vv 1 dS 0
Z
V
gradpdV
(1.3)
となる。ここで、右辺を体積積分に変換し、V として微小体積を考えれば、
@v
@t
+ divvv = 0gradp
(1.4)
が得られる。この式は添字をつけて書くと、
@vi
@t
=
=
@p
i vj
+ @v
= 0 @x
@x
j
i
(1.5)
=
と書ける。もちろん、x1 x 、x2 y 、x3 z を表している。
さて、断熱の場合のエネルギーに関する方程式は、
@e
@t
+ div(e + p)v = 0
(1.6)
ここで e は単位体積あたりの全エネルギーで
e
= jvj2 =2 + (1.7)
第1項は単位体積当たりの運動エネルギー、第2項はおなじく熱エネルギーを表す。理想気体の場合は
p= である。
= ( 0 1)
4
1
第 章 差分解法・数値流束
+ div( ) = 0
質量と同じように、全エネルギーの増減はエネル[email protected]
ev
となる
@t
はずであるが、そうではない。熱力学の第1法則で断熱の場合を考えると、内部エネルギー U と体積 V
は
+ p dV
=0
dt
(1.8)
+ div(v) = 0pdivv
(1.9)
dU
dt
という関係で変化する。これから、単位体積あたりの熱エネルギー は
@
@t
(1.4)
という関係にしたがって変化することが簡単な計算でわかる。これと式
から得られる運動エネル
ギーの変化を表す(この式の右辺が単位体積・単位時間に流体素片になされた仕事を表すことに注意)
2 + div( jvj2 v) = 0v 1 gradp
(1.10)
2
の和をとれば全エネルギーに関する方程式 (1.6) が得られる。
流体力学の基礎方程式は、式 (1.2) 、(1.4) 、(1.6) ということになる。1次元の流体の基礎方程式は、
@jvj2 =
@t
3つの式をベクトルのように並べて書いて、
@Q
@t
+ @E
=0
@x
(1.11)
と書けるただし、Q と E は物理量をならべたベクトルで、
0
Q
=
1
B
C
@ u A ; E
e
0
=
1
u
B
C
@ p u2 A ;
e pu
+
( + )
(1.12)
= a2 (1.13)
である。
1.2
特性方程式
ガスが断熱であると仮定すると、
p
=
@p
@
s
であるから、基礎方程式は、オイラーの運動方程式が
のように書き直せる。
であるから、
@v
@t
@v
@
+ v @x
= 0a2 @x
(1.14)
a2
= p = K01
(1.15)
2
d
=
01
と従属変数 を a に書き直すことができる。式( 1.14 )は
@v
@v
2a @a = 0;
+
v
+
@t
@x 0 1 @x
da
a
(1.16)
(1.17)
1.2.
5
特性方程式
t
C
F=const
x
図
1.1: 特性線とリーマンの不変量
のように書き直せる。同じように従属変数 を a に書き直すと、連続の式、
@
@t
は、 で割って a を掛けることによって、
2
0
1
@a
@t
+ @v
= 0;
@x
(1.18)
@v
@a
+ a @x
+ 20v 1 @x
= 0;
(1.19)
となる。
式
(1.17) と式 (1.19) を加え合わせると、
2 a + (v + a) @ x + 2 a = 0;
@
x+
@t
01
@x
01
が、式 (1.17) から式 (1.19) を引くと、
@
2 a + (v 0 a) @ x 0 2 a = 0;
v0
@t
01
@x
01
(x; t) 上の曲線 C に沿って、x と t の関数 F (x; t) が一定値をとるという条件は、
dF
=
=
@F
@F
dx;
@x
@t
dx @F
@F
dt
dt
@t
dt @x
(1.20)
(1.21)
+
+
=0
(1.22)
すなわち、
@F
@t
@F
+ dx
=0
dt @x
(1.23)
6
1
第 章 差分解法・数値流束
が成り立つことであることがわかる。ここで曲線 C の勾配が dx
dt である。したがって、式
は
(1.21)
dx
dt
( )
という x; t 上の曲線に沿って、
= v 6 a;
(1.20) と式
(1.24)
= v 6 02 1 a;
(1.25)
が一定であることを示している。この J6 をリーマン不変量と呼ぶ。流体力学の基礎方程式が式 (1.20)
と式 (1.21) のように波動の伝搬を表す式に書き直せることは重要である。現在の流体力学の数値計算法
J6
はこのことに基礎をおいて作られているといっても過言ではない。
1.3
差分法
微分の定義式
df dx x=x0
1
f (x0 + 1x) 0 f (x0 )
= 1lim
x!0
1x
(1.26)
' f (x0 + 11xx) 0 f (x0 )
(1.27)
で、 x を小さくとって、近似した式
df dx x=x0
のことを差分と呼ぶ。差分法とは、微分方程式をこのような形に差分化し、その解を求める方法をいう。
1.3.1
差分
空間微分の差分化を考える。空間を格子に分割しその格子点上の物理量 Qj を計算する。ここで j は
格子点につけた番号を表す。
@Q @x j
@Q @x j
@Q @x j
= Qj+110x Qj ;
= Qj 01Qx j01 ;
= Qj+1210xQj01 :
(1.28)
(1.29)
(1.30)
2
この3つの表現は x が小さい極限ではおなじ dQ
dx を与えるが、 x のオーダで違う量を与える(詳し
くは後に述べる)
。式
のようにとるのを前進差分、式
のようにとるのを後退差分、式
のようにとるのを中心差分と呼ぶ。前進差分
と後退差分
を使って、2階微分は
1
(1.28)
1
(1.29)
(1.28)
(1.29)
@ 2 Q Qj +1 0 2Qj + Qj 01
=
@x2
1x2
(1.30)
(1.31)
のような3点差分を使うことが多い。
どのような差分をとるかということは、計算方法の安定性、収束性などに影響を与える。これはこの
授業の主要な内容にあたる。
1.4.
1.3.2
式
7
数値流束
陽解法と陰解法
(1.12)
@Q
@t
+ @E
= 0;
@x
(1.32)
を時間について差分化する方法も、前進、後退、中心差分を考えることができる。時間に関する添字を
n とし、n ステップめの物理量から n
ステップめのそれを求める(時間進化を求める)ことを考え
る。これが流体力学の非定常問題の解を求めることに相当している。
時間に関しては前進差分、中心差分そして後退差分を適用すると、それぞれ
+1
Qnj +1
n
= Qnj 0 1t @E
;
@x
(1.33)
n
Qnj +1
;
= Qnj 01 0 21t @E
@x
(1.34)
Qnj +1
= Qnj 0 1t @[email protected]
(1.35)
n
+1
;
のようになる。ここで、前の2つは Qnj +1 を求めるのに、過去の情報 E n E Qn や Qn01 のみを使っ
ている。したがって既に知れている情報で求める1つ時間ステップ進んだ後の物理量の値を決めること
ができる。このようにな方法を陽解法と呼ぶ。それに対して、式
では、E n+1
E Qn+1 だか
ら、1つ時間ステップ進んだ後の物理量の値を決めるのに、それを右辺に含んでいる。このような解法
は陰解法とよばれ、1つ時間ステップ進んだ後の物理量の値を決めるのに、それら自身の間にある関
係を使わなければならないことになっている。具体的には、陰解法では、空間の格子点数 従属変数
の種類だけの変数に関する連立方程式を1ステップごとに解く必要がある。
= ( )
(1.35)
= (
)
2
1.4
数値流束
(1.33) (1.34) (1.35)
前節の式
、
、
では、流束の空間微分に関する部分をどうとるかということには触
れなかった。流体力学の数値解法は、つまるところ、空間微分に関する差分と時間微分に関する差分を
どのようにとるか、ということに集約できる。
式
を
(1.12)
Qnj +1 0 Qnj
t
1
+
~
~
Ejn+1=2 0 Ejn01=2
1x
=0
~
~
(1.36)
このような差分に書き直した時、グリッド の中心での値である Qj に対して、Ej +1=2 や Ej 01=2 はグリッ
ド の境界での値であるべきである。
このときどのようにして Ej +1=2 などの量を決めれば良いのだろうか。
ここで、Ej+1=2 を Ej E Qj 、Ej 01=2 を Ej 01 E Qj 01 とすれば、空間後退差分・時間前進差分
のスキームが Ej +1=2 を Ej+1 E Qj +1 、Ej 01=2 を Ej
E Qj とすれば、空間前進差分・時間前進
E Qj +1 = 、Ej 01=2 を
差分のスキームができる。空間中心差分のスキームには、Ej +1=2 を E Qj
E Qj
E Qj 01 = のようにすれば良いことは容易にわかるだろう。
数値流体力学の計算方法とは、この Ej +1=2 などの流束の値を、どのよう決めればよいかという問題
に帰着する。このような立場に立ってスキームに現れる Ej +1=2 などの流束のことを数値流束と呼ぶ。数
値流速の決め方こそ数値流体力学の計算方法のポイントといえる。数値計算スキームによってこの数値
流束のとりかたは異なるが、いくつかの例について後に見ることになる。
~
( ( )+ (
~
= ( ) ~
= ( ) ~
~
)) 2
~
= (
)
= ( )
~
( ( )+ (
~
)) 2 ~
8
1
第 章 差分解法・数値流束
Ej-1/2
Qj-1
Ej+1/2
Qj
図
Qj+1
1.2: 物理量と流束の関係
1.5 精度
1.2 節で見たように、流体力学の方程式系はある変換によって波動方程式に書き換えることができ、
これが、後に風上差分法(この意味は後述 )を作る上での基礎となっている。そこで、式 (1.20) や式
(1.21) 型の方程式のプロトタイプ
@ @v
+ @x = 0;
(1.37)
@t
で、v = c =const とした
@
@
+ c @x
= 0;
(1.38)
@t
の型の方程式の数値解法を考えよう。
1.5.1 FTCS スキーム
時間前進差分・空間中心差分の
を
は、式
(1.38)
FTCS(Forward in Time and Central Dierence in Space) スキーム
nj +1 0 nj
t
1
のように差分化して
nj +1
+ c j+1210xj01 = 0
n
n
1t (n 0 n );
= nj 0 c 21
j 01
x j+1
1 として解く方法である。Pascal のプログラムでは
2 for j:=0 to jm do
3 w[j]:=u[j];
4 for j:=1 to jm-1 do
5 u[j]:=w[j]-nu/2*(w[j+1]-w[j-1]);
(
)
(1.39)
(1.40)
= 1 1
のよう書けば良い。 プログラム中でnu と書かれている量 はクーラン数と呼ばれ、 c t= x で
定義される量である。1ステップで進む時間の長さと、2つの格子点の間の距離の間の関係をつけてい
る量であるといえる。
1.5.
9
精度
1.3: FTCS スキームで、初期値として、j = 1; : : : ; 50 に対して = 1 、j = 51; : : : ; 100 に対して
= 0 としクーラン数 c1t=1x = 0:25 で 50 ステップ、100 ステップ計算した時の をプロットし
た (左) 。同じ初期条件で、クーラン数 = 0:5 で 25 ステップ、50 ステップ計算したものが右図である。
図
このスキームで、初期値として、
(
j = 1; : : : ; 50
= 10;; :: :: :: for
(1.41)
for j = 51; : : : ; 100
をとり、クーラン数 = c1t=1x を 0.25 と 0.5 の場合について 2 Nstep = 25 まで計算した結果が図
1.3 に示してある。この初期条件に合致する厳密解は、密度が 1 から 0 へ飛ぶ位置が、1ステップごと
に右に だけ移動してゆくものにならなければならないことは容易にわかるだろう。しかし、図 1.3 に
示したように、この FTCS スキームを用いて計算すると、厳密解の回りに余分の振動が発生してしま
うことがわかる。
この差分化にどれくらいの誤差が混入しているかを考える。時間方向にテイラー展開して、
n
@ 2 1t2
@ n
+1
n
= + @t 1t + @t2 2 + O(1t3 );
n
また空間方向にテイラー展開して、
(1.42)
1x2 + @ 3 1x3 + O(1x4);
j +1 = j +
1x +
(1.43)
2 @x3 j 6
@ @ 2 1x2
@ 3 1x3
j 01 = j 0
1x + @x2 2 0 @x3 6 + O(1x4 );
(1.44)
@x j
j
j
となる。O1t3 は、1t3 に比例する項もしくはそれよりも高次の項( 1t3 より小さい)が省略されてい
ることを示す。これを (1.39) に代入すると、(1.37) と完全に同じになるわけではなく、
n
@
@ 3 1x2
@ 2 1t
@
+ c @x + c @x3 6 + @t2 2 = 0;
(1.45)
@t
@ @x j
@ 2 @x2 j
j
となる。このように、空間方向に 1x2 、時間方向に 1t の誤差を含むスキームを、空間2次精度、時間
1次精度のスキームと呼ぶ。
10
1
第 章 差分解法・数値流束
1.4: 1.3 とおなじ。ただし採用した計算スキームは Lax 法で、左はクーラン数 = 0:25 、右はクー
= 0:5 のもの。FTCS スキームに比べて、振動はおさえられているが、密度の飛びは鈍って
図
図
ラン数 いる。
1.5.2 Lax スキーム
FTCS
1.3 に見えるように、正しい解の回りで振動を生じており、実際に使うこと
(1.39) の右辺第1項 j を (j01 + j+1)=2 で置き換えた式
nj01 + nj+1
1t (n 0 n );
nj +1 =
0 c 21
(1.46)
j 01
2
x j +1
1 を解く方法が Lax 法(図 1.4 )である。この方法のプログラムは
2 for j:=0 to jm do
3 w[j]:=u[j];
4 for j:=1 to jm-1 do
5 u[j]:=(w[j-1]+w[j+1])/2-nu/2*(w[j+1]-w[j-1]);
この
スキームは図
はできない。そこで、式
FTCS スキームの数値的振動は押えられているが、密度の飛びは著しく鈍って
のように書けば良い。
しまっている。
1.5.3 Lax-Wendro 法
時間2次精度のスキームとするためには、
n
n
n
n
+1 = n + @ 1t + @ 2 1t2 + @ 3 1t3 + O(1t4);
@t
@t2 2
@t3 6
(1.47)
で、
@
@t
@
= 0c @x
;
(1.48)
1.5.
11
精度
1.5: 1.3
= 0 25 50
Lax-Wendro 法で、左はクーラ
= 0 5 25 ステップ、50 ステップ
図
図
とおなじ。ただし採用した計算スキームはオリジナルの
ン数 : で ステップ、
ステップ後、中央はクーラン数 : で
後、右はクーラン数 : で ステップ、 ステップ計算したもの。
100
= 0 95 13
26
@2
@t2
であることに注すると、
@ +1 = n 0 c
2
@ = +c2 @x
2;
(1.49)
c2 @ 2 @x2 j
1t + 2
1t2 + O(1t3 );
(1.50)
j
となる。この式を差分化するとオリジナルの Lax-Wendro スキームが得られる。
c1t n
c2 1t2 n
n
n
n+1 = n 0
(
j +1 0 nj01 ) +
(1.51)
21x
21x2 (j+1 0 2j + j01 );
これは FTCS スキームに第3項の拡散項を加えたものとなっている。拡散項は、FTCS スキームで生
n
1
2
3
4
5
@x
じる振動を拡散によって鎮める効果を持つのでこのスキームはより安定化されることになる。プログラ
ムは
for j:=0 to jm do
w[j]:=u[j];
for j:=1 to jm-1 do
u[j]:=w[j]-nu/2*(w[j+1]-w[j-1])+nu*nu/2*(w[j+1]-2*w[j]+w[j-1]);
とすれば良い。
1.5.4
空間1次精度のスキーム
nj +1
= nj 0 c11xt (nj 0 nj01);
(1.52)
とすれば、時間については1次精度前進差分、空間については1次精度後退差分をとったスキームがで
きる。プログラムは
12
1
2
3
4
5
1
第 章 差分解法・数値流束
for j:=0 to jm do
w[j]:=u[j];
for j:=1 to jm do
u[j]:=w[j]-nu*(w[j]-w[j-1]);
のようにすれば良い。
この問題では、現象は左から右へ伝わってゆくので、それをとらえられるように、空間については後
退差分をとった。図
から、空間1次精度のスキームは、密度の飛びの部分で余分な振動を起こさな
いことがわかる。しかし同時に、密度の飛びの部分が時間が進むにつれて、ぼやけてゆくこともわかる。
このように、空間1次精度のスキームは余分な振動を起こさない良い性質をもっているが、飛びの部
分がぼやけてゆくという悪い性質を持っていることがわかる。
テイラー展開を用いてこのスキームにど
の程度の誤差が含まれているかを見てみると、式
は元の微分方程式
に2つ余分な項を含
んでおり
!
1.6
(1.52)
@
@t
+ 12t @@t2 = 0c
2
と等しくなることがわかる。ここで、@
@t
@
@t
となる。この式で
@
@x
(1.38)
@ 0 12x @x
2
2
(1.53)
@
= 0c @x
に気をつけて変形すると、これは
2
@
@ + c @x
= 2c (1x 0 c1t) @x
2
1t < 1x=c
(1.54)
(1.55)
のようにとると、右辺は正の拡散係数を持つ拡散項であることがわかる。つまりこのようにして作った
空間1次精度のスキームは正の拡散係数を持つ拡散項を差分化にともなう誤差として含むために余分の
振動を生じなくなっていることがわかる。
ここで式
は早さ c で伝わる現象について時間進化を計算する場合に、情報が格子1つ分以上
伝搬する時間以上に時間ステップをとってはならないという条件に相当しており、数値流体力学で非
定常問題を解く場合の時間ステップに関する条件を与えるものである。これは、
条件(
条件)と呼ばれる。
(1.55)
Lewy(1928)
1.6:
1.3
Courant, Friedrich,
CFL
図
図
とおなじ。ただし採用した計算スキームは空間1次精度後退差分をとったスキームで、
クーラン数は左から : 、
: 、
: 、
: をとっている。
= 0 25
=05
=08
= 0 95
1.6.
1.6
13
安定性
安定性
この節の内容は今後の議論に不可欠というわけではないので、最初にこの資料を勉強する場合には飛
ばして先に進んでも良い。
前節で見たように、空間2次精度のスキームでは、余分な振動が発生した。空間1次精度のスキーム
は余分な振動は発生しなかったが、密度の飛びは時間がすすむとともに鈍ってゆくことがわかった。こ
れらの差はどこからくるのだろうか。差分方程式に含まれる元の微分方程式からのずれ(誤差)に起因
することは前に述べた。
ここでは、
の安定性解析を使って余分な振動や解の鈍りと差分スキームの間の関係を
定量的に見てゆく。解をフーリエ分解して、
von Neumann
nj
=
X
gn
exp(ij)
(1.56)
1
=2 1
と書く。ここで j は格子点の番号で、 は格子間隔 x を単位とした波数つまり k x で
を表す。
(これは x だけ進んだ時の位相の変化を与える)
この式をスキームに代入し
ij に比例する項を取り出す。
スキーム
1
exp( )
FTCS
nj +1
に対しては、
gn+1
となる。これから複素振幅 g が
(0 )
= nj 0 2 (nj+1 0 nj01 );
(1.57)
= gn 0 2 gn (exp(i) 0 exp (0i));
(1.58)
= 1 0 i sin ;
(1.59)
= nj0 を満足していなければならないので、
nj +1 = g n+1 exp(ij) = nj0 = gn exp(i(j 0 )) ;
(1.60)
g
となっていることがわかる。
この問題では厳密解は、nj +1
が成り立たねばならず、これから厳密解に対する複素振幅 g は
gexact
j
= exp(0i) ;
(1.61)
j=1
で与えられる。これは当然 gexact
で、厳密解では位相だけが動いてゆくことになる。この gexact か
らどれだけ異なっているかを調べれば、差分化にともなって、微分方程式の厳密解からずれた部分がわ
かることになる。
複素振幅 g の絶対値はある波の1ステップ後との振幅の比を与えるが、式
に対しては
(1.59)
0
jgj2 = 1 + 2 sin2 ;
(1.62)
0 の波は増幅されることがわかる。
Im(g)
= arcsin
(1.63)
で与えられる。 > に対しては波数 >
つぎに複素振幅 g の位相 は、
jgj
14
1
第 章 差分解法・数値流束
で、厳密解の位相に対する比は
arcsin p1+ sin2 sin 2 =
(1.64)
のように計算できる。これを図に書くと図 1.7 のようになる。jg j と =exact を図示する Mathematica
exact
1 のプログラムは
2 G[nu_,theta_]=Sqrt[1+nu^2 Sin[theta]^2]
3 Plot[G[0.5,theta],{theta,0,Pi},Frame ->
1
2
3
True]
および
Phi[nu_,theta_]=ArcSin[nu Sin[theta]/Sqrt[1+nu^2 Sin[theta]^2]]/(nu theta)
Plot[Phi[0.5,theta],{theta,0.0001,Pi},Frame->True]
である。
1.12
1
1.1
0.8
1.08
0.6
1.06
0.4
1.04
0.2
1.02
1
0
0
0.5
1
1.5
2
2.5
1.7: FTCS
3
0
=05
0.5
1
1.5
2
2.5
3
( )
図
スキームで、クーラン数 : でにする複素振幅 g の絶対値 左 と複素振幅 g の位相
と厳密解の位相に対する比(右)。この図で、それぞれ に近いほど厳密解に近いことになる。
0
1
0
1
ここで、複素振幅 g の絶対値が < < ( と 以外)に対して を越えているので、初期値に含
む周期的な波はかならず増幅されることになることがわかる。また、 の大きな(短波長の)波の方で
で短波長の波が後方に取り残されているように見
厳密解に比べて、位相が遅れることもわかる。図
えることもこのように説明される。
同じようにして、
のスキームでは、
1.3
Lax
= cos 0 i sin ;
(1.65)
jgj2 = cos2 + 2 sin2 ;
(1.66)
g
となり、絶対値と位相の比は、
arcsin (cos + sinsin )
2 2 2 1=2
=
()
で与えられる。これを図示する Mathematica のプログラムは
exact
(1.67)
1.7.
1
2
3
4
5
6
15
数値流束に関するまとめ
G[nu_,theta_]=Sqrt[Cos[theta]^2+nu^2 Sin[theta]^2]
Plot[G[0.5,theta],{theta,0,Pi},Frame -> True]
Phi[nu_,theta_]=ArcCos[Cos[theta]/G[nu,theta]]/(nu theta)
Plot[Phi[0.5,theta],{theta,0.0001,Pi},Frame->True]
のようになる。
1
2.2
0.9
2
1.8
0.8
1.6
0.7
1.4
0.6
1.2
0.5
1
0
0.5
1
1.5
2
2.5
1.8: Lax
図
スキームで、クーラン数 の厳密解の位相に対する比(右)
。
1.8
3
0
0.5
1
1.5
2
2.5
3
= 0:5 に対する複素振幅 g の絶対値 (左) と複素振幅 g の位相 jj 1
図 (左)を見ると、すべての で g < であり、周期的な振動はその振幅が減少してゆくこと
がわかる。特に、 = つまり グリッド で 波長となる波が大きく減衰することがわかる。図
(右)からは、 の大きい波に対して、正の位相誤差とを生じており、
スキームで飛びが鈍ってゆく
につれて前方に崩れてゆく構造がこれが原因となっていると見ることができる。
また、1次のスキームに対しても同様の解析は可能で、
= 2
4
1
1.8
Lax
= 1 0 (1 0 cos ) 0 i sin ;
jgj2 = (1 0 (1 0 cos ))2 + 2 sin2 ;
arcsin sin
j
gj
= ;
(1.68)
(1.69)
(1.70)
g
exact
1.9
で、図
で複素振幅 g の絶対値は、大きい波数に対して減少しており、短波長の波ほど減衰しやすい
ことがわかる。
: の場合は の場合とともに複素振幅 g の位相誤差は になるが(図
右)、
: < < では =exact > で前進位相誤差を、 < < : では =exact < で遅延位相誤差を
持つ。
05
1.7
1
=05
1
=1
0
0
05
1
1.9
数値流束に関するまとめ
1.4
ここまで、いくつかのスキーム(差分法)について見てきた。これを
節で導入した数値流束とい
う面からまとめておく。以降では差分式に出てくる従属変数 ではなく u と書く。
16
1
第 章 差分解法・数値流束
1
2
0.8
1.5
0.6
1
0.4
0.5
0.2
0
0
0
0.5
1.9: 1.7
1
1.5
2
2.5
3
0
0.5
1
1.5
図
図
と同じ、ただし、1次精度空間後退差分スキームで、クーラン数 。
幅 g の絶対値 左 と複素振幅 g の位相 の厳密解の位相に対する比(右)
( )
2
2.5
3
= 0:5 でにする複素振
1.7.1 FTCS
空間2次精度の
FTCS スキームは、
unj +1
であり、これを
unj +1
= unj 0 2
unj+1 0 unj01 ;
= unj 0 11xt Fj+1=2 0 Fj01=2 ;
(1.71)
(1.72)
= c uj +2uj+1 ; Fj01=2 = c uj +2uj+1 ;
(1.73)
であるとわかる。つまり、この方法では、セルの境界である j 6 1=2 での数値流速として、セルの中心
の形に書き直すと、
Fj +1=2
での値の算術平均を取っていることがわかる。
1.7.2 Lax-Wendro スキーム
Lax-Wendro スキームでは
unj +1
2
= unj 0 2 (unj+1 0 unj01 ) + 2 (unj+1 0 2unj + unj01 );
(1.74)
であるから、数値流速は
= c 1 02 uj+1 + 1 +2 uj ; Fj01=2 = c 1 02 uj + 1 +2 uj01 ;
(1.75)
となる。セルの境界での数値流速として、セルの中心での値に (1 0 )=2 と (1 + )=2 の重みをつけて
Fj +1=2
平均を取っていることがわかる。
1.7.
17
数値流束に関するまとめ
1.7.3 Lax スキーム
unj +1
n
n
= uj01 +2 uj+1 0 2 unj+1 0 unj01 ;
(1.76)
であり、
1
0
1
=
1
+
1
=
1
0
1
=
1
+
1
=
(1.77)
Fj +1=2 = c
2 uj+1 + 2 uj ; Fj01=2 = c 2 uj + 2 uj01 ;
となる。セルの境界での数値流束として、セルの中心での値に (1 0 1= )=2 と (1 + 1= )=2 の重みをつ
けて平均を取っていることがわかる。
1.7.4
空間1次精度のスキーム
空間1次精度のスキームでは風上側を含んだ差分を取る必要がある。すなわち、c >
ては後退差分、
をc<
unj +1
= unj 0 unj +1
= unj 0 0( < 0) に対しては前進差分、
0( > 0) に対し
unj 0 unj01 ;
(1.78)
(1.79)
unj+1 0 unj ;
を取る必要がある。このようにとるスキームを1次精度風上差分とよぶ。これは、1つの式で書くこと
もできて、
= unj 0 +2 jj unj 0 unj01 0 02 j j unj+1 0 unj ;
(1.80)
0j j
+j j
2 が、 > 0 の時は と等しく、 < 0 の時は 0 となることは容易にわかるだろう。 2 は、 < 0
の時は と等しく、 > 0 の時は 0 となる。数値流束は
c 0 jcj
c + jcj
c 0 jcj
c + jcj
Fj +1=2 =
(1.81)
2 uj+1 + 2 uj ; Fj01=2 = 2 uj + 2 uj01 ;
unj +1
となるが、この式は
Fj +1=2
= Fj+12+ Fj 0 j2cj (uj+1 0 uj ); Fj01=2 = Fj +2Fj01 0 j2cj (uj 0 uj01 );
のようにも書き直せる。
(1.82)
19
第 2 章 ポリト ロピック・ガスの流体力学数値解法
まず、等温を代表とするポリトロピックガスに関する流体力学を解く解法について見てみる。これはた
だちに一般の気体に関する解法に応用できる。
圧力 p が密度 だけで決まっているような場合、すなわち p p のようになっている性質をバロ
トロピーと呼ぶが、その特別の場合で圧力が p K のように密度 のベキに比例する場合をポリト
2 03 <n< 10 03 星間気体の中では
ロピックと呼ぶ。
は等温の場合に当たり密度の高い
良い近似であることが知られている。
=
=1
2.1
= ()
(10 cm 10 cm )
流束ヤコビアン A
以下では特に述べない限り等温の場合に限って議論を進める。等温音速を
= (p=)1=2;
a
(2.1)
( 1.12) は等温の場合は、
@Q @E
+ @x = 0
@t
と表すと、流体力学の基礎方程式 式
と書ける。ただし、Q と E は物理量を並べたベクトルで、
Q
である。ここで、u
=
u
!
;E
u
=
=
( + u2 )
m
!
=
a2 m
;
(2.3)
(2.4)
;
に書き直すと、E の方はこれを使って
E
!
a2
= m と書いて従属変数を
Q
(2.2)
!
+ m2
;
(2.5)
のように書ける。質量密度、運動量密度、エネルギー密度は(大きな体積で積分すれば保存する量なの
で)保存量と呼ばれる。
式
を
(2.1)
@Q
@t
+ A @Q
=0
@x
(2.6)
X
@Ek @Ql
@Ql @x
(2.7)
のように書き直すことを考える。
@Ek
@x
=
l
l
= Ak;l @Q
@x
20
2
第 章 ポリトロピック・ガスの流体力学数値解法
のように書きなおせて、ここに出てくる
Ak;l
k
= @E
@Q
(2.8)
l
を保存量 Q に対する流束ヤコビアンと呼ぶ。具体的に書き下すと、
=
Ak;l
すなわち、
@
@t
これが、保存量で書いた @Q
@t
2.1.1
m
!
0
2
a2 0 m
2
0
+
2
a2 0 m
2
2m
;
@
@x
m
+ A @Q
@x = 0 型の方程式である。
基本量に対する流束ヤコビアン
(2.9)
2m
!
1
!
1
!
= 0;
(2.10)
M
質量密度、速度という基本量 q で方程式を書くこともできて、このときの基礎方程式は、シンボリッ
クに
@q
@t
@q
+ M @x
= 0;
(2.11)
と書く。M は基本量に関する流束ヤコビアンである。具体的な形では、
u
@
@t
これが、基本量で書いた @q
@t
2.1.2
!
u
a2
+
u
!
@
@x
u
!
= 0;
(2.12)
@q
+ M @x
= 0 型の方程式である。
対角化
(2.10)
(2.12)
式
や式
を流束ヤコビアンが対角成分のみを持つ形に変形することを考える。
そのためには、流束ヤコビアンの固有値と固有ベクトルを求める。簡単な計算(つまり det A I
を満たす を計算する。)から、A 、M のこの固有値は
( 0 )=0
1
2
=
=
(2.13)
(2.14)
u 0 a;
u
+ a;
であることがわかる。それぞれの固有値に属する固有ベクトルは、A では
x A1 =
また、M では、
xM 1 =
1
!
u0a
0a
1
; x A2 =
u+a
!
; xM 2
a
= +
!
(2.15)
;
!
(2.16)
;
となる。これから、固有ベクトルを並べた行列(右固有行列 R )およびその逆行列 R01 は、それぞれ、
A では
!
!
RA
=
1
u0a
1
u+a
; RA 01
u+a
1
= 02ua0a 012a
2a 2a
;
(2.17)
2.2. FDS (Flux Dierence Splitting) 法
また M では、
RM
=
21
!
0a a
; RM 01 =
1
2
0 21
0 21a
!
1
2a
(2.18)
;
となる。
右固有行列 R およびその逆行列 R01 を用いると、流束ヤコビアンは
RA 01 ARA
や
= 3;
(2.19)
= 3;
(2.20)
RM 01 M RM
と対角化される。ここで
3=
である。もちろん、逆変換
1
2
!
(2.21)
A
= RA 3RA01 ;
(2.22)
M
= RM 3RM 01;
(2.23)
や
が成り立つ。
2.2 FDS (Flux Dierence Splitting) 法
波動方程式に関する差分化で、式 (1.82) つまり
Fj +1 + Fj jcj
Fj + Fj 01 jcj
Fj +1=2 =
2 0 2 (uj+1 0 uj ); Fj01=2 = 2 0 2 (uj 0 uj01);
(2.24)
を手本として、流体力学の式を風上差分化することを考える。つまり、
Ej +1=2
01
01
= Ej+12+ Ej 0 RAj32jRA (Qj+1 0 Qj ); Ej01=2 = Ej +2Ej01 0 RA j32jRA (Qj 0 Qj01); (2.25)
1
である。
法では RA R0
= で、Ej 01=2
A の部分は、Ej +1=2 については、j
については、j
= での値で評価する。
= での物理量は、その左側の格子点での値を添字 L をつけて、また、右側の格
この方法では、j
子点での値を添字 R をつけて表すと、
Flux Dierence Splitting
01 2
61 2
j3j
ave
uave
Roe
+1 2
= ppLR ; p
= pLuL ++ pR uR ;
L
R
(2.26)
+1
のように決める。この方法を、
の近似リーマン解法と呼び、時間 n ステップの j および j
のふた
ステップの時に j
=
つのグリッド 内部が格子点上の値に等しい一定の値を持っていた時の時間 n
でとるであろう値を衝撃波管( リーマン)問題として解き、それを数値流束を計算する時に必要な中間
+1
+1 2
22
2
第 章 ポリトロピック・ガスの流体力学数値解法
1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
値とする。このときに厳密に t だけ時間を進めるリーマン問題を解く代わりに、それを近似する量で
置き換えて中間値に用いる方法を近似リーマン解法と呼ぶ。
プログラムを見てみよう。これは n ステップの密度
、速度
の値から、 t だけ時間の進
んだ n
ステップの値を求めるサブルーティンである。 のついた量は
平均をあらわし、他の
、
、
は、固有値、右固有行列、およびその逆行列をあらわす。また、 に対する、
に、運動量 m に対する、数値流束が
に求められて時間推進に使われる。 は等
数値流束が
温音速をあらわす。
+1
Lambda(:) R(:,:) L(:,:)
E1(:)
rho(:)
p
E2(:)
subroutine FDS()
#include "param.h"
#include "dim.h"
#include "var.h"
! Local Variable
real*8 rhop(0:Nx), up(0:Nx), E1(0:Nx), E2(0:Nx)
integer j,n1,n2,k
! Matrix
real*8 Lambda(2),R(2,2),L(2,2),w(2,2)
! Roe Average
do j=0, jm-1
rhop(j)=Sqrt(rho(j)*rho(j+1))
up(j)=(Sqrt(rho(j))*u(j)+Sqrt(rho(j+1))*u(j+1))
*
/(Sqrt(rho(j))+Sqrt(rho(j+1)))
end do
rhop(jm)=rhop(jm-1)
up(jm)=up(jm-1)
do j=0, jm-1
Lambda(1)=Abs(up(j)-cs)
Lambda(2)=Abs(up(j)+cs)
R(1,1)=
R(1,2)=
R(2,1)=
R(2,2)=
1d0
1d0
up(j)-cs
up(j)+cs
L(1,1)= 0.5d0*(up(j)/cs+1d0)
L(1,2)= -0.5/cs
L(2,1)= -0.5d0*(up(j)/cs-1d0)
L(2,2)= 0.5/cs
u(:)
Roe
1
cs
2.2. FDS (Flux Dierence Splitting) 法
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
23
do n1=1, 2
do n2=1, 2
w(n1,n2)=0.0
do k=1, 2
w(n1,n2)=w(n1,n2)+R(n1,k)*Lambda(k)*L(k,n2)
end do
end do
end do
E1(j)=0.5*(m(j+1)+m(j)
*
-w(1,1)*(rho(j+1)-rho(j))
*
-w(1,2)*(m(j+1)-m(j)))
E2(j)=0.5*(m(j+1)**2/rho(j+1)+m(j)**2/rho(j)
*
+cs2*(rho(j+1)+rho(j))
*
-w(2,1)*(rho(j+1)-rho(j))
*
-w(2,2)*(m(j+1)-m(j)))
end do
do j=1, jm-1
rho(j)=rho(j)-dtdx*(E1(j)-E1(j-1))
m(j)=m(j)-dtdx*(E2(j)-E2(j-1))
end do
rho(0)=rho(1)
m(0)=m(1)
rho(jm)=rho(jm-1)
m(jm)=m(jm-1)
do j=0, jm
u(j)=m(j)/rho(j)
end do
end
222
do
w(n1,n2)
j3j
2 2 2 2 2 do
ここで
の配列要素
には、R R01 を格納する。通常は
の
ループを書か
ないで、 ループで行なわれていることを展開してプログラムするが、ここでは定式化と同じになる
ようにプログラミングしてある。
2.2.1 FDS 法の場合の結果
次に数値計算結果を見る。
24
2
第 章 ポリトロピック・ガスの流体力学数値解法
2.1: Roe Flux Dierence Splitting 法による衝撃波管問題の解。初期条件は同じ。左は空間1次精
MUSCL の = 01 の場合である)。太い線は計算値、細い線は解析解
図
の
度、右は空間2次精度(ただし
の値である。
これは初期状態
(
= 10;:1; j =j =1;4140; 100
(2.27)
uj = 0; j = 1; 100
(2.28)
とし、1t = 0:25 、ステップ数 n = 0; 40; 80; 120 すなわち、t = 0; 10; 20; 30 の状態をプロットしたもの
である。1x = 1 なので a1t=1x = 0:25 である。
j
この場合の解析解については、補章を参照せよ。
2.3 MUSCL による空間2次精度化
2.3.1 MUSCL 内挿
セルの境界での物理量をセルの中央でのそれから内挿する方法を考える。図
xj 01=2 x xj +1=2 に対して、物理量 u を内挿する。
"
1 (x 0 xj )j u + 3 (x 0 xj )2 0 1x2
u(x) = uj +
1x
21x2
12
ただし、
= uj+1 02 uj01 ;
j2 u = uj +1 0 2uj + uj 01 ;
j u
#
2.2 に示したように、
j2 u;
(2.29)
(2.30)
(2.31)
2.3. MUSCL による空間2次精度化
25
uj+1
uRj+1/2
uj
uRj-1/2
uLj+1/2
uj-1
uLj-1/2
図
である。右辺第3項で
2.2: MUSCL 内挿法。
1x2 =12 という項がつくのは、
Z x
j +1=2
xj01=2
1x=2
2
x3 2
= 112x
(x 0 xj ) dx = 3 01x=2
(2.32)
()
+1 2
= uj + 21 j u + 4 j2 u;
= uj + 41 (1 0 )(uj 0 uj01 ) + 14 (1 + )(uj+1 0 uj );
であるから、この区間について、u x を積分して得た平均値と、uj を等しくするためである。このよ
= での値は、左つまり j 側から補間すると、
うに区分多項式で近似したときの、j
(uL)j+1=2
+ 1 側から補間すると、
(uR )j+1=2 = uj+1 0 21 j+1 u + 4 j2+1u;
= uj+1 0 41 (1 + )(uj+1 0 uj ) 0 14 (1 0 )(uj+2 0 uj+1 );
この表式は の値によって、異なった精度の内挿公式となっている。 = 01 とすると、
(uL)j+1=2 = uj + 12 (uj 0 uj01);
(uL)j01=2 = uj01 + 21 (uj01 0 uj02);
右つまり j
であるが、これから
+ 21 (uj 0 2uj01 + uj02 );
8
2
!
< @2u
@u
1
' 1x 4 @x + 2 1x : @x2
0
j
j 01
2
! 3
2
3
@u
@ u 5
= 1x 4 @x
0 12x @x
;
3
j
j
(uL)j+1=2 0 (uL)j01=2 =
uj 0 uj 01
!
@ 2u
@x2
j
(2.33)
(2.34)
(2.35)
(2.36)
93
=
5;
;
(2.37)
26
2
第 章 ポリトロピック・ガスの流体力学数値解法
であるから、この表式は空間2次精度になっていることがわかる。従って、左から uL を外挿して得る
のに適した風上の3点( j
、j
、j )を使ってできる2次の完全風上差分
となっていることがわかる。
また、
とすると、
02 01
=0
(2nd-order fully upwind)
+ uj+140 uj 0 uj01 04 uj02 ;
u
0 2uj + uj01 + uj 0 2uj01 + uj02 ;
uj 0 uj 01 + j +1
4 !
!4
2
2
2
2
1x @ u + 1x @ u + 1x2
= 1x @u
0
@x j
2 @x2 j 4 @x2 j01 4
3
2
2 @3u !
1
x
@u
5;
= 1x 4 @x 0 4 @x3
j
j 01=2
(uL)j+1=2 0 (uL)j01=2 =
=
(2.38)
uj 0 uj 01
1
@2u
@x2
!
;
j
(2.39)
(2.38)
となる。これが空間2次精度になっていることは誤差が x2 となっていることからわかり、式
が左側から内挿したときに、風上側3点 j
、j
、j と風下側1点 j
をつかって表した2次精度
であることがわかる。
の風上重点差分
とすると、式
は
02
=1
01
+1
(2nd-order upwind biased)
(2.33)
(uL)j+1=2 = uj +2uj+1 ;
(2.40)
のように境界での値はその中央での値の算術平均に一致する。
2.4 FDS 法の空間2次精度化
前 2.2 節で見た FDS 法は、空間1次精度である。これを2次精度にするには 2.3 節で取り上げた
MUSCL 法を用いれば良い。
FDS 法に MUSCL 法を適用するには、次のようにすればよい。
(1)j + 1=2 の物理量を評価する Roe 平均にとられる QR と QL を MUSCL で内挿されたグリッド 境界
の値をとる (式 (??) 、(??)) 。j + 1=2 の物理量を、左側から内挿して求める場合は、
h
i
(QL)j+1=2 = Qj + 14 (1 0 )10j + (1 + )1+j ;
(2.41)
であり、右側から内挿して求める場合は、
1 =
h
i
(QR )j+1=2 = Qj+1 0 14 (1 0 )1+j+1 + (1 + )10j+1 ;
1 = 0
1
ここで、 + Qj +1 Qj であり、 0 Qj Qj 01 である。ただし、実際には 6 をそのまま使う
のではなく、後述する流束制限関数を用いて「鈍らせた」ものを 6 のかわりに計算に用いる。
内挿をどの物理量に対して行なうかは任意性がある。経験的には、基本量(密度、速度、圧
力)を内挿して、それから他の量を求める方がうまく行くといわれている。
プログラムで変更する部分は、まず、
内挿を行なう部分で
0
MUSCL
1
2
3
(2.42)
Muscle
! Muscle
1
2.5.
27
流束制限関数
4
5
6
7
8
9
call Muscle(rho,rhoL,rhoR)
call Muscle(u,uL,uR)
do j=0,jm
mR(j)=rhoR(j)*uR(j)
mL(j)=rhoL(j)*uL(j)
end do
Muscle(w,wL,wR)
wL(:) ( )
w(:) に格納された物理量の Muscle 内挿し、(QL)j+1=2 の
wR(:) に格納するサブルーティン副プログラムである。
となる。ここで
は、配列
に、 QR j +1=2 の値を配列
値を配列
次に数値流束を
(2)
~
Ej+1=2
から
~
Ej +1=2
= 12
Ej +1
+ Ej 0 Rj3jR01(Qj+1 0 Qj )
= 12 E(QR ) + E (QL) 0 Rj3jR01(QR 0 QL)
(2.43)
(2.44)
1 に高精度化する。この部分は次のようにプログラムに書けて、
2
E1(j)=0.5*(mR(j)+mL(j)
3
*
-w(1,1)*(rhoR(j)-rhoL(j))
4
*
-w(1,2)*(mR(j)-mL(j)))
5
E2(j)=0.5*(mR(j)**2/rhoR(j)+mL(j)**2/rhoL(j)
6
*
+cs2*(rhoR(j)+rhoL(j))
7
*
-w(2,1)*(rhoR(j)-rhoL(j))
8
*
-w(2,2)*(mR(j)-mL(j)))
2.1
となる。図
にこの方法で求めたテスト問題の解を掲げる。空間1次精度の図
されていることがわかる。
2.5
2.1 左に比べて改善
流束制限関数
(2.41)
(2.42)
式
や式
で、右辺が1項目だけであれば、1次精度の差分を与える。第2項目が高次の
精度に寄与しているわけだが、第1章でも見たように、高次精度の差分解法は振動を起こす。これは、
衝撃波面のようなするどい飛びを持つ部分で起こるので、このようなするどい飛びを持つ部分では、自
動的に1次精度に落ち、スムーズな領域では高次精度を保つようにする工夫を考える。これが流束制限
関数と呼ばれる概念である。
内挿にこれを適用する場合には、するどい飛びを持つ部分、ある
6
が になり、スムーズな領域ではそのままの値をとるように
いは、振動を起こしている部分で
MUSCL
1 0
1 6 = B(16)
(2.45)
のようにフィルターを掛けることを考える。
振動を起こしている部分では + と 0 とが異符合であるから、そのような部分では B
また + と 0 とが同符合である場合は、それになるように選択する。
関数を制限関数に使う場合は、式
、
の + および 0 を
1 1
minmod
1
1
(2.41) (2.42) 1
1 +j = minmod(1+j ; b10j );
1
= 0 とする。
(2.46)
28
2
第 章 ポリトロピック・ガスの流体力学数値解法
1 0j = minmod(10j ; b1+j );
で置き換えればよい。ただしここで、b = 310
0 であり、minmod 関数は
8
>
< 0; x 1 y < 0
minmod(x; y) = > x; x 1 y > 0 かつ jxj > jyj
:
y x 1 y > 0 かつ jxj < jy j
(2.47)
(2.48)
のように定義されている。
2.6
時間2次精度化
いままでの計算法は、空間方向の差分の精度では1次と2次のものを見てきたが、時間方向の差分の
精度については1次精度であった。
空間方向と時間方向の差分の精度は独立に議論してはいけないのだが、常微分方程式の修正オイラー
法つまり、dy=dx f x; y を
= ( )
y n+1
=
=
x +1
n
として解くオイラー法に対して、
y n+1
xn+1
=
=
+ f (xn ; yn)1x
xn + 1x
(2.49)
(2.50)
yn
+ f xn + 12x ; yn + f (xn; yn) 12x 1x
xn + 1x
yn
(2.51)
(2.52)
とする修正オイラー法風の方法で2次化する。
Qn から Qn+1 に写す時間1次精度の演算を
Qn+1
とし、時間2次精度の演算を
= Qn + L(Qn)1t;
(2.53)
= Qn + L(Qn) 12t ;
Qn+1 = Qn + L(Qn+1=2 )1t;
(2.54)
Qn+1=2
で実現する。
2.7
多次元問題
ここではポリトロピック・ガスに関する流体力学の多次元問題(2次元以上)をオペレータ分割法を
u; v に対
用いて取り扱う方法をのべる。2次元のポリトロピック・ガスの流体力学の方程式は、v
して、
=( )
@Q
@t
@F
+ @E
+
=0
@x
@y
(2.55)
2.7.
29
多次元問題
2.3: Roe Flux Dierence Splitting 法による衝撃波管問題の解。初期条件は前と同じ、空間、時間
MUSCL の = 0 )の場合。太い線は計算値、細い線は解析解の値である。
図
の
とも2次精度(ただし
0
Q
である。例えば
=
1
B
C
@ u A ; E
v
0
u
1
= B@ a2 + u2 CA ; F =
uv
0
1
v
B
C
uv
@
A;
2
2
a v
+
Lax-Wendro 法について考えれば、O(1t3 ) を無視すると、
!
@Q
1
t2 @ 2 Q
n+1
n
Q
= Q + 1t @t + 2 @t2
(2.56)
(2.57)
で1次元と同じように流束ヤコビアン A と B を使って、
= 0AQx 0 BQy ;
(2.58)
= A(AQx)x + A(BQy )x + B(AQx)y + B (BQy )y ;
(2.59)
Qt
また
Qtt
であることに注意すると、
Qn+1
となることがわかる。
ここで、時間分割法は、これを
と
2
= Qn 0 1t(AQx + BQy ) + 12t [A(AQx)x + A(BQy )x + B(AQx)y + B (BQy )y ] ;
Qn+1=2
Qn+1
(2.60)
= Qn 0 1tExn
(2.61)
= Qn+1=2 0 1tFyn+1=2
(2.62)
30
2
第 章 ポリトロピック・ガスの流体力学数値解法
との連続した2つの1次元の方程式の解法に分解するものである。ただし、このままで式
に代入すると、
Qn+1
式
(2.61) を (2.62)
2
= Qn 0 1t(AQx + BQy ) + 12t [A(AQx)x + 2B(AQx)y + B(BQy )y ] + O(1t3 );
(2.60) と同じにはならない。これは流束ヤコビアン(行列)A と B が可換ではない
AB 6= BA
正確にいえば、
(
A BQy
)x 6= B(AQx)y
(2.63)
(2.64)
(2.65)
であるために発生する現象である。これを避けるには、分ステップを4分割して、その順序を回転させ
れば良い。すなわち、
= Qn 0 12t Exn
1t F n+1=4
Qn+1=2 = Qn+1=4 0
2 y
1t F n+1=2
Qn+3=4 = Qn+1=2 0
2 y
1t E n+3=4
Qn+1 = Qn+3=4 0
2 x
とすればよい。ここで、式 (2.67) と (2.68) は1つの
Qn+3=4 = Qn+1=4 0 1tFyn+1=4
Qn+1=4
(2.66)
(2.67)
(2.68)
(2.69)
(2.70)
とすれば良いことは明らかであろう。
分ステップのオペレータを Lx のように書くと、
Qn+1
= Lx1=2Ly Lx1=2Qn ;
(2.71)
のように書ける。このように書くと、1ステップ進むのに、3分ステップの計算が必要なように見える
1=2
が、もし Qn+1 の値が必要でないなら、直前の Lx 1=2 は Qn+1 の直後に出てくる Lx と一緒に計算して
1 1 1 Lx Ly Lx Ly L1x=2Qn;
(2.72)
のように計算すれば良い。これは3次元にすると、
Qn+1
= Lx1=6Ly 1=6Lz 1=3Ly 1=6Lx1=3Lz 1=6Ly 1=3Lx1=6Lz 1=3Lx1=6Ly 1=3Lz 1=6Lx1=6Qn ;
のようにすれば良いことになる。
(2.73)
2.8. FVS (Flux Vector Splitting)
31
2.8 FVS (Flux Vector Splitting)
流体力学の風上差分解法を作るもう一つの有力な方法が
流束を、式
すなわち
(1.81)
Fj +1=2
Flux Vector Splitting である。これは、数値
= c 02 jcj uj+1 + c +2 jcj uj ; Fj01=2 = c 02 jcj uj + c +2 jcj uj01 ;
(2.74)
これを書き換えると、
Fj +1=2
=
(
for c > 0
for c < 0
cuj ;
cuj +1 ;
Fj 01=2
=
(
for c > 0
for c < 0
cuj 01 ;
cuj ;
(2.75)
を手本に書き直す。つまり、固有値
36 = 3 62 j3j
(2.76)
を考える。3+ は 3 のうち正の固有値のみを持つ行列、30 は負の固有値のみを持つ行列である。この
速度でやりとりされる流束の方も、3+ によって移動する E + と 30 によって移動する E 0 とに分離す
ることを考える。そこで、E + については後退差分を、E 0 については前進差分をとれば自然と風上差
分が実現する。この方法を、
流速分離)法という。空間微分の項は
Flux Vector Splitting(FVS:
~j01=2
E~j +1=2 0 E
@E
=
@x j
1x
+
0
+
0
= (Ej + Ej+1)10x(Ej01 + Ej )
のように書ける。ここで、流束は
E+
= R3+R01Q;
E 0 = R30 R01 Q;
(2.77)
(2.78)
(2.79)
のように計算される。
法にしても
法にしても、1次精度風上差分の自然な拡張になっていることがわかる。
FDS
FVS
33
第 3 章 流体力学の風上差分法
3.1
断熱気体の流体力学の差分解法
これまでの基礎の上に立って、この章では、流体力学の方程式を差分法によって解くことを学ぶ。
3.1.1
オイラー方程式
1次元の流体の基礎方程式は、3つの式をベクトルのように並べて書いて、
@Q
@t
+ @E
=0
@x
(3.1)
と書けるただし、Q と E は物理量を並べたベクトルで、
0
Q
である。ここで、u
=
1
B
C
@ u A ; E
e
0
=
1
u
B
C
@ p u2 A ;
e pu
+
( + )
(3.2)
= m と書いて従属変数を
0
Q
=
1
B
C
@ m A;
e
(3.3)
に書き直すと、E の方はこれを使って
1
0
E
m
B
30 m2 C
B 0 e
2 C
A;
@
10 m23 em
2 = (
1) +
+
(3.4)
のように書ける。質量密度、運動量密度、エネルギー密度は(大きな体積で積分すれば保存する量なの
で)保存量と呼ばれる。
式
を @Q
A @Q
のように書き直すことを考える
@t
@x
(3.1)
+
=0
@Ek
@x
=
X
l
@Ek @Ql
@Ql @x
のように書きなおせて、ここに出てくる
Ak;l
l
= Ak;l @Q
@x
k
= @E
@Q
l
(3.5)
(3.6)
34
3
第 章 流体力学の風上差分法
を流束ヤコビアンと呼ぶ。具体的に書き下すと、
0
Ak;l
=
@m
@
B @ ( 01)e+(30 )m2 =2
B
@
@
@ (10 )m3 =22
@
0
@m
@m
( 1) +(30 )m2 =2
@m
@ (10 )m3 =22
@ 0 e
0
@m
1
0 302 m22
3
0
@m
@e
( 1) +(30 )m2 =2
@e
@ (10 )m3 =22
@ 0 e
1
C
C;
A
@e
1
C
01 C
(3 0 )
A;
3(1
0
) m2
m
m
em
e
( 0 1) 3 0 2 2 2 + 0
1
0
1
0
(3 0 )u 0 1 CCA
= BB@
0 302 u2
( 2 0 1)u3 0 0 1 p u 0 1 p + 3022 u2 u
=
B
B
@
m
保存量 Q(質量密度、運動量密度、エネルギー密度)を使って書いた式の代わりに、基本量 q
(密度、速度、圧力)を用いて非保存形の式を書くこともできて、
0
1
@ B C
@ u A
@t
p
0
u
+ B@ 0
0
u
p
0
1
0
= (; u; p)
1
C @ B
C
= A
@ u A
@x
u
p
1
(3.7)
= 0;
(3.8)
のように書くこともできる。シンボリックな書き方をすれば
@q
@t
@q
+ M @x
= 0;
(3.9)
となる。M は基本量 q に関する流束ヤコビアンである。
3.1.2
式
特性速度
(3.2) や (3.8) を第1章で見た advection の形に書き直すことを考える。すなわち、
@Q01
@t
@Q02
@t
@Q03
@t
0
1 = 0;
+ 1 @Q
@x
0
2 = 0;
+ 2 @Q
@x
0
3 = 0;
+ 3 @Q
@x
(3.10)
の形に変形するのである。これができれば、それぞれの式に第2章でみた高精度風上差分の方法を適応
することによって流体力学の方程式も解くことができる。
式
を上の形に書き直すのは行列 M を対角化するということに相当する。行列 M を対角化する
とは、
(3.8)
Mx
= x;
(3.11)
となる固有ベクトル x と固有値 をもとめることに帰着する。これは、
det (M 0 I ) = 0
(3.12)
3.1.
35
断熱気体の流体力学の差分解法
0
u0
det B@ 0
0
1
0
1
u0
= A
p
u0
C
=0
(3.13)
の特性方程式の解を求めると、
(u 0 )3 0 (u 0 )(p=) = 0;
[ 0 (u 0 c)][ 0 u][ 0 (u + c)] = 0;
であるから、固有値は音速 c = (p=)1=2 を使って、
1 = u 0 c;
2 = u;
3 = u + c;
(3.14)
(3.15)
(3.16)
(3.17)
(3.18)
であることがわかる。
d={{u,rho,0},{0,u,1/rho},{0,g p,u}}
Simplify[Eigenvalues[d]]/.Sqrt[g]Sqrt[p]/Sqrt[rho]-> cs
{u, -cs + u, cs + u}
( )
( )
この速度は、それぞれ、物体の流れる速度 2 、物体に対して左方向に伝わる音波の速度 1 、おなじ
く右方向に伝わる音波の速度 3 である。さらに、固有ベクトルは
( )
= 1 x1 ;
(3.19)
+ = 0;
+ = 0;
+ = 0;
(3.20)
M x1
を解いて、
8
>
>
<
cx1 x2
cx2 1 x3
>
>
: px
2 cx3
から
0
x1 =
となる。同じようにして、
0
x2 = B
@
1
B
C
@ 0c A
c2
1
(3.21)
0
0 CA ; x3 =
0
1
B
C
@ c A;
c2
(3.22)
となることがわかる。この固有ベクトルを横に並べた行列(右固有行列)
0
R
= (x1; x2; x3 ) =
1
B
C
c A
@ 0c
c2
c2
0
0
(3.23)
36
3
第 章 流体力学の風上差分法
および、その逆行列 R01 を考え、M に左から R01 を、右から R を掛けると、
0
R 01 M R
= B@
0 0 1
0 2 0 CA 3
0 0 3
1
(3.24)
となる。
R の逆行列は
R={{rho,rho,rho},{-c,0,c},{rho c^2,0,rho c^2}}
2
2
{{rho, rho, rho}, {-c, 0, c}, {c rho, 0, c rho}}
Inverse[R]//MatrixForm//TeXForm
0
R01
0
01
1 1
2 c2 C
0 c21 A
1
2 c2 2c
= B@ 1 0
0 21c
(3.25)
のようにして求められる。
式
に左から R01 を掛けると、
(3.9)
R 01
@q
@t
@q
+ R01M RR01 @x
=0
(3.26)
となるが、R01 が空間、時間の微小部分ついては一定と考えられるならば、
@R01 q
@t
01
+ 3 @[email protected] q = 0
(3.27)
であり、3つの独立したスカラー方程式に分離できたことになる。もうすこし正確にいうと、
R01
@q
@t
@q
+ 3R01 @x
=0
(3.28)
01 q( R01 は R01 の第 i 行のみでできた行ベクトル)が一定であること
とは、 x= t
i に沿って、Ri
i
01
1
1
; 21c ; 2c12 ; c u;
p
をいっている。たとえば、R101
2c2 ; c; で R1 q 2c2
12 p c u が dx=dt u c に沿って になることを意味している。
c
この独立したスカラー方程式それぞれは R01 q の3つの成分それぞれに対する
の方程式に
u c 、2 u 、3 u c と異なっていることが違ってい
合致している。ただし、その速度が 1
る。それぞれの変数に対して風上差分を考えることは第2章と同じようにできる。これが流体力学に関
する高解像度風上差分の基本である。
なお流束ヤコビアン A から保存量に対する固有行列を求めておくと
1 1 =3
= (1 0 1 )
1
= (0 0
)=
0
= 0
= 0
=
1
1
0
R
= B@
0
R01
= BB@
u0c
H 0 uc
1 0b1 + u 1
2
c
1 0 b1
1 0b1 0 u 1
2
c
1 =
(0 0 1)
b2 u
1 1 0 b2 2 c
0b 2
1 b2
2
1 111 )
advection
= +
1 1
C
u
u+c A
1 u2 H + uc
2
0 12 1c + b2 u 21 b2
(0 1 1
(3.29)
1
C
C
A
(3.30)
3.2. FDS 法
37
である。ただしここで
2
= 02 1 uc2
01
b =
(3.31)
+ AQx = 0;
(3.33)
b1
2
(3.32)
c2
を表す。
3.2 FDS 法
Qt
に左固有行列(右固有行列の逆行列)R01 を左から掛けると、
R01 Qt
よって、R が微小な
+ R01ARR01Qx = 0;
(3.34)
1t 、1x に関して一定とみなせるなら、
@R01 Q
@t
01
+ 3 @[email protected] Q = 0
(3.35)
のように書ける。ここで、風上差分の数値流束を
(3R01Q)j+1=2 = 12 3(R01Q)j+1 + 3(R01 Q)j 0 j3jR01(Qj+1 0 Qj )
左から R を掛けると
R R01 Qj +1=2
3
= 12 R3R01Qj+1 + R3R01Qj 0 Rj3jR01 (Qj+1 0 Qj )
(3.36)
(3.37)
= 12 Ej+1 + Ej 0 Rj3jR01(Qj+1 0 Qj )
(3.38)
と書き直すことができて、FDS( Flux Dierence Spilitting )スキームの数値流速を与える。
ここで Rj3jR01 の項は j + 1=2 についての量なのであるが、FDS スキームでは、Roe の近似リーマ
となるが、この式は
Ej+1=2
ン解法を用いてこの中間値を評価する。
1次元の単純な衝撃波管問題( リーマン問題)の解は解析的に得られる。j および j
のグリッド
の時間 t での物理量がわかっていれば、その後の解を正確にとくことができてこれを j
= の物理量
の評価に用いることができる。これを正確に解くのではなく、近似的に解いて(近似リーマン解法)そ
れを j
= の物理量の評価に用いる。
の近似リーマン解法と呼ばれるものは、境界の左側の物理
量 QL と右側の物理量 QR を用いて、平均の値 Qave を
+1
+1 2
+1 2
Roe
= pLR ;
pLuL + pR uR
uave =
p + p ;
L
R
p H + p H
L L
Have =
p + pR R ;
ave
L
R
(3.39)
(3.40)
(3.41)
38
3
第 章 流体力学の風上差分法
j
j
j+1/2
j+1
j+1
x
3.1:
1
x
図
近似リーマン解法の概念図。物理量がグリッド の内部で空間的に一定だったとして(左の図)、
その t 後の物理量の分布( リーマン問題)を解く。たとえばそれが右の図のようなものであったとし
= の物理量を用いて流束を計算する方法を(1次精度の)
法とい
て、ここで得られた j
う。このようにして得た j
= の物理量を風上差分の数値流束の評価に用いる方法を近似リーマン解
法を用いた差分法と呼ぶ。
+1 2
+1 2
Gonunov
= ( 0 1)(Have 0 12 u2ave);
(3.42)
( H は単位質量あたりの全エンタルピーで H = (e + p)= = u2 =2 + =( 0 1) 1 p= である)のように平
均量( Roe 平均)を得るものである。
c2ave
この評価式は、
流束ヤコビアンに関する @[email protected]
1
2
3
4
5
6
7
8
(1)
= [email protected][email protected] が差分式についてもなり立つこと:
E (QR ) 0 E (QL ) = A(QR ; QL )(AR 0 QL ) = Aave (QR 0 QL ) 。
(2)A(QR ; QL) は A と同じように、実固有値と線形独立な固有ベクトルを持つこと。
(3)A(Q; Q) = A(Q) という対称性を持つこと。
これら( Property U )を満たすものとして Roe により提案されたものである。
プログラムは、第1に
for j:=0 to jm-1 do // Roe Average
begin
rhop[j]:=Sqrt(rho[j]*rho[j+1]);
up[j]:=(Sqrt(rho[j])*u[j]+Sqrt(rho[j+1])*u[j+1])/(Sqrt(rho[j])+Sqrt(rho[j+1]));
Hp[j]:=(Sqrt(rho[j])*H[j]+Sqrt(rho[j+1])*H[j+1])/(Sqrt(rho[j])+Sqrt(rho[j+1]));
cp[j]:=Sqrt(gam1*(Hp[j]-up[j]*up[j]/2))
end;
でrhop[j] などに
1
2
Roe 平均 j +1=2 の評価値を求める。次に、右固有行列 R をR[l,m] に、その逆行列
R01 をL[l,m] に求める。w[n1,n2] は Rj3jR01 が格納される。最後に数値流束をE[j] に求めている。
for j:=0 to jm-1 do
3.2. FDS 法
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
39
begin
Lambda[1]:=Abs(up[j]-cp[j]);
Lambda[2]:=Abs(up[j]);
Lambda[3]:=Abs(up[j]+cp[j]);
R[1,1]:=
R[1,2]:=
R[1,3]:=
R[2,1]:=
R[2,2]:=
R[2,3]:=
R[3,1]:=
R[3,2]:=
R[3,3]:=
1;
1;
1;
up[j]-cp[j];
up[j];
up[j]+cp[j];
Hp[j]-up[j]*cp[j];
Sqr(up[j])/2;
Hp[j]+up[j]*cp[j];
b1:=gam1/2*Sqr(up[j]/cp[j]);
b2:=gam1/Sqr(cp[j]);
L[1,1]:=
L[1,2]:=
L[1,3]:=
L[2,1]:=
L[2,2]:=
L[2,3]:=
L[3,1]:=
L[3,2]:=
L[3,3]:=
0.5*(b1+up[j]/cp[j]);
-0.5*(1/cp[j]+b2*up[j]);
0.5*b2;
1-b1;
b2*up[j];
-b2;
0.5*(b1-up[j]/cp[j]);
0.5*(1/cp[j]-b2*up[j]);
0.5*b2;
for n1:=1 to 3 do
begin
for n2:=1 to 3 do
begin
w[n1,n2]:=0.0;
for k:=1 to 3 do
begin
w[n1,n2]:=w[n1,n2]+R[n1,k]*Lambda[k]*L[k,n2];
end;
end;
end;
E1[j]:=0.5*(m[j+1]+m[j]
-w[1,1]*(rho[j+1]-rho[j])
40
3
第 章 流体力学の風上差分法
図
44
45
46
47
48
49
50
51
52
53
54
55
1
2
3
4
5
6
7
3.2: Roe の Flux Dierence Splitting 法による衝撃波管問題の解。初期条件は同じ。
-w[1,2]*(m[j+1]-m[j])
-w[1,3]*(e[j+1]-e[j]));
E2[j]:=0.5*(gam1*(e[j+1]+e[j])+gam3/2*(Sqr(m[j+1])/rho[j+1]+Sqr(m[j])/rho[j])
-w[2,1]*(rho[j+1]-rho[j])
-w[2,2]*(m[j+1]-m[j])
-w[2,3]*(e[j+1]-e[j]));
E3[j]:=0.5*(gamma*(e[j+1]*m[j+1]/rho[j+1]+e[j]*m[j]/rho[j])
-gam1/2*(Sqr(m[j+1]/rho[j+1])*m[j+1]+Sqr(m[j]/rho[j])*m[j])
-w[3,1]*(rho[j+1]-rho[j])
-w[3,2]*(m[j+1]-m[j])
-w[3,3]*(e[j+1]-e[j]));
end;
この数値流束を用いて時間進化を計算すれば良い。これを使って得られた解を図
for j:=1 to jm-1 do
begin
rho[j]:=rho[j]-nu*(E1[j]-E1[j-1]);
m[j]:=m[j]-nu*(E2[j]-E2[j-1]);
e[j]:=e[j]-nu*(E3[j]-E3[j-1]);
end;
3.2 に掲げた。
3.3. FVS 法
41
3.3 FVS 法
固有値
(3.43)
36 = 3 62 j3j
を考える。3+ は 3 のうち正の固有値のみを持つ行列、30 は負の固有値のみを持つ行列である。この
速度でやりとりされる流束の方も、3+ によって移動する E + と 30 によって移動する E 0 とに分離す
ることを考える。そこで、E + については後退差分を、E 0 については前進差分をとれば自然と風上差
流速分離)法という。
分が実現する。この方法を、
空間微分の項は
Flux Vector Splitting(FVS:
@E
@x
j
のように書ける。ここで、流束は
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
1
2
=
~
~
Ej +1=2 0 Ej01=2
x
+
Ej
Ej0+1 0 Ej+01
x
= (
1
+ ) (
1
+ Ej0)
(3.44)
E+
= R3+R01Q;
E 0 = R30 R01 Q;
(3.45)
(3.46)
のように計算される。
プログラムを見てみよう。右固有行列 R をR[l,m] に、その逆行列 R01 をL[l,m] に求めるところは
と同じである。w[n1,n2] は R + R01 が格納される。最後にこれに伴う数値流束 R + R01 Q を
Ep[j] に求めている。
FDS
3
3
for n1:=1 to 3 do
begin
for n2:=1 to 3 do
begin
w[n1,n2]:=0.0;
for k:=1 to 3 do
begin
w[n1,n2]:=w[n1,n2]
+R[n1,k]*max(0.0,Lambda[k])*L[k,n2];
end;
end;
end;
E1p[j]:=w[1,1]*rho[j]+w[1,2]*m[j]+w[1,3]*e[j];
E2p[j]:=w[2,1]*rho[j]+w[2,2]*m[j]+w[2,3]*e[j];
E3p[j]:=w[3,1]*rho[j]+w[3,2]*m[j]+w[3,3]*e[j];
30 R01 が格納される。最後にこれに伴う数値流束 R30R01Q をEm[j] に求
つぎに、w[n1,n2] は R
めている。
for n1:=1 to 3 do
42
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
1
2
3
4
5
6
7
3
第 章 流体力学の風上差分法
begin
for n2:=1 to 3 do
begin
w[n1,n2]:=0.0;
for k:=1 to 3 do
begin
w[n1,n2]:=w[n1,n2]
+R[n1,k]*min(0.0,Lambda[k])*L[k,n2];
end;
end;
end;
E1m[j]:=w[1,1]*rho[j]+w[1,2]*m[j]+w[1,3]*e[j];
E2m[j]:=w[2,1]*rho[j]+w[2,2]*m[j]+w[2,3]*e[j];
E3m[j]:=w[3,1]*rho[j]+w[3,2]*m[j]+w[3,3]*e[j];
end;
ここで得られた数値流束は E + Ep[j] については後退差分、E 0 Em[j] については前進差分をとっ
て、時間進化を計算することによって、 t だけ進んだ状態を計算することができる。
(
)
1
(
)
for j:=1 to jm-1 do
begin
rho[j]:=rho[j]-nu*(E1p[j]+E1m[j+1]-E1p[j-1]-E1m[j]);
m[j] :=m[j] -nu*(E2p[j]+E2m[j+1]-E2p[j-1]-E2m[j]);
e[j] :=e[j] -nu*(E3p[j]+E3m[j+1]-E3p[j-1]-E3m[j]);
end;
この数値流束を用いて時間進化を計算すれば良い。これを使って得られた解を図
3.3 に掲げた。
3.3. FVS 法
図
43
3.3: Steger Warming の Flux Vector Splitting 法による衝撃波管問題の解。初期条件は同じ。
45
付 録A
A.1
衝撃波管の解析解
等温衝撃波
衝撃波に止まった系で、衝撃波前面の速度を u1 、衝撃波後面の速度を u2 、それぞれの密度を、1 、
2 、とおくと、ランキン・ユゴニオ関係から、
= 2u2 ;
1(u21 + a2 ) = 2(u22 + a2 );
1 u1
ここで、式
(A.1) から
と置くと、式
(A.2) は
となり、これから、
2
1
(
= uu1 = x;
(A.3)
+ 1)x + M12 = 0;
(A.4)
2
x2 0 M12
(A.1)
(A.2)
2
u1
2
2
x = M1 =
=
;
a
1
よって、
u1 u2
= a2
(A.5)
(A.6)
となる。衝撃波の前方は静止していたとき、静止系から見た時の衝撃波後面の速度 U2 は、衝撃波面の
進む速度を Vs u1 として
=
U2
2
= Vs 0 u2 = Vs 0 aV
s
であり、衝撃波後面の密度は
2
= 1
Vs
a
;
2
(A.7)
(A.8)
となる。
A.2
等温の場合の Riemann 不変量
連続の式
を で割ったもの
@
@t
@
@u
@
+ @x
+ u @x
= 0;
ln + @u + u @ ln = 0;
@t
@x
@x
(A.9)
(A.10)
46
付録
と、運動方程式
@u
@t
を等温音速 a で割ったもの、
から、式
@M
@t
2
衝撃波管の解析解
@
+ u @u
+ a @x
= 0;
@x
(A.11)
M + a @ ln = 0;
+ u @@x
@x
(A.12)
(A.10) と式 (A.12) を加えると、
@ ln + M
+ (u + a) @ ln + M = 0;
@t
また、式
A
@x
(A.13)
(A.10) から式 (A.12) を引くと、
@ ln 0 M
+ (u 0 a) @ ln 0 M = 0;
(A.14)
が得られる。すなわち、C+ カーブ dx=dt = u + a では、J+ = ln + M が一定で、C0 カーブ dx=dt = u0 a
では、J0 = ln 0 M が一定となる。この J6 が等温の場合の Rieman 不変量である。
@t
A.3
@x
等温の場合の衝撃波管の解析解
マイナス方向に希薄波が伝搬する。その上では右向きに伝わる C+ に沿って、J+ が一定である。C+
は、もっとも左側の静止している部分につながっていることに注意すると、
J+
= ln + M = ln H ;
(A.15)
ここで H はもっとも左側の静止している部分の密度である。これから、
H
= exp(0M);
(A.16)
が成り立つ。希薄波のもっとも右側の密度と、速度が衝撃波後面のそれ( U2 、2 )に等しいことを用い
ると、
2
H
(A.7 A.8)
= exp 0 Ua2
(A.17)
;
等温衝撃波の条件式
、
から、もっとも右側の静止している衝撃波前面の密度 L と、もっとも
左側の静止している希薄の先の密度 H と、衝撃波の速度の関係が以下のように得られる。
L
これは Vs =a
= とおくと、
Vs
a
2
H
= exp 0
Vs
a
0 Va
s
;
(A.18)
(A.19)
exp 0 1 = H ;
L
となる。この解は数値的に求めると、 H = 10 のときに、 = 1:75194::: で 0 x1 = 1:1811::: のように
L
求められる。
2
A.3.
47
等温の場合の衝撃波管の解析解
すなわち、初期の密度の不連続の位置を原点にとり、a
8
>
>
>
>
<
=>
>
>
>
:
U
L
L 2
H
H
exp(0M)
8
>
>
>
>
<
=>
>
>
>
:
0
0 1
x+t
0
t
となる。これを数値計算の結果と比較する。
= 1 とすると、
x > t
0 1 0 t < x < t
0t < x < 0 1 0 t
x < 0t
(
1)
(
1)
x > t
0 1 0 t < x < t
0t < x < 0 1 0 t
x < 0t
(
1)
(
1)
(A.20)
(A.21)
49
付 録B
計算のこと
B.1 Makele
1
2
一連のプログラムの翻訳(コンパイル)、結合、実行などの手続きは
しておく。そして
Makele というファイルに記述
make ターゲット
make
という
コマンド でターゲットを実現する(ファイルを作ったり、なにかを実行したりする)
。こ
れによって、サブルーティンごと( この例では、プログラムはサブルーティンごと
、
、
、
、
に格納されている)の分割コンパイル時に必要なもののみが再翻訳される。
main.F initia.F
FDS.F FDS2.F Muscle.F
B.1.1
1
2
3
1
2
基本的な Makele を書く方法
基本は、
ターゲット : 元になるもの 1 元になるもの 2 元になるもの 3
<-- タブ -->元になるものからターゲットを作成する方法
という形であり、以下の例では、
main.o: main.F var.h
main.F
1
2
Fortran
#include
var.h
C
CPP
var.h
(CPP)
main.F
<-- タブ -->$(F77) $(CFLAGS) -c -o main.o main.F
main.o
$(F77) $(CFLAGS)
の部分で、
を作成する(翻訳する)コマンド が記述されている。
と
は変
数を引用しており、上の方で定義されている。
ここでは、主プログラム以外のサブルーティンは翻訳した後
というファイルにライブラリ
としてアーカイブしている。
実行ファイルを作成する時には
libFDS.a
1
2
.F
Fortran
.F
の部分で、
をコンパイルする時に(普通の
コンパイラーにとって、拡張子 をもつ
に通されてから
ファイルに格納されたソースプログラムはまず のプ リプロセッサー
コンパイラーに掛けられる。
のような
の制御文が書かれていることを拡張子 は示して
いる。)、
が必要であること(依存関係があり
が変更された場合は、
の再コンパイル
が必要なこと)が示されている。
FDS: main.o
libFDS.a
50
3
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
付録
B
計算のこと
<-- タブ -->$(F77) $(LFLAGS) -o FDS main.o -L. -lFDS
main.o
FDS
-lFDS 引数で libFDS.a から検
で
を入力し、
という実行ファイルを作成する。そのとき、
索して結合することを指定している。
Makefile の例
#!/bin/sh
#
F77 = f77
CFLAGS = -u
LFLAGS =
#
all: FDS
./FDS
main.o: main.F var.h
$(F77) $(CFLAGS) -c -o main.o main.F
#
initia.o: initia.F param.h dim.h var.h
$(F77) $(CFLAGS) -c -o initia.o initia.F
#
FDS.o: FDS.F param.h dim.h var.h
$(F77) $(CFLAGS) -c -o FDS.o FDS.F
#
FDS2.o: FDS2.F param.h dim.h var.h
$(F77) $(CFLAGS) -c -o FDS2.o FDS2.F
#
Muscle.o: Muscle.F param.h
$(F77) $(CFLAGS) -c -o Muscle.o Muscle.F
#
output.o: output.F param.h dim.h var.h
$(F77) $(CFLAGS) -c -o output.o output.F
#
libFDS.a: initia.o output.o FDS.o FDS2.o Muscle.o
ar ru libFDS.a initia.o output.o FDS.o FDS2.o Muscle.o
#
FDS: main.o libFDS.a
$(F77) $(LFLAGS) -o FDS main.o -L. -lFDS
B.2. Gnuplot
51
B.2 Gnuplot
1次元の計算結果を簡単にグラフにする方法のひとつとして
B.2.1
1
2
3
4
5
Gnuplot をあげておく。
データの作成法
データは
do j=0,Nx
write(*,'(2(3x,f10.3))') j,rho(j),u(j)
end do
write(*,*)
のようなプログラムで、1行に配列要素番号、そこの密度、そこの速度がこの順に並び、次の点の
データが次の行に並ぶように作成する。一連のデータが書き終ったら、空の行を1行つけておく。
さらに時間が経った後のデータもこの後に同じようにつけておく。これをファイルに保存しておく。
B.2.2
1
2
3
4
5
6
7
8
1
2
1
2
3
4
5
6
1
プロット 法
まず、プロンプトから
gnuplot コマンド を実行するそうすると、
G N U P L O T
Unix version 3.7
patchlevel 1 (+1.1.9 1999/11/08)
last modified Fri Oct 22 18:00:00 BST 1999
途中省略
Terminal type set to 'x11'
gnuplot>
gnuplo> のプロンプトが出る。データが upwind.d に格納されている場合、
となり、
gnuplot> plot "upwind.d" with linespoints
とすると、別のウインド ウが開き線画が画面に得られる。さらに、印刷するには、
gnuplot> set terminal postscript
Terminal type set to 'postscript'
Options are 'landscape noenhanced monochrome dashed defaultplex "Helvetica-Ryumin" 14'
gnuplot> set output "upwind.ps"
gnuplot> replot
upwind.ps にポストスクリプト形式で図が書き込まれる。
のようにすれば、
52
2
1
2
付録
B
計算のこと
gnuplot> !lpr -Plpn1 upwind.ps
もしくは
% lpr -Plpn1 upwind.ps
のようにすれば、プ リンターに出力することができる。
参考資料
数値流体力学に関する教科書
(1994)
流体力学の数値計算法★★★
東京大学出版会、藤井孝蔵著
はじめての
移流拡散方程式
★★
コロナ社、棚橋降彦
東京大学出版会、保原 充、大久保久明編
数値流体力学 基礎と応用★★
圧縮性流体解析★★
東京大学出版会、数値流体力学編集委員会編
発展方程式の数値解析★★
岩波応用数学叢書、矢嶋信男 野木達夫著
東京大学出版会、荒川忠一著
数値流体工学★
CFD |
(1995)
(1994)
(1992)
|
(1996)
(1977)
,
重要な論文
Roe, P. L.: Approximate Riemann Solvers, Parameter Vectors, and Dierence Schemes, (1985) JCP,
43, 357.
Roe, P. L.: A Survey of Upwind Dierencing Techniques (1989) Lecture Notes in Physics, 323, 69.
van Leer, B.: Flux-Vector Splitting for Euler Equations (1982) Lecture Notes in Physics, 170, 507.
流体力学(圧縮性流体)に関する教科書
(1972)
(1970)
)
(1997)
流体力学★★★
培風館、巽友正著
流体力学★★
東京図書、エリ・ランダウ、イェ・リフシッツ著 竹内均訳 ランダウ リフシッ
ツ理論物理学教程
宇宙流体力学★★
培風館、坂下志郎、池内 了著
;
(
=
数値流体力学に関するホームページ
http://nova.earth.s.kobe-u.ac.jp/shocktube/index-j.html: Java による 1 次元ショックチューブ問題の
アニメーション: 藤原秀和 (神戸大)
http://grape.c.u-tokyo.ac.jp/hachisu/java.shtml: 1-D Numerical simulation by Java: 蜂巣 泉(東大
教養学部)
http://www.camk.edu.pl/tomek/htmls/num meth.html: Numerical methods for Fortran Programmers: Tomasz Plewa ( Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences
)
B.2. Gnuplot
53
http://grape.c.u-tokyo.ac.jp/hachisu/lecture/hydro/hydro.shtml: 数値流体力学の基礎と宇宙気体力
学 (大学院講義ノート ): 蜂巣 泉(東大教養学部)
ダウンロード

Document