第 25 回数値流体力学シンポジウム
C02-4

圧縮性 CFD による低マッハ数流れ計算のための新しい陰的時間積分法
New Implicit Time Integration Method for Compressible CFD in Low Mach number○嶋英志，JAXA/JEDI，神奈川県相模原市由野台 3-1-1，shima.eiji@jaxa.jp
Eiji Shima, JAXA/JEDI, 3-1-1, Yoshinodai, Sagamihara, Kanagawa, Japan
By using all speed numerical flux schemes, such as SLAU or SD-SLAU, in MUSCL approach for compressible CFD, low Mach
number flows and sound propagations can be solved at the same time without loss of accuracy nor parameter tuning. In order for
efficient computation, this paper deals with implicit time integration methods, those are combined with SLAU, and theirtemporal accuracy for low Mach number flows and sound propagations. The influence of parameters in the implicit scheme,
such as the number of iterations in a time step and the cutoff Mach number, on the solution accuracy of low speed flows andsound propagations is studied.
1. はじめに
本研究は，非構造格子を含む，圧縮性 MUSCL (Monotone
Upstream-centered Schemes for Conservation Laws) 型（以下 MUSCL）の FVM（有限体積法）の枠組みの中で，低 Mach(<<1)数流れと音波の同時計算が可能な手法の構築を目的とする． MUSCL はいく
つかの，ほぼ独立の構成要素からなっており，その中の流束関数と陰解法（陰的時間積分法）の改良による実現を目指す． MUSCL
は複雑な現実的形状の解析手法として幅広く用いられており，この改良は，容易に様々な応用計算に適用可能と期待される．
また，これにより，音波が誘起する流れ場に起因する現象を直接計算することが可能となる．このような現象の理解は，液体ロケットエ
ンジン内のレゾネータでの消音機構の解明などに必要である．
圧縮性 CFD スキームを低 Mach 数流れに適用する場合，過大な
数値散逸による解像度の低下および，移流速度と音速の大きな隔たりによるスティフネスが問題となる．
（圧力等のスカラー量の変動が微小になることによる丸め誤差の問題は，基準量+変動量（例
えば p∞+p’）の形式で変数を記憶し，差分を変動量から計算する
ことで比較的容易に回避可能で，ここでは取り上げない．
）
筆者らは AUSM(Advection Upstream Splitting Method)族数値流束
スキームの改良によりパラメータ調整なしに低 Mach 数から高Mach 数まで安定かつ高解像度に計算できる全速度数値流束スキ
ーム SLAU(Simple Low-dissipation AUSM)[1-3]および超音速での特性を更に改良した SD(Shock Detecting)-SLAU[4,5]を提案した．また，これらの手法は時間微分前処理（以下 PC）[6-9]との併用により，低 Mach 数での定常解に関してはスティフネスも解決可能である
ことが分かっている．
PC は圧力場更新に局所時間を使うのと同等なので，そのままでは非定常計算に適用することが出来ない．DTS（二重時間法）の併用が知られているが，筆者らは，PC は陰解法における数値散逸の改良と等価であることを見出し，より計算量の少ない
TC-PGS(Time Consisten Preconditioned Gauss Seidel)を提案し[10,11]，
非定常計算のみならず，定常計算にも有効であることを示した．
音の伝播に関しては，時空間の離散化サイズが音波のスケールより十分小さければ，陰解法を用いても，十分な精度で解析可能である[10,11]が，TC-PGS を含む，PC を併用した方法は，PC に特有のパラメータであるカットオフ Mach 数<<1 の場合，音波の減衰が大きく，実用的ではないことも分かっている．
本研究では，TC-PGS について，より簡潔でパラメータ依存性の少ない新しい定式化を示し，更に，音響伝搬の問題を解決する
新しい手法を示し，各々の特性を数値例で示す．

2. 基礎方程式
保存系の圧縮性 NS(Navier-Stokes)方程式は積分形式で

 Q dv   (Eˆ  Rˆ ) ds  0
t

(2.1)

と書ける．ここで， Q ， Ê ， R̂ は保存変数，非粘性流束，粘性
流束のベクトルである． セル（検査体積）を多面体（２次元では多角形）とすると構造/非構造格子に共通な FVM の基礎式が得ら
れる．

1
1
~
~
Q i   (E i , j  R i , j ) s i , j  0
t
Vi j

(2.2)

Q i  Q n 1 i  Q n i(2.3)

ここで添え字 i，j は i 番目のセルの j 番目の境界面(辺)もしくは，その面で接する隣のセル，n は時間ステップ，V，s はセルの体積~ ~及び境界の面積，Δt は時間刻み幅を示す．E ，R は非粘性項と粘性項の数値流束である．なお，構造格子の場合は，単に，4 角形
あるいは 6 面体の検査体積に適用し，非構造格子と区別しない．
非定常解析のために DTS，3 点後退時間差分，記憶変数ベクト
ル S を導入する．Δτ,k は疑似時間の巾およびステップを示す．
k 1

k

 Si

 S k 1   2S i n  1   2 S i n1
{ 1 i
t
Si

k 1



1  S Vi  Q  i

~

 (E

k 1
i, j

(2.4)

~ k 1
 R i , j ) si , j }  0j

