第 24 回数値流体力学シンポジウム
C3-3非構造格子を用いた all-speed 流れ解析
Numerical simulation of all-speed flows using unstructured grids○ 矢田 和之, IHI 解析技術部, 横浜市磯子区新中原町 1 番地, E-mail:kazuyuki_yada@ihi.co.jp
市東 素明, IHI 解析技術部, 横浜市磯子区新中原町 1 番地, E-mail:motoaki_shito@ihi.co.jp
孫 明宇, 東北大, 仙台市青葉区荒巻字青葉 6-3, E-mail:sun@cir.tohoku.ac.jp
Kazuyuki Yada, IHI, 1, Shin-Nakahara-Cho, Isogo-ku, Yokohama
Motoaki Shito, IHI, 1, Shin-Nakahara-Cho, Isogo-ku, Yokohama
Mingyu Sun, Tohoku University, 6-3, Aoba, Aramaki, Aoba-ku, Sendai
An unsteady compressible/incompressible flow solver has been developed for structured and unstructured grid basedon the Lagrangian-remapping approach that can be readily extended for the interface-tracking method for two-phaseflow simulations. In order to apply the present Riemann-solver-based compressible solver for the simulation of thelow-Mach number flows, the low-Mach fix developed in the AUSM scheme has been adopted in solving theLagrangian step. The solver gains a good convergence rate for flow Mach number ranging from M=0.003 and over.
In this paper, efficiency of several popular iteration methods is evaluated in solving the implicit Lagrangian step forunstructured grids.

１． はじめに
著者らは、混相流を対象として、沸騰・凝縮やキャビテーション現象等を精度よく解析できるような流体解析手法の確立を目指している。いままでに圧縮性 2 相流問題を陽解法で解析可能にしたが(1)、低マッハ数に適用する際に、計算時間の増大が問題となっていた。そこで低マッハ数から高マッハ数まで扱える all-speed 流体解析手法の開発を行っている。
昨今、圧縮性流体と非圧縮性流体を統一的に扱う種々の流体解析ソルバーが開発されている。例えば AUSM スキームの改良型(2)（AUSUM+-up スキーム）やプリコンディショニング法(3)などが圧縮性流体からのアプローチで開発されてきた。一方、非圧縮性流体からのアプローチとしては低マッハ数近似手法が用いられてきた。
本研究では、semi-Lagrangian 法(4)に基づく圧縮性ベースの解法に対し、AUSUM+-up スキームで用いられる all-speed 化の手法を適用する新しい手法を開発した。Lagrange 解法の部分には陰解法を用いるため、非構造格子では特に効率的な行列解法が必要となる。本報では、マッハ数 1.0、0.1、0.01 の場合の翼型まわりの流れを題材として、行列解法（CG 法(5)、BiCG 法(5)、GMRES 法(6)）の効率を検討した。

方程式の陰解法と④マッピングの CFL 条件を組み合わせることから semi-implicit となる。
start
初期値設定
①Euler座標のu,v,pをLagrange座標のu,v,pに代入②Lagrange方程式を陰解法で解く③Lagrange座標の保存量を算出④Euler座標に保存量をマッピング

time<time_max?

YN
end

Fig.1 Flow chart of present semi-implicit method
2.2 基礎方程式（Lagrange 方程式）
2 次元非粘性流体に対する Lagrange 座標での質量保存式、運動量保存式を以下に示す。

２．計算手法rDρ= − ρ∇ ⋅ uDtr
Du
∇p=−Dtρ

2.1 semi-Lagrangian 法と semi-impicit 法
semi-Lagrangian 法は、Lagrange 座標で基礎方程式を解き、Euler 座標にマッピングする手法であり、混相流への拡張を念頭においた場合、異種媒体界面の計算等に高い精度を期待できる。Fig.1 に semi-Lagrangian 法のフローチャートを示す。①
まず最初に Eluer 座標に初期値 u , v, p を与え、これを同時にLagurange 方程式の初期値とする。ここで u, v, p は x 方向の速度、y 方向の速度、圧力である。②Lagrange 方程式を acoustic
linearized Riemann solver で近似し AUSUM+-up スキームの手法を用いて all-speed 化したものを陰解法で解く。③Lagrange座標での保存量を算出する。④Lagrange 座標での保存量をEluer 座標にマッピングし Euler 座標での値とする。マッピン
グに対しては CFL 条件を課す。以上で 1 つの時間ステップの計算は終了し、次ステップへ進む。時間積分は、②Lagrange

