第 17 章 圧縮性流れの解法|Navier-Stokes 方
程式
前章では，2 独立変数の双曲型偏微分方程式の特性の理論を一通り説明した後に，この理論に基づく１次元
オイラー方程式の初期値問題の３つの解法，すなわち流束ベクトル分離法，流束差分分離法，有限体積法に
ついて詳しく述べた．本章では，流束差分分離法と有限体積法を多次元に，また拡散項を含め Navier-Stokes
方程式の解法に拡張する．この拡張に伴って方程式系の式の数と項数は増えるが双曲型の性質は変わらな
いので，この拡張は困難なく実行できる．ここでは３次元のデカルト座標系と一般曲線座標系に対し，また
一般座標系に対しては反変速度または質量流束の運動方程式についてこれらの解法の説明を行う．最後に計
算の高効率化についても述べる．
17.1
3 独立変数以上の双曲型微分方程式
はじめに 非定常２次元または３次元流れの Euler 方程式にも適用できる
3 独立変数以上の双曲型偏微分
方程式の理論について述べる．この理論に関する定番ともいえる解説書は Courant の著書である1 ．しかし
ながら この本を読破するには数学的素養と相当の時間が必要である．以下にはこの本の中から必要最低限
の事項を選びあえて不完全ともいえる図を補い平易に解説する．この解説は入門書で，諸式の導出と解釈
についてはやや詳しく述べるが，定理の証明はほとんど 省略されているので関心のある読者は
Courant の
原著を参照されたい．
双曲型偏微分方程式の理論は，この方程式の初期値問題の解がいかに構成されるのかを数学的に論じ，ま
た双曲型方程式で支配される非定常２次元または３次元空間の波動現象，すなわち特性面に沿って伝播す
る波の挙動を明らかにするものである．この多次元の特性の理論に基づいて作られた解法は物理現象に則
したもので適正な解を与えるものと思われる．しかしながら現時点での解法の主流は，各座標軸方向に１
次元流れの特性の理論に基づく流束分離法を適用して多次元問題に対処しようとするもので，本章の内容
もこの範囲を出るものではない．このような解法は衝撃波の捕獲に関し問題があるとの指摘もあり，このこ
とを可能ならば本章の終わりに検証したいと考えている．
というわけで，本章に述べる二つの解法を理解するには本節を読み飛ばしても特に支障はない．とはい
え，この分野の研究者には本節の内容程度のことは理解されることを期待する一方，大学の数学の講義で
も，コンピュータの時代にあわせて不要不急のものは切捨て，
「 微分方程式論」で式の意味や解の構成，ま
た「数値解法」で解法の選び方や解き方を教授されることを望むものである．
1
R. Courant, \Methods of Mathematical Physics, Vol.II Partial Di erential Equations," 1962, Interscience Publishers.
日本語訳，R. クーラン・D. ヒルベルト，数理物理学の方法，4, 麻嶋訳，東京図書．
1
2
圧縮性流れの解法|NS 方程式 17.1.1
3 独立変数以上の双曲型方程式の Cauchy 問題
3 独立変数以上の 1 階準線形偏微分方程式系 ( rst-order quasi-linear partial di erential equations in more
than two independent variables) は一般に次のように書くことができる．
L u]
ただし
n
X
i=0
Ai Di u = f
(17.1a)
u(x) は k 成分の未知変数ベクトル，Ai は k k 成分の係数行列で，
0 1
0i i
1
0 1
u1
a11 a12 : : : ai1k
f1
B
C
B
C
B
B
BBai21 ai22 : : : ai2k CC
BBf2 CCC
u2 C
B
C
i
u=B
f =B
. C A =B
.. C
B
B@ ... ...
B@ ... CCA
. C
@ .. CA
A
aik1 aik2 : : : aikk
uk
(17.1b)
fk
また Ai と f は一般に x0
x1 : : : xn と u の関数である．時間依存問題では，x0 = t は時間，また x1 : : : xn
は n 次元空間座標である．Di @=@xi ．A0 は単位行列 I とすることができる．式 (17.1) は，1 階微分 Di u
の 1 次式で，1 階準線形偏微分方程式系である．
;1 次元，線は 1 次元，2 次元から m;2 次元の部分空間は多様体 (manifold)
一般に m 次元空間内の面は m
と呼ばれる．当然のことながら，k 次元の多様体の中には
k 次元の座標系を取ることができる．m 次元空
間の座標 x1 : : : xm に１つの拘束条件をかせば面を作ることができる．例えば (x) = 0 は一つの面の式
である． (x) = c と置き，パラメータ c の値を変えれば ，1 パラメータ族の面を発生させることができる．
したがって は面 (x) = 0 のそとの座標になる．上記の n +1 次元の x 空間内に一つの面 (x) = 0 を取
) である．この n 次元の面の中に座標 1 : : : n ，外に 0 = を取る．今独
る．この面の法線は r (
立変数を x から に変えれば ，
n @
X
Di u = i u + @xj u j
j =1 i
ただし
@ =@xi = Di
i
n
X
i=0
である．この関係を用いれば式 (17.1) は次のように書換えられる．
Ai i u = g
(17.2a)
あるいは簡単に
Au = g
ただし
A=
また g は
n
X
i=0
(17.2b)
Ai i
u u
1
(17.2c)
::: u
n
の関数である．
式 (17.1) がどのような条件を満足するときに双曲型になるのかという議論は先送りにし，ここでは式 (17.1)
が双曲型であるとしてその Cauchy 問題 (初期値問題) を考える．この問題は基本的には与えられた初期曲
3 独立変数以上の双曲型微分方程式 3
= 0 上に初期値 u が与えられるときに，この初期値を満足する式 (17.1) の解 u(x) を初期曲面の近傍で
面
求める問題である．まず初期曲面上のいたるところで
Q( ) jAj 6= 0
の場合には，初期曲面上に u の値が与えられれば ，u
: : : u n は既知で，したがって，A と g も既知に
なる (下図参照)．式 (17.2) は u の連立 1 次方程式で，jAj 6= 0 であるから u の値は一意的に決定できる．
このようにして初期曲面上で関数値 u とその法線微分値 u が分かれば ，初期曲面 = 0 のご く近傍の面
= c 上で u の値が分かることになり， = c を初期曲面と見て同様のことを繰返せば ，解 u(x) を初期曲
面 = 0 の近傍に延長することができる．この場合には初期値 u を任意に与えることができ，解 u(x) を一
1
意的に求めることができる．このような初期曲面を無拘束 (free) であるという．
他方 初期曲面上のいたるところで特性条件
jAj = Q( ) = Q( 0
: : : n) = 0
(17.3)
が満足される場合にはこの面は L の特性面と呼ばれる．この場合には u の連立 1 次方程式 (17.2) を構成
する k 個の式は 1 次独立でなくなる．式 (17.2) の g はこの式が正則になるように与えなければならず，初
期値 u は任意にではなく g がそのようになるように与えなければならない2 ．この場合にはまた，式 (17.2)
から u の値は一意的には定まらず，初期曲面上の一つの初期値 u から無数の解を延長することができる．
17.1.2
1
法線錐，射線錐，特性面の構成
時間 t(= x0 ) を別扱いにして特性面を次のように表すことにする (上図参照)．
(x t) = (x) ; t = 0
(17.4)
(x) = const: = t で x 空間内を伝播する．この面の法線方向距離を n とすれば，特性面の
法線速度は v = dn=dt，また の勾配は = r = d =dn = dt=dn となる．ベクトル と v は同じ方向で
この特性面は，
n
X
i=1
i vi =
v=1
= v=v2 ，v = = 2 である．
式 (17.4) から 0 = @ =@t = ;1 i = @ =@xi = @ =@xi (i = 1 2 : : : n)，したがってこの場合の特性の
条件は式 (17.3) から
また
Q(;1
1
: : : n) = 0
2 この意味が納得できないときには
( i = Di )
16.1 節を参照してください．
(17.5a)
4
圧縮性流れの解法|NS 方程式 図
17.1:
法線錐と法線面
あるいは成分で書けば
a111 1 + + an11 n ; 1 a112 1 + + an12 n
a121 1 + + an21 n a122 1 + + an22 n ; 1
..
.
..
.
...
a11k 1 + + an1k n
a12k 1 + + an2k n
=0
..
.
a1kk 1 + + ankk n ; 1
a1k1 1 + + ank1 n
a1k2 1 + + ank2 n
この式は i の k 次の同次方程式である．特性面 = t の単位法線ベクトルを
式 (17.5a) から特性面の速度 v の方程式は次のように導かれる．
Q(;jvj
1
:::
とすれば ，
i = i =j j，
n) = 0
(17.5b)
次に，特性条件の式 (17.3) の幾何学的解釈について述べる．この式は，n +1 次元の xt 空間内のある点
O を通る特性面の点 O における法線 (= r ) を拘束する式である．例えば ３次元空間 xyt では，ある点
O における は，あらゆる方向 (2 パラメータ族の法線) を取り得るのではなく，Q( 0 1 2 ) = 0 により
１つ次元が下がり点 O に頂点を持つ一つの錐面 (１パラメータ族の法線からなる) 内に制限される．一般に
式 (17.3) はこのような n ; 1 パラメータ族の法線からなる錐面の式で，この錐面は法線錐 (normal cone) と
呼ばれる (図 17.1)．また時間 t を特別扱いした式 (17.5a) はこの法線錐の 0 = ;1 面による切り口を表し ，
これは 1 : : : n 空間内の法線面 (normal surface) と呼ばれる．上述のように，特性面の法線速度 v は
と同じ 方向を持ち， v = 1 であるから，法線面を単位球面 2 = 1 に関して反転し ，v1 : : : vn 空間内
に法線速度面 (normal velocity surface) を作ることができる．この面はまた相反法線面 (reciprocal normal
surface) とも呼ばれる．
特性方程式 Q( ) = 0 は ( ) の 3 独立変数以上の 1 階非線形微分方程式である．この項の議論を先に進
めるには 1 階非線形微分方程式の理論を用いる必要がある．以下には読者がこの理論を知らないものとし
てまず簡単な 2 独立変数の 1 階準線形微分方程式，次いでこれとの比較において 3 独立変数以上の 1 階非
線形微分方程式の理論について略述する．
2 独立変数 1 階準線形微分方程式は次のように書くことができる．
aux + buy = c
ただし
u(x y) はスカラーの未知変数，a b c は一般に x y u の関数で，a2 + b2 6= 0．式 (17.6) は，
(a b c) (ux uy ; 1) = 0
(17.6)
3 独立変数以上の双曲型微分方程式 (a)
図
6= 0
(b)
5
=0
17.2: 2 独立変数 1 階準線形微分方程式の初期値問題
のように書換えることができるので，xyu 空間内のベクトル (a
b c) と解曲面 (integral surface) u = u(x y)
の法線ベクトル (ux uy ; 1) が直交することを示している3 ．ベクトル (a b c) を Monge 軸，またいたる
ところで接線が
Monge 軸の方向をとる曲線を特性曲線 (characteristic curve) と呼ぶ4 ．特性曲線に沿って
パラメータ s を適切に取れば，特性曲線は次の常微分方程式で与えられる．
dx=ds = a dy=ds = b du=ds = c
(17.7)
特性曲線は図 17.2(a) に示すように一つの解曲面上に乗っている．
微分方程式 (17.6) は初期値問題として解かれる．初期曲線上に初期値 u を与え曲線 C を
x = x(t) y =
y(t) u = u(t) xt + yt 6= 0 で定義する．式 (17.6) の初期値問題は C 上の各点を通る式 (17.7) の特性曲線族
2
2
x = x(s t) y = y(s t) u = u(s t)
の作る解曲面 u(x
y) を C の近傍で求める問題である (C 上で s = 0)．この面が形成されるためには，上式
の第 1 式と第 2 式を s または t で微分したものが独立でなければならない．すなわち
xs ys
a b
=
xt yt xt yt
で定義される
がゼロでないという条件が満足されなければならない5 ．まず初期曲線が特性でない
の場合には，初期値 u(t) は任意に与えることができ，また解 u = u(x
6= 0
y) を図 17.2(a) に示すようにこの初
期値を通る 1 パラメータ族の特性曲線を引くことによって一意的に求めることができる．特性の理論では
このように偏微分方程式が常微分方程式系に置換えられ，特性の方法ではこの常微分方程式系が初期値問
題として数値的に解かれる．
= 0 の場合には，初期値 u(t) が，式 (17.6) すなわち dx : dy : du =
a : b : c なる条件を満足するように与えられるときにのみ初期値問題は成立し ，図 17.2(b) に示すようにこ
他方 初期曲線が一つの特性曲線
の初期値から無数の解を延長することができる．
F (x y u) u(x y) ; u = 0
Monge
C
s t
rF = (ux uy ; 1) となる．
xs yt ; xtys = 0 となる．
3 念のため，解曲面は
のように表されその法線ベクトルは
4
軸は
空間内の方向場で，特性曲線は
軸の包絡線であるともいえる．
5 逆に 曲線 が特性であれば， 上で と の方向は同じで，
，したがって
Monge
C
xyu
xs : xt = ys : yt
6
圧縮性流れの解法|NS 方程式 図
17.3:
法線錐と Monge 錐
次に 3 独立変数以上の 1 階非線形微分方程式を考える．
F (x1 : : : xn u p1 : : : pn ) = 0
ただし
pi = @u=@xi
結果 ２次元の面は
(17.8)
P F 2 =6 0 である．ここでは上記の理論の独立変数 x y は x : : : x となり，その
pi
1
n
n 次元の面に，線は線または n ; 1 次元の多様体になる．また準線形から非線形への一
般化に伴い Monge 軸は Monge 錐に変身する．式 (17.8) は，式 (17.3) と同様に幾何学的に解釈することが
できる．すなわち n +1 次元の xu 空間内で，一つの点 O を通るすべての可能な解曲面 u = u(x) の点 O に
おける法線ベクトル (p1
pn ;1) は点 O に頂点を持つひとつの錐面になる．この錐面は法線錐と呼ば
れる．敷衍すればこの法線ベクトルは，あらゆる方向を取りうるのではなく，式 (17.8) の条件により次元
が一つ下がり n = 2 の場合には 1 パラメータ族の法線からなる法線錐，n = n の場合には n
; 1 パラメータ
族の法線からなる法線錐内に制約されるということである．
点 O を通るすべての可能な解曲面の点 O における接平面の包絡面も点 O に頂点を持つ錐面になる．この
錐面は Monge 錐または特性錐 (Monge
cone, characteristic cone) と呼ばれる．特性曲線はいたるところで
その接線が特性錐の母線の方向を取る曲線である．次に特性曲線上で成立する式を導く．ここでは，図 17.3
に示すような幾何学的にイメージし易い n = 2 の場合から始める6 ．Monge 錐の母線 (generating line) に
沿ってその頂点からの距離を とすれば，母線方向のベクトル (dx=d dy=d du=d ) は解曲面の法線ベク
トル (p q ; 1) と直交するから次式が得られる．
p( ) ddx + q( ) ddy = ddu
は法線錐上に取られたパラメータである．ベクトル (p
ただし
接平面 (tangent
q ;1) は法線 (normal) と同時に解曲面の
plain) を表すものでもある．Monge 錐は接平面の包絡面であるからその母線上では次式が
成立する7 ．
p ( ) ddx + q ( ) ddy = 0
0
0
他方，この場合の式 (17.8) を
で微分すれば
Fp p ( ) + Fq q ( ) = 0
0
0
u
Monge 錐を先に与えその法線錐を計算で求め示したものである．法線錐は Monge
錐と頂点を共有する 1 対の楕円錐で， = 1 面による切口は楕円または双曲線になる．
7 Monge 錐が接平面の包絡面ということは， に対応する Monge 錐の母線が と +
の二つの接平面の交線の，
0の
ものということである．したがって包絡面の母線の式は交線が によらないものとして上式を で微分することによって得られる．
6 この図は， 一定面による切口が円形になる
u ;
!
3 独立変数以上の双曲型微分方程式 以上の 3 式から
7
Monge 錐の母線の式が次のように導かれる．第 2，3 式から dx : dy = Fp : Fq ，この関係
を第 1 式に用いれば
dx : dy : du = Fp : Fq : pFp + qFq
特性曲線はいたるところでその接線が特性錐の母線の方向を取る曲線であるから，特性曲線の式は
dx = F dy = F du = pF + qF
(17.9)
p
q
ds p ds q ds
ただし s は特性曲線に沿って適切に取られたパラメータである．式 (17.9) の最後の式はこの曲線に面を付
加するもので成帯条件 (strip
condition) といわれれる．
この場合の式 (17.8) を x で微分すれば
Fp px + Fq qx + Fu ux + Fx = 0
式 (17.9) の第 1，2 式と qx
= py を用いれば
dy + F p + F = 0
px dx
+
p
y
ds ds u x
ひとつの積分曲面 u
= u(x y) 上では p q は x y の与えられた関数で，dp=ds = pxdx=ds + py dy=ds である
から
dp + F p + F = 0
ds u x
となる．同様に q の式も導かれる．上記の常微分方程式はまとめて式 (17.8) の特性微分方程式系 (characteristic
system of di erential quations) と呼ばれる．
dx = F dy = F du = pF + qF
p
q
ds p ds q ds
(17.10)
dq = ;;qF + F
dp = ;;pF + F
u x
u y
ds
ds
これら５つの式は５つの未知変数 xyupq を含み決定系をなす．式 (17.10) の解は一般に xyupq 空間内の 4
パラメータ族の曲線になる．１つの曲線に沿って F を s で微分し 式 (17.10) を用いれば dF=ds = 0 なる関
係が得られ，これより１つの曲線に沿って F = const. となる．したがって微分方程式 F = 0 を満足する曲
線は 3 パラメータ族存在することになる8 ．この曲線を特性帯 (characteritic
strip) という．
３独立変数以上の場合にも同様に理論を展開することができ，式 (17.8) の特性微分方程式系は下記のよ
うに n = 2 の場合の式 (17.10) を直接的に多次元に拡張したものになる．
n
dxi = F
du = X
p
i
ds
ds i=1 pi Fpi
dpi = ;;p F + F
i u xi
ds
(i = 1 2 : : : n)
xyu
(17.11)
8 見方をかえれば，2 独立変数準線形の場合には
空間内に Monge 軸の焦曲線である特性曲線は 2 パラメータ族存在した．つ
まり一つの特性曲線をとり，１つのパラメータを変えればこの曲線が面をなし，それから残りのパラメータを変えれば特性曲線が
空間を覆った．非線形の場合には第 3 のパラメータが存在し これを変えれば特性曲線が
空間内の各点でこの点を頂点に持つ一
対の錐曲線体 (conoid) を描くことになり，第 1，2 のパラメータを変えればこの錐曲線体が空間を覆うことになる．
xyu
xyu
8
圧縮性流れの解法|NS 方程式 (a) 6= 0
(b) = 0
図 17.4: 2 独立変数 1 階非線形微分方程式の初期値問題
2 式は，特性曲線に接平面を付加し 特性帯を作る成帯条件 (strip condition) である．特性帯は式
(17.11) の条件 F = 0 を満足する解で，特性帯は一つの解曲面 u = u(xi ) に乗っている．
式 (17.8) は初期値問題として解かれる．n ; 1 次元の初期多様体 C 上にパラメータ t1 : : : tn 1 を取り，
C を ti の関数 x1 : : : xn u で与える．更に C 上に微分方程式 (17.8) と成帯条件
なお第
;
n @x
@u = X
i
( = 1 2 : : : n ; 1)
(17.12)
@t i=1 pi @t
を満足するように ti の関数 p1 : : : pn を与える．このとき C は帯多様体 (strip manifold) C1 になり，積分
多様体 (integral manifold) と呼ばれる．この場合の初期値問題は C1 を含む式 (17.8) の解 u = u(xi ) を見出
す問題になる．C1 に沿って
@xn
@x1
@s
@s
@x1
@xn
@t1
@t1
:::::::::::::::::::::::
@x1
@xn
@tn 1
@tn 1
;
で定義される
=
;
Fp1
Fpn
@x1
@xn
@t1
@t1
:::::::::::::::::::::::
@xn
@x1
@tn 1
@tn 1
;
(17.13)
;
がゼロでない場合には，特性微分方程式系 (17.11) の初期値問題を解くことによって，積
分多様体 C1 の近傍に積分曲面
xi (s t ) u(s t ) pi (s t )
s は特性曲線に沿って取られるパラメータで C1 上で 0 である．
このようにして求めた解がもとの微分方程式 (17.8) を満足し，また pi が @u=@xi になることは容易に証明で
きる．この場合の初期多様体，解曲面と特性錐の関係は，３次元の xyu 空間では，図 17.4(a) に示すように
を一意的に延長することができる．ただし
なり，初期積分帯 (integral
strip) は特性錐に接するように与えられ，また解曲面も特性錐に接するように
求められることになる．
次に微分方程式 (17.8) と成帯条件 (17.12) を満足する初期多様体上のいたるところで特性条件
成立する場合を考える．
nX1
;
=1
@xi = F
pi
@t
= 0 ということは線形の関係
(i = 1 2 : : : n)
= 0が
3 独立変数以上の双曲型微分方程式 が成立するということである．ただし
nX1
;
=1
nX1
;
=1
n
@u = X
@t i=1 pi Fpi
@pi
@t = ;(pi Fu + Fxi )
9
(t1 : : : tn 1 ) は係数である．初期値 u と pi が
;
(i = 1 2 : : : n)
なる関係を満足するならば ，この初期多様体は n
れる．初期多様体のいたるところで
;2 パラメータ族の特性帯からなり，特性帯多様体と呼ば
= 0 の場合に，初期値問題が成立するための必要十分条件は，初期
多様体が特性帯多様体であることで，このとき初期値問題の解は無数に存在することになる．３次元の xyu
空間では，初期特性帯，解曲面，特性錐の関係は，図 17.4(b) のようになる．
以上述べた 2 独立変数準線形方程式の理論と
3 独立変数以上の非線形方程式の理論の間には次の対応関
係がある．
2 独立変数準線形方程式 (17.6)
２次元法線面
Monge 軸
初期曲線
特性曲線 式 (17.7)
n 独立変数非線形方程式 (17.8)
n ; 1 パラメータを持つ法線錐
n ; 1 パラメータを持つ特性錐
n ; 1 次元初期多様体
特性帯 式 (17.11)
(x0 x1 : : : xn ) の 1
= const: = c はもとの微分方程式系 (17.1) の特性面である．n +1 次元
ここで，3 独立変数以上の双曲型微分方程式の説明に戻る．特性方程式 (17.3) は
階非線形微分方程式で，その解
の空間内で，式 (17.3) の解は n
;1 パラメータ族の特性曲線によって作られる．この特性曲線は，演算子 L
に対して従特性 (bicharacteristic) または射線 (ray) と呼ばれる．従特性に
i = Di の付加された従特性帯
(bicharacteristic strip) は，上記の 1 階非線形微分方程式の特性微分方程式系 (17.11) から導かれる次の常
微分方程式系を満足する．式 (17.3) には u が陽に現れないから
x_ i = Q i _i = ;Qxi
(i = 0 1 : : : n)
(17.14)
` _ ' は射線に沿って適切に取られたパラメータに関する微分を意味する．
xt 空間内の一つの点 O を頂点に持つ射線錐 (ray cone) を上記により作る．この射線錐は法線錐と双対的
関係にある (図 17.5(a))．ここで再び x0 = t を別扱する．前方の射線錐の t = 1 面での切り口は射線面 (ray
surface) と呼ばれる．この射線面は，点 O の擾乱から t = 0 に出た波の t = 1 における波面を表している．
ただし
このような波の一部分を平面波
(x t) = vt ; 1x1 ;
で近似する．ただし
とき，ベクトル
v
; nxn = vt ; x = 0
v はこの波の伝播の速さ，
は波面の単位法線ベクトルである． が単位球面上を動く
は法線速度面を作り，また平面波 v
=
x の包絡面は射線面を作る．射線面はまた法
線面の接平面の単位球面に対する極の作る面でもある．２次元の場合の射線面，法線面，法線速度面の関係
を図 17.5(b) に示す．
次に 2 次元の場合の具体例を２つ示す．第 1 の例は射線面がひとつの点 (a
面はひとつの円
xv = a(1 + cos )=2 yv = a(sin )=2 (0
0) の場合である．法線速度
2 ) になり，法線面は y 軸に平行な線
10
圧縮性流れの解法|NS 方程式 r:v= :x v =1
(a)
図
(b)
17.5:
法線面，射線面，法線速度面
x = 1=a になる．図中の ABCD は対応する点を示す．第 2 の例は射線面が円でその一部が単位円の外に出
yr = r sin (0
2 )．法線速度面は xv = (r + a cos ) cos yv =
る場合である．xr = a + r cos
(r + a cos ) sin ，また法線面は xyt 空間内の法線錐 (一対のだ円錐) を t = ;1 平面で切断したときに現れる
= sin =(r + a cos ) になる9 ．擾乱は一般に xt 空間内の特性面に沿って伝
双曲線 = cos =(r + a cos )
播する．今３次元 x 空間を考えれば，ある瞬間の波面は２次元の面で，単位時間後の波面はその瞬間の波面
上の各点で射線面を作り，これらの射線面の包絡面として求めることができる (図 17.6，波面の Huyghens
構成) 10 ．
9 この双曲線の式を通常ように
(x ; C )2 ; y 2 = 1
A2
B2
で表せば係数 A B C と漸近線の式は
A = a2 ;r r2 B = p 21 2 C = a2 ;a r2
a ;r
y=
p
a2 ; r2 (x ; C )
r
Huyghens 構成はこれによる窪みを接平面で補正した凸殻錐 (convex hull cone) の
10 射線錐の一部が凹面になっている場合には，
射線面に対して行われる．
3 独立変数以上の双曲型微分方程式 図
17.1.3
17.6:
11
波面の Huyghens 構成
演算子 lL
微分方程式 (17.1) の演算子 L は，式 (17.2c) の行列 A が対称のときには，A が正定値 (positive
de nite)
ならば双曲型である．ここでは L が双曲型であるための条件を幾何学的に法線錐を使って説明する．図 17.7
に示すような，xt 空間内の点 O に頂点を持つ法線錐 Q( ) = 0 と点 O を通るベクトル を考える．この
を含むすべての２次元平面
が法線錐と
存在すれば 式 (17.1) の演算子 L は点 O
k 個の異なる線で交わるものとする．このようなベクトル
で双曲型である11 ．このようなベクトルを擬似時間方向
が
(time-like
direction) という．この双曲型の定義を代数的に定義し直せば次のようになる．点 O を通り法線錐に交叉
しない n 次元の面を擬似空間面 (sapce-like surface) という．この面内にベクトル を取る．L が双曲型の
ときには と平行な線 = + ( はパラメータ) は法線錐と一般に k 個の異なる点で交叉する (図参照)．
したがって
の式
Q( + ) = 0
が
k 個の異なる実根
を持つときに L は双曲型である (脚注参照)．
図
17.7:
双曲型微分方程式の幾何学的定義
11 法線錐が図のように交叉しているところ，接しているところ，完全に重なっている場合には交線の数はその分だけ少なくなる．
12
圧縮性流れの解法|NS 方程式 式 (17.3) の成立する特性面上では，行列 A に対し
lA = 0 Ar = 0
(17.15)
を満足する左ナルベクトル l と右ナルベクトル r が存在する12 ．演算子 lL は式 (17.2c) を参照すれば
n
n
n
X
X
X
(17.16)
Ai i @@ + @@xj @@ = l Ai @@xj @@
j =1 i j
i=0 j =1 i j
i=0
i=0
のように書換えることができる． 1 : : :
n は特性面内の座標であるから，演算子 lL は特性面の内微分で
あることが分かる．微分方程式 L u] = 0 に l を左から演算すれば ，特性面の内微分のみを含む見通しの良
lL = l
n
X
Ai Di = l
n
X
い方程式が得られる．
本節の始めに述べたように n +1 次元空間の特性条件は
A = I (I は単位行列) で，式 (17.1) が
Q( 0
1
0
L u] = ut +
n
X
i=1
: : : n ) = 0，また x0 = t
Ai Di u = f
;1
0 である．ただしここでの
1
= ;1
(17.17)
のように書かれる n 次元空間では，特性条件は Q(
=j j に対しては Q(;jvj
0
: : : n ) = 0 である．更に単位法線ベクトル =
n ) = 0，また法線速度 v = = 2 に対しては Q(;1 v1 =v 2 : : : vn =v 2 ) =
:::
は n 次元空間のもので
2
1
=
1
2
+ + n 2 ．単位法線ベクトル
で表された特
性条件は具体的に書けば次のようになる．
A~ ; jvjI =
a111 1 + + an11 n ;jvj a112 1 + + an12 n : : : : : :
a121 1 + + an21 n a122 1 + + an22 n ;jvj : : : : : :
..
..
..
a1k1 1 + + ank1 n
ただし
..
..
..
..
=0
...
.
1
: : : : : : akk 1 + + ankk n ;jvj
...
a1k2 1 + + ank2 n
a11k 1 + + an1k n
a12k 1 + + an2k n
...
(17.18)
P
A~ = ni=1 Ai i である．したがって jvj は行列 A~ の固有値，またナルベクトル l と r は A~ の固有ベ
クトルである．
xt 空間内の初期面 t = 0 内に n;1 次元の初期多様体を取り，この多様体を含む n 次元の特性面を初期面
の近傍に延長する (下図参照)．ここでは初期多様体上の擾乱がこれらの特性上を波となって伝播する問題を
(x) = 0 とすれば ，その法線ベクトルは = r ，単位法線ベクトルは i = i =j j
~ が式 (17.18) より，k 個の固有値 jv1 j : : : jvk j(重
となる．双曲型演算子 L に対しては，この多様体上で A
根含む) とこれらに対応する k 個の 1 次独立の固有ベクトル l1 : : : lk を持つことになる．初期多様体上の
擾乱は，一般に k 個の特性面上を，法線速度 v の波となって伝播する．各波の式は微分方程式 (17.17) の
考える．初期多様体を
左から固有ベクトル
jAj
l
を演算することによって求められる．これらの式は特性面内の微分すなわち通常
k
u
A
k
12 特性面上で
= 0 ということは，幾何学的に解釈すれば ， 次元の 空間において を構成する 個の列ベクトルまたは行
ベクトルが 次元空間を張るのではなく
1 次元の面内にあるということである．それゆえこの面に直交するナルベクトル，すな
わち を構成するすべての列ベクトルとのスカラー積がゼロになる行ベクトル ，すべての行ベクトルとのスカラー積がゼロになる
次元多様体内にあるときには，こ
列ベクトル が存在することになる．なお当然ながら更に を構成する 箇のベクトルが
の多様体に直交する
個の独立のベクトルを多様体の外に取ることができる．
A
k
r
k;
m
A
k
l
k ;m
3 独立変数以上の双曲型微分方程式 番めの従特性に沿う微分と
13
方向の局所平面波に沿う (初期多様体方向) 微分を含み，従特性に沿って初期
m 個の因子が同じで v が m 重根のと
きには，この m 重固有値に対応して 1 次独立の固有ベクトル l : : : lm が存在し，m 重特性面に沿って m
個の波が同じ法線速度 v で伝播することになる．
まず，行列 A のランクが k ; 1 の特性面の場合には，ナルベクトル l と r は大きさは任意に取れるが１つ
ずつ決まる．lAr = 0 を従特性に沿って微分すれば ，l_Ar = lAr_ = 0，また式 (17.2c) から @A=@ = A で
値問題として解かれる．なお，式 (17.18) を因子に分解したときに
1
あるから
_ =
lAr
他方，Q(
n
X
=0
n
X
=0
lA r _ = 0
) = 0 を従特性に沿って微分すれば ，
Q _ =0
以上の二つの関係は
_
によらず成立するから，lA
r と @Q=@
は比例関係にある．式 (17.14) の従特性を
求める式を用いれば次の関係が導かれる．
x_ = lA r
( = 0 1 : : : n)
(17.19)
ただし ，` _ ' は従特性に沿って適切に取られたパラメータに関する微分を示す．
k ; m の m 重特性面の場合には，1 次独立のナルベクトル l1 : : : lm r1 : : : rm が
存在する．上と同様に lj Ari = 0 を従特性に沿って微分すれば
次に A のランクが
lj A0 ri _0 +
他方，Q(
n
X
j
=1
l A ri _ = 0
) = 0 + Q~ (
1
_0 + X Q~ _ = 0
n
: : : n ) = 0 を従特性に沿って微分すれば ，
=1
これらの関係も
_
x_ = lj A ri
x_ 0 lj A0 ri
によらず成立するから，式 (17.14) を用いれば次の関係が導かれる．
(i j = 1 : : : m
= 1 2 : : : n)
(17.20)
14
圧縮性流れの解法|NS 方程式 17.1.4
特性面に沿う不連続の伝播
この節の始めの項で述べたように，無拘束の面
(x) = 0 上で u
が微分方程式から一意的に決まること
は，この面を横切って u が常に連続になることを示している．これに対し ，特性面
が一意的に決まらないので，この面を横切って u は連続になるが，u 更には u
(x) = 0 上では，u
等の微分は必ずしも連続
にはならない．以下には特性面を横切って u が不連続になる場合について述べる．微分方程式は準線形で
あるからここでは右辺 f
せば ，
L u] =
n
X
i=0
Ai
= 0 として議論を進める．ある特性面を横切っての a の跳躍 (jump) を a ] で表
Di u ] =
他方 行列 A のランクが k
n
X
i=0
Ai i u ] = A u ] = 0
;1 のときには，Ar = 0 を満足するナルベクトル r は一意的に決まるので，上式
と比較すれば ，
u ] g= r
(17.21)
のように置くことができる．ここに
は跳躍の大きさを示すスカラー量である．
= 0 を横切って関数の微分 u
ことができる (下図参照13 )．
特性面
が不連続になる場合には，その近傍で
u(x t) は次のように置く
u(x t) = h( )g +~u(x t)
(17.22)
ただし ，u
~ は連続関数，また h は Heaviside 関数で
h( ) =
(
1 ( 0)
0 ( < 0)
次にこの不連続の特性面に沿っての伝播について考える．u に演算子 lL を作用させれば ，
lL u] = hlL g] + hlL ]g + lL u~] = 0
この式の中辺第 2 項は lL が
= 0 面内の演算子であるから 0，また第 3 項も u~ が = 0 を横切って連続で
あるから無視される14 ．したがって
lL g] = lL ]r + lL r] = 0
u(x)
13
はベクトル量で図はそのふるまいを象徴的に示したものである．
14 ここでは不連続面
であるが，不連続面
0，
0
この議論は成立する．
=0
=
h=1(
) = 0 ( < 0 ) u = hg( ; 0 )+~u の一般の場合にも
3 独立変数以上の双曲型微分方程式 右辺第 1 項に関しては，式 (17.1a) より L
これより跳躍の大きさ
]=
15
P Ai D
P
i , 式 (17.19) より lAi r = x_ i ，また x_ i Di = _ ．
は従特性に沿って次式により伝播することが分かる．
lL g] = _ + lL r] = 0
(17.23)
= 0 を横切って関数 u が自身を含めて不連続になる より一般的な場合を考え，u(x t) を次
次に特性面
のように置く．
u(x t) = DW ( ) g0 (x t) + W ( ) g1 (x t) + D 1 W ( ) g2 (x t) +
+ u~(x t)
;
ただし
DW ( ) = h( ) W ( ) = h( ) D 1 W ( ) = h( )
;
2
(17.24)
: : : ．また g0 g1 g2 : : : u~ は連続関数の
ベクトルである．
L u] = Ag0 + h(L g0 ]+ Ag1) + h (L g1 ]+ Ag2) +
+ L u~] = 0
( ) = D2 W ( ) は Dirac のデルタ関数．この式が成立するためには次の一連の関係が満足されなければな
らない．
Ag0 = 0
L g0] + Ag1 = 0
L g1] + Ag2 = 0
(17.25a)
(17.25b)
(17.25c)
式 (17.25a) から特性上で式 (17.21) と同様の関係 g 0
r が成立することになり，また式 (17.25b) から
lL g ] = 0，したがって式 (17.23) と同様の 0 の従特性に沿っての伝播の式 lL g0] = _ 0 + lL r] 0 = 0 が導
;
かれる．また式 (17.25b) は A g 1 + A 1 L g 0 ] = 0 のように書くことができるから g 1 は次のように置くこ
=
0
0
;
とができる．
g 1 = 1 r + h1
ただし h1
(17.26)
= ;A 1 L g0 ]．また式 (17.25c) から
;
lL g1] = _ 1 + lL r] 1 + lL h1] = 0
以下同様にして，従特性に沿って伝播する不連続の大きさ 0
らの式を積分すれば跳躍量 g
また，
g=
0
(17.27)
1
2
: : : の式を導くことができる．これ
g g : : : を次々に決定することができる．
1
2
= 0 が m 重特性面の場合には，この面を横切っての不連続量は
m
X
i=1
i ri
(17.28)
と置くことができる．ただし r1
示す量である．
lj L g ] = lj
m nX
n
X
i=1
=0
: : : rm は 1 次独立な右ナルベクトル，
o
A ri D i + L ri ] i = 0
1
:::
m は不連続の大きさを
16
圧縮性流れの解法|NS 方程式 式 (17.20) を用いれば ，
m
X
i=1
( _ i + lj L ri ] i ) = 0
(j = 1 2 : : : m)
(17.29)
i (i = 1 2
: : : m) はこの連立常微分方程式から求めることができる．
以上この節では，xt 空間内の初期面 t = 0 上にある n ; 1 次元多様体 (x 0) = 0 に u の不連続 (擾乱) が
あるときに，この不連続は k 個の波となって初期多様体を含む k 個の特性面 (x t) = 0 に沿って伝播す
ること，ただしこれらの各特性面上の一部は重なっていることもあることを述べた．それから各特性面上を
伝播する不連続量 g
(x t) の式 (17.21)(17.23)(17.26){(17.29) を示した．以下にはこれらの方程式の初期値
問題を完成させるために，初期多様体内の不連続量 g (x 0) が，特性面上の不連続量の初期値 g (x 0) にど
のように分割されるのかについて述べる．初期面上の u(x 0) を式 (17.24) にならって一般性を持たせ次の
ように置く．
u(x 0) = DW ( ) g0 (x 0) + W ( ) g1 (x 0) + D 1 W ( ) g2 (x 0) +
;
+ u~(x 0)
(17.30)
DW ( ) = h( ) W ( ) = h( ) D 1 W ( ) = h( ) 2 : : : ．これらの不連続の各レベルに対し ，
g (x 0) の和は g(x 0) に等しくなるべきであるから，次式が導かれる．
;
ただし
k
X
=1
k
X
=1
0
(x 0)r (x 0) = g0 (x 0)
1
(x 0)r (x 0) = g1 (x 0) ;
k
X
=1
h1 (x 0)
(17.31)
(x 0) = ;A 1 L 0 (x 0)r (x 0)]．たとえば第 1 式は，k 個の未知変数 0 (x 0) を持つ k 個の式
からなる連立 1 次方程式で，これから各特性面に沿って伝播する不連続量の大きさの初期値 0 (x 0) を一
ただし h1
;
意的に決定することができる．
双曲型微分方程式 L u] = f で支配される場においては，xt 空間内の１つの点 O の関数 u の値は点 O に
頂点を持ち過去に広がる射線錐体 (ray
conoid) 内のデータのみによって決まる．この錐体内の領域は依存
域 (domain of dependence) と呼ばれる．また点 O の u の値は点 O に頂点を持ち未来に広がる射線錐体内
の u にのみ影響を及ぼす．この錐体内の領域は影響域 (domain of in uence) と呼ばれる．ある時間 t にお
ける領域 G 内の関数 u の値は，この G の境界 ; 上の各点から過去に広がる射線錐体の外側包絡面内のデー
タのみによって決まる．領域 G の影響域についても同様のことがいえる．流れの亜音速域では境界の形状
や条件の不連続性は数学的には領域の内部に及ばない．これに対し 超音速域ではこれらの不連続性によっ
て領域内部の流れも不連続性を持つことになる．不連続性を持つ流れにあってはテイラー級数の収束を前
提とする差分式の精度は意味をなさなくなる．流れの数値解析では一般に
2 次ないし 3 次精度の差分式で
十分で，より高次の差分式の使用は弱い不連続の計算など 特別の場合に限られる．
3 独立変数以上の双曲型微分方程式 17.1.5
17
非定常圧縮性等エント ロピー流れへの応用
この流れの基礎方程式は，非保存形では次のように書かれる (保存形方程式も，通常 特性の理論を適用
する前に非保存形に書換えられる)．
d =dt + r u = 0
duu=dt + rp = 0
ds=dt = 0
ただし
d=dt @=@t +u r は流跡線に沿う内微分である．第 3 式は流跡線に沿ってエントロピーが保存され
ることを表し ，これより等エントロピー流れの関係式
d ; 1 dp = 0
dt c2 dt
(17.32)
が流跡線に沿って成立することになり，この流れの方程式は次の u
p を未知変数とする方程式系になる．
duu=dt + rp = 0
dp=dt + c2 r u = 0
(17.33a)
(17.33b)
あるいは
0
d=dt
B
B
0
B
B
B
@ 0
c2 @=@x
ここで面
V
と置く．
0
0
d=dt
0
c2 @=@y
(17.33c)
(x t) = 0 を考え，
t +u
r = +u +v +w
= 0 が特性面ならば，特性条件が式 (17.3) と (17.33) から次のように求められる．
V
Q(
10 1
@=@xC B u C
@=@y C
CC BBB v CCC = 0
d=dt @=@z C
A B@ w CA
c2 @=@z @=@t
p
0
)=
0
0
c2
0
V
0
c2
0
0
V
c2
= V 2 V 2 ; c2 ( 2 + 2 + 2 ) = 0
V
これより次の２つの波の式が得られる．
(i)
V =0
(ii) V 2 ; c2 2 = 0
V = 0 の法線錐は，
(17.34a)
(17.34b)
空間の原点を通る平面である．これに対応する射線錐は式 (17.14) から
dx : dy : dz : dt = u : v : w : 1
18
圧縮性流れの解法|NS 方程式 となり，x t 空間内の流跡線になる15 ．すなわち法線錐の式 (17.34a) は流跡線に沿って伝播するエントロピー
波のものになる．
式 (17.34b) は 2 次の法線錐を表している．これに対応する射線錐の式は
dx ; u 2 + dy ; v 2 + dz ; w 2 = c2
(17.35)
dt
dt
dt
となる16 ．これより t = 1 の射線面の式は (x ; u)2 +(y ; v )2 +(z ; w)2 = c2 となり，u に中心を持つ半径 c
の球になる．射線錐は流跡線上の擾乱から音速で広がる圧力波 (acoustic
waves) のものになる．流線上のあ
る点からこのような圧力波が連続に放出されるときに圧力波の包絡面はマッハ錐 (Mach cone) を形成する．
次に 2 次元の場合を考える．式 (17.33c) から式 (17.18) に相当の式は次のようになる．
U;
0
2
c cos
ただし
流速 u の
0
U;
c2 sin
1
;
1
;
cos
sin
U;
= (U ; ) (U ; )2 ; c2 = 0
(17.36)
sin ) は波面の単位法線ベクトルで = r cos
= r sin ，また U = u cos +v sin は
方向成分である (上図参照)．式 (17.36) より一つの波面の固有値とそれに対応する左固有ベク
= (cos
15 式
(17.14) と (17.34a) から
dx=dt = Q =Q = V =V = u
16 式
(17.14) と (17.34b) から
dx = Q = u ; c2
dx ; u 2 = c4 2 = c2 2
2
dt Q
V
dt
V2
3 独立変数以上の双曲型微分方程式 19
トルは 17
1
=U
l1 = (sin
; cos
0)
したがって
1 p sin ; dv + 1 p cos = 0
l1L u] = du
+
x
y
dt
dt
これを書換えれば
du sin ; dv cos = ; 1 dp
(17.37)
dt
dt
d
ただし d=dt = @=@t + u@=@x + v@=@y は流跡線に沿う微分，d=d = sin @=@x ; cos @=@y は t = const. 面内の
波面に沿う微分である．
残り二つの波の固有値と左固有ベクトルは
i=U
c : li = (cos
;
ただし ，
sin
1= c)
(i = 2 3)
に取れば ，複合は+のみで良いことになる．
1
dv 1
1 dp
li L u] = du
dt + px cos ; dt + py sin + c dt + c(ux + vy ) = 0
これを書換えれば
~ du
~
~
1 dp
dv
du
dv
+
cos
+
(17.38)
c dt dt
dt sin = ;c d sin ; d cos
~ = @=@t +(u + c cos )@=@x +(v + c sin )@=@y は圧力波の射線錐の従特性に沿う微分である．d=d
ただし d=dt
は上と同じ定義である．
式 (17.32) と (17.37)(17.38) は適合方程式 (compatibility
equations) とも言われ，前章の 2 独立変数の場
合の式 (16.24) に相当の式である．式 (16.24) は常微分方程式で，図に示す格子特性法などにより解を求め
3 独立変数の式 (17.38) は，従特性に沿う微分演算子のほかに波面に沿う微分 d=d
も含み，また波面の方向 を変えることによって無数の式を作ることができる．そのため この式は 2 独立
変数の式 (16.24) ほど 決定的意味を持つものではないとも言われてきた．3 独立変数の特性の方法も数多く
提案されている．その中で図 17.8 に示すように，点 O を通る流跡線と射線錐体上の４つの適合方程式を用
ることができた．しかし
いる従特性法は標準的なものである．
(17.36) より
;1 cos 1
0
;1 sin A = 0
U;
2
2
c cos
c sin
U;
1
この式に固有値 =
= u cos + v sin を入れ，(l1 l2 l3 ) = (l11 l21 l31 ) と置けば
0 0
;1 cos 1
0
;1 sin A = 0
(l11 l21 l31 ) @ 0
0
c2 cos
c2 sin
0
左固有ベクトルと行列の第 1 列，第 2 列のスカラー積が 0 になるためには固有ベクトルの第 3 成分 l31 = 0，第 3 列から l11 = sin
l21 =
; cos ．固有ベクトルはある任意性を持ち簡単なものをとる．
17 左固有ベクトルを求める式は式
0 U;
(l1 l2 l3 ) @ 0
20
圧縮性流れの解法|NS 方程式 図
17.1.6
17.8:
従特性法
電磁流体力学への応用
圧縮性流れにおいて，流跡線に沿ってエントロピーが保存され圧力波がマッハ錐を形成することは，特性
の理論をまつまでもなく周知のことである．電磁流体力学は本章の範囲外であるが，特性の理論によっては
じめてその全貌が明らかになる例としてここに取上げることにする．磁場の作用するプラズマ流を扱う電
磁流体力学 (MHD,
magnetohydrodynamics) の基礎方程式は次の Maxwell の式と Euler 方程式である．
rD= e
rB =0
r H = J + Dt
r E = ;B t
t + r ( u) = 0
duu=dt = ;rp + j B
ただし記号は通常使われているもので D
密度である．なお
e
(17.39a)
(17.39b)
(17.39c)
(17.39d)
(17.39e)
(17.39f)
B は電束密度と磁束密度，H B は電場と磁場の強さ，j は電流
D t ;B t j B は電荷密度，変位電流，Faraday 起電力，Lorentz 力である．また
Ohm の法則，電束密度，磁束密度は次のようになる．
j = (E + u B )
D= E
B= H
ただし
(17.40a)
(17.40b)
(17.40c)
は導電率 (conductivity)，誘電率 (permittivity)，透磁率 (permeability) である．式 (17.39)
は上から，電場の保存，磁場の保存，Faraday の法則，Ampere の法則，質量の保存，Euler の式である．な
お上記の Maxwell の式はもともとの式とはやや異なり光速度に関する項が省かれている．
次に電磁流体力学の方程式系を導く．式
式
(17.40) の第２式を用いれば ，
j =r B
(17.39c) において変位電流は十分小さいものとし て無視し ，
3 独立変数以上の双曲型微分方程式 21
また Ohm の法則で導電率が十分大きいものとすれば，
E = ;u B
u B の方程式系が導かれる．
これらの MHD 近似を用いれば，式 (17.39) から次の
rB =0
t + r ( u) = 0
duu=dt + rp ; 1 (r B ) B = 0
(17.41a)
(17.41b)
B t ; r (u B ) = 0
(17.41d)
(17.41c)
なお第１式は初期条件として用いられるものである18 ．
まず非圧縮性流れの場合について述べる．上式に
(r B ) B = (B r)B ; 12 rB 2
r (u B ) = (B r)u ; (u r)B ; B r u + u r B
なる関係を用いれば次の非圧縮性 MHD の方程式系が得られる．
r u=0
duu + r p + 1 B 2 ; 1 (B r)B = 0
dt
2
B
dB
dt ; (B r) u = 0
(17.42b)
(17.42c)
(x t) = 0 を考え，
ここで特性面
V
(17.42a)
t +u
r = + u1 1 + u2 2 + u3
3
のように置く．特性条件は式 (17.42) から次のようになる．
0
1
1
2
Q( ) =
3
0
0
0
a
V
0
0
;B
0
0
r r a
2
3
0
0
0
V
0
0
;B
0
V
0
0
;B
0
B1 1 ;B
B1 2
B1 3
V
0
0
B2
0
1
B2 2 ;B
B2 3
0
V
0
B3
B3
1
2
B3 3 ; B
=0
0
0
V
@ r B )=@t = 0 となり，磁場の初期値が条件 (17.41a) を満足す
に対し
(
) = 0 であるから，最後の式は (
るように与えられればその磁場は常にこの条件を満足することになる．
18 任意ベクトル
0
22
圧縮性流れの解法|NS 方程式 この行列式を計算すれば 19 ，
6= 0 であるから
V 2 ; (B )2 2 = 0
V
(17.43)
ここで流れに乗って移動する座標系を用いることにすれば
Alfven 速度を導入すれば上式は次のようになる．
2
V =
，また
b = jB j( ) 1=2 で定義され る
;
; (b )2 2 = 0
これより
(i) エントロピー波：
=0
(ii) Alfven 波： 2 ; (b )2 = 0
(17.44a)
(17.44b)
1 式はこの座標系の流跡線に沿って伝播するエントロピ ー波を表し ている．また第 2 式はアルヴェン
(Alfven) 波と呼ばれるもので，以下これについて簡単に説明する．
Alfven 波の法線面，法線速度面，射線面は，B 方向の座標 1 に関し 軸対称である．法線面は，極座標
r を用いれば式 (17.44b) から 1;(br cos )2 = 0 1 = r cos = 1=b となり， 1 軸に垂直な２つの平面
になる．法線速度面は r = b cos となり，座標の原点に接し 1 軸上に中心を持つ２つの球面になる．ま
た射線面は 1 = b となり 1 軸上の２つの点になる．Alfven 波は磁場の擾乱が磁力線に沿って伝播するも
ので，天文学的スケールのプラズマ (例えば太陽の黒点) では磁場が流体に凍結 (frozen in) されているので
第
重要であるが，実験室レベルでプラズマが磁場内を通りぬける場合には無視できよう．
19 この行列式は次のように計算することができる．まず行列式の各行，各列に適当な係数を掛ければ
0
1
; V =BB
2
3
0 0 0
0
0
0
; V 2 =BB
0
この部分は同じ
Q = V (B )2
0
0
; V 2 =BB
0
1
0
0
1 0 0
0
0
1
0
0 1 0
0
0
0
1
0 0 1
次に第 2 列から第 5 列，第 3 列から第 6 列，第 4 列から第 7 列を引算すれば行列式は小さくなり
0
1
2
3
;B1 1 + A ;B2 1
;B3 1
= V (B )2 1
;B1 2 ;B2 2 + A ;B3 2
2
;B1 3
;B2 3 ;B3 3 + A
3
B ．更に計算すれば
ただし A B ; V 2=B
0
1 =B1
2 =B2
3 =B3
;1+ A=B1 1
;1
;1
= V (B )2 B1 B2 B3 1 2 3 1
1
;1
;1+ A=B2 2
;1
1
;1
;1
;1+ A=B3 3
2
2
2
0
; 2 12 22 32
1
2
3
1 0 0 = V (B )2 A2 0
1 0 0
= V (B )2 A2 1
1
0 1 0
0
0 1 0
1
0 0 1
0
0 0 1
1
2
3
2
3 独立変数以上の双曲型微分方程式 23
17.9: Alfven 波の法線面，法線速度面，射線
図
次に圧縮性流れの場合について述べる．式 (17.41) は
rp = c2 r (1=2)rB 2 = (B r)B などの関係を用
いれば次のようになる．
d + r u=0
dt
duu + c2 r + 1 B rB ; 1 (B r)B = 0
dt
B
Bd
dB
dt ; (B r)u ; dt = 0
(17.45a)
(17.45b)
(17.45c)
特性の条件は次のようになる．
V
c2
c2
c2
Q( ) =
1
2
V
1
0
3
;B1 V=
;B2 V=
;B3 V=
V
0
;B
0
0
1
B2 2 ;B
B2 3
0
V
0
0
;B
0
B2
V
0
0
;B
0
0
0
B1 1 ;B
B1 2
B1 3
0
0
V
0
0
2
0
3
0
B3
B3
1
2
B3 3 ;B
=0
0
0
V
この行列式を計算すれば 20 ．
V V 2 ; (b )2 V 4 ; (c2 + b2 ) 2 V 2 + c2 (b )2
2
=0
(17.46)
更に流れに乗って移動する座標系を用いれば
2
; (b )2
4
; (c2 + b2)
2 2
+ c2(b )2
2
=0
20 上と同様に計算すれば，
0
;c2 1 ; V 2 B1 =(B )
;
;cc22 2 ;; VV 22 BB2 ==((BB ))
Q = V (B )2
0
= c V (B )
2
2
1
2
3
3
3
1
2
1
3
;B1 1 + A ;B2 1
;B3 1
;B1 2 ;B2 2 + A ;B3 2
;B1 3
;B2 3 ;B3 3 + A
3
; V B
3
上と同じ
2
0
B1
B2
B3
1
2
上と同じ
3
24
圧縮性流れの解法|NS 方程式 これより
(i) エントロピー波：
=0
2
(ii) Alfven 波：
; (b )2 = 0
4
; (c2 + b2) 2 2 + c2(b )2 2 = 0
(iii) 圧力波：
(17.47a)
(17.47b)
(17.47c)
エントロピー波と Alfven 波については既に述べたので以下には圧力波についてのみ説明する．
圧力波の法線面の式は式 (17.47c) から
1 ; (c2 + b2) 2 + c2 (b )2 2 = 0
を極座標 (r
) で表せば
1 ; (c2 + b2)r2 + c2 b2r4 cos2 = 0
この r の 4 次式は D
r = p1
1
(c2 + b2)2 ; 4c2b2 cos2 = (c2 ; b2 )2 +4c2b2 sin2
1 ;c2 + b2 pD
p1
12 =
2
2
0 であるから 4 実根を持つ．
特に
r = 1=c 1=b
p
r = 1= c2 + b2
=0 :
= =2 :
漸近線は
1
!
=2 1=r ! 0 と置けば得られる．
= r cos =
p2
c + b2 cb
CG で描いた圧縮性 MHD 流れの法線面，法線速度面，射線面をそれぞれ図 17.10，17.11，
17.12 に示す．Alfven 波は非圧縮性流れのものと同じである．またエントロピー波も同じで座標の原点にあ
以上を基に
1 項は非圧縮性流れと同じで，第 2 項の行列式 q2 は例えば次のように計算できる．
0
1 =B1
2 =B2
3 =B3
;1+ A=B1 1
;1
;1
q2 = B1 B2 B3 1 2 3 B1 = 1
B2 = 2
;1
;1+ A=B2 2
;1
B3 = 3
;1
;1
;1+ A=B3 3
1
0
1
1
1
0
0
1 =B1
2 =B2
3 =B3
= B1 B2 B3 1 2 3 1
B1 = 1
A=B1 1
0
0
1
B2 = 2
0
A=B2 2
0
1
B3 = 3
0
0
A=B3 3
A
0
B1 1 B2 2 B3 3
2
2
2
0
0
1
2
3
A ;B
;B2
=A 1
=
A
= A V 2 ; B2 2
B1 = 1
1
0
0
; 2 ;B
1
B2 = 2
0
1
0
1
B3 = 3
0
0
1
第 2 式は第 1 式の上と左に要素を補って計算したものである．また第 4 式は第 3 式の第 1 列の第 3 要素以下を 0 にし ，それから第
2 列の第 3 要素以下を 0 にすれば得られる．
この式の第
3 独立変数以上の双曲型微分方程式 25
る．しかし圧力波は通常の圧縮性流れのものとは全く異なり，法線速度面上でいえば
る速い波 (fast
Alfven 波の外側にあ
waves) と内側にある遅い波 (slow waves) の２つになる．なおこれら 3 つの図では波面の方
向の同じ波は同じ色で表わされている．遅い圧力波の法線面には変曲点があり，これらの点は射線面の尖点
(cusp) に対応している．磁場の作用するプラズマ中を磁場の方向に伝播する擾乱は一般に１つのエントロ
ピー波，２つの Alfven 波，4 つの圧力波からなる．
c<b
図
c<b
図
17.10:
17.11:
c<b
図
c=b
c>b
c=b
c>b
c=b
c>b
圧縮性 MHD 流れの法線面
圧縮性 MHD 流れの法線速度面
17.12:
圧縮性 MHD 流れの射線面
26
圧縮性流れの解法|NS 方程式 17.2
圧縮性流れの基礎方程式
圧縮性流れの基礎方程式は質量，運動量，エネルギーの保存の式で次のように書かれる．
@ + u=0
(17.48a)
@t
@
+f
(17.48b)
@t u + ( uu + p) =
@e + Huu = ( u) q +f u
(17.48c)
@t
ただし は密度，u は流速，e = ( + u 2=2) は単位体積当たりの岐点内部エネルギー，p は静圧，H =
h+u2=2 = (e+p)= は岐点エンタルピー， は粘性応力テンソル，q は熱流束，f は外力である．式 (17.48)
r
r
r
r
r
;r
は圧縮性ナヴ ィエ・ストークス方程式 (the
compressible Navier-Stokes equations) と呼ばれ，層流または
乱流の直接数値シミュレーション (DNS, direct numerical simulation) に用いられる．乱流の実用計算につ
いては次章に述べる．式 (17.48) から拡散項すなわち と q を含む項を除いた式は圧縮性オイラー方程式
(the compressible Euler equations) と呼ばれる．
h = +p= は静温
ここで圧縮性流れの数式によく出てくる記号をあらかじめまとめて説明しておく．T
度，比内部エネルギー，比エンタルピーで，完全気体を仮定すれば，p = RT d = cv dT dh = cp dT であ
R は定積比熱，定圧比熱，比熱比，気体定数でこれら 4 定数の間には = cp =cv R = cp cv
る．cv cp
なる関係があり，２つが分かれば他の２つは決まる．c =
RT は音速，M = u =c はマッハ数である．ま
た粘性応力テンソル (viscous stress tensor) と熱流束 (heat ux) q の成分は次のようになる．
;
2
(i j = 1 2 3)
(17.49)
ij = uj i + ui j 3 ij uk k
2
@c
qi = T i = 1 1 P @x
(i = 1 2 3)
(17.50)
;
p
j
j
;
;
ただし
;
r
;
i
uj i = @uj=@xi ， ij はクロネッカー
ル数 (Prandtl
number) である．
関数， は粘性係数， は熱伝導率，Pr
= cp =
はプラント
17.2.1 デカルト 座標系の圧縮性 Navier-Stokes 方程式
圧縮性ナヴ ィエ・ストークス方程式 (17.48) はデカルト座標系 (cartesian
coordinates) x では次のように
なる．
@q + @Fi + D + g = 0
@t @xi
0 1
B
C
B
u1 C
B
CC
q =B
B
u
2C
B
C
B
@ u3CA
e
0
1
ui
BB u u + pCC
BB 1 i 1i CC
Fi = B
BB u2ui + 2ipCCC
@ u3ui + 3ipA
Hui
(17.51a)
(i = 1 2 3)
(17.51b)
圧縮性流れの基礎方程式 D=
;
0
BB
@ B
B
BB
@xi B
@
0
1i
2i
3i
ij uj ; qi
1
CC
CC
CC
CA
g=
27
0 1
BB f0 CC
BB 1 CC
BB f2 CC
B@ f3 CA
;
fi ui
ただし 上式の q は未知変数，Fi は流束，D は拡散項，g は外力項である．
17.2.2
一般曲線座標系の圧縮性 Navier-Stokes 方程式
デカルト座標系を x ，一般曲線座標系 (general
の変換の測度 (metrics)
x ::: と
0
10
x
x
x
B
CCBB x
B
y
y
y
@
A@ x
z
z
z
1 0
で表せば ，これらの座標間
1
(17.52)
(jacobian)J は
また変換のヤコビアン
x x
J= y y
z z
: : : の間には次の関係が成立する．
zC B 1 0 0 C
zC
A = B@ 0 1 0 CA
0 0 1
z
y
y
y
x
x
curvilinear coordinates) を
x
y
z
(17.53)
となる．これらは２次元の場合には次のようになる．
x =J
y
x = J
x に関する微分を
@ @j @
@xi = @xi @ j
;
y = J
y
;
x
y =J
x
J =x y x y
;
(17.54)
に関する微分に書き換えれば次のようになる21 ．
r
= ( j ) @@
j
(17.55)
r
具体的には，例えば
y
@ =1 z
@x J
@=@
y
z
@=@
y
z
@=@
また２次元では
@
@ 1 @
@x = J y @ y @
;
21 ここでは Einstein 総和規約
次元では次のようになる．
(Einstein cnovention)，a1 b1 +a2 b2 +a3 b3 = a b
i i
@ @ @ @ @ @
@
@x = @x @ + @x @ + @x @
i
i
i
i
が用いられる．例えば式 (17.55) の第１式は３
28
圧縮性流れの解法|NS 方程式 また一般に @ (J@
j=@xi )=@ j = 0
@ (J j )=@ j = 0 で次の関係が成立する．
@ =J@ j @ = @ J@ j
J @x
@xi @ j @ j @xi
i
J = @@ J ( j )
j
@
J 2 = @ Jgjk @@
j
k
ただし (g ij ) は測度テンソル (metric tensor) で，その成分は次のように定義される．
@i@j
gij = i j = @x
k @xk
r
r
r
r
r
r
流線に沿う内微分 (interior
u
r
(17.56a)
(17.56b)
(17.56c)
(17.57)
derivative) は次のように書き換えられる．
@ =u @ j @ =U @ =U ~
= ui @x
i @x @
j@
i
i j
j
(17.58)
r
Uj は反変速度 (contravariant velocity) の j 方向成分である．
写像空間の反変速度 U と物理空間の流速 u の間には次の関係がある．
@ ju
iU
Uj = @x
ui = @x
i
j
@
i
j
u は成分で表せば
流れの渦度 と流速 u の関係 =
@uk
i = ijk @x
j
ただし
(17.59)
r
ただし
( ijk ) は Eddington の
である22 ．写像空間の反変渦度 (contravaiant
(17.60)
vorticity)Z は物理空間の渦度
に対し次のように定義される．
@j
Zj = @x
i
i=
i
@xi Z
@j j
(17.61)
このとき反変渦度と反変速度の関係は次のようになる．
Zl =
r
l
=
r
l
r
m
= J1 lmn @ @ (gnj Uj )
m
ただし
@uu 1
@m=J
@xi @ui = 1
n@ m J
lmn @
lmn @
@
m
@xx
@nu
(17.62)
(gij ) は測度テンソルで，(gij );1 = (gij )，その成分は次式で定義される．
k @xk
gij = @@xx @@xx = @x
@
i j
i @ j
(17.63)
なおデカルト座標系から一般曲線座標系への変換に関して不明の点は 13.6 節を参照いただきたい．
Eddington の は次のように定義される．
8< 1 ((ijk) = (123) (231) または (312))
= : ;1 ((ijk) = (321) (213) または (132))
0 (その他の場合)
例えば， 1 = 1 @u =@x = @u3=@x2 ; @u2=@x3 である．
22
ijk
jk
k
j
29
圧縮性流れの基礎方程式 一般曲線座標系 の圧縮性 Navier-Stokes 方程式は，式 (17.48) に上述の関係を用いれば比較的容易に導
出できる．すなわち式 (17.48) にヤコビアン
0
@ JB
B
@
@t
J を乗じ，式 (17.56b) を用い，反変速度 Ui を導入すれば
1
0
1
U
i
C+ @ J BB
C @
uC
A @ i @ uUi + ipCA = @ i (J
r
e
0
B
i) B
@
r
HUi
0
u q
1 0
CC+ J BB 0
A @g
1
CC
A
gu
;
これより次の一般曲線座標系の圧縮性ナヴ ィエ・ストークス方程式が導かれる．
@ q^ + @ F^i + D^ +^g = 0
@t @ i
0 1
BB u CC
BB 1CC
q^ = J B
BB u2CCC
@ u3 A
e
D^ =
;
0
BB
@J B
B
BB
@ i i jB
@
(17.64a)
0
1
Ui
BB u U + pCC
BB 1 i i 1 CC
^
Fi = J B
BB u2Ui + i 2pCCC
@ u3Ui + i 3pA
0
j1
j2
j3
u
jk k ; qj
1
CC
CC
CC
CA
HUi
g^ =
(i = 1 2 3)
0 1
BB f0 CC
BB 1 CC
JB
BB f2 CCC
@ f3 A
(17.64b)
;
fk uk
q^ は未知変数，F^i は流束，D^ は拡散項，g^ は外力項である．また J = @ (x y z )=@ (
)， i j =
@ i =@xj は変換のヤコビアンと測度である．式 (17.64) は，式 (17.51) の x の微分を の微分に変え，曲線
ただし
座標格子を用いる数値計算に適するようにしたものである．
式 (17.64) は，数学的に出てきたものでその解釈は不要との考えもあろうが，格子形成や境界条件の設定
のためにはその意味をある程度把握しておくべきであろう．まずヤコビアン
J の意味を考えよう．図 17.13
を参照すれば ，
x x
J= y y
z z
Jは
x
y
z
(x 100 x 000 ) (x 010 x000 ) (x 001 x000 )
;
;
;
空間内の単位立方体要素の x 空間における体積，すなわち 物理空間の局所体積] : 写像空間の相当体
積] の比を表している．次に p にかかる係数，例えば
J
J
r
は
r
= Jx Jy Jz
y y
=
z z
z
x
J
z
x
r
の意味を考えよう．下図を参照すれば ，
x
y
x
y
!
(x 010 x000 ) (x001 x 000 )
;
;
= const. 面に相当する x 空間の面の面積ベクトルである．またその
= const. 面の局所面積] : 写像空間の相当面積] の比を表している．
空間内の単位立方体要素の
大きさは 物理空間の
ここで物理空間 x の流れと写像空間 の流れの対応を考えれば ，x 空間で ある流体素片が１つの流跡線
に沿って流速 u で移動するときに，その流体素片は
空間では相当の流跡線に沿って反変速度 U で移動す
30
圧縮性流れの解法|NS 方程式 H
HH W6 HHHH
HH
U H
HHj V :
:
U
V
H
HHj
HH W6 HHH
HHH
001
001
000
6:
HH
j
010
100
000
010
000
000
100
図
17.13:
空間内の単位立方体格子
る23 ．また 空間の密度，運動量，比内部エネルギー，比エンタルピー，物体力はそれぞれ x 空間のものに
J をかけたものになり，
つまり
空間の静圧，応力テンソル，熱流束は x 空間のものに J r をかけたものになる．
空間における諸量は x 空間の諸量に対し ，流速は写像に伴う流線方向の長さの変化に比例し ，密
度，運動量などは写像による体積の変化に反比例し ，静圧などは写像による面積の変化に反比例して変化
する．例えば密度に関しては２つの流れ場の対応部分の質量は同じになり，また静圧に関しては対応する局
所面積に作用する力の大きさは同じになるということである．式 (17.64) はこのように定義された諸量に対
する質量，運動量，エネルギーの保存の式と解釈することができる．
式 (17.64) には x1
x2 x3 方向成分の運動方程式が含まれるが，以下では 1 2 3 方向成分の運動方程
式を導く．反変速度成分 Ul は流速 u に
l をかけたものであるから，この方程式は式 (17.48b) に J l
r
r
をかけることによって導かれる．なおその過程では
u = J Ul
J l ( uu + p) =
J
r
r
l
r
r
l
@
@ i J ( Uiu +
r
i p) =
@
@
il
@ i J ( Ui Ul + g p) J ( Uiu + i p) @ i
;
r
r
l
のように計算される．このようにして得られた運動方程式に上記の連続方程式とエネルギー方程式を補っ
たナヴ ィエ・ストークス方程式は次のようになる．
^
~
B @@tq^ + @@Fi + D^ +^g = @@tq~ + @@Fi + R~ + D~ +~g = 0
0i 1
0 i
1
Ui
BB U CC
BB U U + gi1pCC
BB 1CC
BB 1 i i2 CC
~
q~ = J B
F
=
J
BB U2Ui + g pCC
i
BB U2CCC
B@ U3Ui + gi3pCA
@ U3A
e
HUi
(17.65a)
(i = 1 2 3)
x の変化，実質微分 @=@t+u r はその間における流体素片に関する
U u r は単位時間における流体素片の の変化すなわち写像空間の流速
u
23 流速ベクトル は物理空間内の流体素片の単位時間における
量の変化を表す微分演算子である．したがって反変速度 =
である．
y
6
u
-
x
* +u r
6
U
-
1 +U
31
線 形 化 と 対 角 化 R~ = J ( Uiu +
;
0
BB
@B
B
i p) @ B
iB
B@
0
r
r
r
r
D~ =
0
1
C
1C
CC
2C
CC
3A
(17.65b)
0 1
0 1
0
BB CC
BB 00 CC
1 i0 C
CC
@ Jg k @T B
@j @ B
B
C
B
J @x
B
C
B
0
ik 2 i0 C @
ij @ B C
k@ j B
jB C
B@ 3 i0 CA i
@ 0 CA
g~ =
;
;
0
BB 10
B
B=B
BB 0
B@ 0
0
ui
0
0
0
11
12
13
21
22
23
31
32
33
0
0
0
1
0
C
0C
C
CC
0C
0C
A
0
BB
B
J fB
BB
B@
0
r
;
r
r
1
u
1
C
1C
CC
2C
CC
3A
(17.65c)
1
~ は付加項，D~ は拡散項，g~ は外力項で，式 (17.65c) の B
ただし式 (17.65) の q~ は未知変数，F~i は流束，R
は流速 u を反変速度 U に変換する行列である．また g ij
付いている量は微分しないものとする．@ (a0 b)=@
= i k j k は反変測度テンソルの成分で，添字 0 の
= a0 @b=@
．
ux)J U の運動方程式である．この質量流束の i 成分 J Ui
は 空間の単位立方体格子を x 空間に写像した格子セルの i = const. 面を通る質量流量である．この運動
方程式の 3 成分の式は，式 (17.64) の保存形運動方程式の 3 成分の線形結合を取ったもので，非保存形であ
~ が生じている．この
る．式 (17.65a) は，対流項を無理遣り 保存形に書換えたもので，そのために付加項 R
式 (17.65) の運動方程式は質量流束 (mass
付加項の大きさは滑らかな格子では小さくなる．なお下記の解法は
形近似因子法によるもので，その左
辺の演算子には式 (17.65a) の @ q~=@t + @ F~i=@ i ，右辺には式 (17.65a) の左辺からの式が用いられる．
17.3
線形化と対角化
圧縮性ナヴィエ・ストークス方程式は，拡散項を含むので正しくは双曲型ではないが，この方程式の拡散
項は流れの大部分の領域でほかの項に比べ十分小さいので，この方程式は，オイラー方程式に拡散項がた
だ単に加わったものとして，双曲型のオイラー方程式の解法で解かれている．後述の流束差分分離法または
有限体積法では，オイラー方程式の流束の微分が流束分離法によって離散近似される．流束分離法は，各座
標方向に伝播するエントロピー波，圧力波，剪断波などの波を，それぞれの向きよって個別に上流化するも
ので，物理現象に則した妥当な解を安定に得るものである．本節ではそのために，各座標方向の１次元オイ
ラー方程式に対し線形化 (linearization)，非保存形化，更に対角化 (diagonalization) を行うが，これは前章
に述べた 1 次元双曲型微分方程式の特性の理論によって，１次元オイラー方程式をこれらの波の常微分方
程式の系に変換することにほかならない．最近は，衝撃波の捕獲に保存性が重要という観点から，基礎方程
式として保存形 (conservative
form) で書かれたナヴ ィエ・ストークス方程式 (17.48) が用いられているが，
32
圧縮性流れの解法|NS 方程式 以前には簡単な次の非保存形 (nonconservative
form) ナヴ ィエ・ストークス方程式も用いられていた24 ．
@ +u + u = 0
@t
@uu + u u + 1 p = 1 (
+ f)
@t
@p + u p + c2 u = ~ q
@t
r
r
r
(17.66a)
r
r
(17.66b)
r
r
;
(17.66c)
r
17.3.1 デカルト 座標系のオイラー方程式
はじめに デカルト座標系のオイラー方程式を考える．１つの空間次元の部分のみを取り出せば
@q + @Fi = 0
(i = 1 2 3)
(17.67a)
@t @xi
ただし 添字 i は 1 または 2 または 3 で，ここではアインシュタイン総和規約には従わないものとする．未
知変数ベクトル q と流束ベクトル Fi は次のようになる．
0 1
B
C
B
u1 C
B
CC
q =B
B
u
2C
B
C
B
@ u3CA
e
0
1
ui
BB u u + pC
BB 1 i 1i C
C
C
Fi = B
u
2 ui + 2i pC
BB
C
@ u3ui + 3ipC
A
(i = 1 2 3)
(17.67b)
Hui
流束 Fi は q の成分の 1 次の同次式であるから，オイラーの同次関係を用いれば
i
dFi = @F
@q dq = Ai dq Fi = Ai q dAi q = 0
(17.68)
となり，式 (17.67) は次のように線形化することができる．
@q + A @q = 0
@t i @xi
ただし ヤコビ行列 Ai = @Fi =@q は
0
0
BB u u +
BB 1 i i1
Ai = B
BB u2ui + 2i
@ u3ui + 3i
;
;
;
i1
Di1
i1 u2 ; i2 ~ u1
i1 u3 ; i3 ~ u1
i1 H ; ~u1 ui
2
2
2
Hui + 2 ui
;
(i = 1 2 3)
i2
i2 u1 ; i1 ~u2
Di2
i2 u3 ; i3 ~u2
i2 H ; ~u2 ui
i3
i3 u1 ; i1 ~u3
i3 u2 ; i2 ~u3
Di3
i3 H ; ~u3 ui
(17.69a)
1
C
i1 ~C
CC
~
i2 C
CC
i3 ~A
0
(17.69b)
ui
24 ナヴィエ・ストークス方程式は周知のように，熱力学でいう流体要素のシステムに対し質量・運動量．エネルギーの保存則を適用
することによって導かれた保存形方程式が先にあり，非保存形方程式はその線形結合を取ることによって導かれるものである．圧力
の輸送方程式は 2 (17.48a) ~ (17.48b)+~ (17.48c) と置くことによって導かれる．
2
; u
; ~u u + ~ e = p
2
r u ; ~u r( uu + p)+~ r Huu = u rp + c2r u
線 形 化 と 対 角 化 ~=
;
1
2
33
= ~ u2=2 Dij = ui + ij (1 ~)uj である25 ．
;
デカルト座標系の非保存形オイラー方程式の１つの空間次元の式は，式 (17.66) から次のようになる．
@q + A @q = 0
@t i @xi
0 1
BBu CC
B 1CC
q =B
BBu2CC
B@u3CA
(i = 1 2 3)
0
BB u0i
B
Ai = B
BB 0
B@ 0
p
i1
i2
0
0
ui
ui
0
i1
i3
0
c
2
0
0
0
c
2
i2
ui
i3 c2
1
C
i1 = C
CC
i2 = C
CC
i3 = A
(17.70a)
0
(17.70b)
ui
式 (17.70a) は，
dq = @q
@q dq = Ndq
で定義される行列 N を式 (17.69a) の左から演算したものである．なお
(17.71)
Ai = NAi N ;1
である26 ．
保存形から非保存形への変換の行列 N
0
1
B
B
u1 =
B
B
N =B
u2 =
B
B
@ u3 =
;
;
;
~u2=2
0
1
;
BB u
B 1
N ;1 = B
BB u2
B@ u3
u2=2
25
0
0
1=
0
0
1=
0
0
~ u 1 ~ u2
0 0 0
0 0
0
0
0 0
1
0 0
C
0 0C
C
CC
0 0C
1= 0C
A
~u3 ~
1
0
C
0C
C
CC
0C
0C
A
u3 1 = ~
;
u1
= @q =@q とその逆行列 N ;1 は次のようになる．
u2
(17.72a)
;
(17.72b)
A の要素は，まず F の成分の式を q の成分で表し ，
F 2 = u1 u + 1 p = u1 u + 1 ~e ; 1 2~ ( u1 )2 +( u2 )2 +( u3 )2
F = Hu = e u ; ~ u ( u )2 +( u )2 +( u )2
i
i5
i
i
i
i
i
2
i
2
i
i
1
2
3
q の成分で偏微分すれば容易に求めることができる．なお
2
2
2 = 1 ~ (u 2 + u 2 + u 2 )
H = e+p = c + = e ; 2
それから
~
q
2
q
1
2
3
は の成分の 1 次の同次式ではないから，上記のような関係
の同次関係については代数の教科書あるいは 16.2 節の脚注参照．
26 蛇足ではあるが，
q = Nq と dNq = 0 は成立しない．オイラー
34
圧縮性流れの解法|NS 方程式 これらの行列の導出については脚注参照27 ．
ヤコビ行列 Ai は直接対角化することもできるが，その計算はかなり面倒になるので，ここではまず
を対角化しそれから Ai を次式によって対角化することにする．
Ai = N ;1Ai N = N ;1 Ri i Li N = Ri i Li
Ai
(17.73)
i は行列 Ai の固有値 ik からなる対角行列，Li は Ai の左固有ベクトル lik からなる行列，Ri = Li ;1 は
右固有ベクトル rik からなる行列である．ここでは次に Ai の固有値 ik と左固有ベクトル lik の求め方に
ついて簡単に説明する．固有値 ik は i の 5 次代数方程式
ui
Ai
j
iI j =
;
の根で， i1
i
;
0
0
0
0
ui
= ui (3 重根)
;
;
0
0
i1
2
i1 =
ui ; i
i2 = = 0
0
ui ; i i3 =
2
i2 c
i3 c2 ui ; i
0
i
c2
0
0
i3 = ui ; c となる28 ．左固有ベクトル lik を求める式は固有値
i1
0
0
0
0
27
0
i3
0
BB00
B
lik5 )B
BB0
B@0
i I ) = (lik lik
1
i2
i2 = ui + c
i = ui に対しては
li (Ai
i1
i1
i2
0
0
0
c2
i2
1
C
i1 = C
CC
i2 = C
CC = 0
i3 = A
0
i3
0
0
0
c2
i3
c2
0
N の要素は，まず q の成分の式を
q2 = u1 = u1 =
q5 = p = e ; 2~ ( u1 )2 +( u2 )2 +( u3 )2
q
q
のように の成分で表し ，それから の成分で偏微分すれば容易に求めることができる．
逆行列の計算の仕方は代数の教科書にあると思うが，その概略を説明すれば次のようになる．行列
の要素と単位行列 の要素
を横に並べた 5 行 10 列の行列を紙に間をあけて書く． が単位行列 になるように，次の３つの演算を適宜行う．
１つの行のすべての要素に同じ 0 でない数または変数をかける．
１つの行を別の行と入替える．
１つの行のすべての要素に同じ 0 でない数または変数をかけ別の行に加える．
の場合にはまず第 2, 3, 4 行に をかけ，第 5 行に 1 ~ をかける．次に第 1 列が の第 1 列になるように第 1 行に 1 をかけ第 2
行に加える，などする．第 2 行に 1 ，第 3 行に 2 ，第 4 行に 3 をかけ第 5 行に加える．この時点ではじめ
の部分は になり，
の部分に逆行列 ;1 が得られる．古いものを / で消し新しいものを余白に書入れる．
28 上三角行列または下三角行列の行列式の値は対角要素の積になるという性質がある．また不完全な三角行列の場合にも同様の性
質が利用できる．
N
N
I
u
N
u;
=
u
I
u
=
= =0
jA ; I j = (u ; )
0
u;
3=
2
2
2
u;
1 c
2 c
3 c
次にこの行列式の第 1, 2, 3 行に 1=(u ; ) をかけ，その後で第 1 行に ;
i
i
i
i
i
0
0
0
i
u;
i
0
0
i
i
i
i
i
i
i
4 行に加え，上三角行列のものにする．
jA ; I j = (u ; )3 ;(u ; )2 ; c2 = 0
れらを第
i
i
i
i
i
i
N
I
I
N
u
I
i1
i2
i
i
i
i
i1
c2 ，第 2 行に ;
i2
c2 ，第 3 行に ;
i3
c2 をかけこ
35
線 形 化 と 対 角 化 固有値
i = ui
li (Ai
;
c に対しては
0 c
BB 0
B
lik5 )B
BB 0
B@ 0
i I ) = (lik lik
1
2
i1
i2
0
0
c
c
0
i1
0
0
c
2
i2
1
C
i1 = C
CC
i2 = C
CC = 0
i3 = A
0
i3
0
0
c
i3 c2
c
2
c
となる．これらの式から 1 次独立の 5 個の左固有ベクトルが得られる29 ．
i の固有値の並べ方は，任意性を持つが，ここでは Li したがって Ri
の対角要素の大きさが小さくならないようにしている．固有値 ui に属する 3 個の固有ベクトルの選び方も
ここで多少補足すれば，対角行列
任意性を持つが，ここでは常微分方程式が流跡線に沿って伝播する１つのエントロピー波，２つの剪断波の
式になるようにしている．なお固有値 ui
c に属する固有ベクトルの方向は一意的に決まり，その常微分方
程式は２つの圧力波の式になる．固有値と固有ベクトルは対応するところに配さなければならない．行列
Ai の固有値も i である．これは固有値 k が物理的に波の位相速度を表していることを考えれば当然のこ
とである．Ai の左固有ベクトルは Li = Li N ，右固有ベクトルは Ri = N ;1 Ri で，これらの行列は次のよ
うに表される．
0
0
BB u0i u +0 c 00
0
i i1
BB
0
ui + i2 c
0
i=B
BB 0
0
0
ui + i3 c
@0
0
0
BB 10
B
Li = B
BB 0
B@ 0
0
0
0
0
1
0
0
0
0
1
0
0
0
0
1
i1
i2
i3
0
1
CC
CC
CC
CA
(17.74a)
ui c
;
;
(17.74b)
1= c
;
= u1 に対する意味のある式は l1 1 + c2 l1 5 = 0 l1 2 = 0 の２つである．u1 は 3 重根であ
るから，1 次独立の左固有ベクトル 3 個が存在する可能性がある．その１つは第 1 式を満足する
l11 = (1 0 0 0 ; 1=c2 )
29
i = 1 の場合を考える．固有値
1
1=c2
C
i1 = c C
CC
i2 = c C
CC N
i3 = c A
0
0
0
0
1
k
k
k
である．残りは
l13 = (0 0 1 0 0) l14 = (0 0 0 1 0)
n
n
である．一般に固有ベクトルは大きさと向きは任意で方向のみが特定される．また行列の次元が縮退し固有値が 重根の場合に 個のベクト
ルの方向は一義的には定まらない．上記のベクトルは左固有ベクトルの行列の対角要素がゼロにならないように，また最もシンプルになるよ
1
1
2
2 5 =0
3
4
2
5
うに選んだものである．固有値 1 = 1
に対しては
1k = 0
1k
1k+
1k
1k = 0
1k = 0 1k
1k = 0
5 =0
の５つの式が得られるが独立のものは４つである．これより 1k1 = 1k3 = 1k4 = 0，また 1
に対してそれぞれ 1k2
1k
が得られ，左固有ベクトルは次のようになる．
u c
cl
l12 = (0 1 0 0 1= c) l15 = (0 1 0 0 ; 1= c)
i = 2 3 の場合も同様になる．
l
l
l
cl
l
cl
u
cl
c
cl
l
l
cl
cl
36
圧縮性流れの解法|NS 方程式 0
BB 10
B
Ri = N ;1 B
BB 0
B@ 0
i1 =2c
1 ; i 1 =2
i2 =2c
0
i3 =2c
0
0
1 i2=2
0
1 i3=2
i1 c=2
i2 c=2
i3 c=2
0
0
0
;
;
1
=2c
C
i1 =2 C
CC
i2 =2 C
CC
i3 =2 A
c=2
;
(17.74c)
;
非保存形の式 (17.70) と等価な常微分方程式系は，この式の左から左固有ベクトル Li を乗じることによっ
て得られるが，これは保存形の式 (17.67) の左から左固有ベクトル Li を乗じたものでもある．
@Fi
@q
@q
Li @q
@t + @xi = Li @t + iLi @xi = 0
(17.75)
この常微分方程式系は具体的に書けば次のようになる．
; @ +u @
@t i @xi
n@
;
1 ; @ +u @ p = 0
c2 @t i @x
i
@o
ij n @
o
@
@t +(ui + ij c) @xi uj + c @t +(ui + c) @xi p = 0
n@
@ ou 1 n @ +(u c) @ op = 0
+(
u
c
)
i
@t
@x i c @t i @x
;
;
i
;
(j = 1 2 3)
(17.76)
i
これらの式は特性曲線 dxi =dt = ui
dxi =dt = ui +c dxi =dt = ui c に沿う内微分 @=@t+ui@=@xi @=@t+(ui +
c)@=@xi または @=@t+(ui c)@=@xi のいずれか１つを含むので常微分方程式で，上から dxi =dt = ui に沿って
伝播するエントロピー波，dxi =dt = ui に沿う剪断波または dxi =dt = ui + c に沿う圧力波，dxi =dt = ui c
;
;
;
に沿う圧力波を表している．
17.3.2
一般曲線座標系のオイラー方程式
一般曲線座標系の Euler 方程式の場合にも，同様に線形化，対角化することができる．
@ q^ + @ F^i = 0
@t @ i
0 1
BB u CC
BB 1CC
q^ = J B
BB u2CCC
@ u3 A
e
0
1
Ui
BB u U + pCC
BB 1 i i 1 CC
^
Fi = J B
BB u2Ui + i 2pCCC
@ u3Ui + i 3pA
流束 F^i は q^ の成分だけでなく測度
(i = 1 2 3)
(i = 1 2 3)
(17.77a)
(17.77b)
HUi
i j の関数であるが，時空間内の固定された点を考えれば ，F^i は q^ の成
分だけの 1 次の同次式で Euler の同次関係から
^
dF^i = @@Fq^i dq^ = A^i dq^ F^i = A^i q^ dA^i q^ = 0
(17.78)
37
線 形 化 と 対 角 化 となり30 ，式 (17.77) は次のように線形化される．
@ q^ + A^ @ q^ = 0
@t i @ i
^i
ただしヤコビ行列 A
(i = 1 2 3)
= @ F^i =@ q^ は
0
BB u U 0+
BB 1 i i 1
^
Ai = B
BB u2Ui + i 2
@ u3Ui + i 3
i1
^
Di1
2
;
2
;
2
;
HUi + 2 Ui
;
i2
i 1 u2 ; ~ i 2 u1
i 1 u3 ; ~ i 3 u1
i 1 H ; ~Ui u1
D^ ij = Ui +(1 ~) i j uj である31 ．
= Ui @=@
非保存形方程式 (17.66) は，u
i 2 u1 ; ~ i 1 u2
D^ i2
i 2 u3 ; ~ i 3 u2
i 2 H ; ~Ui u2
(17.79a)
1
C
i 3 u1 ~ i 1 u3 ~ i 1 C
CC
i 3 u2 ~ i 2 u3 ~ i 2 C
C
D^ i3
~ i 3C
A
i3
0
;
;
i 3 H ; ~Ui u3
(17.79b)
Ui
;
r
i
=
r
r
i @=@ i であるから簡単に一般曲線座標系の式に書換
えることができ，その１つの空間次元の式は次のようになる．
@ q^ + A^ @ q^ = 0
@t i @ i
(17.80a)
^i は次のようになる．
ただし未知変数ベクトル q^ と行列 A
0
BBU0i
B
A^i = B
BB 0
B@ 0
q^ = q
0
i1
Ui
0
0
i1
i2
0
0
0
Ui
c2
0
i2
c2
Ui
i3
1
C
i 1= C
CC
=
i2 C
CC
i 3= A
0
i3
c2
(17.80b)
Ui
式 (17.80a) は，
dq^ = @@q^q^ dq^ = NJ dq^
であるから，保存形方程式 (17.79a) の左から N=J を演算したものである．ただし
N は式 (17.72) で定義さ
れる非保存形への変換の行列である．
^i は式が簡単で容易に対角化することができ，これをもとに A^i も対角化することができる．
行列 A
A^i = N ;1 A^i N = N ;1 R^i ^i L^ i N = R^i ^i L^ i
A^
の要素は， ^i の成分を次のよう
F
q^ の成分で表してから，q^ の成分で微分して求める．
J U = (J u )
J ( u1 U + 1 p) = (J uJ1 )(J u ) + 1 ~Je ; 2J1 ~ (J u)2
(J u ) n Je ; ~ (J u )2o
J
2J
2
なお (1=2J )@ (J u ) =@J u1 = u1 .
31 例えば D
^ 1 = U +(1 ; ~ ) 1 u1 で，ここにはアインシュタインの総和側は適用されないことに注意．
30
i
i
i j
i
i j
i
i j
i
j
j
i
i
i
i
i
(17.81)
38
圧縮性流れの解法|NS 方程式 A^i の固有値と固有ベクトルは，デカルト座標系の Ai の固有値，固有ベクトルとほぼ同様の計算で求める
ことができるので結果のみ示すことにする32 ．
0
Ui
0
0
0
p
B
11
B
0
0
B 0 Ui + 1i g c
p
^i = B
22
B
0
0
Ui + 2i g c
0
B
p
B
0
0
Ui + 3i g33 c
@0
0
01
BB0
B
L^ i = B
BB0
B@0
0
0
0
0
0
ii
1i 1 2 ; 2i 2 1
1i 1 3 ; 3i 3 1
2i 2 1 ; 1i 1 2
ii
2i 2 3 ; 3i 3 2
3i 3 1 ; 1i 1 3
3i 3 2 ; 2i 2 3
i1
i2
ii
i3
0
BB1
BB
BB0
BB
R^i = N ;1 B
BB0
BB
BB0
B@
0
p1i
2 g11 c
0
0
p2i
2 g22 c
3
i1 i2
2 2i 1 gii i i
D i1
;
3
i2 i1
1i 1
2
gii i i
3
i3 i1
2 1i 1 gii i i
D i2
;
;
1i c
p
2 g11
2i c
p
2 g22
2
= 23 ji 1 giii j +(1 ji) 1
ii
ii
式 (17.80) と等価な常微分方程式系は次のようになる．
L^ i @ q^ + @ F^i = L^ @ q^ + ^ L^ @ q^ = 0
D ij
;
J @t @
;
i
i
i i
@t
@
(17.82a)
;
;
;
p3i
2 g33 c
3
i1 i3
2 3i 1 gii i i
3
i2 i3
3i 1
2
gii i i
;
;
3
i3 i2
2 2i 1 gii i i
;
1
CC
CC
CC
C
pgii cA
Ui
1
1=c2
p 11 C
1i g = cC
pg22= cCCCN
2i
p 33 CC
3i g = cA
pgii = c
0
0
0
0
D i3
3i c
p
2 g33
;
1
p
2 gii c C
C
i1 C
C
2gii C
CC
i2 C
C
2gii C
C
i3 C
C
2gii C
CC
c A
p
2 gii
(17.82b)
(17.82c)
;
(j = 1 2 3)
(17.83)
i
あるいは
@ +U @
1 @ +U @ p = 0
i
@t @ i
c2 @t i @ i
n @ ; p ii @ o pgii n @ ; p ii
@t + Ui + g c @ i Ui + c @t + Ui + g c
@
@
@
@
i i @t + Ui @ uj i j @t + Ui @ ui = 0
i
i
n @ ; p ii @ o pgii n @ ; p ii
@t + Ui g c @ i Ui
c @t Ui g c
;
@ op = 0
@i
(j = i)
;
;
;
;
;
32 特性条件は
6
@ op = 0
@i
jA^ ; ^ I j = (U ; ^ )3 f(U ; ^ )2 ; g c2 g = 0
p
p
となり，固有値は ^ 1 = U (3 重根) ^ 2 = U + g c ^ 3 = U ; g c となる．
i
i
i
i
i
i
ii
i
i
ii
i
ii
(j = i)
(17.84)
39
線 形 化 と 対 角 化 p
p
i =dt = Ui d i =dt = Ui + g ii c d i =dt = Ui ; g ii c に沿う内微分
p
p
i @=@t +(Ui + g ii c) @=@ i @=@t +(Ui ; g ii c) @=@ i のいずれか１つのみを含む常微分方程式
空間内の特性曲線 d
これらの式は
@=@t + Ui@=@
で，写像空間の特性曲線に沿って伝播するエントロピー波，圧力波または剪断波を表している．
17.3.3
一般曲線座標系の質量流束のオイラー方程式
一般曲線座標系の質量流束の運動方程式から構成される Euler 方程式の場合にも，流束を線形化し，その
ヤコビ行列を対角化することができる．
@ q~ + @ F~i = 0
@t @ i
1
0
Ui
BB U U + g1ipCC
BB 1 i 2i CC
~
Fi = J B
BB U2Ui + g3ipCCC
@ U3Ui + g pA
0 1
BB U CC
BB 1CC
q~ = J B
BB U2CCC
@ U3A
e
(i = 1 2 3)
(17.85a)
(i = 1 2 3)
(17.85b)
HUi
F~i は q~ の成分の 1 次の同次式であるから
~
dF~i = @@Fq~i dq~ = A~i dq~ F~i = A~i q~ dA~i q~ = 0
(17.86)
となり33 ，式 (17.85) は次のように線形化することができる．
@ q~ + A~ @ q~ = 0
@t i @ i
~i = @ F~i =@ q~ はヤコビ行列で
ただし A
0
0
BB
BB U1Ui + g1i
A~i = B
BB U2Ui + g2i
BB U3Ui + g3i
@
;
;
;
1i
2
2
2
HUi + 2 Ui
;
33
F~
i
(17.87a)
D~ i1
2i
1i U2 ~g
1
3i
1i U3 ~g
1
1i H ~ 1 Ui
;
;
;
2i
3i
U ~g1i 2
D~ i2
3i
2i U3 ~ g
2
2i H ~ 2 Ui
U ~g1i 3
2i
3i U2 ~ g
3
D~ i3
3i H ~ 3 Ui
2i 1 ;
3i 1 ;
;
;
;
;
1
C
~g1i C
CC
CC
~g2i C
C
~g3i C
A
0
Ui
q~ の成分で表せば
n
o
J ( U1 U + g1 p) = (J U1J)(J U ) + g1 ~Je ; 2J~ (J u)2
n
o
J HU = JJU Je ; 2~J (J u )2
の成分を
i
i
i
i
i
i
ただし
(J u)2 = g11 (J U1 )2 + g22 (J U2 )2 + g33 (J U3 )2 +2g23 (J U2 )(J U3 )+2g31 (J U3 )(J U1 )+2g12 (J U1 )(J U2 )
また
@xx @xx U = @xx u
1 @ (J u )2
2J @J U1 = g11 U1 + g12 U2 + g13 U3 = @ 1 @
@1
j
j
1
(17.87b)
40
圧縮性流れの解法|NS 方程式 D~ ij = Ui + ji Uj ~gji
i = (@xk=@ i )uk = gij Uj である．
j
;
式 (17.87) に対応する非保存形方程式の未知変数ベクトル q~ は，下記の式 (17.89b) のように置くのが適
~
当であろう．このとき保存形から非保存形への変換の行列 N
なる34 ．
0
BB U1 =
B 1
N~ = B
BB U2=
B@ U3=
0
1=
0
0
~
;
;
;
2
;
0
BB U1
BB 1
;
1
~
N =B
BB U2
@ U3
1
0
u =2
;
0
0
0
0
2
0
0
1=
0
~
1
3
;
0
0
0
0
1
0
0
0
1=
~
2
2
0
0
0
0
1=~
3
= J@ q~ =@ q~ とその逆行列 N~ ;1 は次のように
0
C
0C
CC
0C
C
0C
A
~
(17.88a)
1
CC
CC
CC
CA
(17.88b)
~ を演算したもので次のようになる．
この質量流束の非保存形 Euler 方程式は，式 (17.87) の左から N=J
@ q~ + A~ @ q~ = 0
@t i @ i
0 1
0
BBU CC
BBU0i
B 1CC
B
q~ = B
BBU2CC A~i = BBB 0
B@U3CA
B@ 0
p
q~
1i
2i
0
0
Ui
Ui
0
1i
0
c
2
2i
0
3i
0
0
0
c
2
1
C
g1i= C
CC
2
i
g =C
C
g3i= C
A
(17.89a)
Ui
2
3i c
(17.89b)
Ui
~i
はこの場合の未知変数ベクトル，行列 A
= N~ A~i N~ ;1 である35 ．
~i は簡単で容易に対角化でき，これより A~i も対角化できる．
行列 A
A~i = N~ ;1 A~i N~ = N~ ;1 R~i ~i L~ i N~ = R~i ~i L~ i
q~ の成分で表せば
= J =J U1 = J U1 =J p = ~ Je=J ; ~ (J u)2=2JJ
(J u)2 については前注参照．逆行列は，まず j N j I j の第 2, 3, 4 行に
34
q~
の成分を
容易に求めることができる．
j
U = @@xx uU = u2
j
2
2
; ~ + U = u2
j
j
j
35
(17.90)
~H = c2 +
2
j
g =U
ji
i
j
を掛け，第 5 行を
~ で割り，それから次の関係を用いれば
41
形 陰 解 法 A~i の特性条件と固有値は上記の A^i のものと同じになり， ~i ，L~ i ，R~i は次のようになる．
~i = ^i
0
BB1
BB0
B
L~ i = B
BB0
BB0
@
(17.91a)
1
p 11 CC
12 22
13 33
1
2i g =g
3i g =g
1i g = cC
CC
p
21
11
23
33
C~
22
1
1i g =g
3i g =g
2i g = cCN
pg33= cCCC
31=g 11
32=g 22
g
g
1
1i
2i
3i
p ii A
0
g =c
1i
2i
3i
0
p 1
p
p
p
=2 gii c
1 1i =2 g11 c 2i =2 g22 c 3i =2 g33 c
C
BB
12
22
13
33
CC
BB0 1 1i=2
g1i=2gii C
2i g =2g
3i g =2g
C
B
23
33
R~i = N~ ;1 B
g2i=2gii C
3i g =2g
CC
BB0 1i g21=2g11 1 2i=2
BB0 1i g31=2g11 2ig32=2g22 1 3i=2
3
i
ii
C
g =2g C
@
p ii A
p 22
p 33
p 11
0
0
1=c2
0
;
;
;
;
(17.91b)
;
;
;
;
;
;
(17.91c)
;
;
0
c=2 g
1i
2i
c=2 g
3i
c=2 g
;
c=2 g
~ i を演算することに
非保存形の式 (17.89) と等価な常微分方程式系は，この式の左から左固有ベクトル L
~ i =J を演算したものでもある．
よって求められるが，この常微分方程式系は保存形の式 (17.85) の左から L
L~ i @ q~ + @ F~i = L~ @ q~ + ~ L~ @ q~ = 0
i i@
i @t
J @t @ i
i
(17.92)
なおこの常微分方程式系の個々の式は，当然のことながら式 (17.84) すなわち写像空間の特性曲線に沿って
伝播するエントロピー波，圧力波，剪断波の式になる．
式
(17.77)
- 式 (17.80) L^ - 式 (17.83)
N=J
i
B
~
?
N=J
- 式 (17.89) L~ - 式 (17.92)
式 (17.85)
=
i
17.4
形陰解法
圧縮性流れの計算には
形陰解法 (delta-form implicit method) が広く用いられている．その理由はいくつ
かあるが， 形にすることにより対流項が精度に無関係に 1 次上流差分で近似でき計算が安定になること，ま
た
形では近似因子法などの適用が可能になり計算量を大幅に削減できることである．圧縮性流れのクーラン
数は各座標方向に伝播するエントロピー波，圧力波，剪断波などに対して定義され，CFL = maxj
または max ^ i t= i
j
j
i
t= xi
j
のようになる．圧縮性流れでは一般にクーラン数の制約が非常にきびしいが，それ
は衝撃波の捕獲のために細かい格子を用いることや，速度
i が壁近くでも音速の大きさであることによる．
形陰解法に後述の近似因子法を組み合わせればクーラン数はかなり大きく取ることができ，計算量は更に
削減されることになる．
デカルト座標系における Navier-Stokes 方程式 (17.51) の
形陰解法の式は，定常流れと非定常流れに対
42
圧縮性流れの解法|NS 方程式 しそれぞれ次のようになる．
n
o
I + t @x@ (Ani + i ) qn = rhsn
i
i
rhs = t @F
@x + D + g
(17.93)
D
;
qn+1 = qn +
i
n
q
o
;
I + t @x@ (Ani + i ) q(m) = (q(m;1) qn)+ 21 rhsn + rhs(m;1)
i
q(m) = q(m;1) + q(m)
n
D
;
(17.94)
;
また一般曲線座標系における Navier-Stokes 方程式 (17.65) の
形陰解法の式は，定常流れと非定常流れに
対しそれぞれ次のようになる．
n
o
n
I + t @@ (A~ni + ~ i ) q~n = rg
hs
i
^
rg
hs = t B @@Fi + D^ +^g
;
n
(17.95)
D
q~n+1 = q~n + q~n
i
o
; hsn + rg
(m;1)
I + t @@ (A~ni + ~ i ) q~(m) = (~q (m;1) q~n )+ 12 rg
hs
i
(m)
(m;1)
q~ = q~
+ q~(m)
ただし
D
@(
i
D
q)=@xi @ ( ~i q~)=@
D
;
(17.96)
;
i は拡散項の主要部，また左辺の演算子 @=@xi は
q に，また @=@
i は q~
にも作用していることに注意すべきである．
上式に関しては 16.3.2 項に述べた１次元流れの場合とほぼ同じことが言える．定常流れの式 (17.93)(17.95)
は，反復法によって解かれる．その右辺は式 (17.51)(17.65) の残差，
qn
q~ n はその修正値を意味し，解
が定常流れに漸近すれば，残差は零に近づき，同時に修正値も零に近づく．したがって左辺の演算子は，解
の精度には無関係で，安定性と計算の容易さから
1 次上流差分で近似され，拡散項が適当に省略され，ま
た近似因子法が適用される．他方 右辺は解の精度に直接関係するので正確に計算しなければならない．非
定常流れの式 (17.94)(17.96) も 非定常流れの各時間ステップごとに 反復計算によって解かれる．この式は
t における修正値 qn を (q(m;1) qn)+ q(m) のように分け，その既知の部分 q(m;1) qn を陽
的に，未知の部分 q (m) を陰的に求めるものである．曲線座標格子の q~ n についても同様である．この時
時間間隔
;
;
間ステップごとの反復計算は数回で十分収束する．
これらの
形陰解法の式は，次節に述べる流束差分分離法または
17.5 節に述べる有限体積法を適用し
て解かれる．これらの節には 3 次の Chakravarthy-Osher 型スキームとそれを高次にしたものについて述べ
る．衝撃波の捕獲という観点からは 2 次ないし
3 次スキームで十分である．他方 滑り面や境界面などのい
わゆる弱い不連続の鮮明な捕獲には高次スキームの使用が現状では不可欠である．弱い不連続は上流の不
連続が下流に伝播または過去の不連続が未来に伝播するものでいったん拡散してしまうと修復不能であり，
格子を細かくしても数値的拡散を減らすことは難しく，スキームを高次化して流れの急激な変化に柔軟に対
応させることが必要になる．
流 束 差 分 分 離 法 17.5
43
流束差分分離法
本節では第 16 章に述べた１次元流れの流束差分分離法を多次元に拡張する．まず
Chakravarthy-Osher
型 TVD スキームの流束差分分離法の式を示し，この式中にあらわれる分離流束差分の計算式をデカルト座
標格子と曲線座標格子に対し導出する．また流束差分分離法の高次化についても述べる．以下の記述に疑
問が生じたときには 16 章を参照してください．
17.5.1
Chakravarthy-Osher 型 TVD スキーム
一般に保存形差分スキームは次のように書くことができる．
@F = 1 ;H
(17.97)
@x l
x l+1=2 Hl;1=2
ここに Hl+1=2 は数値流束で，3 次の Chakravarthy-Osher 型 TVD スキームでは次のようになる．
;
1 ~+
1 ~+
1 ~;
1 ~;
Hl+1=2 = Hl(1)
+1=2 + 6 F l;1=2 + 3 Fl+1=2 6 Fl+3=2 3 F l+1=2
;
F~j+1=2 = minmod Fj+1=2 b Fj;1=2
(j = l l +1)
;
;
;
F~ j+1=2 = minmod Fj+1=2 b Fj+3=2
1; +
1;
;
Hl(1)
+1=2 = 2 Fl + Fl+1
2 Fl+1=2 Fl+1=2
;
(j = l 1 l )
(17.98a)
(17.98b)
;
(17.98c)
;
minmod は minmod(x y) = sign(x) max 0 min x sign(x)y ] のように定義される制限関数であ
る．式 (17.98) の数値流束 H は TVD 安定な 1 次上流差分の H (1) とそれを 3 次にする補正項からなり，補
正項の大きさは解が TVD 安定になるように minmod 関数によって制限されている．また 式 (17.101a) で
は，流束ベクトル分離法 ( ux-vector splitting) により，流束 F が x の正または負の向きに伝播する波の
ただし
fj
j
g
F + と F ; に分離され ，Fl+1=2 の値が上流側に片寄った点の F の値すなわち正方向の波に対しては
Fl+;1 Fl+ Fl++1 の値，また負方向の波に対しては Fl; Fl;+1 Fl;+2 の値をもとに補間式のようなものによっ
流束
て求められると解釈できる．
F
F = (R L q) と置けば，周知のように計算結果が不自然に振舞
うことがある．それは例えば位相速度 u c の圧力波は，音速点でその符号が変わるため，ここで式 (17.101a)
内の補間式が壊れてしまうからである．あるいは翼面上に発達する境界層中の横断方向の流速 v の符号の
式 (17.98) の流束差分
は単純に
;
変わるところでも同様の不自然な振舞いの起きることがある．ここではこのような振舞いが起きないよう
に，次のいずれかによって線形化対角化を行うことにする．
Fj+1=2 = (R L)l+1=2 qj+1=2
Fj+1=2 = Rj+1=2 sign l+1=2 ( L q )j+1=2
j
第 1 式は，A
j
(17.99a)
(17.99b)
= R L を l +1=2 の値で代表しようというものである．この近似により精度は形式的に 2 次
に落ちるが実質的にはあまり変わらない．それは基の方程式が準線形 (quasilinear) で A がその係数である
ことによる．第 2 式は精度が落ちないように l+1=2 の符号のみを用いるものである．これら両者の計算結
果は通常ほとんど 違わない．
44
圧縮性流れの解法|NS 方程式 一般曲線座標系では保存形差分スキームは次のように書かれる．
(@ F=@x)l = H^ l+1=2 H^ l;1=2
(17.100)
;
F F^ ，H^ l+1=2 は数値流束で 3 次の Chakravarthy-Osher 型 TVD スキームでは次のようになる．
1 ~+
1 ~+
1 ~;
1 ~;
H^ l+1=2 = H^ l(1)
(17.101a)
+1=2 + 6 Fl;1=2 + 3 Fl+1=2 6 Fl+3=2 3 Fl+1=2
;
F~ j+1=2 = minmod Fj+1=2 b Fj;1=2
(j = l l +1)
(17.101b)
;
F~j+1=2 = minmod Fj+1=2 b Fj+3=2
(j = l 1 l)
1 ; F+
1;
;
(17.101c)
H^ l(1)
+1=2 = 2 Fl + Fl+1
2 l+1=2 Fl+1=2
ただし
;
;
;
;
17.5.2
;
流束差分の計算
始めに波の伝播方向によって分離された流束差分，式 (17.99a) の値を求める計算式を導出する．非保存
形の変数
q = Nl+1=2 qj+2=1 の差分は次のようになる．
0
BB u1 =
BB 1
qj+2=1 = B
BB u2=
@ u3=
0
1=
0
0
~u1
;
;
;
2
0
0
1=
0
~u2
;
1 0
0
0
0
1=
~u3
;
1
0
1
0
BB M = CC
C BB ( u )CC
0C
BB 1 CC
CC BB 1 CC
0C
CC BBB ( u2)CCC = BBB M2= CCC
0 A @ ( u3 )A
@ M3= A
P
~ l+1=2
e j+1=2
;
(17.102)
ただし
Mi = ui
1i
;
P=
ここでは差分は
Ai = Ri
~u1
2
;
3i
~u2
;
0 l+1=2 qj+1=2
~u3 ~ l+1=2 qj+1=2
;
j +1=2 のもの，その他の変数は l +1=2 のものを取ることにする．また非保存形の行列
i Li は次のようになる．
0
BB 0i1
B
Ai = B
BB 0
B@ 0
0
ただし
2i
i1 = ui
1i
ia
=c
i1 + 1i ib
0
0
1i
ia
i2 = ui+c
c
2i
ia
=c
0
i1 + 2i ib
0
2i
ia
i3 = ui;c
c
3i
=c
ia
0
0
i1 + 3i ib
3i
ij = ( ij
c
ia
j
である．この式はまず A1 を計算し ，次にこれより A2
1
ib =c C
1i ia = cC
CC
2i ia = cC
CC
3i ia = cA
i1 + ib
2
ij j)=2
A3
ia = ( i2 ; i3 )=2
(17.103)
ib = ( i2 + i3 )=2; i1
を推定し，これらを Kronecker のデルタ関
流 束 差 分 分 離 法 45
数を用いて重ね合わせれば容易に導くことができる．これらの式から
0
BB
BB
BB
B
qj+1=2 = B
BB
BB
BB
B@
(Ai )l+1=2
更に
1 M+ 1 P 1
+
i
ia c
i1
ib c2
CC
1 M + 1i P + 1i M C
C
1
ia c
1C
i1
ib
CC
1 M + 2i P + 2i M C
C
i1
ia
2
ib
c
1 M + 3i P + 3i
3
ia c
ib
i1 P + ia c Mi + ib P
i1
CC
C
M3 C
CC
A
2
( Fi )j+1=2 = (N ;1 Ai )l+1=2 qj+1=2 を計算すれば ，デカルト座標格子の場合の分離流束差分の計
算式が次のように求められる36 37 ．
( Fi )j+1=2 = (Ri i Li )l+1=2 qj+1=2
n
o n
o
= i1 q + ib Mi + ia cP qai + ia Mi + ib cP qb
(i = 1 2 3 j = l l 1)
(17.104a)
ただし
i1 = ui
i2 = ui + c
ij = ( ij
ij j)=2
j
Mi = ui
1i
;
P=
ia = ( i2 ; i3 )=2
2i
~u1
2
i3 = ui ; c
3i
~u2
;
;
ib = ( i2 + i3 )=2 ; i1
0 l+1=2 qj+1=2
~u3 ~ l+1=2 qj+1=2
0 1
0 1
0 1
0
1
B
C
B
C
B
B
BB 1i CC
BBu1CCC
u1 C
B
C
C qai = BB 2i CC qb = 1 BBu2CC
q=B
B
u
2C
B
C
BB CC
cB
B
B@u3CCA
@ u3CA
@ 3i A
e
(17.104b)
;
ui
(i = 1 2 3)
H
36 この計算には次の関係が用いられた．
~ = ;1
2
= ~ u2
2
u2
2
c2 + u2 = H
~ 2
+ u M + 1~ P = e
j
j
(17.98) の流束差分は i = 1 の場合に次のようになる．
~ )q 1 + ( + M~ 1 + + P=c
~ )q
( F~ +1 ) ;1 2 = +11 q~ + ( +1 M~ 1 + +1 P=c
1
1
M~ 1 = ;u ~+ ~u
P~ = 2 ~ ; ~(u ~u + v ~v + w ~w) + ~ e~
~u ;1 2 = minmod; ( u) ;1 2 b ( u) +1 2
37 蛇足を加えれば，式
l
=
l
=
a
a
b
l
=
l
a
=
b
b
46
圧縮性流れの解法|NS 方程式 分離流束差分の式 (17.99b) の計算式は，この式に対しては
Mi =
ui
P= p N q
q
であ
るから次のようになる．
( Fi )j+1=2 = (Ri )j+1=2 sign i l+1=2 ( i Li q )j+1=2
= sign i1 l+1=2 i1 q
uiqai ( p=c)qb j+1=2
h 1
i
+ sign i2 l+1=2 i2 2 ( ui + p=c)qai +( ui + p=c)qb
ij+1=2
h 1
+ sign i3 l+1=2 i3 2 ( ui p=c)qai ( ui p=c)qb
j
j
j
j
j
j
j
j
;
;
;
;
;
j +1=2
(17.105)
次に曲線座標格子の場合の分離流束差分の式 (17.99a) の値を求める計算式を導出する．なおここでは ^ の
= F q^ = q A^ = A ^ = k．非保存形の
変数 q はデカルト座標系のものと同じである． q = q ．また行列 Ai = Ri Ki Li は次のようになる．
1
0
1
i
1
i
2
i
3
k
k p
k p
k p
k
付いている文字の代わりに Fraktur 体文字を用いることにする．F^
BB i1
BB
BB 0
BB
Ai = B
BB 0
BB
BB 0
B@
0
これより
(Ai )l+1=2
CC
C
ki1 + kib giii1 kib i g1iii 2 kib i g1iii 3 kia p iii1 C
C
g cC
CC
2
kib i g2iii 1 ki1 + kib gi ii2 kib i g2iii 3 kia p iii2 C
g cC
CC
2
i
3 i2
i
3
i
3 C
i
3 i1
kib gii
ki1 + kib gii kia p ii C
kib gii
g cC
C
i
2 c
i
3 c
i
1 c
kia p ii
kia p ii
kia p ii
ki1 + kib A
ia
gii c2
ia
g
gii c
g
ia
ib c2
gii c
(17.106)
g
0
1 M +k 1 P 1
p
k
+
k
i ib c2
ia g ii c
CC
BB i1
C
BB 1
BBki1 M1 + kib giii 1 Mi + kia pgiii1 c P CCC
CC
BB 1
i
2
i
2
B
qj+1=2 = Bki1 M2 + kib gii Mi + kia p ii P C
g c C
C
BB
BBki1 1 M3 + kib i 3 Mi + kia p i 3 P CCC
BB
gii
gii c C
@ k P + k pc Mi + k P CA
i1
ia g ii
ib
ただし
Mi = i k Mk =
更に (
Ui
;
i1
i2
i 3 0 l+1=2
qj+1=2
Fi )j+1=2 = (JN ;1 Ai )l+1=2 qj+1=2 を計算すれば ，曲線座標格子の場合の分離流束差分の計算式
が次のように求められる．
( Fi )j+1=2 = (J Ai )l+1=2 qj+1=2
n
o
n
o
= J ki1 q + J kib p1 ii Mi + kia 1c P qai + J kia p1 ii Mi + kib 1c P qb
g
g
(i = 1 2 3 j = l l 1)
(17.107a)
47
流 束 差 分 分 離 法 ただし
p
p
ki1 = Ui ki2 = Ui + gii c ki3 = Ui gii c
kij = (kij kij )=2 kia = (ki2 ki3 )=2 kib = (ki2 + ki3 )=2 ki1
Mi = i k Mk = Ui i 1 i 2 i 3 0 l+1=2 qj+1=2
j
;
j
;
;
;
P=
2
~u3 ~ l+1=2 qj+1=2
0 1
0 1
0
1
~ u2
~u1
;
(17.107b)
;
;
0 1
B
C
BBu CC
BB CC
B
u1 C
1C
i 1C
B
C
B
1 B C q = q = 1B
C
B
C
p
q=B
q
=
B
C
B
B
C
u
u
ai
b
b
2
2C
ii B i 2 C
B
C
B
c
g
B
B@u3CCA
B@ i 3 CA
@ u3CA
e
(i = 1 2 3)
H
Ui
なおここでも差分は j +1=2 のもの，その他の変数は l +1=2 のものを取ることにする．この場合の分離流束
Fi )j+1=2 = (Ai )l+1=2 qj+1=2 ではなく式 (17.107a) のように置かれたのは，Ai の格子依存性
が大きいのに対し，J Ai は i 方向に格子依存性がほとんどなく，点 ( i )l+1=2 の値で代表させても誤差があ
差分が，(
まり生じないことによる．
分離流束差分の式 (17.99b) の計算式は，
Mi = i k uk
Ui
P = p であるから次のように
なる．
( Fi )j+1=2 = (Ri )j+1=2 signKi l+1=2 (J Ki Li q )j+1=2
oi
h n
= signki1 l+1=2 J ki1 q p ii Ui qai 1c p qb
j
j
j
h
1n
+ signki2 l+1=2 J ki2 2
h n
+ signki3 l+1=2 J ki3 12
j
j
j
j
式 (17.98) のスキームが
らない．b の値は，
j
;
g
;
j +1=2
oi
pgii Ui + 1c p qai + pgii Ui + 1c p qb j+1=2
oi
pgii Ui 1c p qai pgii Ui 1c p qb j+1=2
;
;
;
(17.108)
TVD 安定になるためには，b の値とクーラン数が次の条件を満足しなければな
1<b 4
によって制限され，その上限
j
j
t= x は
(1
;
)b1CFL 1
b = 4 に取れば制限関数の働く範囲は最も狭くなる．またクーラン数 CFL =
b1 = 1+ 16 + 31 b
によって制限される．完全陰解法 (
= 1) ではクーラン数 CFL は一応自由に選べることになる．
48
圧縮性流れの解法|NS 方程式 次に前節に示した
形陰解法の式の左辺の計算について述べる．長方形格子の場合には，近似因子化す
る前の左辺は１次上流差分を用いることにすれば次のようになる．
lhsijk = qijk +
x
+ y
+ z
t= x
x =
性速度 ( )ijk
ただし
，jA
;
;
;
(A+1 )ijk
(A+2 )ijk
(A+3 )ijk
= A+
j
;
qi;1 j k + A1 ijk qijk + (A;1 )ijk qi+1 j k
qi j;1 k + A2 ijk qijk + (A;2 )ijk qi j+1 k
(17.109)
qi j k;1 + A3 ijk qijk + (A;3 )ijk qi j k+1
A; である．この左辺の計算では，上流化は点 xi yj zk における特
j
j
j
j
j
j
= 1 2 3 の符号によって行われ る．特性速度の符号によって分離されたヤコビ 行列
Ai = N ;1Ai N の値は式 (17.103) の Ai と式 (17.72) の N N ;1 の積をコンピュータで計算することに
よって求められる．なお行列 Ai の式は次のようになりかなり複雑である．
0
i1 a1
BB
BB (a1u1 bc 1i)
B
Ai = B
BB (a1u2 bc 2i)
BB (a1u3 bc 3i)
@
;
;
;
;
;
;
;
ia + b ib
a12
ia +(b ; 1i )u1 ib
ia +(b ; 2i )u2 ib
(a12 u2 a13 2i ) ia a14u2 ib
(a12 u3 a13 3i ) ia a14u3 ib
(a12 H a13ui ) ia (a14 H 1i u1 ) ib
ia +(b ; 3i )u3 ib
;
a22
;
a32
2
(a32 u1
(a22 u1 a23 1i ) ia a24u1 ib
2
2
2
(a32 u2
i1 +(a2 u2 a3 2i ) ia (a4 u2 2i ) ib
3
(a22 u3 a23 3i ) ia a24u3 ib
i1 +(a2 u3
(a22 H a23 ui ) ia (a24 H 2i u2 ) ib (a32 H
;
;
;
;
;
;
;
;
;
a33 1i )
a23 2i )
a33 3i )
a33 ui )
;
;
;
3
ia ; a4 u1 ib
3
ia ; a4 u2 ib
3
;
;
;
ia ; a4 ib
;
;
;
;
ia ; a4 ib
;
1
1
1
1
i1 +(a2 u1 ; a3 1i ) ia ; (a4 u1 ; 1i ) ib
(a1 H bcui ) ia +(bH ui2 ) ib
;
ia ; a4 ib
3
ia ; (a4 u3 ; 3i ) ib
3
ia ; (a4 H ; 3i u3 ) ib
;
a6 ib
a5 1i ia + a6u1
a5 2i ia + a6u2
a5 3i ia + a6u3
i1+ a5 ui ia+ a6 H
1
CC
ib C
CC
C
ib C
CC
ib C
A
ib
a1 ui =c ak2 ki =c ak3 ~uk =c ak4 ~uk =c2 a5 ~=c a6 ~=c2
ただし b = ~u 2=2c2 ，また H = c2= ~ +u 2=2 は岐点エンタルピーである．
他方 前節の 形陰解法の式の両辺に行列 N を乗じ 形陰解法の非保存形の式を作れば その左辺は次の
ようになる．
qi;1 j k + A1 ijk qijk + (A1; )ijk qi+1 j k
qi j;1 k + A2 ijk qijk + (A2; )ijk qi j+1 k
qi j k;1 + A3 ijk qijk + (A3; )ijk qi j k+1
ただし Ai は式 (17.103) で定義されたもの，また q = N q は次のようになる．
(N lhs)ijk = qijk + x
+ y
+ z
0 1
BB u CC
BB 1CC
q =B
BB u2CCC
@ u3A
p
;
;
;
(A1+ )ijk
(A2+ )ijk
(A3+ )ijk
j
j
j
j
j
j
(17.110)
(17.111)
49
流 束 差 分 分 離 法 次に曲線座標格子の場合の
f = q~ijk + t
lhs
;
;
;
形陰解法の式の左辺を示せば 次のようになる．
q~i;1 j k + A~1 ijk q~ijk + (A~;1 )ijk q~i+1 j k
q~i j;1 k + A~2 ijk q~ijk + (A~;2 )ijk q~i j+1 k
q~i j k;1 + A~3 ijk q~ijk + (A~;3 )ijk q~i j k+1
(A~+1 )ijk
(A~+2 )ijk
(A~+3 )ijk
j
j
j
j
j
j
(17.112)
i j k における特性速度 ( ~ )ijk
= 1 2 3 の符号によって行われ る．
;
1
~
~
~
~
特性速度の符号により分離されたヤコビ 行列 Ai = N Ai N を求めるために，まず非保存形の行列
A~i = R~i ~ L~ i を式 (17.91a)(17.91c)(17.91b) から計算すれば次のようになる．
この左辺の計算では上流化は点
0~
BB i1
BB 0
BB
A~i = B
BB 0
BB
BB 0
@
1i
~ia =ci
~i1 + 1i ~ib
2i
~ g
1i ib ii
g
3i
~ g
1i ib ii
g
2
~
1i ia c =ci
2i
~ia =ci
2i
1i
~ib g ii
g
~i1 + 2i ~ib
3i
~ g
2i ib ii
3i
1
C
1i C
g
~ia C
C
ci C
2i C
g
~ia C
C
ci C
C
3i C
g
~ia C
ci C
A
~ia =ci
3i
~ib =c2
1i
~ib g ii
g
2i
g
~
3i ib ii
g
(17.113)
~i1 + 3i ~ib
~ 2
~i1 + ~ib
0
3i ia c =ci
p ii ~ ~ ~
ただし ~ i1 = Ui ~ i2 = Ui + ci ~ i3 = Ui ci ci
g c ij = ( ij ij )=2 ~ia = (~i2 ~i3 )=2 ~ib =
(~i2 + ~i3 )=2 ~i1 である．
~i = N~ ;1A~i N~ の値は式 (17.113) の A~i と式 (17.88)
特性速度の符号によって分離されたヤコビ行列 A
~ N~ ;1 の積をコンピュータで計算することによって求められる．なお行列 A~i の式は次のようになり
のN
g
2
~
2i ia c =ci
;
j
j
;
;
一層複雑である．
0
~i1 f1 ~ia + b~ib
BB
BB (f1U1 f2g1i)~ia +(bU1 f3g1i)~ib
B (f U f g2i)~ +(bU f g2i)~
A~i = B
BB 1 2 2 ia 2 3 ib
BB (f1U3 f2g3i)~ia +(bU3 f3g3i)~ib
@
;
;
;
;
;
;
;
;
;
;
f41 ~ia f61 ~ib
~i1 +(f41U1 f51g1i )~ia (f61U1 f71g1i )~ib
(f41 U2 f51g2i )~ia (f61 U2 f71g2i )~ib
(f41 U3 f51g3i )~ia (f61 U3 f71g3i )~ib
(f41 H f51Ui )~ia (f61 H f71Ui )~ib
;
;
;
;
;
;
;
;
;
;
(f1H f2 Ui )~ia +(bH f3Ui )~ib
1
f42 ~ia f62 ~ib
f43 ~ia f63 ~ib
f9 ~ib
C
(f42 U1 f52g1i )~ia (f62U1 f72g2i )~ib (f43 U1 f53g1i )~ia (f63U1 f73g1i )~ib f8 g1i ~ia + f9U1 ~ib C
CC
~i1+(f42U2 f52g2i )~ia (f62 U2 f72g2i )~ib (f43 U2 f53g2i )~ia (f63U2 f73g2i )~ib f8 g2i ~ia + f9U2 ~ib C
CC
C
(f42 U3 f52g3i )~ia (f62U3 f72g3i )~ib ~i1+(f43U3 f53g3i )~ia (f63U3 f73g3i )~ib f8 g3i ~ia + f9U3 ~ib C
CA
(f42 H f52 Ui )~ia (f62H f72 Ui )~ib
(f43 H f53 Ui )~ia (f63H f73 Ui )~ib ~i1+ f8Ui ~ia+ f9H ~ib
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
f1 Ui =ci f2 ci q=gii f3 Ui =gii f4k ki =ci
f6k ~ k =c2 f7k ki =gii f8 ~=ci f9 ~=c2
;
;
f5k ~ k =ci
50
圧縮性流れの解法|NS 方程式 b = ~u2=2c2 H = c2=~ +u2=2 である38 ．
~ を乗じた非保存形の式の左辺は次のようになる．
他方 この場合の 形陰解法の式の両辺に N=J
ただし
N~ lhs
f
J ijk = q~ijk +
x
+ y
+ z
ただし
A~i
;
;
;
(A~1+ )ijk
(A~2+ )ijk
(A~3+ )ijk
q~i;1 j k + A~1 ijk q~ijk + (A~1; )ijk q~i+1 j k
j
j
q~i j;1 k + A~2 ijk q~ijk + (A~2; )ijk q~i j+1 k
q~i j k;1 + A~3 ijk q~ijk + (A~3; )ijk q~i j k+1
は式 (17.113) で定義されたもの，また
0 1
BB U CC
BB 1CC
~
N
q~ = J q~ = B
BB U2CCC
@ U3A
q~
j
j
j
j
(17.114)
は次のようになる39 ．
(17.115)
p
17.5.3
流束差分分離法の高次コンパクト TVD スキーム
16.7.1 節に述べた高次コンパクト TVD スキームは多次元にも適用できる．ここにはその概要について
述べる．3 次 Chakravarthy-Osher 型 TVD スキームの数値流束の式 (17.98) は，その補正項に含まれるす
F を次の DF で置き換えることによって，高次化することができる．
べての
1 ~+
1 ~+
1 ~;
1 ~;
Hl+1=2 = Hl(1)
+1=2 + 6 D F l;1=2 + 3 D Fl+1=2 6 D Fl+3=2 3 D F l+1=2
;
DF~j+1=2 = minmod DFj+1=2 bDFj;1=2
(j = l l +1)
;
;
;
DF~ j+1=2 = minmod DFj+1=2 bDFj+3=2
(j = l 1 l )
(17.116a)
(17.116b)
;
38
g U
U = u2
@ @
@ @x @ u = @ u = U
g = @x
@ u @x @x = @x @
@x
@x
なお最後の式で (@ =@x )(@x =@ ) は k = l のときに 1，k 6= l のときに 0 になる．
i
ij
i
i
j
i
k
ij
k
i
i
l
k
i
j
i
k
j
l
l
l
i
l
k
j
j
i
39 なおここでの
は空間差分の演算子ではなく同一空間点における修正値を意味するものであるから，ヤコビアン
に対し定数として扱われこれを自由に出入りできる．
テンソルは
@x
1 2
U = g U U = @x
@ U @ U =u u = 2 u
p = ~e ; 2
p = ~ e; 2 ; 2
これらの式から q~ の最後の要素は次のように求められる．
2
;~
U +~ e
2
=
; ~u ( u )+ p + 2 + 2
= p +~(u 2 + u u ; u u )
= p
i
k
k
i
ij
j
k
i
k
j
j
i
i
i
k
k
k
k
k
k
i
k
k
J ，測度，測度
流 束 差 分 分 離 法 51
この式の minmod 制限関数のすべての DF を
DFj+1=2 = Fj+1=2 16 3Fj+1=2
;
4 次スキームになり，また
3 3F
DFl+1=2 = Fl+1=2 20
l+1=2
1
DFj+1=2 = Fj+1=2 5 3Fj+1=2
(j = l l 1)
(17.117)
のように与えれば
;
;
ように与えれば
(j = l 1)
(17.118)
5 次スキームになる．この 4 次スキームは 5 次に近いもので弱い不連続を鮮明に捕獲する
目的にも適っている．
この高次化に伴って，上記のスキームを TVD 安定化する２つの条件は修正を要することになり，また付
F も局所的不安定を招く虞があるので次のようにその大きさを制限する必要がある．
3
加項
Fj+1=2 = 2F~j+1 2F~ j
2F
~j+1 = minmod; 2Fj+1 b2 2Fj
~ = minmod; 2F b2 2F
2F
j
j
j +1
2F = F
Fk;1=2
k
k+1=2
3
;
(17.119)
;
ただし b2 は DF が極値にならないように 1 < b2
4 の範囲に取られる．3 次の TVD スキームでは 3 次の
補正項の大きさを，この補正項のために新たな F の極値が発生しないように，勾配 F の大きさを minmod
制限関数によって制限したが，この高次スキームでは高次の補正項の大きさを，この補正項のために新たな
F の極値 (F の変曲点) が発生しないように曲率 2F の大きさを制限することになる．
この高次スキームの TVD 安定条件は厳格には次のようになる．
0 < b 6C;1 2
;
(1
ただし
;
)b1CFL 1
b1 = 1+ 16 C(1+2b)
C = max(Dqk = qk ) である．制限補正勾配 DF に変曲点が現れないためには倍率 b2 の上限は次式
によって制限される．
1 < b2 4
ここで b2
= 4 に選べば ，C の大きさは 1=2 C 1+1=8 となり，これよりこの高次スキームが TVD 安定
になる条件は次のようになる．
0 < b 10=3
(1 )b1CFL 1
;
b1 = 39=16
(17.120a)
(17.120b)
52
圧縮性流れの解法|NS 方程式 17.6
有限体積法
本節では第 16 章に述べた１次元流れの有限体積法を多次元流れに拡張する．有限体積法で解を
る計算は３つの段階，すなわち再構成 (reconstruction)，発展 (evolution)，セル平均 (cell
t 進め
averaging) から
構成される．ここではまず区分的 2 次多項式による再構成段階，次に Roe の近似リーマン解による発展段階
について，それから有限体積法の高次化についても述べる．以下の記述に疑義が生じたときにはまず 16 章
を参照してください．
17.6.1
再構成段階
始めにデカルト座標系 xyz における有限体積法について述べる．その基礎微分方程式は Euler 方程式
@q + @F1 + @F2 + @F3 = 0
@t @x @y @z
である．流れ場を格子間隔 x y z の直方体格子セルに分割する．上式を ijk 番目の格子セルにわたって
積分すれば次式が得られる．
t
n+1 = q n ; t (H1 )
qijk
i+1=2 j k ; (H1 )i;1=2 j k ; y (H2 )i j +1=2 k ; (H2 )i j ;1=2 k
ijk
x
; zt (H3 )i j k+1=2 ; (H3 )i j k;1=2
(17.121)
ただし qijk は q (x y z ) の ijk 番目のセルにわたる平均値である．また (H1 )i+1=2 j k は数値流束 F1 のセル
境界面 xi+1=2 にわたる平均値 H1 (t qiL+1=2 j k qiR+1=2 j k ) の時間間隔 t における平均値として定義される．
qijk = x 1y z
Z xi+1=2 Z yj+1=2 Z zk+1=2
dx
xi;1=2
n+1
1 Zt
yj;1=2
dy
zk;1=2
q(x y z ) dz
(17.122a)
(H1 )i+1=2 j k = t
H1 (t qiL+1=2 j k qiR+1=2 j k ) dt
(17.122b)
tn
q(x y z ) は各セル内で連続な x y z の区分的多項式で，セル境界面上では一般に不連続である．添字 L R
はセル境界の左右の値を示す．
再構成段階では，セル平均値 qijk を基に区分的 2 次多項式が再構成され，そのセル境界値が求められる．
ijk 番目セルの境界 x = xi+1=2 = xi;1=2 における q のセル断面平均値 qiL+1=2 qiR;1=2 は次式で与えられる．
qiL+1=2 = qi + 16 minmod( qi;1=2 b qi+1=2 )+ 31 minmod( qi+1=2 b qi;1=2 )
qiR;1=2 = qi ; 61 minmod( qi+1=2 b qi;1=2 ) ; 31 minmod( qi;1=2 b qi+1=2 )
これらは，q の x 一定セル断面平均値に対する区分的 2 次多項式 qjk (x)
= q(x) を 16.5.4 項にならって再構
成しそれから求めたものである．q の y または z に関する区分的 2 次式のセル境界値も同様に求められる．
セル境界面の左右の qlL+1=2 ，qlR+1=2 は 3 次の有限体積法 TVD スキームでは次のようになる．
qlL+1=2 = ql + 61 q~l;1=2 + 13 q~l+1=2
qlR+1=2 = ql+1 ; 61 q~l+3=2 ; 31 q~l+1=2
(17.123a)
53
有 限 体 積 法 ;
q~j+1=2 = minmod qj+1=2 b qj;1=2
;
q~j+1=2 = minmod qj+1=2 b qj+3=2
(j = l l +1)
(17.123b)
(j = l ; 1 l)
区分的 2 次式のセル中心の値は minmod 関数が働かないときには
1 2q
q(xi ) = qi ; 24
i
(17.124)
となる．この値はセル断面 x = xi にわたる平均値で，q (yj )
境界 xi+1=2 における変数 q の跳躍量 (jump) は
q(zk ) とは一般に完全に一致しない．またセル
qi+1=2 = qiR+1=2 ; qiL+1=2 = ; 61 3 qi+1=2
となる．この跳躍量はもとの関数 q (x
(17.125)
y z ) が滑らかに変化するところでは通常無視できる大きさになる．
次に曲線座標格子セルの場合について述べる．この場合の Navier-Stokes 方程式 (17.64a) に相当の Euler
方程式は次のようになる40 ．
@ q^ + @ F^1 + @ F^2 + @ F^3 = 0
@t @
@
@
ここでは物理空間 xyz 内の曲線座標格子セルが写像空間
内の単位立方体格子セル (
に対応付けられるものとする．上式を ijk 番目の格子セルにわたって積分すれば
qnijk+1 = qnijk ; t (H^ 1 )i+1=2 j k ; (H^ 1)i;1=2 j k + (H^ 2 )i j+1=2 k ; (H^ 2 )i j;1=2 k
+ (H^ 3 )i j k+1=2 ; (H^ 3 )i j k;1=2
=
=
= 1)
(17.126)
qijk = q^ijk は q^(
) の点 i j k における値ではなく q^(
) の ijk 番目のセルにわたる平均
^ 1 )i+1=2 j k は数値流束 F1 = F^1 のセル境界面 i+1=2 にわたる平均値 H^ 1 (t qLi+1=2 j k qRi+1=2 j k )
値，また (H
の時間間隔 t における平均値として次のように定義される．
ただし
qijk =
q(
Z
i+1=2
i;1=2
d
Z
j +1=2
d
j ;1=2
tn+1
Z
(H^ 1 )i+1=2 j k = 1t
tn
Z
k+1=2
k;1=2
q(
)d
H^ 1 (t qLi+1=2 j k qRi+1=2 j k ) dt
(17.127a)
(17.127b)
の区分的多項式で，セル境界面上では一般に不連続である．
) は各セル内で連続な
あるいは Navier-Stokes 方程式 (17.65a) に相当の Euler 方程式
@ q~ + B @ F^1 + @ F^2 + @ F^3 = 0
@t
@
@
@
を写像空間内の単位立方体格子セルにわたって積分すれば次式が得られる．
n+1 = q~n ; t (BH^ 1 )
q~ijk
i+1=2 j k ; (BH^ 1 )i;1=2 j k + (BH^ 2 )i j +1=2 k ; (BH^ 2 )i j ;1=2 k
ijk
+ (BH^ 3 )i j k+1=2 ; (BH^ 3 )i j k;1=2
+ t (Bi+1=2 j k ; Bi;1=2 j k )(F1 )ijk + (Bi j+1=2 k ; Bi j;1=2 k )(F2 )ijk
+ (Bi j k+1=2 ; Bi j k;1=2 )(F3 )ijk
40 記号については
17.1.2 項参照
(17.128)
54
圧縮性流れの解法|NS 方程式 ただし
q~ijk は q~(
) の ijk 番目のセルにわたる平均値，また上式の最終項は部分積分を取る際に出てく
る項で次のようなものである41 ．
0
BBP3
BP3
B F1 + B F2 + B F3 = B
BB
B@P3
1
^ C
=1 (@ x =@ )f +1 C
CC
^
CC
(
@
=@
)
f
x
+1
=1
C
^
=1 (@ x =@ )f +1 A
0
0
(@ x=@ )ijk = ( x )i+1=2 j k ; ( x )i;1=2 j k ，また f^12 は列ベクトル F1 の第 2 成分である．
再構成段階では，セル平均値 qijk を基に区分的 2 次多項式が再構成されそのセル境界値が求められる．ijk
R
番目のセルの境界 = i+1=2 = i;1=2 における q のセル断面平均値 qL
i+1=2 qi;1=2 は次式で与えられる．
qLi+1=2 = qi + 16 minmod( qi;1=2 b qi+1=2 )+ 13 minmod( qi+1=2 b qi;1=2 )
qRi;1=2 = qi ; 16 minmod( qi+1=2 b qi;1=2 ) ; 31 minmod( qi;1=2 b qi+1=2 )
これらの境界値は，q の 一定セル断面平均値に対する区分的 2 次多項式 qjk ( ) = q( ) を 16.5.4 項にな
らって再構成しそれから求めたものである．q の または に関する区分的 2 次式のセル境界値も同様に求
R
められる．セル境界面の左右の qL
l+1=2 ，ql+1=2 は 3 次の有限体積法 TVD スキームでは次のようになる．
qLl+1=2 = ql + 61 ~ql;1=2 + 13 ~ql+1=2
qRl+1=2 = ql+1 ; 16 ~ql+3=2 ; 31 ~ql+1=2
;
~qj+1=2 = minmod qj+1=2 b qj;1=2
~qj+1=2 = minmod; qj+1=2 b qj+3=2
(17.129a)
(j = l l +1)
(j = l ; 1 l)
(17.129b)
区分的 2 次式のセル中心の値は minmod 関数が働かないときには
1 2q
q( i ) = qi ; 24
i
(17.130)
となる．この値はセル断面
= i にわたる平均値で，q( j ) q( k ) とは完全に一致しない．またセル境界
i+1=2 における変数 q の跳躍量は
qi+1=2 = qRi+1=2 ; qLi+1=2 = ; 16 3 qi+1=2
となる．この跳躍量ももとの関数 q(
41
Z
B
@ F1
d
@
=
B F1
;
Z
B
F1 d
(17.131)
) が滑らかに変化するところでは通常無視できる大きさになる．
有 限 体 積 法 17.6.2
55
発展段階|Roe の近似リーマン解
セル境界面 xl+1=2 における数値流束 H は次式から求められる．
(Hi )l+1=2 = 12 Fi (qlL+1=2 )+ Fi (qlR+1=2 ) ; 21 jAi (qlL+1=2 qlR+1=2 )j q
(17.132)
q(x) は各セルごとに独立の x の区分的多項式で，その値はセル境界面で一般に連続にはならない．qlL+1=2
と qlR+1=2 はセル境界面の左右の q の値で，この式の右辺第 1 項は境界面の左右の流束の値の平均を取るこ
とを意味する．また第 2 項はこの不連続 q のために境界から発生する波の項である．セル境界面における
;
+
;
流束の跳躍量は Fi = Ai q = Ri +
i Li q + Ri i Li q で，Ri i Li q Ri i Li q はそれぞれ x の正または
負の向きに伝播する波を横切っての流束の跳躍量である．TVD 安定な解は，(Hi )l+1=2 をこれらすべての
波の上流側に取ることによって得られる．この流束の値は境界面左側の流束に負の向きの波を横切っての流
束の跳躍量を加えたもの Fi (qlL+1=2 )+ Ri
跳躍量を引いたもの Fi (qlR+1=2 )
;Li q ，または右側の流束から正の向きの波を横切っての流束の
i
;Ri +i Li q である．式 (17.132) はこれらの式の平均を取り対称性を持たせ
たものである．
ここではセル境界の不連続に起因する波を Roe の近似リーマン解で求めることにする．ヤコビ行列 A は
次の条件を満足するように決定できれば ，スキームは保存形になる．
F = A q FR ; FL = A(qR ; qL)
(17.133)
A は qL と qR の関数で，qL = qR = q ならば A(q q) = A(q) = @F=@q ，また A は実固有値と 1 次独
立の固有ベクトルを持つものとする．Roe はこれらの 3 条件を満足する A を巧妙な手段で見出している42 ．
未知変数ベクトル q と流束ベクトル F は
ただし
0 1
0 1
z1
1
C
B
B
B
BB u CCC
z2 C
C
B
z=B
B
CC = p BBB v CCC
z3 C
B
B
B@ w CA
@z4CA
z5
(17.134)
H
で定義されるベクトル z の成分の 2 次式で表すことができる．すなわち
0 1 0 z12
B
C BB z1z2
B
uC
B
CC BB
q=B
B
CC = BB z1z3
v
B
B
B
@ wCA B@ z1z5 z1z~4
e
+ 2 zk zk
1
CC
CC
CC
CC
A
(17.135a)
Roe, P. L., Approximate Riemann solvers, parameter vectors and di erence schemes, J. Comput. Phys., 43(1981),
357{72.
42
56
圧縮性流れの解法|NS 方程式 0
z1 zi+1
0
1 BB
ui
BBz2zi+1 + 1i ~ ;z1z5 ; 1 zk zk
B
C
B
2
uui + 1i p C
B
CC BB
;
B
~
B
Fi = B
= Bz3 zi+1 + 2i z1 z5 ; 1 zk zk
vui + 2i p C
B
C
2
B
B
@ wui + 3ipCA BBBz z + ~ ;z z ; 1 z z
Hui
@ 4 i+1 3i 1 5 2 k k
z5 zi+1
1
CC
CC
CC
CC
CC
CC
A
(i = 1 2 3)
(17.135b)
zk zk = z22 + z32 + z42 ~ = ; 1 である．q と Fi の中の 2 次式 z12 z1 z2 : : : の差分は， (ab) =
(ab)R ; (ab)L = a b + b a のように書き換えることができる．ただし u = uR ; uL u = (uL + uR)=2 であ
る． q Fi に対しこの書換えを実行し整理すれば次式が得られる．
ただし
1
0
2z1 0
0
0
0
BB z z 0 0 0 CC
CC
BB 2 1
CC
z
0
z
0
0
3
1
q=B z B=B
BB
B@ z4 0 0 z1 0 CCA
1z ~z ~z ~z 1z
5
2
3
4
1
0
1i z1
2i z1
BB zi+1
~z
1 z +z
BB 1i 5 1i 2 i+1 2i z2 ; 1i ~ z3
B ~
Fi = Ci z Ci = B
BB 2i z5 1i z3 ; 2i ~ z2 2i 1 z3 + zi+1
BB ~
B@ 3i z5 1i z4 ; 3i ~ z2 2i z4 ; 3i ~ z3
0
これより
1i z5
(17.136a)
3i z1
3i z2 ; 1i
3i z3 ; 2i
~z
~z
4
4
1 z +z
4 i+1
3i
2i z5
3i z5
1
CC
1i 1 C
C
~z C
C
2i 1 C
C
~z C
C
3i 1 C
A
0
~z
(17.136b)
zi+1
Fi = Ci z = Ci B ;1 q ，式 (17.133) を満足する Ai = Ci B ;1 は次のように求められる．
0
0
1i
B
~
z
z
z
z
B
k
k
2
i
+1
B
; z 2 + 1i 2 z 2 1i (1 ; ~) zz2 + zzi+1
B
1
1
1
1
B
z
z
B
~
z
z
z
z
3
i
+1
k
k
3
2
; z 2 + 2i 2 z 2
1i z ; 2i ~ z
Ai = B
B
1
1
1
1
B
B
~
z
z
z
z
z
z
k
k
4
2
4
i
+1
B
; z 2 + 3i 2 z 2
1i z ; 3i ~ z
B
1
1
1
1
B
@; z5zi+1 + ~ zk zk zi+1
z5 ; ~ z2 zi+1
1i z
2
2
z12
1
0 z1 2 z1 z1
0
1i
2i
B
~
2
B
;uui + 1i 2 u
D i1
B
2i u ; 1i ~v
B
B
D i2
;vui + 2i ~2 u2 1i v ; 2i ~u
=B
B
B
B
B
;wui + 3i 2~ u2 1i w ; 3i ~u 2i w ; 3i ~v
B
@
~ 2
;H ui + 2 u ui
1i H ; ~uui
2i H ; ~v ui
2i
z2 ; ~ z3
1i z
1
1
z
z
3
i+1
2i (1 ; ~) z + z
1
1
z4 ; ~ z3
2i z
3i z
1
1
z5 ; ~ z3 zi+1
2i z
z2
2i z
1
z2
z4
3i z ; 1i ~ z
1
1
z3 ; ~ z4
3i z
2i z
1
1
z
z
4
i+1
3i (1 ; ~) z + z
1
1
z5 ; ~ z4 zi+1
3i z
z12
1
1
0
3i
CC
3i u ; 1i ~ w 1i ~C
CC
3i v ; 2i ~w 2i ~C
CC
C
D i3
3i ~C
CA
1
3i H ; ~w ui
ui
3i
1
CC
1i ~ C
CC
C
2i ~ C
CC
C
3i ~ C
CC
zi+1 A
0
z1
(17.137)
57
有 限 体 積 法 ただし u
2 = u 2 +v 2 +w 2 D
ij
ui または H の平均値である．
p L uiL + p R uiR
ui = p + p
L
= ji (1; ~)uj +ui ，また ui = zi+1 =z1 H = z5=z1 は次のように定義される
R
p L H L +p R H R
H = p +p
L
R
この 行列Ai は，Roe の第 1 条件 (17.133) を満足するように導かれたものであり，また ui
り第 2 条件
H の定義式によ
qL = qR = q のとき Ai (q q) = Ai (q) も満足されることになる．また Ai は保存形方程式のヤ
コビ行列 Ai の式 (17.69b) と同形で，5 個の実固有値と 5 個の 1 次独立の固有ベクトルを持つもので第 3 条
件も満足されることになる．
保存形の Ai を非保存形のものにする行列 N とその逆行列 N ;1 は式 (17.72) に倣えば次のようになる．
0
1
BB ;u=
B
N =B
BB ;v=
B@;w=
0
0
0
0
0
0
0
BB u1
B
N ;1 = B
BB v
B@ w
1
0
C
1=
0C
C
CC
0 1=
0C
0
0 1= 0 C
A
2
~ u =2 ;~u ;~v ;~w ~
0
u2=2
0
0
u
0
0
0
v
0
0
0
0
0
0
0
w 1= ~
1
CC
CC
CC
CA
(17.138)
またこの N が実際の計算に便利に利用されるには
0
1
BB ;u=
B
N q=B
BB ;v=
B@;w=
0
0
0
1=
0
0
0 1=
0
0
0 1=
2
~ u =2 ; ~ u ; ~ v ; ~ w
10
0
CBB
0C
CBB
CCBB
0C
0C
AB@
~
1 0 1
C BB u CC
( u) C
C BB CC
CC = BB v CC = q
( v) C
( w)C
A B@ wCA
e
(17.139)
p
でなければならない．すなわち次の２つの式が満足されなければならない．
;
ui = 1 ; ui + ( ui )
(i = 1 2 3)
p = ~ e ; 12 ~ ( u 2 ) = 12 ~ u2 ; ~uk ( uk )+~ e
ところでこれらの式の中の ( ui ) と ( u 2 ) は次のように書き換えることができる43 ．
;
( ui ) = ( ui )R ; ( ui )L = (p ui )R +(p ui )L (p R ; p L )+ p L R (ui )R ; (ui)L )
= ui + p L R ui
( u2 )= ( u2 )R ; ( u2 )L = (p uk )R +(p uk )L (p uk )R ; (p uk )L
= uk (p R + p L ) (p uk )R ; (p uk )L
= u ( u )+ p
u u
k
k
R L k k
したがって上記の２つの式の上式は
=p L R
43 ニュートンの総和則を用いる．
( u) = u + u すなわち密度の平均値
を
58
圧縮性流れの解法|NS 方程式 のように定義すれば成立する．また下式もこの
の定義のもとで成立する．
非保存形の係数行列 Ai = NAi N ;1 は次のようになる．
0
BBu0i
B
Ai = B
BB 0
B@ 0
2i
0
0
ui
ui
0
ただし
1i
1i
0
0
c2
2i
1
C
1i = C
CC
=
2i C
CC
3i = A
0
3i
0
0
ui
3i c 2
c2
(17.140)
ui
c は次式で定義される平均音速である．
;
c 2 = ~ H ; 12 u2
この Ai も形式的に式 (17.70b) の Ai と同形で，その固有値
ui ;
0
0
0
0
i
1i
ui ;
0
0
1i
2i
i
ui ;
0
c2
2i
i
c2
jAi ; i I j = 0 すなわち
0
3i
0
iは
1i =
2i =
3i =
0
0
ui ; i
3i c 2 ui ;
= (ui ; i )3 (ui ; i )2 ; ( 1i + 2i + 3i )c 2 = 0
i
から計算することができ次のように求められる．
0
ui
B
B
ui + 1i c
B
B
=
ui + 2i c
i B
B
B
ui + 3ic
@
0
0
1
CC
CC
CC
CA
(17.141)
ui ; c
(k = 1 2 : : : 5) は lik (Ai ; ik I ) = 0 から求められる．その行列を Li ，そ
の逆行列を Ri とすれば ，これらの固有ベクトルの行列は次のように式 (17.74b) または (17.74c) と形式的
また Ai の左固有ベクトル lik
に同じものになる．
0
BB10
B
Li = B
BB0
B@0
0
0
1
0
0
0
0
1
0
1i
2i
1
0 ;1=c 2
C
0 1i = c C
C
CC
0 2i = c C
1 3i = c C
A
3i ;1= c
1
0
1 1i =2c 2i =2c 3i =2c ; =2c
C
BB0 1 ; =2 0
0
1i
1i =2 C
BB
CC
Ri = B
CC
0
0
1
;
=
2
0
=
2
2
i
2
i
BB
0
1 ; 3i=2 3i =2 C
@0 0
A
0
1i
c=2
2i
c=2
3i
(17.142)
c=2 ; c=2
Ai の固有値は Ai のものと同じになり，Ai の固有ベクトルの行列は Li = Li N Ri = N ;1Ri となる．こ
れより Ai は確かに 5 個の実固有値 (重根を含む) と 5 個の 1 次独立の実固有ベクトルを持つことが分かる．
実際の計算に必要な特性速度，特性変数の差分，右固有ベクトルすなわち i @ W i = Li q Ri の成分
有 限 体 積 法 59
はまとめて書けば次のようになる．
i1 = ui
i5 = ui ; c
i k+1 = ui + ki c
(k = 1 2 3)
(17.143a)
@ w1i = ; p=c 2 @ wki +1 = uk +
@ w5i = ui ; p= c
ki
p= c
0 1
0
1
1
ki
B
C
BB u + (2 ; )cCC
B
u1 C
B
C
BB ki 1 k1 1i CC
B
C
k
+1
1
ri = B
u2 C
B
CC ri = 2c BBB ki u2 + k2(2 ; 2i)cCCC
B
@ u3 A
@ ki u3 + k3(2 ; 3i)cA
u2=2
0
1
1
BBu ; cCC
B 1 ki CC
ri5 = ; 2c B
BBu2 ; kicCC
B@u3 ; kicCA
(k = 1 2 3)
(17.143b)
(k = 1 2 3)
ki H +(2 ; ki )c uk
(17.143c)
H ; c ui
= p L R = R L ui = uiL1++RRuiR H = HL1++RRHR
(17.143d)
p
c 2 = ~ H ; 12 u2 R
=
R L
以上述べた Roe の近似 Riemann 解を用いれば ，セル境界における数値流束は次のようになる．
5 ;
X
(Hi )l+1=2 = 12 Fi (qlL+1=2 )+ Fi (qlR+1=2 ) ; 21
j ik j @ wki rik l+1=2
k=1
(17.144)
ik @ wki rik はセル境界面から発生する k 番目の波を横切っての流束の跳躍量である．この近似リーマン解
を導入しリーマン問題を特性の理論によっ
は，物理変数に代わるものとしていわゆる Roe 平均 ui H
て陽的に解くもので，リーマン問題を厳密に解くときに現れる膨張扇と衝撃波は圧力波として求められる
ことになる．Roe の近似リーマン解において Ai
q は保存性を持つ量で，したがってこのスキームは保存性
を有する．またこのスキームは不連続をかなりの精度で捕獲する能力を秘めている．セル境界面における q
S
S
の成分 S
l+1=2 ( ui )l+1=2 el+1=2 (S = L R) の値は前項の式 (17.123) から求められ，セル境界面における
流束 Fi (qlS+1=2 ) (S = L R) はこれらの値を用い次式から求められる．
0
1
ui
BB( u )( u )= + pCC
1i C
BB 1 i
C
Fi (q) = B
BB( u2)( ui)= + 2ipCCC
@( u3)( ui)= + 3ipA
(e + p)( ui)=
ただし
p = ~fe ; ( uk )( uk )=2 g である．
(17.145)
60
圧縮性流れの解法|NS 方程式 17.6.3
発展段階|曲線座標格子における Roe の近似リーマン解
^ を求めることができる．
曲線座標格子に対しても 同様にして セル境界面における数値流束 H
(17.146)
(H^ i )l+1=2 = 12 Fi (qLl+1=2 )+ Fi(qRl+1=2 ) ; 21 jAi (qLl+1=2 qRl+1=2 )j q
q( ) は の区分的多項式，qLl+1=2 qRl+1=2 はセル境界面の左右の q の値，また最後の項はセル境界面におけ
る不連続による波の項でここでは Roe の近似リーマン解によって求められる．
Roe の近似リーマン解では，まずヤコビ 行列 A が次の 3 条件を満足する A で近似される．すなわち qL
と qR の関数 A は (i) 保存の条件
F=A q
FR ; FL = A(qR ; qL)
(17.147)
qL = qR = q ならば A(q q) = A(q) = @ F=@ q，(iii) また特性の理論を適用すべく実固有値と
1 次独立の固有ベクトルを持つように決定される．その導出過程は上記と同様である．未知変数ベクトル q
と流束ベクトル F も式 (17.134) のベクトル z の成分の 2 次式で表すことができる．すなわち
1
0 1 0
z12
BB u CC BB z1z2 CC
CC
B CC BB
(17.148a)
q = JB
BB v CC = J BB z1z3 CC
B@ wCA BB z1z4 CC
A
@1
z1z5 + 2~ zk zk
e
を満足し ，(ii)
0
z1Zi
0
1 BB
Ui
BB uU + p CC BBz2Zi + i 1 ~ ;z1z5 ; 12 zk zk
BB i i 1 CC BB
BBz3Zi + i 2 ~ ;z1z5 ; 1 zk zk
Fi = J B
=
C
vU
+
p
i
i
2
BB
C B
2
@ wUi + i 3pCA BBBz Z + ~ ;z z ; 1 z z
HUi
@ 4 i i3 1 5 2 k k
1
CC
CC
CC
CC
CC
CC
A
(17.148b)
z5Zi
ただし Zi = i 1 z2 + i 2 z3 + i 3 z4 zk zk = z22 +z32 +z42 である．q と Fi の中の 2 次式 z12 z1 z2 : : : の差分は，
(ab) = (ab)R;(ab)L = a b+b a のように書き換えることができる．ただし u = uR;uL u = (uL+uR )=2
である． q Fi に対しこの書換えを実行し整理すれば次式が得られる．
q = BJ z
(17.149a)
Fi = Ci J z
0
1
Z
0
i
i
1 z1
i
2 z1
i
3 z1
BB ~
C
BB i 1 z5 1 i 1 z2 + Zi i 2 z2 ; ~ i 1 z3 i 3 z2 ; ~ i 1z4 ~ i 1 z1CCC
BB ~
Ci = B
BB i 2 z5
BB ~
B@ i 3 z5
0
~
i 1 z3 ; i 2 z2
z ;~ z
i1 4
i3 2
i 1 z5
1
i 2 z3 + Zi
i 2 z4 ;
~
i 2 z5
i 3 z3
C
~
~ zC
i 3 z3 ; i 2 z4
i 2 1C
CC
1 z +Z
~ zC
i3 4 i
i 3 1C
C
i 3 z5
Zi
A
(17.149b)
61
有 限 体 積 法 B は式 (17.136a) で定義された行列である．これより Fi = Ci J z = Ci B ;1 q，式 (17.147) を満足
する Ai = Ci B ;1 は次のように求められる．
ただし
0
0
B
B
B
; zz2Z2 i + 2~ i 1 zzk z2k
B
1
1
B
B
z
Z
~
z
3
i
k
B
; z 2 + 2 i 2 z z2k
Ai = B
B
1
1
B
z
Z
~
z
B
4
i
k
B
; z 2 + 2 i 3 z z2k
B
B
@ z 1Z ~ z z Z1
i1
i2
i3
(1 ; ~) i 1 zz2 + Zz i i 2 zz2 ; ~ i 1 zz3 i 3 zz2 ; ~ i 1 zz4
1
1
1
1
1
1
z3 ; ~ z2 (1 ; ~) z3 + Zi
z3 ; ~ z4
i 1z
i 2z
i 2z z
i 3z
i 2z
1
1
1
1
1
1
z4 ; ~ z2
z4 ; ~ z3 (1 ; ~) z4 + Zi
i 1z
i 3z
i 2z
i 3z
i 3z z
1
1
1
1
1
1
z
z
z
z
z
z
Z
Z
Z
; z5 2 i + 2 zk 2k z i i 1 z5 ; ~ z2 2 i i 2 z5 ; ~ z3 2 i i 3 z5 ; ~ z4 2 i
1
1
1
1
1 1
1
1
1
0
0
B
B
;
uU i + 2~ i 1u2
B
B
B
;vU i + ~2 i 2u2
=B
B
B
B
B
;wU i + ~2 i 3u2
B
@
~ 2
;H U i + 2 u U i
ただし u
2
1
C
D^ i1
i 2u; ~ i 1v i 3u; ~ i 1w ~ i 1C
CC
C
D^ i2
i 1v ; ~ i 2u
i 3v ; ~ i 2w ~ i 2C
CC
C
D^i3
~ i 3C
i 1w ; ~ i 3u i 2w ; ~ i 3v
CA
i1
i 1 H ; ~uU i
i2
i
z1
0
i3
i 2 H ; ~v U i
1
CC
~ i 1C
CC
C
~ i 2C
CC
C
CC
~ i 3C
C
ZA
0
i 3 H ; ~w U i
(17.150)
Ui
= u 2 + v 2 + w 2 D^ij = (1 ; ~) i j uj + U i ，また ui = zi+1 =z1 U i = Zi =z1 H = z5 =z1 は次のよ
うに定義される ui
Ui H の平均値である．
p
p L uiL + p R uiR
p R UiR
U i = LpUiL +
ui = p + p
p
L
R
L+ R
H=
p L H L +p R H R
p L +p R
この行列 Ai は，Roe の第 1 条件 (17.147) を満足するように導かれたもので，その第 2 条件も ui
Ui H の
定義式によって満足されることになる．また Ai は保存形方程式のヤコビ行列 Ai の式 (17.79b) と同形で 5
個の実固有値と 5 個の 1 次独立の固有ベクトルを持つもので第 3 条件も満足されることになる．
保存形の Ai を非保存形のものにする行列は式 (17.138) の N である．またこの N が実際の計算に便利に
利用されるためには
0
1
BB ;u=
B
N q=B
BB ;v=
B@;w=
0
0
0
1=
0
0
0 1=
0
0
0 1=
2
~ u =2 ; ~ u ; ~ v ; ~ w
10
0
C BB
0C
C BB
CCJ BB
0C
C B@
0A
~
1 0 1
C BB u CC
( u) C
C BB CC
CC = J BB v CC = J q
( v) C
( w )C
A B@ wCA
e
(17.151)
p
なる関係が成立しなければならない．この式は直角座標の場合の式 (17.139) と等価なもので，密度の平均
値
を
=p L R
のように定義すれば成立することになる．
62
圧縮性流れの解法|NS 方程式 非保存形の係数行列 Ai
0
Ui
B
B
0
B
B
B
Ai = B
0
B
B
B
@0
i1
i2
0
0
Ui
Ui
1
CC
=
i1 C
CC
i 2= C
CC
i 3= C
A
0
i3
0
0
0
Ui
i 3 c2
0
i 1 c2
0
ただし
= N Ai N ;1 は次のようになる．
i 2 c2
(17.152)
Ui
c は次式で定義される平均音速である．
;
c 2 = ~ H ; 12 u2
この Ai も形式的に式 (17.80b) の Ai と同形で，その固有値 ki は jAi ; ki I j = 0 すなわち
Ui ; ki
0
0
0
0
i1
i2
0
0
Ui ; ki
Ui ; ki
i 1 c2
0
0
i 2 c2
0
i3
i 1=
0
0
3
2 ii 2
i 2 = = (Ui ; ki ) (Ui ; ki ) ; g c = 0
Ui ; ki i 3 =
i 3 c 2 Ui ; ki
から計算することができ次のように求められる．
0
BB Ui
BB Ui + 1ipg11c
B
p
Ki = B
Ui + 2i g22 c
BB
pg33c
BB
U
+
i
3
i
@
0
0
p
1
CC
CC
CC
CC
CC
A
(17.153)
Ui ; gii c
(k = 1 2 : : : 5) は lik (Ai ;kik I ) = 0 から求められる．その行列を Li ，そ
の逆行列を Ri とすれば ，これらの固有ベクトルの行列は次のように式 (17.82b) または (17.82c) と形式的
また Ai の左固有ベクトル lik
に同じものになる．
0
BB 1
BB 0
B
Li = B
BB 0
BB 0
@
0
0
0
0
ii
1i 1 2 ; 2i 2 1
1i 1 3 ; 3i 3 1
2i 2 1 ; 1i 1 2
ii
2i 2 3 ; 3i 3 2
3i 3 1 ; 1i 1 3
3i 3 2 ; 2i 2 3
ii
i1
i2
i3
1
;1=c 2
CC
p
1i g 11 = cC
pg22= cCCC
2i
C
pg33= cCCC
3i
p ii A
; g =c
(17.154a)
63
有 限 体 積 法 0
BB1
BB
BB0
BB
Ri = B
BB0
BB
BB0
BB
@0
D ij
p1i
2 g11 c
p2i
2 g22 c
3 ;1 i 1 i 2
2 2i
gii i i
D i1
3 ;1 i 2 i 1
2 1i
gii i i
3
i3 i1
2 1i ; 1 gii i i
D i2
3 ;1 i 3 i 2
2 2i
gii i i
2i c
p
2 g22
1i c
p
2 g11
2
= 32 ji ; 1 giii j +(1 ; ji) 1
ii
ii
1
p3i
2 g33 c
3
i1 i3
2 3i ; 1 gii i i
3 ;1 i 2 i 3
2 3i
gii i i
; p ii C
2 g cC
i1
2gii
i2
2gii
i3
D i3
2gii
3i c
p
2 g33
; pc ii
2 g
CC
CC
CC
CC
CC
CC
CC
CA
(17.154b)
(j = 1 2 3)
Ai の固有値は Ai のものと同じになり，Ai の固有ベクトルの行列は Li = Li N Ri = N ;1Ri となる．こ
れより Ai も確かに 5 個の実固有値 (重根を含む) と 5 個の 1 次独立の実固有ベクトルを持つことが分かる．
実際の計算に必要な特性速度，特性変数の差分，右固有ベクトルすなわち Ki @ Wi = Li q Ri の成分
はまとめて書けば次のようになる．
p
ki1 = U i kik+1 = U i + ik gii c
p
ki5 = U i ; gii c
(k = 1 2 3)
(17.155a)
p
w1i = ; p=c 2 wki +1 = ki Ui + i i uk ; i k ui + ki gii p= c
p
w5i = Ui ; gii p= c
(k = 1 2 3)
0
1
B
B
u1
B
B
1
ri = B
u2
B
B
@ u3
1
0 1
CC
BBu1 CC
CC k+1
BB 1CC
ki
CC ri = 2pgii c BBu2CC +
CA
B@u3CA
u 2=2
H
0 1
0 1
1
BB 0 CC
BBu CC
BB i 1 CC
B 1CC
ri5 = ; p ii B
+
B i 2C
B
C
u2
2 g cB
B@u3CCA 2gii BB@ i 3 CCA
Ui =
=p
c2 = ~
Ui
0 1
BB 0 CC
BB 1k CC
BB 2k CC
B@ 3k CA
(k = 1 2 3)
uk
(17.155c)
Ui
H
ただし
3
2
0 1
BB 0 CC
BB i 1 CC (1 ; ki)
i
k
ki ; 1 g ii B
i 2C+
i iB
ii
B@ i 3 CCA
(17.155b)
i k uk ，
L R=R L
H ; 12 u2
ui = uiL1++RRuiR
R
p
R= L
U i = UiL1++RRUiR
H = HL1++RRHR
(17.155d)
64
圧縮性流れの解法|NS 方程式 以上述べた Roe の近似 Riemann 解を用いれば ，セル境界における数値流束は次のようになる．
5 ;
X
jkik j@ wki rik l+1=2
(17.156)
(H^ i )l+1=2 = 12 Fi (qLl+1=2 )+ Fi(qRl+1=2 ) ; 21
k=1
この Roe の近似リーマン解に関しても前項と同様のことが言える．セル境界面における q の成分の (J )S
l+1=2
(J ui )Sl+1=2 (Je)Sl+1=2 (S = L R) の値は式 (17.129) から求められ，セル境界面における流束 Fi (qSl+1=2 ) (S =
L R) はこれらの値を用い次式から求められる．
0
1
J ui
BB(J u )(J U )=J + JpCC
i1 C
B 1 i
Fi (q) = B
BB(J u2)(J Ui )=J + i 2JpCCC
B@(J u3)(J Ui )=J + i 3JpCA
(17.157)
(Je + Jp)(J Ui)=J
ただし
J Ui =
i k J uk
Jp = ~fJe ; (J uk )(J uk )=2J g である．
この項の終わりにひとこと加えれば ，曲線座標格子のセル境界面における流束分離法で上流化された数
値流束は，セル境界面の左右の物理空間における諸量に対し前項の Roe の近似リーマン解を適用し求める
こともできようが，本項では保存性をより厳密に追求すべく写像空間における諸量に対し
Roe の近似リー
マン解を導出したのである．
17.6.4
発展段階|時間積分
時間積分は 16.5 節に述べた 1 次元流れの陽的 2 段階法，Crank-Nicholson 型時間積分法，Roe の波動分
離法を多次元に拡張して行うこともできる．ここには陽的 Runge-Kutta 法と
型陰解法について述べる．
3 次 Runge-Kutta 法の式は長方形格子の場合には次のようになる．
n+1 = q n + 1 ;k 1 +4k 2 + k 3 ; t(Dn + g n )
qijk
ijk 6 ijk
ijk ijk
ijk ijk
n
1 ;
1 = ; t 1 ;H (q n
n
n
n
kijk
x 1 i+1=2 j k ) ; H1 (qi;1=2 j k ) + y H2 (qi j+1=2 k ) ; H2(qi j;1=2 k )
o
;
) ; H (q n
)
+ 1 H (q n
2
また kijk
z
3 i j k+1=2
3 i j k;1=2
3 は k 1 の式中の q n をそれぞれ q n + k =2 または q n + k 2=2 で置換えたもの，k は q n + k 1=4
kijk
ijk
ijk
で置換えたものである．
4 次 Runge-Kutta 法の式は曲線座標格子の場合には次のようになる．
; 1 +2k^2 +2k^3 + k^4 ; t(D^ n +^gn )
qnijk+1 = qnijk + 16 k^ijk
ijk
ijk ijk
ijk ijk
n
1 =; t H
^ 1 (qni+1=2 j k ) ; H^ 1(qni;1=2 j k ) + H^ 2 (qnij+1=2 k ) ; H^ 2(qnij;1=2 k )
k^ijk
o
+ H^ 3 (qnij k+1=2 ) ; H^ 3(qnij k;1=2 )
3 k^4 は k^1 の式中の qn をそれぞれ qn +k^1=2 qn +k^ =2 または qn +k^
k^ijk
ijk
ijk
n
1
2
^
^
は q +(k + k )=4 で置換えたものである．
2
^ijk
また k
^ijk
で置換えたもの，k
有 限 体 積 法 あるいは
;
n+1 = q~n + 1 Bijk k^1 +2k^2 +2k^3 + k^4 ; tBijk (D
n +^g n )
^ ijk
q~ijk
ijk
ijk
ijk ijk
ijk
ijk 6
+ t (B F1 )ijk + (B F2 )ijk + (B F3 )ijk
n
(B k^1 )ijk = ; t Bi+1=2 j k H^ 1 (qni+1=2 j k ) ; Bi;1=2 j k H^ 1 (qni;1=2 j k )
+Bi j+1=2 k H^ 2 (qnij+1=2 k ) ; Bi j;1=2 k H^ 2 (qnij;1=2 k )
o
+Bi j k+1=2 H^ 3 (qnij k+1=2 ) ; Bi j k;1=2 H^ 3 (qnij k;1=2 )
次に
型陰解法の式は定常流れを長方形格子を用いて計算する場合には次のようになる．
n + t
qijk
n 1 ; +n n
;n
n
n
n
x (A1 )ijk qi;1 j k + jA1jijk qijk +(A1 )ijk qi+1 j k
;
n +(A; )n q n
+ 1y (A+2 )nijk qinj;1 k + jA2jnijk qijk
2 ijk i j +1 k
o
;
n +(A; )n q n
n
+ 1z (A+3 )nijk qinj k;1 + jA3jnijk qijk
3 ijk i j k+1 = rhsijk
n1;
rhsijk = ; t
x (H1 )i+1=2 j k ; (H1 )i;1=2 j k
;
+ 1y (H2 )i j+1=2 k ; (H2 )i j;1=2 k
o
;
+ 1z (H3 )i j k+1=2 ; (H3 )i j k;1=2 + Dijk + gijk
n+1 = q n + q n
qijk
ijk
ijk
また非定常流れを曲線座標格子を用いて計算する場合には次のようになる．
(m) + t (A+ )n q (m) + jA jn q (m) +(A; )n q (m) )
qijk
1 ijk ijk
1 ijk i;1 j k
1 ijk i+1 j k
(
m
)
(
m
)
) )
+
n
n
+ (A2 )ijk qi j;1 k + jA2jijk qijk +(A;2 )nijk qi(mj+1
k
(
m
)
(
m
)
(
m
)
+
;
n
n
n
+ (A3 )ijk qi j k;1 + jA3jijk qijk +(A3 )ijk qi j k+1 )
n
(m;1)
(m;1) ; q n )+ 1 ;rd
hsijk
= ;(qijk
ijk 2 hsijk + rd
rd
hsijk = ; t (H^ 1 )i+1=2 j k ; (H^ 1 )i;1=2 j k + (H^ 2 )i j+1=2 k ; (H^ 2)i j;1=2 k
+(H^ 3 )i j k+1=2 ; (H^ 3 )i j k;1=2 + D^ ijk + g^ijk
(m) = q (m;1) + q (m)
qijk
ijk
ijk
65
66
圧縮性流れの解法|NS 方程式 17.6.5
有限体積法の高次コンパクト TVD スキーム
式 (17.123) はそのすべての
q を以下に示すように Dq で置き換えることによって 4 次または 5 次にす
ることができる．
qlL+1=2 = ql + 61 Dq~l;1=2 + 31 Dq~l+1=2
(17.158a)
qlR+1=2 = ql+1 ; 61 Dq~l+3=2 ; 13 Dq~l+1=2
;
Dq~j+1=2 = minmod Dqj+1=2 bDqj;1=2
;
Dq~j+1=2 = minmod Dqj+1=2 bDqj+3=2
(j = l l +1)
(j = l ; 1 l )
4 次スキームの場合には，Dq は次のように１つの式で統一的に与えることもできる．
(j = l l 1)
Dqj+1=2 = qj+1=2 ; 16 3 qj+1=2
また 5 次スキームでは，Dq は次のようになる．
3 3q
Dql+1=2 = ql+1=2 ; 20
l+1=2
Dqj+1=2 = qj+1=2 ; 15 3 qj+1=2
(j = l 1)
この高次化に際しての不安定性を除くために，3 階差分
3q
j +1=2 =
j +1 ; 2 q~j
2 q~
2 q~ = minmod; 2 q
j +1
j +1 b0 2 qj
(17.158b)
(17.159)
(17.160)
3 q は次のように置かれる．
2 q~ = minmod; 2 q b0 2 q
j
j
j +1
(17.161)
式 (17.158b) の b，式 (17.161) の b0 は通常 4 に取られる．
q のセル境界面における不連続 (跳躍量) q の値は，q(x) が滑らかで制限関数の働かないところでは
q
qlR+1=2 ; qlL+1=2 =
; 16 3 ql+1=2
(3 次スキーム)
1 (1 ; ) 5 q
l+1=2 (4 次または 5 次スキーム)
4!
(17.162)
となり一般に非常に小さい．しかしながら 流れの不連続や制限関数の働くところではこの跳躍量は小さく
なく前述のようにここに近似リーマン解が補われる．流れの不連続の近傍を除けば スキームの精度は
2次
でもよく 3 次精度あれば十分である．上記の 4 次以上の高次スキームを用いてもスキームの精度そのもの
4 次以上になるわけではないが，近似多項式の次数が上がれば 式の適応性ないし 柔軟性が増し滑り面な
どの弱い不連続の呆け (smearing 不鮮明化) を防止することができることになる．
曲線座標格子の場合の式 (17.129) の qL qR に対しても上記と全く同様に高次化することができる．
が
67
計 算 の 高 効 率 化 17.7
計算の高効率化
形陰解法の式は，以前には 近似因子法を適用し１次元問題に分解し ，更に対角化することによって計
算量の軽減が図られていたが，その後 Jameson らによる
LU-SGS 法の適用によってなお一層の計算量の削
減が図られた．ここには袁新らによるその改良版についても述べる．
17.7.1
近似因子法
３次元構造格子の未知変数の定義される格子点の数を I
J K 個とする．もともとの
形陰解法は，未
5IJK の非常に大きい連立 1 次方程式を一つ解くものであった．これに Yanenko の近似因子法
(approximate-factorization) 44 を適用し１次元問題に分解すれば ，未知数 5I の連立 1 次方程式を JK 個，
5J のものを KI 個，5K のものを IJ 個解く問題になる．更に特性の理論によって対角化 (diagonalization)
すれば ，未知数 I の連立 1 次方程式を 5JK 個，J のものを 5KI 個，K のものを 5IJ 個解く問題になり計
知数の数
算量はかなり軽減される．
解 法
未知数の数
方程式の数
係数行列の性質
5IJK
1
帯行列 (バンド 幅 2JK +1)
近似因子法を適用
5I
5J
5K
対角化を併用
I
J
K
LU-SGS-GE 法
5
JK
KI
IJ
5JK
5KI
5IJ
2IJK
形陰解法のみ
長方形格子，定常流れの場合の
5 5 ブロック 3 重対角行列
スカラー 3 重対角行列
半対角化された行列
形陰解法の式 (17.93) は近似因子法を適用し１次元問題に分解し更に対
角化すれば次のようになる．
R1 I + x ( +1 rx + ;1 x ) M21;1
I + y ( +2 ry + ;2 y ) M32;1
I + z ( +3 rz + ;3 z ) L3 qn = rhsn
ただし
xi =
;
i xi + i xi
+r
(17.163)
t= xi M21 = L2 R1 M32 = L3R2 ，また rxi xi は xi 方向の後退または前進差分で，
は 1 次上流差分の演算子である．M21 M32 の成分の式は，R1 L2 R2 L3 に含まれる熱
44 Yanenko, N.N., On the implicit di erence computing methods for solving multidimensional heat conduction equations(Russian), Izv. Uchebn. Zaved., Mathematika, 4(1961), 148{157.
Yanenko, N.N., The Method of Fractional Steps. Springer-Verlag, Berlin, 1971.
Beam, R.M. and Warming, R.F., An implicit factored scheme for the compressible Navier-Stokes equations. AIAA J.,
16(1978), 393{402.
68
圧縮性流れの解法|NS 方程式 力学的状態量はすべて消え次のようになる45 ．
0
BB10
B
M21 = B
BB0
B@0
0
1=2
1=2
0
0 ;1=2
1
0
0
1
0
1
0 0
C
0 1=2 C
C
CC
0 ;1=2C
1 0 C
A
0 1=2
(17.164a)
0 0
1 0
0 1=2
0 1=2
0 0 ;1=2
0 0
C
0 0 C
C
CC
0 1=2 C
1 ;1=2C
A
1 1=2
(17.164b)
0
BB10
B
M32 = B
BB0
B@0
1
式 (17.163) は次のように 7 段階に分けて解かれる．
q1 = L1 rhsn
I + x ( +1 rx + ;1 x ) q2 = q1
q3 = M21 q2
I + y ( +2 ry + ;2 y ) q4 = q3
q5 = M32 q4
I + z ( +3 rz + ;3 z ) q6 = q5
qn = R3 q6
(17.165)
この 7 段階のうち奇数段階は各格子点ごとのベクトルと行列の積の計算である．また偶数段階は未知変数
q2 q4 または q6 の連立 1 次方程式の計算であるが，その係数行列は帯行列で 5 組のスカラー３重対角行列
の連立 1 次方程式を解く計算になる．例えば第 2 段階は x 方向に並ぶ I 個の点列に対し 各 5 個の連立 1 次
方程式，都合 5JK 個の連立 1 次方程式を解く計算になる．
非保存形の式 (17.110) に対しては， 形陰解法の式は次のように因子化できる．
R1 I + x ( +1 rx + ;1 x ) M21;1 I + y ( +2 ry + ;2 y ) M32;1
I + z ( +3 rz + ;3 z ) L3 q n = N rhsn
(17.166)
この式も上記のように 7 段階に分けて解くことができる．
残りの
形陰解法の式も同様に近似因子化，対角化して解くことができる．曲線座標格子，非定常流れ
の場合の式 (17.96) は次のようになる．
R~1 I +
I+
45 M
21
( ~+1 r + ~;1 ) M~ 21;1 I + ( ~+2 r + ~;2 ) M~ 32;1
; hsn + rg
(m;1)
( ~+3 r + ~;3 ) L~ 3 q~(m) = ;(~q(m;1) ; q~n )+ 12 rg
hs
= L2 R1 = L2 R1
(17.167)
69
計 算 の 高 効 率 化 =
ただし
t M~ 21 = L~ 2R~1 M~ 32 = L~ 3R~2 また ~+i r i + ~;i
i
は
i 方向の 1 次上流差分演算子である．
M~ 21 M~ 32 は成分の式で表せば次のように測度テンソル成分のみの関数になる46 ．
0
0
BB1 1 n g021 g12 o
12
BB0 2 1 ; g11 g22
; gg22
BB 1 n g21 p
o
B
22=g 11
0
+
g
1
~
B
M21 = B 2 g11
BB 1 n g31 g21 g32 o g32
BB0 2 g11 ; g11 g22 ; g22
B@ n 21 p
o
0 1 g ; g22=g11
1
0
BB1
BB0
BB
B0
M~ 32 = B
BB
BB0
BB
@
0
0
0
1
0
(17.168a)
0 12 gg11 + g22=g11
2 g11
0
1
1
1 n1 ; g21 g12 o C
CC
2
g11 g22 C
C
1 n g21 ; pg22=g11 oC
CC
2 g11
C
1 n g31 ; g21 g32 o C
C
C
2 g11 g11 g22 C
n 21 p
oCA
0
0
1
1 n g12 ; g32 g13 o C
CC
2 g22 g22 g33 C
C
1 n1 ; g32 g23 o C
C
C
2
g22 g33 C
CC
1 n g32 ; pg33=g22 oC
22
C
2 g
oCA
n 32 p
0
1 n g12 ; g32 g13 o ; g13
2 g22 g22 g33
g33
1 n1 ; g32 g23 o ; g23
0
22 33
33
n2g32 gp 33g 22o g
1
0 2 g22 + g =g
1
o
n 32 p
1 12 gg22 + g33=g22
0 0 12 gg22 ; g33=g22
式 (17.167) も 7 段階に分けて同様に解かれる．
非保存形の式 (17.114) に対しては， 形陰解法の式は次のように因子化できる．
(17.168b)
R~1 I + ( ~+1 r + ~;1 ) M~ 21;1 I + ( ~+2 r + ~;2 ) M~ 32;1
~ ; n g(m;1)
I + ( ~+3 r + ~;3 ) L~ 3 q~ (m) = ; (~q (m;1) ; q~ n )+ 2NJ rg
hs + rhs
ただし
q~
17.7.2
(17.169)
は式 (17.89b) で定義されたものである．この式も 7 段階に分けて同様に解かれる．
LU-SGS 法
symmetric Gauss-Seidel method) 47 について述べる．この方法では
陰解法の式 (17.93) に流束分離法を適用した式 (17.109) は次のように書き換えられる．
次に LU-SGS 法 (lower-upper
; x(A+1 )ijk qin;1 j k ; y (A+2 )ijk qinj;1 k ; z (A+3 )ijk qinj k;1
+ (I + x jA1 j + y jA2 j + z jA3 j) qn ijk
+ x(A;1 )ijk qin+1 j k + y (A;2 )ijk qinj+1 k + z (A;3 )ijk qinj k+1 = rhsnijk
形
(17.170)
~ = L~ 2 R~ 1 = L~ 2 R~ 1
Jameson, A. and Yoon, S., Lower-upper implicit schemes with multiple grids for the Euler equations. AIAA J.,
25(1987), 929{935.
Yoon, S. and Kwak, D., Implicit Navier-Stokes solver for three-dimensional compressible ows. AIAA J., 30(1992), 2653{59.
46 M
21
47
70
圧縮性流れの解法|NS 方程式 ただし
xi =
t= xi である．いま
; x(A+1 )ijk qin;1 j k ; y (A+2 )ijk qinj;1 k ; z (A+3 )ijk qinj k;1 (L qn )ijk
(I + x jA1 j + y jA2 j + z jA3 j) qn ijk (D qn )ijk
;
;
n
n
n
n
x (A;
1 )ijk qi+1 j k + y (A2 )ijk qi j +1 k + z (A3 )ijk qi j k+1 (U q )ijk
と置くことにすれば上式は次のようになる．
(U + D + L) qn = rhsn
ただし
L D U は演算子で
L=
X
U=
rxi
xi A+
i (rxi ; 1)
i
X
X
i
D=I+
なお
(17.171a)
qin+1 j k ;
i
xi jAi j
(17.171b)
xi A;
i ( xi +1)
xi は後退差分または前進差分の演算子で ，rx (
n である．
qijk
n ) = qn ; qn
qijk
ijk
i;1 j k
x(
n ) =
qijk
LU-SGS 法では因子化が
(D + L)D;1(D + U) qn = rhsn
(17.172a)
のように行われ，この式は 2 段階に分けて解かれる48 ．
q = D;1 (rhsn ; Lq )
qn = q ; D;1 U qn
(17.172b)
(17.172c)
式 (17.172b) の計算は，i + j + k
= 3 の点から始め，i + j + k = const. 面上の点の式をできるだけ並列処理
するように進め，i + j + k = I +J +K の点まで掃引される．また式 (17.172c) の計算は，i + j + k = I +J +K
の点から逆向きに i + j + k = 3 の点まで掃引される．この前半の計算では，２次元の場合には図 17.14 に
示すように， 印で示す i + j = 8 線上の点の q の値は 印で示す i + j = 7 線上の点の q の既知の値から
求めることができる．３次元の場合も同様のことが言える．
式 (17.172a) の左辺の演算子は (D+L)D;1 (D+U) = L+D+U+LD;1 U のように書換えられるから，基の
式 (17.171a) のものと比べれば ，この式の因子化に伴う (左辺の演算子の) 誤差は LD;1 U となる．L
D U
の大きさはクーラン数の大きさで，この誤差 LD;1 U の大きさもクーラン数の大きさである．一方，前記の
在来型の近似因子法では相当の誤差は２次元ではクーラン数の 2 乗，３次元ではクーラン数の 3 乗になる．
このゆえに LU-SGS 型解法は，大きい Courant 数に対しても解が不安定になりにくく，Courant 数を大き
く取って計算を加速することができる．
48
D;1 (D + U) qn = q と置けば次式が得られる．
(D + L)q = rhsn (D + U) qn = Dq
71
計 算 の 高 効 率 化 y
6
x
@
i+j
=2
@
i+j
=8
図 17.14: ２次元問題の LU-SGS 法
Jameson らの LU-SGS 法では，式 (17.172) の D;1 の計算を容易にするために，行列 D が次のように対
角化近似される．
Ai = (Ai
ただし
i = 1:01
i I )=2
j
jAi j = A+i ; A;i = i I
(17.173)
i max j， i max は Ai の最大固有値である．この解法は簡明で優れた方法である．
この対角化近似はしかしながら
形陰解法の反復計算の収束性を悪くする．袁らによる LU-SGS-GE 法
は，式 (17.172) を対角化近似せずにガウス消去法で正確に解くものである．LU-SGS-GE 法では LU-SGS
法の D を固有ベクトルの行列 Li
Ri を用い D = Ri BLi のように置く．このとき B は後で示すように半
対角化された行列になる49 ．今 q = q 3 と書くことにすれば式 (17.172) は
Ri BLi q3 = rhsn ; Lq3
Ri BLi qn = Dq3 ; U qn
(17.174a)
(17.174b)
= q1 Li q3 = q2 BLi qn = q4 Li qn = q5 と置き，それぞれ 3 段階に
分けて解くことができる．結局 LU-SGS-GE 法では，式 (17.172a) が次のように 6 段階に分けて解かれるこ
となる．これらの式は，BLi q 3
とになる．
q1 = Li (rhsn ; L q3)
Bq2 = q1
q3 = Ri q2
q4 = q1 ; LiU qn
Bq5 = q4
qn = Ri q5
(17.175a)
ただし
B = Li DRi
(i = 1 or 2 or 3)
(17.175b)
これらの式の第 1, 3, 4, 6 段階は各計算点ごとのベクトルと行列の積の計算である．また第 2 段階と第 5 段
階は連立 1 次方程式の計算で Gauss 消去法で解かれる．
49 半対角化の効果は，乱流の計算で乱れ量の輸送方程式が加わる場合にはより顕著である．
72
圧縮性流れの解法|NS 方程式 B は長方形格子の場合には
B = L1DR1 = L1 (I + x jA1 j + y jA2 j + z jA3 j)R1
= I + x j 1 j + y M12 j 2 jM21 + z M13 j 3 jM31
ただし
(17.176a)
M12 = L1 R2 M21 = L2 R1 M13 = L1 R3 M31 = L3 R1 である．
01
BB0
B
M12 = B
BB0
B@0
0 0
1 1=2
0 1=2
0 0
0 1 ;1=2
1
0 0
C
0 ;1=2C
C
CC
0 1=2 C
1 0 C
A
0 1=2
01
BB0
B
M21 = B
BB0
B@0
0
0
1=2 0
1=2 1
0 0
0 ;1=2 1
1
0 0
C
0 1=2 C
C
CC
0 ;1=2C
1 0 C
A
0 1=2
1
0 1
j 2j
0
0
0
0
BB
1
2
3
2
3
1
2
3 C
BB 0 j 22 j + j 42 j + j 42j j 22j ; j 22 j 0 j 22 j ; j 42 j ; j 42 j CCC
B
j 22 j ; j 32 j
j 22 j + j 32 j 0
j 22 j + j 32 j C
CC
M12j 2 jM21 = B
;
BB 0
CC
4
4
2
2
4
4
BB 0
1j
0
0
j
0
CA
2
@
j 12 j ; j 22 j ; j 32 j ; j 22 j + j 32 j
2
4
4
2
2
2 = ui + c 3 = ui ; c である．同様に
i
i
0
ただし 1i
= ui
0
j 12 j + j 22 j + j 32 j
2
4
4
0 1
1
j 3j
0
0
0
0
BB
1
2
3
2
3
1
2
3 C
BB 0 j 23 j + j 43 j + j 43j 0 j 23 j ; j 23 j j 23 j ; j 43 j ; j 43 j CCC
B0
CC
0
j 13 j
0
0
M13j 3 jM31 = B
BB
CC
j 23 j ; j 33 j
j 23 j + j 33 j
j 23 j + j 33 j C
BB 0
0
;
4
4
2
2
4
4 C
@
A
0
j 13 j ; j 23 j ; j 33 j
2
4
4
式 (17.176b) の行列が優対角 (diagonally
0
(17.176c)
2
3
1
2
3
; j 23 j + j 23 j j 23 j + j 43 j + j 43 j
dominant) であることは，第 1, 3, 4 行に対しては明らかで，残
る第 2 行と第 5 行についても下図を参照すれば次式が成立するので明らかである．
2juj + ju + cj + ju ; cj
(17.176b)
2ju + cj; 2ju ; cj + 2juj;ju + cj;ju ; cj
また式 (17.176c) の行列も同様に優対角で，行列 B は優対角行列である．
H 2juj + ju + cj + ju ; cj
A
A
4
c
A
@
;HH
AA
2ju + cj; 2ju ; cj
;
A@
A 2@;
c;
@
A
ju + cj H
juj
ju ; cj
;
H
HH HH A H
H@
H
@
H; H
A
H
c @
u
;c
2juj;ju + cj;ju ; cj
73
計 算 の 高 効 率 化 非保存形方程式は，保存形に比べかなり簡単になるが，在来の近似因子法では前記のように，両者の最
終的に解かれる式は実質的に同等のものであった．次に LU-SGS 法ではどのようになるのかを見ていこう．
非保存形の式 (17.110) は次のように書き換えられる．
(U + D + L ) q n = N rhsn
ただし演算子 L
L =
X
X
X
i
i
(17.177a)
は次のようになる．
xi Ai + (rxi ; 1)
i
D =I+
U =
D U
rhs n
xi jAi j
(17.177b)
xi Ai ; ( xi +1)
この方法では因子化が
(D + L )D ;1 (D + U ) q n = rhs n
(17.178a)
のように行われ，この式は 2 段階に分けて解かれる．
q = D ;1 (rhs n ; L q )
q n = q ; D ;1U q n
(17.178b)
(17.178c)
Jameson らの LU-SGS 法では，式 (17.172) の D ;1 の計算を容易にするために，行列 D
は次のように
対角化近似される．
Ai = (Ai
i I )=2
jAi j = Ai + ; Ai ; = i I
(17.179)
演算子 L と U は L
U に比べ簡単になるが問題の D は D と同じである．
また LUSGS-GE 法では，D = Ri B Li と置かれるが，
B = Li D Ri = Li DRi = B
となり，式 (17.178) は次のように置かれる．
Ri BLi q 3 = rhs n ; L q 3
Ri BLi q n = D q 3 ; U q n
(17.180a)
(17.180b)
3 段階に分けて解くことができる．それらの式の中の q1 q2 q4 q5 は
保存形のものと同じ 値になり，ガウス消去法による連立 1 次方程式の計算量も保存形の場合と同じである
これらの式は上記のようにそれぞれ
が，全体の計算量は多少少なくなる．
次に曲線座標格子の場合の LU-SGS-GE 法について述べる．その諸式は長方形格子の場合の式の A
L D U B に上添え字 \ ~ "を付けたものになる．また x y z はすべて =
q rhs
t と置かれる．それらの
74
圧縮性流れの解法|NS 方程式 ~ の式のみ示す
式の部分段階法による解き方は上記と同じになるのでここでは省略し ，LU-SGS-GE 法の B
ことにする．
B~ = L~ 1D~ R~1 = L~ 1 I + (jA~1 j + jA~2j + jA~3 j) R~1
;
= I + j ^ 1 j + M~ 12j ^ 2 jM~ 21 + M~ 13j ^ 3 jM~ 31
0
0
BB1 0
1
BB0 1 2a (1+ b)
BB0 ; g12
c
~
M12 = B
11
g
2
BB g13 1
BB0 ; 11 ; (1 ; b2)
2
@ g
1
0
0
0
0
1
; 2a (1 ; b) 0 2a (1+ b)
1
0
0
0
BB1
12
1
g
2)
BB0
(1
;
b
;
2
g22
BB
1
1
M~ 21 = B
BB0 1 2 a(1+ b)
BB0 ; g13=J 2g11g22 ; g2232
g
@ 21
0
(17.181a)
1
0
C
; 21a (1 ; b)C
CC
c C
C
2 C
CC
1
2
; 2 (1 ; b )C
CA
1
; 2 a(1 ; b)
1
0
0
0
1
0
1
1 (1 ; b2) C
CC
CC
2
1
; 2 a(1 ; b) C
CC
C
1
; 2 g13 =J 2 g11 g22 C
CA
1
0
2 a(1+ b)
0 1
0
BBj^2j
1
^
BB 0 (1 ; b2) j 2 j +(1+ b)2 j^22 j +(1 ; b)2 j^32j
BB
2
4
4
n
2j
1j
^
^
^3 o
BB
j
j
j
0 a(1 ; b2) ; b 22 +(1+ b) 42 ; (1 ; b) 42j
M~ 12j ^ 2 jM~ 21 = B
BB
n j^12j
BB 0
j^22j +(1 ; b) j^32jo
ac
b
;
(1+
b
)
BB
2
4
4
1
2
3
@
^
^
^
2 j 2j j 2j j 2j
(1 ; b ) 2 ; 4 ; 4
1
0
0
0
CC
^12 j j^22 j j^32 j
j
1 n; bj^1 j +(1+ b) j^22j ; (1 ; b) j^32j o 0
CC
2
(1 ; b ) 2 ; 4 ; 4
2
C
a
2
2
n
2j
3j
1j
2j
3joC
^
^
^
^
^
C
j
j
j
j
j
b2 j^12 j +(1 ; b2) 22 +(1 ; b2) 22
0 a(1 ; b2) ; b 22 ; (1 ; b) 42 +(1+ b) 42 C
CC
CC
n
o
^2 ^3
^1
^2
^3
CC
c j^12 j; j 22 j ; j 22 j
j^12 j
ac b j 22 j +(1 ; b) j 42j ; (1+ b) j 42j
C
^12 j
^22 j
^32 j A
1 n; bj^1 j; (1 ; b) j^22j +(1+ b) j^32j o 0
j
j
j
2
2
2
(1 ; b ) 2 +(1 ; b) 4 +(1+ b) 4
2
a
2
2
0
(17.181b)
75
計 算 の 高 効 率 化 0 1
0
0
BBj^3j
1
2
3
^
^
^
BB 0 (1 ; e2) j 3 j +(1+ e)2 j 3j +(1 ; e)2 j 3 j
0
BB
2
4
4
n ^1
^2
^3 o
B
j^13 j
d f e j 23 j ; (1+ e) j 43j +(1 ; e) j 43j
M~ 13j ^ 3 jM~ 31 = B
BB 0
BB
n 1
2
3o
BB 0 d(1 ; e2) ; e j^23j +(1+ e) j^43j ; (1 ; e) j^43j 0
B@
^1 ^2 ^3
2 j 3j j 3j j 3j
(1 ; e ) 2 ; 4 ; 4
0
0
1
CC
^13 j j^23 j j^33 j
j
1 n; ej^1 j +(1+ e) j^23j ; (1 ; e) j^33j o
CC
2
(1 ; e ) 2 ; 4 ; 4
3
C
d
2
2
n
2 j j^3 j
1j
2j
3jo C
^
^
^
^
CC
j
j
j
j
f j^13 j; 23 ; 23
d f e 23 +(1 ; e) 43 ; (1+ e) 43
CC
n
o
2
3
1
2
3
^
^
^
^
^ C
CC
e2 j^13 j +(1 ; e2) j 23 j +(1 ; e2) j 23 j d(1 ; e2) ; e j 23j ; (1 ; e) j 43j +(1+ e) j 43j C
C
^1
^2
^3 A
j^23j
j^33j o
1 n ^1
2 j 3j
2 j 3j
2 j 3j
0
0
d ; ej 3 j; (1 ; e) 2 +(1+ e) 2
ただし
^1i = Ui
p
^2i = Ui + pgii c
p
(17.181c)
(1 ; e ) 2 +(1 ; e) 4 +(1+ e) 4
^3i = Ui ; pgii c
a = g22 =g11 b = g12= g11g22 c = g23 =J 2 g11 g22
p
p
d = g33 =g11 e = g13= g11 g33 f = g23 =J 2 g11 g33
である50 ．曲線座標格子の場合の行列
B^ は半対角されるが優対角にならないので，このことに留意して
Gauss 消去法で解かなければならない．この計算法では対角化近似をしていないので収束性が良く，計算時
間はもとの LU-SGS 法に比べて一般にかなり少なくなる．
50
= g12 g23 ; g22 g13
g23
= g12 g13 ; g11 g23
2
g33 =J = g11 g22 ; (g12 )2
g11 g13 + g12 g23 + g13 g33 = 0
g13 =J 2
=J 2
76
圧縮性流れの解法|NS 方程式 主な 記号
添字 ^ または Fraktur 体文字
一般曲線座標系の方程式
添字 ~
一般曲線座標系の質量流束の方程式
添字
非保存形方程式
添字
Roe の近似リーマン解
式 (17.143)(17.155) 等
添字 + 正の固有値，正方向に伝播する波
添字
;
負の固有値，負方向に伝播する波
Ai ヤコビ行列 式 (17.69b)
Ai A^i ヤコビ行列 式 (17.79b)
A~i ヤコビ行列 式 (17.87b)
Ai 式 (17.70b)
Ai A^i 式 (17.80b)
A~i 式 (17.89b)
Ai 分離された Ai 式 (17.103)
Ai 分離された Ai 式 (17.106)
A~i 分離された A~i 式 (17.113)
B 式 (17.65c)
c = p RT 音速
cp cv 定圧比熱，定積比熱
D 拡散項 式 (17.51b)
D^ 拡散項 式 (17.64b)
D~ 拡散項 式 (17.65b)
DFj+1=2 F の高次化された差分 式 (17.117)(17.118)
Dqj+1=2 q の高次化された差分 式 (17.159)(17.160)
e = ( +u2=2) 単位体積当たりの岐点内部エネルギー
Fi 流束ベクトル 式 (17.51b)
Fi F^i 流束ベクトル 式 (17.64b)
F~i 流束ベクトル 式 (17.65b)
g 外力項 式 (17.51b)
g^ 外力項 式 (17.64b)
g~ 外力項 式 (17.65b)
gij 測度テンソル 式 (17.57)
gij 測度テンソル 式 (17.63)
H = h +u2=2 = (e + p)= 岐点エンタルピー
Hl+1=2 数値流束 式 (17.98)(17.116)
Hl+1=2 数値流束 式 (17.101)
(H1 )i+1=2 j k 数値流束 F1 のセル境界面 xi+1=2 の平均値
式 (17.144)
付 録 (H^ 1 )i+1=2 j k 数値流束 F1 のセル境界面 i+1=2 の平均値
h = + p= 比エンタルピー
I 単位行列
J ヤコビアン 式 (17.53)
Ki ^ i 固有値 ki の対角行列
ki ^i 固有値，圧力波などの位相速度の i 成分
Li = Li N 左固有ベクトルの行列 式 (17.74b)
L^ i = L^ i N 左固有ベクトルの行列 式 (17.82b)
L~ i = L~ i N~ 左固有ベクトルの行列 式 (17.91b)
M = ju j=c マッハ数
Mij = Li Rj = Li Rj 式 (17.164)(17.176)
M~ ij = L~ i R~j = L~ i R~j 式 (17.168)(17.181)
N 保存形から非保存形への変換の行列 式 (17.72)
N~ 保存形から非保存形への変換の行列 式 (17.88)
Pr = cp = プラントル数
p 静圧
q 未知変数ベクトル 式 (17.51b)
q q^ 未知変数ベクトル 式 (17.64b)
q~ 未知変数ベクトル 式 (17.65b)
qiL+1=2 セル境界 xi+1=2 の左側の q の平均値
qiR+1=2 セル境界 xi+1=2 の右側の q の平均値
qLi+1=2 セル境界 i+1=2 の左側の q の平均値
qRi+1=2 セル境界 i+1=2 の右側の q の平均値
q 未知変数ベクトル 式 (17.70b)
q q^ = q 未知変数ベクトル q~ 未知変数ベクトル 式 (17.89b)
q 熱流束
R 気体定数
Ri = N ;1 Ri 右固有ベクトルの行列 式 (17.74c)
R^i = N ;1 R^i 右固有ベクトルの行列 式 (17.82c)
R~i = N~ ;1 R~i 右固有ベクトルの行列 式 (17.91c)
R~ 付加項 式 (17.65b)
rhs 式 (17.93)
rg
hs 式 (17.95)
T 静温度
U 反変速度 式 (17.59)
u 流速
i = (@xk=@ i )uk = gij Uj
式 (17.156)
77
78
圧縮性流れの解法|NS 方程式 比熱比
~ = ;1
( Fi )j+1=2
( Fi )j+1=2
分離された流束差分 式 (17.104)
分離された流束差分 式 (17.107)
Mi 運動量の差分 式 (17.104b)
Mi 運動量の差分 式 (17.107b)
P 静圧の差分 式 (17.104b)
qj+1=2 q の空間差分
qj+1=2 q の空間差分
q q の修正値 式 (17.111)
q~ q~ の修正値 式 (17.115)
qj+1=2 q の空間差分 式 (17.102)
qj+1=2 = qj+1=2 q の空間差分
ij
クロネッカー
関数
比内部エネルギー
ijk
Eddington の
脚注 2
渦度
熱伝導率
i
固有値
i の対角行列 式 (17.74a)
^ i 固有値 ^i の対角行列
~ i = ^ i 固有値の行列
xi = t= xi
= t
式 (17.82a)
粘性係数
i j = @ i=@xj
粘性応力テンソル
密度
2
= ~u 2=2
rxi r i xi または
xi または
xi
i
i 方向の後退差分
i 方向の前進差分

新 しい単クロ ン抗体を用 いてのIL -2受容 - 東北大学

Document

Download

QRSK=AGHANSTAUVWEXYCZFGHA[O[OHJG=ATSQO=N

第8章 次元縮約 pLSIについて

4月号 - 玉野市教育委員会

粘性流体力学講義プリント-0

報告書_6/6 （PDF:1640KB）

スライド 1

A LOOPY BELIEF PROPAGATION APPROACH TO THE SHAPE