係数θは，時間刻み変化を伴う 3 点後退差分法（時間 2 次精度）の場合，下記であらわせる．

 r  2 r 1(1 ,  2 )  , r 1 r 
t n 1  rt n(2.5)

また，後退 Euler 法（時間一次精度）は次で与えられる．

(1 ,  2 )  (1,1)(2.6)

S は，丸め誤差の回避や，コーデイング上の都合で任意に選ぶことができるが，ここでは，非定常計算でも保存則を満足する，

Copyright © 2009 by JSFM

第 25 回数値流体力学シンポジウム
C02-4
S=Q，∂S/∂Q=I を選択する．
また，tとは独立にτは自由に選択でき，τに関して収束すれば，正しい時間発展が達成されるが，陰解法に，SGS(Symmetric
Gauss-Seidel)反復などの手法を用いると，Δτ が大なるほど収束が良い．そこで Δτ→∞とすると次式を得る．

  2Q i  1   2 Q i
t1
~ k 1 ~ k 1  ( E i , j  R i , j ) si , j  0
Vi j

1Q i

k 1

また， B̂ は次で定義される流束ベクトルの Jacobian である．

 Vn

 xn


ˆ   yn
B

 zn

 0


n 1

n

(2.7)

3. 線形近似並びに作業変数の導入
式(2.7)に対する陰解法は，左辺に一次風上法を用い線形化を行って，次のように書ける．

 

k

Hi 


c 2 z n

Vn

0

0

0

Vn

0

0

0

Vn

0

0

0

0

0


0

0

Vn 

(3.9)

(3.10)

ここで，xn,yn,zn は外向き法線単位ベクトルである．
陰解法全体の手続きは，次のように書ける．

 

1.
2.

k

At n physical time step, set k:=0, Qk:=Qn
Repeat until prescribed convergence level or iteration count;2.1. Compute RHS=  M i Hi k using Qk

1Q i k   2Q i n  1   2 Q i n 1
t

(3.2)

1
~ k ~ k (E i , j  R i , j ) s i , j
Vi j

3.

~
~
~
A  E / Q  R / Q
Q i  Q i

k 1

 Qi

k

(3.4)

 

~
はi
に関して収束すれば，Qn+1=Qk+1 が得られる．ここで， Ai, j

番目のセルの j 番目の境界面(辺)での正方向（i から j への方向）の数値流束の Jacobian である．
次に，Jacobian や陰的散逸の簡略化による計算量減少を狙って，
次式のエントロピー変数を導入する．

W  p, u , v , w, p  c 2  ) T

4. 数値散逸項の検討
計算量の観点から，数値散逸 D’を次のような対角行列とする．

(3.5)

dp

0
D'   000

ここで，ρ,p,u,v,w,c は密度，圧力，デカルト座標の速度，音速である．4 節で議論する非粘性項に対応するマトリックス数値散逸 D’を導入し，また，エントロピー変数に式(3.1)を変換すると次を得る．

1I 1
 D' ~     si , j   I   Wi
2
 i , j  t Vi j

(3.6)

ただし，M は次で定義される変換マトリックスである．
(3.7)

0

dv

0

0

0
0

dv
0

0
dv

0

0

0

0

0
0

0
d s 

(4.1)

d p | Vn | c(4.2)

d v  d s | Vn |  (1   )c

(4.3)



(4.4)



u v w
M  min 1.0,c2

ものを用いる近似を導入している．̂ は粘性項に対応する，次で2定義されるスペクトル半径である．

 iVi

0

  (1 M )2

また，式(3.1)から(3.6)の変形に当たっては，M を隣接するセルの(    T ) si , j

0

SLAU（Appendix 参照），SD-SLAU の数値散逸は近似的に次のよ
うにあらわせる．

ˆ  D'
B

1
k
j ,i
  si , j  ~I Wj  Μ i H i


2
Vi j


W  M Q

2.2. Solve Eq.(3.6) approximately to obtain ΔW
2.3. Compute Qk+1:=Qk+M-1ΔW
2.4. Set k=k+1
Next Time Step, Qn+1:=Qk+1, n:=n+1

陰解法に現れる物理量に関し定義場所に自由度があり，ここでは，セル i の更新に関わる係数には全てセル i のものを用いる．
時間精度や保存則の成立性は，Newton 反復の収束次第であるが，収束の程度に関わらず，物理的な時間発展とは無矛盾である．また，線型方程式式(3.6)は既に様々な近似を含んでいるので，仮に
厳密解が得られたとしても Newton 反復は必須である．したがって，全体的な計算時間短縮には，線形方程式解法の精度と計算量のバランスをとることが重要である．この手法の特徴はエントロ
ピー変数への変換と数値散逸の選択にある．次節以降，ステップ
2.2 の式(3.6)の解法に関して述べる．

(3.3)

これは，右辺の非線形方程式の残差を左辺の線形方程式で修正してゆく手順，あるいは Newton 反復とみなせる．上記反復が kˆ c 2 y nVn  uxn  vy n  wzn

1I 1~ 
1~    si , j A i , j  Q i   si , j A j ,i Q j (3.1)
Vi j t Vi j H ic 2 x n

2

2






(4.5)

一方，式(3.6)の行列全体が優対角となるために，係数行列が常に半正定値となる必要十分条件は次であらわされる．

(3.8)
2