（１）

（２）

(1)式に等エントロピ関係式を用いると(3)式となる。

r
Dp= − ρc 2 ∇ ⋅ uDt（３）r
ここで t は時間、 ρ は密度、u は速度ベクトル（x 方向の速度は u, y 方向の速度は v）、p は圧力、c は音速である。(1)式は(2)式、(3)式を導出する際に用いているので、(2)式、(3)式の u、v、p の 3 変数を連立して解く。有限体積法を用い、i1

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C3-3番目の検査体積内の質量を mi,、Lagrange 座標での粒子軌跡の界面（こでは粒子界面と呼ぶ。2.4 節参照）での垂直ベクト
ルを(Sxij, Syij）とすると、(2)式、(3)式は以下の式となる。

mi

Dui*= −∑ ( pi* ) n+1 S xij
Dt

（４）

mi

Dvi*
= −∑ ( pi* ) n+1 S yij
Dt（５）(Dpi*I2= − i ∑ (ui* ) n+1 S xij + (vi* ) n+1 S yij
Dtmi

2.4 Lagrange 座標の保存量への変換
Lagrange 座標で解いた u*、v*、p*を用いることで粒子界面での流束を算出し、有限体積法により、Lagrange 座標での保存量（密度、運動量、エネルギー）を算出する。これを 1 次元でみた概略
図を Fig.2 に示す。Fig.2 の緑色領域が Lagrange 座標での検査体積となる。粒子界面とは Fig.2 の赤の点線部分となる。密度に対しては Lagrange 座標では質量は時間変化しないので、体積変化の流束から求まる。Lagrange 座標での運動量は圧力変化の流束から求まる。エネルギーは p*u*の流束から求まる。ここで Lagrange 座標での保存量は Euler 座標重心の値と等しいと仮定した。

) （６）u*i-1, v*i-1, p*i-1

ここで n は時間ステップを表し、n+1 は次の時間ステップ

u*i, v*i, p*ixを表す。また、Ii は音響インピーダンス ρc で、前のステップの値(n ステップ）で近似する。添え字*は Lagrange 座標を表す。(4)-(6)式右辺は、粒子界面に対して有限体積法を適用することで、Lagrange 座標での積分となる。実質上、Lagrange解法での体積変化は Geometric Conservation Law(7)を満たしている。

Lagrangian Interface
Lagrangian Interface

Δtxxi-1

2.3 Lagrange 方程式の all-speed 化
Lagrangian ステップでは、(4)-(6)式を陰解法で解く。ここでLagrange 座標での u*、v*、p*は、Riemann 問題における接触不連続面での状態量である。有限振幅の等エントロピ波の関係式（式(7.1)、(7.2)）を用いると式(8)の acoustic linearized
Riemann solver(8,9)が得られる。添え字 L は波の界面の左側の物理量、R は右側の物理量である。U は粒子界面に垂直な速度である。

2.5 Lagrange 座標の保存量から Euler 座標へのマッピング
その後、Euler 座標にマッピングする。Lagrange 座標の速度は求まっているので、検査体積界面から移流する流束は求まる。これ
を 1 次元で見た概略図を Fig.3 に示す。青の領域は Lagrange 座標の速度 ui で検査体積 i から流出する保存量 Qi+1/2 であり、赤の領域は Lagrange 座標の速度 ui-1 で検査体積 i に流入する保存量 Qi-1/2 である。Qi-1/2、Qi+1/2 は検査体積界面での物理量を表し、検査体積重心の保存量 Q より補間する。その際、Q は移流方向により検査体積界面の右側と左側のどちらの Q を採択するかが決定される。これで検査体積界面での流束が求まり、再び有限体積法により、
Euler 座標での保存量は求めることが出来る。マッピングを隣接するセルに行う必要性があるため、流速基準の CFL が妥当な判定であると考えた。求まった Euler 座標での保存量は、次の時間ステップの Lagrange 座標で解く初期値として使用し、非定常問題を解くことになる。ここで仮定として、Lagrange 座標と Euler 座標の勾配は等しいとし、minmod limiter を採用し、MUSCL 法により空間 2 次精度を確保した。本解法の詳細は別途発表を予定しているのでこれを参考にしていただきたい。

（7.1）

後退波

dp + ρcdu = dp + Idu = 0（7.2）(U *) n = U + β1 ( pL − p R ) /( I L + I R )( p*) n = p + β 2 s (U L − U R )U ≡ (I LU L + I RU R ) / (I L + I R )（８）p ≡ (I L p R + I R p L ) / (I L + I R )s ≡ (I L I R ) / (I L + I R )ui-1ΔtQi-1/2

ここで、AUSM+-up スキーム(2)において低マッハ数に拡張する all-speed 化の手法を用い、修正係数β1、β2（式(9)）を導入する。M は上流マッハ数である。

β1 = max(1, 1 / β 2 )
β 2 = min(1, M (2 − M ) )Eulerian Interface
xi+1/2

Fig.2 Flux integral of Lagrangian coordinate

前進波

dp − ρcdu = dp − Idu = 0

Eulerian Interface
xi-1/2xi

Qi+1/2

Qi

uiΔt

Qi-1（９）x

修正係数β1、β2 は、圧縮性解法を低マッハ数に適用する際に圧力振動などの不安定性を回避する修正係数である。

xi-1

β1 = 1, β 2 = 1 のとき、厳密な acoustic linearized RiemannInterfacexi-1/2xiInterface
xi+1/2

Fig.3 Mapping from Lagrangian coordinate to Eulerian coordinatesolver となる。
2Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C3-32.6 陰解法
陰解法は反復法として、共役勾配法(Conjugate Gradient method：CG 法)、双共役勾配法（Bi-Conjugate Gradient method：BiCG 法）、
Generalized Minimal Residual algorithm：GMRES 法のいずれかを採用し、時間積分は Backward Euler 法とした。

( )r r*r r*となり、しかも ri , ri ≠ 0 および (Ap , p ) ≠ 0 の直交条件を満iirr*
たす 2 つのベクトル p k ,および p k （双対共役基底）をつくり、CG 法と同様に計算する。

rrv
rk +1 = rk − α k Apk ,2.7 データ保存点
検査体積はセル節点法で、検査体積重心に保存量をもつ手法（検査体積中心法）で計算を行った。この手法の略図を Fig.4 に示す。
Fig.4 で青の領域はセル節点法と同様の方法で算出する検査体積であり、セル重心とエッジの中点を結んだ領域である。このことで壁面での境界条件の扱いが、セル中心法と同様に扱うことができ、かつ勾配ベクトルを精度よく算出可能である。

(rrrrk*+1 = rk* − α k AT pk*)r r*
直交条件 rk +1 , rk = 0 より上式を満たすαk を決定する。さらr*に、 p k +1 に共役勾配法を適用して次式を定義し、β k を

r r*(Ap , pData storage point（present control volume center method）

kk +1) = 0 （双対直交性）から求める。

rrr
pk*+1 = rk +1 + β k pk*r
収束判定は残差ベクトルのノルムが閾値 ε とベクトル b のノルムの積よりも小さいときに収束とする。
iii) GMRES 法
r rr
Ax = b の近似解はベクトル z を用いて以下の式とする。

r r rx = x0 + z
Data storage point（conventional cell vertex method）

GMRES 法は非対称行列 A に対して、以下のノルムが最小になrるようにベクトル z を決定する。ここで残差ベクトルvrrrr0 = b − Ax0 であり x0 は解ベクトルの初期値である。

Fig.4 Data storage pointrr rr
min b − A( x0 + z ) = min r0 − Az

2.8 反復解法(5, 6)i) CG 法r r
Ax = b の解を求める際に、以下の汎関数の最小値問題に帰着（12）r
ベ ク ト ル z は 、 ク リ ロ フ 部 分 空 間 Kk に 属 す る 。

rz∈K ,させて解く。

r rrr rF ( x ) = ( x , Ax ) − 2(b , x )

（10）

rr
r
xk +1 = xk + α k pkk{r r 2rm − 1r
r , Ar , A r , K , A
r
000r r rz = Vk y k ,（11）rT r
H k ≡ V AVk ,rr(p , Ap ) = 0 （直交条件）を満たすように定める。

kjr−1 r ry k = H k r e10r
Te1 ≡ (1,0,0, K ,0)r
このVk と Hk には以下の関係が成立する。

r rr
収束判定は残差ベクトル r = b − Ax がある閾値より小さい場合に収束とする。

AVk = Vk +1H kii) BiCG 法
対称でない行列 A を係数とする連立一次方程式に対して、CG法を修正したものである。本解析では、行列 A を不完全 LU 分解（13）式(13)の関係式を用いると、式(12)の最小値は以下のようになる。

rrrmin( J ( y )) ≡ min βe1 − H k y（14）rここでβ= r0 である。

r r*r r*で前処理を行った。i≠j のとき、 ri , r j = 0 および (Ap , p ) = 0i}r
ここで Hk（ヘッセンベルグ行列）、 e1 は以下の定義である。

rr
F( x )が最小となるように α k を求める。また p k は i≠j のとき( )0rrr
ここではベクトル z を直接解かず、 z の正規直交ベクトルVkrを用いる。これを用いるとベクトル z は以下のように表される。