Copyright © 2011 by JSFM

d p  Vn c

第 25 回数値流体力学シンポジウム
C02-4
用すると，非粘性の場合，式(3.6)への，次のような陰的数値散逸の適用と同等であることが導ける．[10,11]

(4.6)



d v  Vn   c(4.7)

d S  Vn

(4.8)

d v  d s | Vn |  (1   ' )c

(5.2)



'  (1 M' )2
M '  min 1.0, M 2  Mc 2



(6.3)

(6.4)

1   '2  1   'Vn / c
1  1   'Vn / c(6.5)

式(5.3-5)の定義を用いれば，式(6.5)もM2 のオーダーであり，前項の手法は，PC の特別なケースと見ることもできるのが，TC-PGS1 と呼ぶ理由である．ただし，普通，PC の定式化では εがセル毎に一定なのに対し，式(6.5)はセル境界ごとに異なるので完全に一致する訳ではない．
TC-PGS1 の以前のバージョン TC-PGS に対する改良は，数値散逸の定義が簡潔なこと，
8.1 で示すように，広い計算条件で Mc=M∞が最適に近く最適値の選択が容易なこと，があげられる．
7. 手法２：FGMRES の利用
よどみ点での無散逸を避け，SLAU の数値散逸と近いものとして次を用いる．

ここで Mc は零割と無散逸を防ぐためのパラメータで，M∞ (一様流 Mach 数)と同じオーダーとなるように選択する．これにより優対角性が維持されるので次の SGS 反復を用いることが出来る．

Wi
1d p | Vn | c

(7.1)