ここで括弧は内積を表す。行列 A は対称行列である。本解析でrは、不完全 LU 分解で前処理を実施した。解ベクトル x は、スカr
ラー値 α k とベクトル p k を用いて以下の式を反復させて解く。

iK ≡ spankjr
したがって、GMRES 法は、J ( y ) を最小にすることで解ベクトル

3Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C3-3r rrrx = x + V y を求める手法である。
k0kk

本解析では、GMRES 法にも前処理として不完全 LU 分解を実施した。
本報では、上記 3 つの反復解法を用いて、計算精度や残差の収束性および計算時間を比較して、各反復解法が semi-implicit 法に対してどのような特性を持つかを調査する。
2.9 計算条件
翼型は NACA4412 とし、迎え角は 0[deg], 翼弦長は 1[m]とした。四角形格子で格子点数は 15120 点である。計算格子を Fig.5、6 に示す。格子データは非構造形式となっている。
流体は空気とし、流入条件の圧力は大気圧、温度は 300[K]である。ほぼ定常解に達するまで非定常計算を実施した。計算はマッ
ハ数 1.0 では 0.2[s]、0.1 では 2[s]、0.01 では 30[s]まで計算した。
今回は 2 次元非粘性計算を実施した。境界条件は、翼壁面は非粘性壁とし、外部境界は、低マッハ数領域では速度規定条件、出口圧力規定条件を設定し、M=1.0 では勾配ゼロ条件を課した。

Fig.7 Pressure contour at M=1.0
3.2 局所マッハ数コンター図（M=1.0）
流入マッハ数1.0 のときの局所マッハ数コンター図をFig.8 に示す。圧力低下に伴いマッハ数が上昇しているのが分かる。

Fig.8 Mach number contour at M=1.03.3 圧力係数
マッハ数 1.0 のときの 0.2[s]後（無次元時間 10）の圧力係数をFig.9、マッハ数 0.1 のときの 2[s]後（無次元時間 100）の圧力係数
を Fig.10、マッハ数 0.01 のときの 30[s]後（無次元時間 1500）の圧力係数を Fig.11 に示す。従来手法との比較として、市販ソフト（fluent）を用いた結果も示す。従来法は二次精度 Roe スキームを用いた preconditioning 法であり、本解析と同じ時間まで非定常計算を行った。
いずれの結果も、従来法との良い一致を示している。
以上より本解析手法により、all-speed 解析が可能であることを示すことが出来た。

Fig.5 Computational grid

1.5

Fig. 6 Computational grid（detail）

present(GMRES)

1
Cp [-]

３．計算結果マッハ数 M=1.0、0.1、0.01 において、GMRES 法、BiCG 法、CG 法の計算結果は、同様であった。したがって 3.1、3.2、3.3 節
は GMRES 法を代表して計算結果を示す。

commercial software

0.5000.2

0.40.60.8

1-0.5

3.1 圧力コンター図（M=1.0）
マッハ数 1.0 のときの圧力コンター図を Fig.7 に示す。単位は[Pa]である。翼前縁付近と後縁から弱い衝撃波がでていることが確認できる。

-1x [m]Fig.9 Pressure coefficient at M=1.04

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C3-3
1.00E-06

1.5

BiCG
GMRES

present(GMRES)
commercial software

residual of density

Cp [-]

10.5000.2

0.4

0.60.8

1-0.51.00E-07-1

100x [m]100010000
number of iteration

100000

Fig.10 Pressure coefficient at M=0.1
1.5

Fig.14 Residual history at M=0.01

present(GMRES)Cp [-]

1commercial software

Fig.12 は、高マッハ数ではいずれの反復解法でも同様の収束性であることを示している。Fig.13 より CG 法と BiCG 法はほぼ同様の結果を示しているが、GMRES 法が最も収束が早くかつ収束性が良いことを示している。Fig.14 より GMRES 法の方が BiCG法よりも早く収束し、かつ収束性が良いことが分かる。M=0.01
では CG 法は反復回数の上限として設定した 500 回を超えてしまったので、現実的な計算が不可と判断した。従って、反復解法の収束性はマッハ数に依存することが分かった。高マッハ数では反復解法の違いによる収束性に差異はないが、低マッハ数に近づく
につれて GMRES の収束性の優位性が見られ、逆に CG 法の収束性の悪化が見られた。

0.5000.2

0.4

0.60.8

1-0.5-1x [m]Fig.11 Pressure coefficient at M=0.013.4 各反復解法の収束性
つぎに、M=1.0、0.1、0.01 のときの密度の収束性（ノルム）を調査した結果をそれぞれ Fig.12、Fig.13、Fig.14 に示す。

3.5 各反復解法の計算時間計算時間の比較を行った結果を Fig.15 に示す。
図中 M はマッハ数で 1.0、0.1、0.01 の結果を示した。M=1.0 のときは 0.2[s]までの計算、M=0.1 のときは 2[s]までの計算、そして M=0.01 のときは30[s]までの計算時間である。いずれもほぼ定常状態までの計算時間とみなせる。Fig.15 より、高マッハ数では各反復解法に大きな差は見られない。しかし CG 法は低マッハ数になるにつれて計算時間が急激に増加する。BiCG 法と GMRES 法はどのマッハ数の領域でもほぼ同等の計算時間であることが分かる。

residual of density

1.00E-03
CGBiCGGMRES
1.00E-04

1.00E-05
1000
10000
number of iteration

10000040

Fig.12 Residual history at M=1.0

residual of density

1.00E-04

M=1.0
M=0.1
M=0.0135
calclation time [hr]

100

30
25
20
15
105CG
BiCG
GMRES

0
GMRES

BiCG

CG

Fig.15 Total calculation time of various iteration methods
1.00E-05

つぎに反復解法の一回の反復にかかる平均計算時間を Fig.16 に示す。Fig.16 は、GMRES 法が一回の反復にかかる時間が最も少ないことを示した。ただし Fig.15 から GMRES 法は、全体の計算時間が BiCG 法と同等であることから、一回の時間ステップに必要な反復回数が BiCG 法より多いことが分かる。

1.00E-06
100100010000
number of iteration100000