d v  d s | Vn |  (1   ' )c(7.2)

これを用いると係数行列が非優対角となる．
この線形方程式（非対称疎行列）の反復解法として，ここでは Saad による
FGMRES(Flexible Generalized Minimum RESidual method)[12]を用いる．リスタートパラメータを k，近似解のために内部反復を m にとどめた手法をここでは FGMRES(k,m)と書く．これを用いたス
テップ 2.2 は次のようにまとめられる．

(5.6)

右辺の逆行列は対角行列なのでブロック行列の反転は生じない．
SGS 反復を用いたステップ 2.2 は次のようにまとめられる．

2.2 Linear solver using FGMRES(k,m)
2.2.1 Define b≡[-MiHi], x≡[ΔWi]

2.2 Linear solver using SGS
2.2.1 SetΔWi:=0
2.22 Repeat until prescribed convergence level or iteration count;
2.2.2.1 Forward sweep of Eq.(5.6)
2.2.2.1 Backward sweep of Eq.(5.6)

 I 1 B j ,i  D' ~ 
1
 D'
 
 I W j 
Ax   1   s i , j   ~  Wi   s i , j 



t
V
V
2
2




j
j
i, j 
i
i





2.2.2 x0:=0
2.2.3 For l=1,…,m do
2.2.3.1 r0=b-Ax0, β:=|| r0||2 v1 := r0 /β
2.2.3.2 For j=1,…,k do;zj:= Pj-1vj

w:=Azj
For i=1,…j do {hi,j:=(w,vi), w:=w-hi,jvi}
hj+1,j=|| w ||2, vj+1=w/hj+1,j2.2.3.3 Define Zm:=[ z1…, zk]
where
ym=argminy||βe1-Hmy||2
2.2.3.4
x=x0+Zmym
e1=[1,0,…,0]T and Hm is the Hessenberg matrix consist of h
2.2.3.5 x0:=x

全体的な計算効率からは，下記で表わされる Nopt 回程度の往復反復が最適であることが分かっている．

Nopt  max( 1, min( 10 , N 1 / 10 ))

一方，式(5.1,2)による定義と式(6.1-3)は次の ε の定義の元で等価となることが導ける．

(5.5)

 I 1 
 D'
  1   si , j   ~I  i , j 
2 t Vi jˆ  D'
1Bk  si , j  j ,i
 ~I Wj  Μ i H i 2Vi j






1
2
(1   ) | Vn |  (  1) 2 Vn  4c 2
2

  Min(1, Max( M 2 , Mc 2 ))(5.4)

M  (u 2  v 2  v 2 ) / c 2(6.2)

ε は M2 のオーダーの関数であり通常，次のように定義される．

(5.3)



dv  dS  

5. 手法１：数値散逸の変更による SGS 法：TC-PGS1
優対角を維持し，よどみ点(u=v=w=0)での零割と無散逸を逃れるために，非粘性の数値散逸項を次のように修正する．
(5.1)

(6.1)

ここで σ は PC された系での次のようなスペクトル半径である．

ただし，α>0 である．式(4.2-5)で定義される数値散逸は優対角性を満足しないので，古典的な SGS 反復等を用いることはできない．
（優対角性がないとThomas法も不安定なのでAF（近似因数分解）法等でも同様の問題がある）式(3.6)の解は近似的で良いので，陰解法には，次の二つの選択が考えられる．
・手法 1：優対角となるよう数値散逸を変更し SGS を用いる．
・手法 2：式(4.2-5)の数値散逸のまま線型解法を工夫する．
次節以降で各々の手法を示す．

d p | Vn | c /(1   ' )d p   /(5.7)

ただし，N1/10 は内部反復の相対 L2 残差が 1/10 となる内部反復回数である．
6. TC-PGS1 と時間微分前処理法との関係
ここでは，前項で示した数値散逸は，時間微分前処理法と密接な関係があることを示す．Weiss & Smith の前処理法を陰解法に適
3Copyright © 2011 by JSFM,

第 25 回数値流体力学シンポジウム
C02-4
間は各々のケースで TC-PGS1 の Mc= M∞を基準としている．
この結果からは次のように言える．
(1) どちらの手法でも最適な Mc が選択できれば計算時間に大きな違いはない．
(2) TC-PGS1 では 4/6 ケースで Mc=M∞が最適で残りのケースでも最適に近い．したがって，Mc=M∞に固定するのが妥当な選択である．
(3) TC-PGS では最適値は Mc=2.0 M∞~6.0 M∞に分散している．
したがって，今後の TC-PGS1 の計算例では，Mc=M∞に固定する．

ここで P-1 は全体行列の近似逆行列である前処理行列であるが，PC での前処理とは全く別のものである．FGMRES の特徴は，
GMRES の各ステップで異なる前処理行列を用いることが出来ることにある．本研究では，5 節で定義した TC-PGS1 を前処理行列として用いる．また，FGMRES の特徴を用いて各 GMRES ステッ
プで異なる SGS 反復回数を用いる．ここでの，FGMRES 利用に関する留意点を記す．
(1) SGS 反復を用いることにより全体行列の定義なしに前処理が実施可能である．
(2)上記 2.2.3.2 では行列・ベクトル積 Ax が必要であるが，演算内
容はGS の1 スィープと殆ど同じなので，全体行列も不要である．
(3)いくつかの作業ベクトルの他に，Krylov 部分空間の定義のために 2k 組のベクトルが必要，即ち 3 次元 NS 計算では少なくとも，セル毎に 10k 変数が増加する．
(4)SGS 反復回数を固定すれば通常の GMRES が利用可能で，その場合，部分空間に必要なベクトルは k 組に減少する．
(5)局所最小値に落ち込むことにより，GMRES の収束が，停止することがある．これは全ての ym が 0 に近くなることで検知できるので，そのような場合には GMRES 探索を，TC-PGS1 に切り替える．これには， |ym |<<β の時，y1=β, yk=0(k≠1)とすることで実現できるので，計算量の増加はわずかである．
3.0

8.2 翼型周りの定常流れ
NACA0015 翼型周りの低 Mach 数（M∞=0.01）非粘性定常流
れを例に，FGMRES(k,m)の収束性に関する，パラメータ(k,m)の影
響および，TC-PGS1 との比較を示す．Fig.3 に FGMRES(1,1), (2,1),(4,1), (8,1), (16,1), (4,4)の運動量に関する L2 ノルムの収束を横軸時間ステップで示す．
k および m が増大するほど，収束も早まるが，(8,1)あたりで頭打ちになっていることが分かる．また，(16,1)と(4,4)は同程度の計算コストであり，収束特性も同等である．(k,m)
が増加するにつれ，
時間ステップあたりの計算量も増加するので，計算量の最適値が存在する．
計算量の目安として CPU 時間を用い，運動量と密度の収束履歴を Fig.4,5 に示す．運動量と密度で多少傾向は異なるが，(2,1),(4,1)あたりで最速の収束が得られる．このような小さな k 値が最適なのは，線形方程式が近似的なものであり線形方程式の収束と Newton 反復の収束が一致しないことによると考えられる．GMRES 系の手法では Krylov 部分空間のための記憶容量増加が大規模問題では負担となるが，この場合では，k≤4が最適なので，容量増加は最小限である．
次に，TC-PGS1 との比較を示す．Fig.6,7,8 に CPU 時間で見た運動，密度，抵抗係数の収束履歴を示す．図中の tcpgs1-Mc=1 は，PC を使わない高速用のスキームに対応する TC-PGS1-Mc=1 は明らかに収束が遅く，抵抗の収束では図からはみ出るので示されていない．変数によって状況は多少異なるが，このケースでは
TC-PGS1 と FGMRES(4,1)は，ほぼ同程度の収束特性で，両者は同程度の性能であると云える．

TC‐PGS1

Relative CPU Time

2.5
CIRC‐M=0.001

2.0

CIRC‐M=0.01

1.5

CIRC‐M=0.1

1.0

Rect‐0deg
Rect‐45deg

0.5

Rect‐30deg
0.0
0.0

2.0

4.0

6.0

8.0

10.0

12.0

Mc/M∞

Fig.1 Effect of cutoff Mach number to convergence for new
TC-PGS1
3.0

8.3 円柱周りの非定常剥離流れ
低 Mach 数非定常流れとして Kármán 渦列が発生する Reynolds
数=100，M∞=0.01 の円柱周りの流れを取り上げる．瞬時の圧力分
布の例をFig.9 に示す．
格子形状はO 型，
セル数は100(周方向)x150（半径方向），半径方向の最小格子間隔は 10-2R，計算領域半径は
外部境界の影響を避けるため 10000R と大きくとった．時間刻み
は移流速度に基づく Courant 数が約１となるように定め，その結果，音速ベースの Courant 数は，700 程度となっている．計算初期のインパルシブスタートが初期擾乱となるので各ケースで初期
300 ステップは共通のものを用いた．
時間精度の定性的な傾向をみるために，TC-PGS1 とFGMRES(4,1)で Newton 反復回数を変えた場合の揚力の履歴を示
す．Table.1 に各ケースの手法と相対的計算時間，Fig.10 に参照解と一致するケース，Fig.11 に一致しないケースを示す．なお参照
解はTC-PGS1を用い運動量残差が10-4 に低下するまで反復することで求めた．
ここで着目する流れ場の精度に関する限り，
TC-PGS1
で 4 回，FGMRES(4,1)で 2 回以上の反復が必要であり，その場合にはFGMRES(4,1)は約1.7倍程度のCPU時間を要していることが分かる．これは，揚力の履歴の見た目の比較と云う，定性的なものなので，細かい数値に意味はないが，TC-PGS1 と FGMRES(4,1)
が同程度の性能であるとはいえる．

TC‐PGS

Relative CPU Time

2.5
CIRC‐M=0.001

2.0

CIRC‐M=0.01

1.5

CIRC‐M=0.1

1.0

Rect‐0deg
Rect‐45deg

0.5

Rect‐30deg
0.0
0.0

2.0

4.0

6.0

8.0

10.0

12.0

Mc/M∞

Fig.2 Effect of cutoff Mach number to convergence for old
TC-PGS
8. 数値例
8.1 TC-PGS1 の最適 Mc の選択
PCではMcの選択が収束性に大きな影響を及ぼすことが知られているが，TC-PGS1 においても，状況は同じである．(ただし，
SLAU を用いる限り Mc の選択は右辺の数値散逸には影響しな
い ． )2 次 元 円 柱 (Re=100,M=0.001,0.01,0.1) と 2 次 元 角 注（Re=100,M=0.01,α=0°,30°,45°）周りの非定常流れに関して，
TC-PGS と TC-PGS1 への Mc の影響を Fig1,2 にまとめる．CPU 時
4

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C02-4
Steady Flow round NACA0015 M=0.01
0.1

FGMRES(1,1)
FGMRES(2,1)
FGMRES(4,1)
FGMRES(8,1)
FGMRES(16,1)
FGMRES(4,4)

tcpgs1
tcpgs1-Mc=1
FGMRES(4,1)

0.1
0.01

Momentum Residual

0.01

Momentum Residual

Steady Flow round NACA0015 M=0.01
1

0.001

0.0001

1e-05

0.001
0.0001
1e-05

1e-06

1e-06

1e-07
0

1e-07

200 400 600 800 1000 1200 1400 1600 1800 2000

0

20

40

STEP

Fig.3 Convergence history of momentum in physical time step forsteady inviscid flow around NACA0015 at M=0.01.

Momentum Residual

100

120

Steady Flow round NACA0015 M=0.01
1

FGMRES(1,1)
FGMRES(2,1)
FGMRES(4,1)
FGMRES(8,1)
FGMRES(16,1)
FGMRES(4,4)

0.01

80

Fig.6 Convergence history of momentum in CPU time for steadyinviscid flow around NACA0015 at M=0.01.

Steady Flow round NACA0015 M=0.01
0.1

60
CPU(sec)

tcpgs1
tcpgs1-Mc=1
FGMRES(4,1)

0.1
0.01
Density Residual

0.001

0.00010.001
0.0001
1e-05

1e-05

1e-06
1e-06

1e-07
1e-08

1e-07
0

200

400

600

800

1000

0

1200

20

40

CPU(sec)

Fig.4 Convergence history of momentum in CPU time for steadyinviscid flow around NACA0015 at M=0.01.

120

tcpgs1
tcpgs1-Mc=1
FGMRES(4,1)

0.00070.00060.001

0.0005CD

Density Residual

0.01

100

Steady Flow round NACA0015 M=0.01
0.0008

FGMRES(1,1)
FGMRES(2,1)
FGMRES(4,1)
FGMRES(8,1)
FGMRES(16,1)
FGMRES(4,4)

0.1

80

Fig.7 Convergence history of density in CPU time for steadyinviscid flow around NACA0015 at M=0.01.

Steady Flow round NACA0015 M=0.01
160
CPU(sec)

0.00010.00041e-05

0.00031e-06

0.00021e-07

0.0001

1e-08
0

200

400

600

800

1000

0

1200

0

CPU(sec)

10

20

30

40

50

CPU(sec)

Fig.5 Convergence history of density in CPU time for steadyinviscid flow around NACA0015 at M=0.01.

Fig.8 Convergence history of drag coefficient in CPU time forsteady inviscid flow around NACA0015 at M=0.01.
5

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C02-4
Unsteady Flow around cylynder M=0.01 Re=100
0.8

Refference
TCPGS 1
TCPGS 2
FGMRES 1

0.6

Lift Coefficient

0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
40

Fig.9 Instantaneous pressure distribution around a cylinder.

METHOD

1
2
4
1
2
4

0.21
0.44
1.10
1.00
1.71
3.53

Refference
TCPGS 3
FGMRES 2
FGMRES 3

0.6

Lift Coefficient

0.4
0.2
0

-0.2
-0.4
-0.6
-0.8
40

50

60

70

80

90

100

110

80

90

100

110

120

8.4 1 次元音波の伝搬
低 Mach 数流れ中での音の伝播の基礎として，陰解法で Courant
数>>1 を用いた場合の静止媒質中の 1 次元音波の解析例を示す．
以前の研究で次のようなことが分かっている．
① Newton 反復が収束すれば，Courant 数>>1 の陰解法でも十分な精度で計算可能．
② 2 次精度法では空間時間の離散化幅が音波の波長，周期の1/40 以下であることが必要．
③ Mc=1 の高速用の陰解法であれば，4 回程度の内部反復で実用的な収束が得られ，ほぼ，無視できる振幅・位相誤差で計算可能である．
④ TC-PGS の Mc<0.03 では 20Newton 反復でも減衰が大きく，
より多数回の反復が必要で計算負荷が大きい．
ここでは，Mc=0.01 とし，
条件②（Courant 数1 で40PPW，
Courant
数 10 で 400PPW）を満足するように，400PPW，Courant 数=10，とする． x=0 で音波の圧力・密度・速度を与え準定常状態になるまで計算を行った．x>20 ではメッシュを徐々に拡大することで音波を減衰させ境界は十分遠くに取り反射の影響から逃れている．
Table.2 に各ケースの手法及び，相対的なCPU 時間を示す．
Fig.12に TC-PGS1 各ケース，
Fig.13 に FGMRES(4,1)各ケースの圧力分布を示す．厳密解とほぼ一致する解を得るには，TC-PGS1 では 300
回以上の反復を要するのに対し，FGMRES(4,1)では 8 回で十分である．CPU 時間で見ると，このケースでは，FGMRES は TC-PGS1に比較して 10 倍以上高速であると云える．前項で示したように，
FGMRES(4,1)では流れ場に関しても，2 回程度の反復でほぼ正確な時間発展を捉えられるので，同時に音波を考慮しても，8 回程度の反復で済むことが期待される．したがって，低 Mach 数流れ
場と同時に音の伝搬を計算する場合には，FGMRES(4,1)はTC-PGS1 より大幅に(10 倍程度)計算効率が良いことが期待される．

Unsteady Flow around cylynder M=0.01 Re=100
0.8

70

Fig.11 Histories of lift coefficient for poorly agreed cases.

Non-linear Iteration Relative CPU TimeTCPGS 1 TCPGS1TCPGS 2 TCPGS1TCPGS 3 TCPGS1
FGMRES 1 FGMRES(4,1)
FGMRES 2 FGMRES(4,1)
FGMRES 3 FGMRES(4,1)

60

Nondimensional Time

Table.1 Method and CPU time for each case
CASE

50

120

Nondimensional Time

Fig.10 Histories of lift coefficient for well agreed cases.

Table.2 Method and CPU time for each case
CASE
METHOD
Non-linear Iteration Relative CPU TimeTCPGS 1 TCPGS18
1.28
TCPGS 2 TCPGS1
100
16.66
TCPGS 3 TCPGS1
200
33.69
TCPGS 4 TCPGS1
300
50.89
FGMRES 1 FGMRES(4,1)
2
1.00
FGMRES 2 FGMRES(4,1)
4
1.95
FGMRES 3 FGMRES(4,1)
8
3.91

6

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C02-4
Unsteady Flow around cylynder M=0.01 Re=100

1D sound by TCPGS1 400PPW CFL=10 Mcutoff=0.01
0.0015

Reference
TCPGS 1
TCPGS 2
TCPGS 3
TCPGS 4

tcpgs1
FGMRES(4,1)
FGMRES(16,1)
FGMRES(4,4)

1e-06
1e-07

Density Residual

0.001

Pressure

1e-05

0.0005

0

1e-08
1e-09
1e-10
1e-11
1e-12

-0.0005

1e-13
1e-14

-0.001

0

0

5

10

15

20

30

40

50

60

70

80

90

100

Fig.15 Residual histories of density are shown in every non-linear
evaluation step.

Fig.12 Pressure distribution of 1D sound by TC-PGS1
1D sound by FGMRES 400PPW CFL=10 Mcutoff=0.01

Unsteady Flow around cylynder M=0.01 Re=100
1e-05

Reference
FGMRES 1
2
FGMRES 3
FGMRES 3

0.001

20

STEP

X

0.001510

tcpgs1
FGMRES(4,1)
FGMRES(16,1)
FGMRES(4,4)

1e-06
1e-07
Momentum Residual

1e-08

Pressure

0.0005

1e-09
1e-10

0

1e-11
1e-12

-0.0005

1e-13-0.001

1e-14

0

5

10

15

0

20

0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
CPU(sec)

X

Fig.13 Pressure distribution of 1D sound by FGMRES(4,1)

Fig.16 Residual histories of momentum are shown in CPU time.

Unsteady Flow around cylynder M=0.01 Re=100
1e-05

Unsteady Flow around cylynder M=0.01 Re=100

tcpgs1
FGMRES(4,1)
FGMRES(16,1)
FGMRES(4,4)

1e-06
1e-07

1e-05

tcpgs1
FGMRES(4,1)
FGMRES(16,1)
FGMRES(4,4)

1e-06

Momentum Residual

1e-07

1e-08

Density Residual

1e-08

1e-09

1e-09

1e-10

1e-10

1e-11

1e-11

1e-12
1e-12

1e-13
1e-13

1e-14
0

10

20

30

40

50

60

70

80

90

100

1e-14
0

STEP

Fig.14 Residual histories of momentum are shown in every
non-linear evaluation step.

0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
CPU(sec)

Fig.17 Residual histories of densiyu are shown in CPU time.
7

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C02-4
流れ場だけに興味がある場合には，簡単でメモリ消費も少ない
TC-PGS1 で十分である．
一方，音波との同時解析においては，FGMRES の有効性は明らかであり，TC-PGS1 による流れ場だけの解析と，あまり変わらな
い計算時間で同時に音波を解析可能で，このような応用には有望である．
今後の課題・展開として以下が挙げられる．
(1) 非対称疎行列の反復解法として FGMRES を用いたが，その他にも様々な選択肢があり，計算効率の更なる向上の可能性がある．
(2) 反復法に用いる前処理行列として，TC-PGS1 を，そのまま用い，充分有効ではあったが，改良の余地が予測される．
(3) 大規模な計算には並列計算が必須であるが，前処理は領域分
割により性能が劣化することがあり調査の必要がある．
(4) 本研究では空間離散化手法として2 次精度MUSCL を用いたが，陰的数値散逸の改良のアイデアは，保存系の圧縮性解法
であれば，その他の方法（高次精度差分法など）にも，ほとんどそのまま適用できるはずである．また，一旦エントロピ
ー変数へ変換してしまえば，元の独立変数は，なんでも良いので，適用範囲は保存系の定式化に限定されない．
(5) 低マッハ数流れの他の計算手法としては，非圧縮性解法を源流とする，圧力ベースソルバーが存在する．本研究で提案したような密度ベースとも呼ばれる圧縮性ソルバーとの関係を整理し，特性を比較する．
(6) 本研究の目的である空力音響問題に適用し有効性や適用限界を調査する．

8.5 FGMRES と TC-PGS1 の特性の分析
前項，8.3,8.4 で示したように流れ場計算については，
FGMRES(4,1)と TC-PGS1 は，同程度の効率なのに対し，音波の計算では大幅に FGMRES(4,1)が早い．この分析のために，8.3 の円
柱周り流れのケースについて，非定常計算での残差収束の特性を示す．１物理時間ステップ内の Newton 反復を 20 回に固定した場
合の，運動量，密度の残差履歴を，Fig.10-13 に示す．Fig.10,11 では横軸を反復回数とし，Fig.12,13 では CPU 時間としている．
Fig.10,11 から(k,m)を増すほどステップ毎の収束特性は改善され，
TC-PGS1 との比較では，密度残差の改善が顕著であることが分かる．また(16,1)と(4,4)は同程度の特性を示している．
Fig.12,13 の CPU 基準でみると，運動量と密度では傾向が異なることが明らかである．運動量では TC-PGS1 が最速で 10-3 程度まで低下している．一方，FGMRES では 1 ステップの計算量との関係で，(4,1)が最良の収束特性を示すが，TC-PGS1 に比較するとやや遅い．一方，.Fig.13 に示される密度残差では，TC-PGS1 の残差が落ちていないことが分かる．
この性質は定式化からも説明できる．Eq.(5,1)(6,1)の比較から
TC-PGS1 では圧力に関する陰的数値散逸がほぼ 1/(1-χ)倍で，Mc=0.01 では約 50 倍と極めて大きくなっているので収束が遅く，圧力と密接な関係がある密度の収束も遅い．
（ここには示さないが密度と圧力の収束はほぼ同じである．
）
一方，10-3 程度の運動量の収束が得られれば非定常な流れ場は，
ほぼ正確に計算できるので，流れ場に関しては TC-PGS1 が早い．
また，FGMRES は，速度，圧力の双方に適切な大きさの数値散逸を使うことが出来るので，両変数とも類似の良好な収束特性を示す．

謝辞
反復解法に関し有益なご助言をいただいた，九州大学情報基盤
研究開発センターの藤野清次先生に感謝いたします．また，JAXA
宇宙科学研究所の野々村拓氏，河合宗司氏，東京大学生産技術研
究所の飯塚宣行氏との議論に感謝します．

9. まとめと今後の課題
本稿の内容を次のようにまとめる．
(1) MUSCL 型の圧縮性 FVM 手法において，全速度数値流束関
数 SLAU との組み合わせによる，低 Mach 数流れと音波の同時解析を目指して陰解法の検討を行った．
(2) 陰的時間積分に現れる数値散逸を検討し，SLAU に対応する
数値散逸では係数行列に優対角性がなく SGS 等の古典的な線形解法は利用不可能であることを示した．
(3) これを解決する一方法が，陰的数値散逸の修正によって優対角性を回復する手法であり，この手法を TC-PGS1 として定式化した．エントロピー変数を作業変数として導入することと特殊な数値散逸を用いる点が特徴である．
TC-PGS1 は時間微分前処理法との類似性を示すことが出来るが，定式化はより簡潔で最適カットオフ Mach 数の選択が容易という利点を有する．
(4) もうひとつの方法は，非優対角のまま，近代的な反復線形解法を用いる手法で，ここではSaad によるFGMRES(k)を用い，収束加速のための前処理行列としては，TC-PGS1 の SGS 反復を用いた．
(5) FGMRES(k)は Krylov 部分空間での探索ベクトルのために，余分の記憶領域が必要であるが，本研究の範囲では，k=4 程度が全体的な計算時間最小化に最適であり，記憶領域の増加は最小限であることが分かった．
(6) 低 Mach 数流れ場との同時計算を意識した，
カットオフ Mach
数=0.01 での 1 次元音波の計算では，全体の計算時間で
FGMRES(4)は TC-PGS1 に比較して 10 倍以上高速であった．
(7) 定常，非定常の流れ場の計算に関しては，FGMRES(4)とTC-PGS1 は同程度の計算効率であった．
したがって，音波を含まない，低 Mach 数流れ場だけの解析に関しては，FGMRES(4)と TC-PGS1 は同程度の性能であるので，参考文献[1] 嶋英志，「AUSM 族全速度スキームと前処理付き陰解法」，
第 22 回数値流体力学シンポジウム講演論文集(2008)[2] E.Shima and K.Kitamura, “On New Simple Low-Dissipation
Scheme of AUSM-Family for All Speeds”, AIAA-Paper
2009-136, 2009[3] E.Shima and K.Kitamura, ” Parameter-Free Simple
Low-Dissipation AUSM-Family Scheme for All Speeds”,
AIAA Journal 2011, 0001-1452 vol.49 no.8 (1693-1709)(2011)[4] E.Shima and K.Kitamura, ”On AUSM-Family Scheme for All
Speeds with Shock Detection for Carbuncle-Fix”,
AIAA-Paper 2009-3544,2009[5] 嶋英志，北村圭一，「全速度数値流束スキームの改良について」，流体力学講演会講演論文集(2010)[6] E.Turkel, “Preconditioned methods for solving theincompressible and low speed compressible equations”, J.
Comp.Phys., Vol 72, Issue 2, (1987)[7] J.M. Weiss and W.A. Smith， Preconditioning Applied to
Variable and Constant Density Flows， AIAA J.， Vol. 33，
No.11， pp.2050-2057 (1995)

8

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C02-4[8] D. Unrau and D.W. Zingg, “Viscous Airfoil Computations using
Local Preconditioning， AIAA Paper 90-2027 (1997)Vn [9] X. Lv, Y. Zhao, X.Y. Huang, G.H. Xia, Z.J. Wang, “An
efficient parallel/unstructured-multigrid preconditioned implicit
method for simulating 3D unsteady compressible flows withmoving objects”, Journal of Computational Physics 215 (2006)  Vn     Vn q  q   q (A.9)

g  g   g   0, 1;   g  min max M ,0 ,1[10] 嶋英志，北村圭一，「全速度スキームと前処理付陰解法を用いた空力音解析手法について」，第 24 回数値流体力学
シンポジウム講演論文集(2010)g    max min M  , 0 ,1[11] E.Shima and K.Kitamura, ”CFD Method for Aero-acoustics
using All-speed Numerical Flux and Preconditioned Implicit
Time Integration”, AIAA-2011-3045(2011)  (1  M ) 2







(A.11)2222

1 u   v   u   v  M  min 1.0,


2
c



(A.12)

ここで c は左右音速の算術平均である．上記の関数 g は強い膨張波への対策なので低 Mach 数では，簡単に g=0 としてもよい．また，平均圧力は Mach 数を M として下記で与えられる．

Appendix A. 右辺の評価
粘性項は中心差分によって評価する．非粘性項は MUSCL(Monotone Upwind Scheme for Conservation Laws) の流儀に従い，該当セルとその周辺のセルの値を用いて．中心での基本変数の勾配を計算し，それを用いて計算されたセル境界での値を求める．この一般には跳びを含む，セル境界の左(+)右(-)の物理量（Q+,Q-）と左から右方向の法線ベクトル N を用い，AUSM (Advection
Upstream Splitting Method) 族の全速度流束関数である SLAU
(Simple Low-dissipation AUSM) により求める．なお，低速では
SLAU と SD-SLAU (Shock Detecting SLAU) の挙動は共通である．

(A.10)

[12] Y.Saad, “A flexible inner-outer preconditioned GMRES
algorithm”, SIAM Journal on Scientific Computing, (1993)~ ~F  F Q , Q , N(A.8)

  

p   p     ~p(p  p )22p  p 1         1
2







2
 1



 4 2  M M 1 , M 1 1
 (1  sign (  M  )) , otherwise
2


(A.1)

M 

~
ここで F 保存変数空間でのセル表面垂直外向きの数値流速であ(A.13)(A.14)

Vnu  xn  v  yncc

(A.15)

り，AUSM 族の数値流束関数は次のように書ける．

~ m  m  m  m  ~FΦ Φ  pN22Φ  1, u,v, w, h 

(A.2)

T

(A.3)

Ν  0, xn , y n , z n ,0

T

(A.4)

h  (e  p ) / (A.5)
ここで，  u,v, w はデカルト座標の速度，  xn , y n , z n  はセル表
面の外向き法線ベクトル，e,p, ρ,h は密度，体積当たりの全エネル
ギー，圧力，密度，エンタルピー，である．SLAU の質量流束は次で与えられる．

1m  {( Vn )   ( Vn )   Vn }(1  g )
2




2c

(A.6)

p

Vn  xn u  y n v  z n w

(A.7)

9

Copyright © 2011 by JSFM