Fig.13 Residual history at M=0.15Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C3-3５．まとめ
解析手法の妥当性を示し、不完全 LU 分解を前処理として用いて、各種行列計算法を比較した。
その結果以下のことが分かった。

time per iteration [s]

0.25
M=1.0
M=0.1
M=0.01

0.20.15・本解析手法により all-speed 流れ解析が現実的な計算時間で可能となった。
・GMRES 法が今回比較した中で最も収束性のよい反復解法である。
・GMRES 法と BiCG 法は、ほぼ同等な計算効率である。
・CG 法は低速流には適さない。

0.1
0.05
0GMRES

BiCG

CGFig.16 Calculation time per iteration
参考文献(1) 矢田, 市東, 孫 "ハイブリッド格子を用いた気液２相流解析",
第 23 回数値流体力学シンポジウム, A5-3 (2009).
(2) M.S.Liou, “A Sequel to AUSM, Part II: AUSM+-up for AllSpeeds”, J.Comp.Phys., 214(2006), pp.127-170.
(3) J.M.Weiss and W.A.Smith, “Preconditioning Applied to Variable
and Constant Density Flows”, AIAA J., 33(1995), pp.2050-2056.
(4) S.Y.Yoon and T.Yabe, “The Unified Simulation for Incompressibleand Compressible Flow by the Predictor-Corrector Scheme Basedon the CIP Method” Comp. Phys. Comm., 119(1999), pp. 149-158.
(5) 小国, "行列演算ソフトウェア WS、スーパーコン、並列計算機", 丸善, (1991).
(6) Y.Saad and M. H. Schultz., “GMRES: A Generalized Minimal
Residual Algorithm for Solving Nonsymmetric Linear Systems”,
SIAM J. Sci. Stat. Comput., 7 (1986), pp. 856-869.
(7) P.D.Thomas and C.K.Lombard, “Geometric Conservation Law andits Application to Flow Computations on Moving grids”, AIAA J.,
17(1979), pp.1030-1037.
(8) R.Saurel, F.Petitpas and R.A.Berry, “Simple and Efficient
Relaxation Methods for Interfaces Separating Compressible Fluids,Cavitating Flows and Shocks in Multiphase Mixture”,
J.Comp.Phys., 228(2009), pp.1678-1712.
(9) E.F.Toro, “Riemann Solvers and Numerical Methods for Fluids
Dynamics”, Springer-Verlag., Berlin, (1997).
(10) H.Luo, J.D.Baum and R.Loehner, “A Fast, Matrix-free Implicit
Method for Compressible Flows on Unstructured Grids”,
J.Comp.Phys., 146(1998), pp.664-690.

４．考察
r r
Ax = b の解を求める際の行列 A の構造は以下のようになっている。

o o o• • •v • • •0o o o0p • • •o o o• • •• • •0K
A=• • •• • • o o o• • • o o o• • • o o o0Mu：自分自身のブロックデータ（対角要素）：隣接点のブロックデータ（非対角要素）：各方程式（u、v、p）の要素
ここで、黒丸、白丸（黒丸の隣接点のブロックデータ）がデーr
タのある部分であり、それ以外は 0 である。ベクトル x は以下のようになっている。

x = [u1 v1p1 u 2v2p2 L]

T

解析手法では u、v、p の 3 変数について方程式 3 つを解いているので、1 つの検査体積あたり 3×3 行列のブロックデータが存在する。行列 A の行には自分自身の点のブロックデータ（黒丸）が一つと自分自身の点と隣接する点のブロックデータ（白丸）が複数存在する。
このように係数行列 A は対称行列とはなっておらず、効率的な計算を行うには前処理が必要である。現段階では不完全 LU 分解を行っているが、CG 法のような対称行列を得意とする計算手法
では、不完全 LU 分解を前処理として用いても本解析には不向きであるといえる。また GMRES 法は残差を最小にするアルゴリズムであるため、収束性が良い特性があり、今回の結果はそれを裏づけている。より高速化を求めるならば、流体解析用の前処理として LU-SGS を用いるとより効果的であるとの報告もある(10)。

6

Copyright © 2010 by JSFM