第 24 回数値流体力学シンポジウム
C7-3

全速度スキームと前処理付陰解法を用いた空力音解析手法について
CFD method for aero-acoustics using all-speed numerical flux and preconditioned
implicit time integration○嶋英志，JAXA/JEDI，神奈川県相模原市由野台 3-1-1，shima.eiji@jaxa.jp
北村圭一，JAXA/JEDI，神奈川県相模原市由野台 3-1-1，kitamura.keiichi@jaxa.jp
Eiji Shima, JAXA/JEDI, 3-1-1, Yoshinodai, Sagamihara, Kanagawa, Japan
Keiichi Kitamura, JAXA/JEDI, 3-1-1, Yoshinodai, Sagamihara, Kanagawa, Japan
By using all speed numerical flux schemes in MUSCL approach for compressible CFD like SLAU or SD-SLAU, which thepresent authors developed, 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 preconditioned time integration method,which is combined with SLAU, and its accuracy for low Mach number flows and sound propagations. The influence ofparameters 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 and sound propagations is studied.

1. はじめに
圧縮性 CFD スキームを低マッハ数流れ(M<0.1)に適用する場合，過大な数値散逸による解像度の低下および，移流速度と音速の大きな隔たりによるスティフネスが問題となる．
（密度，圧力等のス
カラー量の変動が微小になることによる丸め誤差の問題もあるが，，基準量+変動量（例えば p∞+p’）の形式で変数を記憶し，差分を変動量から計算することで比較的容易に回避可能である．
）
筆者らは AUSM 族数値流束スキームの改良によりパラメータ
調整なしに低マッハ数から高マッハ数まで安定かつ高解像度に計算できる全速度数値流束スキームSLAU[1,2]および超音速での特性を更に改良した SD-SLAU[3,4]を提案した．また，これらの手法は
時間微分前処理（以後，PC と表記）[5-12]との併用により低マッハ数での収束解を短時間で算出できることも分かっており，少なくとも，定常解に関してはスティフネスも解決可能である．
SLAU 等は低マッハ数流れを高解像度で解析できる圧縮性ソルバーなので，低マッハ数流れの空力音響問題について，流れと音
波を同時に解析する直接解法への利用が考えられる．一方，低マ
ッハ数流れに適した非圧縮 CFD では音波を計算できないので，
他の方法で扱う（分離解法）必要がある．分離解法での音波の計算において，ポテンシャル法を用いると，流れ場の影響を正しく考慮できない問題点がある．また，線形オイラー(LEE)を用いれば，流れ場の影響を考慮できるが，LEE の計算の手間は Euler や NSと殆ど変わらず，更に，平均流れ場が定義できないような状況では適用は難しい．したがって，直接解法のメリットは大きい．
陽解法を用いれば，SLAU 等で音波が計算できることは，ほぼ自明であるし，実証もされている．しかし，スティフネスの解消のために必須の陰解法と大クーラン数や，音速を変化させてしまう PC の利用など，低マッハ数流れ解析との両立にはいくつもの疑問と問題点がある．このような，低マッハ数流れにおける音波の取り扱いは，流体機械の空力騒音をはじめ，ロケットエンジン
内で燃焼安定性向上のために用いられる音響減衰器の設計・解析などにも重要である．
PC はそのままでは時間発展を壊してしまうので，工夫が必要である．そのための，DTS（二重時間法）と併用した非定常計算法
複数の市販ソルバーにも導入さには多くの先行研究があり[11,12 等]，れている．しかし， PC 付きの時間積分法の音波との関連は必ずしも明らかではない．また，筆者らが提案した SLAU と MFGS

陰解法の組み合わせは，低マッハ数定常流れに対する有効性は確認されているが，非定常流れや空力音に関しては，十分に調査されていない．
本研究の目的は，MFGS をベースに低マッハ数の空力音響問題に適用可能な新しい時間積分スキーム(TC-PGS)を構築し，その適用性を調査することである．
本稿では，まず，陰解法による流体音直接解法の可能性を簡単な推定で示す．次に，陰解法を低マッハ数非定常流れ解析に適した手法に拡張し，更に計算量を削減した陰解法(TC-PGS)を導く．
この際に，前処理法ではしばしば，保存系とは異なる変数を用いるが，記憶変数と作業変数の概念を導入して統一的で柔軟な手法の記述を導入する．最後に数値例で，低マッハ数定常/非定常流れ
での計算時間短縮での有効性や音の伝播に関する性質などを示す．
2. 陰解法による空力音響問題解析に関する簡易推定
数値スキームの検討のため，一次元線形移流方程式を考える．

qt  aqx  0, a  0(2.1)

どの計算手法でも傾向は同じであるが，簡単な例として，空間を中心差分，時間を一次 Euler 陰解法で差分化すると，フォン・ノ
イマン解析による増幅係数は次のようにあらわされる．

A

1
1  i sin kx

(2.2)

(2.3)
  at / x : クーラン数
ここで，i は虚数単位，Δt，Δx は時間空間の刻み幅である．また
簡単のため k>0 とする．陰解法を用いることでνによらず無条件
安定となる．一方，厳密解の増幅係数は，次式となる．

AExact  exp(ikx)  coskx  i sinkx(2.4)

したがって，式(2.2)が式(2.4)の良い近似であるには，次の二つの条件が必要である．
(2.5)
kx  2x / l  1
(2.6)
kx  akt  2t / T  1
ここで，l は波長，T は周期である．つまり，考慮すべき波の波長よりも十分小さい格子を使い，周期よりも十分小さい時間刻みを使えば，陰解法の波の解析精度はクーラン数νに関係しないこと

Copyright © 2009 by JSFM

第 24 回数値流体力学シンポジウム
C7-3
係数θは，時間刻み変化を伴う 3 点後退差分の場合，下記であらわせる．

を示している．
次に，流体音を考える．ストローハル数が St の非定常流体現象
により発生する音の波数 k は，U,C, L M=U/C を代表的な速度，音速，長さ，マッハ数，として次のように表わされる．
(2.7)
k  2  St  M / L
これを式(2.5)に代入するとΔx に関し次の条件を得る．

x 

 r  2 r 1(1 ,  2 )  , r 1 r 
t n 1  rt n
また，時間一次精度陰解法の場合；

L
2  St  M(2.8)(1 ,  2 )  (1,1)

したがって，マッハ数の低下とともに必要な格子幅の条件は緩やかになる． また式(2.6)に式(2.7)を適用すると次を得る．

akt  U  C 2St  M  t / L  1

   1  1 / M 

(2.6)

(2.8)

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

QP{ 1 i

この場合，式(2.6)は次のように書きなおせる．

L

2  St    1  M 

(3.6)

となる．
S は記憶変数ベクトルであり，丸め誤差の回避や，コーデイン
グ上の都合で任意に選ぶことができるが，ここでは，非定常計算
でも保存則を満足する，S=Q，∂S/∂Q=I を選択する．
また，物理時間 t とは独立に疑似物理時間τは自由に選択でき，τに関して収束すれば，正しい時間発展が達成される．
陰解法に，Gauss-Seidel 反復などの近似因数分解を用いない手法
を用いると，Δτが大なるほど収束が早まるので，Δτ→∞とす
ると次式を得る．
（同様の手法は文献[11]で DTS と LU-SGS を併用する場合にも用いられている．
）

スティフネスの問題を避けるため時間刻みを移流速度と格子サイズだけから決定して，
Δt=αΔx/U(2.7)とすると，クーラン数は下記であらわされ，M の低下とともに増
大する．

x 

(3.5)

(2.9)

k 1

n 1

n

(3.7)

したがって，α，St が一定であれば，低マッハ数での流体現象による音の伝播解析に必要な条件は，ほとんどマッハ数に依存せず，
α<1 であれば StΔx/L<<1 のみであることが分かる．この条件は物体近傍では流体運動解析のための必要性は自明である．
つまり，空力音の伝播解析は，大クーラン数を用いた陰解法でも十分可能であることを示している．

数学的には，PC 行列 P の乗算は無意味であるが，数値解法としては，収束加速に有意である．少し意味合いは異なるが，線型方
程式 Ax=b を反復解法で解く時に前処理行列を作用した PAx=Pb
を代わりに解くことと類似している．また，式(3.7)は，数学的には PC を導入していない元の方程式と同じであるから，物理的な時間発展を追うことができる．

3. 基礎方程式
保存系の圧縮性 NS 方程式は積分形式で次の様に書ける．

4. Gauss-Seidel 反復及び作業変数の導入
式(3.7)の陰解法は左辺に一次風上法を用い線形化を行って，次のように書ける．

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

(3.1)



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

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

(3.2)

k

Q i  Q n 1 i  Q n i

Hi 

Si

 Si

k 1



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

~

 (E

k 1
i, j

(4.1)

k

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

Q i  Q ik
 1S i k 1   2 S i n   1   2 S i n 1
 P{
t



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

(3.3)

ここで添え字i，j はｉ番目のセルのｊ番目の境界面(辺)もしくは，その面で接する隣のセル，n，n+1 は時間ステップ，V，s はセル~ ~，R は非粘性項と粘性項の数値の体積及び境界の面積を示す．E流束である．
低マッハ数の非定常解析のために DTS，PC，3 点後退時間差分，
記憶変数を導入する．k, k+1 は疑似時間反復を示す．
k 1

1Pi 1~   si , j Pi Ai , j  Qi t Vi j1~   si , j Pi A j ,i Q j
Vi j

k 1

 Qi

k

(4.2)

(4.3)
(4.4)

 ~ 



ここで， Pi A i , j 等はｉ番目のセルのｊ番目の境界面(辺)での正(3.4)

方向（i から j への方向）の固有値成分のみを分離した流束の PC
後のジャコビアンである．なお，構造格子の場合は，単に，4 角形あるいは 6 面体の検査体積に適用すれば良いので，非構造格子との区別は必要ない．また，一般には，~ k 1
 R i, j )si, j }  0j2

Copyright © 2010 by JSFM

P A~   P A~ 


i



i, j

i

第 24 回数値流体力学シンポジウム
C7-3でσは PC 後の非粘性流束に関するスペクトル半径（後述）であ
る．

(4.5)

i, j

であり，
式(3.7)ではPC 行列の有無は数学的に無意味であったが，PC 行列の追加で数値解法としては変化していることが分かる．
な
お，収束すれば，所定の時間精度が得られることに違いはない．
風上法を適用することで優対角となるので連立一次方程式(4.1)の解法として Gauss-Seidel 反復を利用できる．保存変数のままでも GS 反復を導入できるが，前処理行列が後述のエントロピー変数の空間で定義されることを考慮すると，エントロピー変数など
を作業変数として導入するのが便利である．
作業変数との変換行列等について， PC 付の 2 次元 Euler 方程式を用いて説明する．なお，下記具体例では簡単のため 2 次元で示すが 3 次元への拡張は容易である．まず，保存形式では次のよ
うに書ける．

Q t  PE x  Fy   0

~i , j   i , j  s Bˆ
i, j

(4.9)


P  M 1PM

(4.10)(4.11)

Bˆ  xn B  ynC(4.12)

Step0：k=0, Qk=Qn
Step1：外部非線形反復：

これらの関係を用いると，式(4.1)はΔW の変換の誤差を除いて，
厳密に次のように変形できる．

 ~ 
1Pi 1  si , j Pi Bi , j  Wi t Vi j~11  si , j Pi B j ,i Μ i M j Wj
Vi j 1 k Pi Μ i Hi







Step1.1：Qk を用いて Pi M i H i k を更新

Step1.2：内部線形反復
Step1.2.1：ΔW=0
Step1.2.2：式(4.19)の 10 回程度の往復反復
End Step1.2
Step1.3：Qk+1=Qk+M-1ΔW, k=k+1
Step1.4：収束すれば終了，でなければ繰り返し
End Step1(4.13)

ここで，下記の近似を導入する．

Mi  M j(4.14)

5. エントロピー変数および前処理行列の導入
Turkel の PC 手法[5]の一種である Weiss&Smith[6,7]の PC 行列は，次のようなエントロピー変数の空間で定義される．

この近似の導入により，式(4.13)は下記のように変形できる．

 ~ 
1Pi 1  si , j Pi Bi , j  Wi t Vi j1~   si , j Pi B j ,i Wj
Vi j 1 k Pi Μ i H i











W  p, u , v, p  c 2 )



T

(5.1)

PC 行列は，この空間では次のような対角行列であらわされる．

(4.15)

P  diag ( ,1,1,1)(5.2)

0   1

セル面に垂直な流束ジャコビアンの固有値λは PC の導入により次のように変化する．

これは，変数が作業変数に変わり，ジャコビアン等がその変数空間のものに変わっただけで,式(4.1)と同じ形になっている．
更に， LU-SGS におけるのと同様の下記の近似を導入する．

~~  Pi Bˆ i , j   i , j
Bi , j 
2(4.19)

右辺のΔW はその時点での最新のものを用いる．また，亜音速流れに適用する場合，1 スィープごとに前進後退する往復反復が有効である．
作業変数を導入することの長所短所は，具体的な選択次第ではあるが，次のようにまとめられる．
×保存変数⇔作業変数の変換の計算時間の増加○流束ジャコビアンがシンプルになり計算時間が減る○PC 行列が対角行列になりスカラーの乗算で済む○ジャコビアンに複雑な熱力学的変数が含まれないので，GS
反復は状態方程式にかかわらず不変
最終的に，この GS 反復を用いた式(3.7)の陰解法の計算手順は次のようにあらわされる．なお，どのような作業変数を用いても外部非線形反復が収束すれば，保存則が満足される．

また，保存形式の流束ジャコビアンとは次の関係がある．

ˆ  Eˆ / Q  M 1Bˆ M
A(4.18)


~
P I
new
[ 1 i   si , j i , j )]Wi 
t Vi j2
Pi Bˆ j ,i  ~i , j
1k[ { si , j
Wj }  Pi M i H i ]Vi j2

このとき，変換行列等は次のように定義される．

M  W / Q

0

これらの関係を用いて，GS 反復を下記のように定義できる．

(4.7)

(4.8)

i, j

j

任意の作業変数ベクトル W を用いると次のように表わされる．

W  MQ ， Q  M 1 W(4.17)

iVi

また，保存系と同様に下記の関係が成立する．

(4.6)


Wt  PBWx  CWy   02(   T ) si , j1  2  Vn  xn u  y n v(4.16)

3 ,  4 

粘性項を含めて，スペクトル半径~ は次のように表される．ここ
3



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



(5.3)

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C7-3
MFGS よりも計算量が少ない．PC を導入すると利点は更に大きくなる．
(3) 式(5.7,8)に示されるように，PC 導入による計算量の増加がほ
とんどない．
(4) エントロピー変数に変換後は，GS 反復が状態方程式によらないので理想気体以外への適用に向く．

したがって，粘性項を含めた，スペクトル半径σは次のように表
される．



1
2
(1   ) | Vn i , j |  (  1) 2Vn i , j  4ci2
22(   T ) si , j

~i , j 




(5.4)

iVi

εは 1 の時，PC 前の固有値を回復し，(マッハ数)2 の時，全固有値が移流速度のオーダーとなりスティフネスが解消される．εは6. MFGS および pMFGS式(4.19)で作業変数=保存変数とすることで，保存変数に対するGS 反復が容易に定義され，更に PC 無しの MFGS(Matrix Free
Gauss-Seidel)陰解法[13]に対するのと同様の近似を導入することにより，pMFGS（preconditioned MFGS）陰解法が次のように定義される．

次のように定義される．

  min(1, max( M 2 , M cutoff 2 ))(5.5)

Mcutoff の選択が PC の安定性と収束性に大きな影響を与えるが，本


P 1
new
[ 1 i   si , j i , j )]Qi 
t Vi j2
ˆ (Q k  Q )
 Q j (6.1)
E1j  si , j i , j[ {Pi  si , j j ,i i
}
Vi
2
2
j
j

研究では 8.3 項での検討を除いて下記のように設定した．

M cutoff  3M 

(5.6)

なお，本研究で右辺の数値流束に用いている SLAU には Mcutoff は含
まれないので，Mcutoff の選択は空間離散化には影響を与えない．
この PC によって，固有値が密集するので，式(4.16)のような近似を用いても誤差が少なくなり，収束性が改善される．
エントロピー変数空間での流束ジャコビアンは疎行列であり，PC 行列は単位行列なので式(4.19)中の行列は簡単なものとなる．
まず，左辺の行列は下記のような対角行列である．

~

 1 Pi Γ i1
t


1
t

ここで，右辺のブロック行列の乗算を流束差分による逆線形化を用いて計算量を減らしている．Ê は次式で定義され，SLAU 等の数値流束ではない．

 i, jI  s i, j
Vi j
2diag  1 1 1 k

 Pi H i ]ˆ  V , uV  px , vV  py , e  p V t
E
n
n
n
n
n
n

(5.7)

 i, j
I
si, j

Vi j
2

また，Weiss-Smith の PC 行列がエントロピー変数空間で定義されることを背景に，保存変数からエントロピー変数への変換行列
を M として，次のように，PC 行列 P に起因するブロック行列の反転を回避できる．

また，右辺のジャコビアンは下記のような疎行列となる．


Pi Bˆ j ,i  ~I Vn  ~ c 2 xn c 2 yn
0 


x
c
 n
0
0 
Vn  ~


 

yn c
~
0
0 
Vn  


 
 0
0
0
Vn  ~ 


(6.2)

1

M [



 1 Pit i, j1
news i, j)]MQ i  R.H .S . (6.3)

Vi j
2

Weiss-Smith の PC 行列を用いた場合，粘性項を含めて，スペクト
ル半径σは TC-PGS と同様に式(5.4)で与えられる．また左辺[ ]内の定義は式(5.7)で与えられる．
ε=1，
P=M=Iの時，PCなしのMFGS
となる．

(5.8)

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

流束はセル境界で定義されるので，式の各要素の定義にどこの物
理量を用いるかは任意性があるが，ここではセルセンター”i”の値を用いる．式(5.8)での法線ベクトル(Vn の決定にも必要)は境界面の面積ベクトルから決定し，方向は，周辺セル”j”から中心セル ”i”の向きとする．
本手法の時間精度は，時間積分スキーム（θの選択）次第でもあり，また外部非線形反復の収束度合いにもよるが，例え収束しない場合でも，（反復を含まない LU-SGS 等と同程度に）ある程
度の精度で物理的時間発展を追うことができる．したがって，精度はともかく，時間無矛盾とは云える．また，エントロピー変数
では，MFGS におけるような流束の差分による逆線形化を用いたMatrix Free 化は不可能なので，本手法を TC-PGS(Time Consistent
Preconditioned Gauss-Seidel)と呼ぶことにする．
既存の手法と比較した本手法の特徴は次のようにまとめられる．
(1) PC の効果を陰解法の中だけに導入することでDTS を用いずに時間無矛盾性と計算効率を両立できる．また陰解法の部分だけのコード改修で導入が可能である．
(2) 保存変数の残差をエントロピー変数に変換した後に GS 反復を行うことで計算量が削減され，後述のように PC なしでも~ ~F  F Q , Q , N



(7.1)

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

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

(7.2)

T

(7.3)

Ν  0, xn , yn 0

T

4

(7.4)
Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C7-3

h  (e  p ) / 

8. 数値例
8.1 最適内部線形反復回数の検討
外部非線形反復（以後，外部反復）の収束は時間精度や保存則のために必要であるが内部線形反復（以後，内部反復）は多くの近似を含んでいるので，完全に収束させる必要はなく，計算時間が短くなるように設定すればよい．ここで，TC-PGS 等の 1 時間
ステップの計算量は次のようにあらわせる．
Total=A+(B+C*Nin) *Nout(8.1)
ここで，Nin,Nout は内部・外部反復回数，C は１内部反復の計算，
B は右辺の流束計算など，A は 1 時間ステップ内のその他の計算である．TC-PGS 等では，Nin が増加するほど Nout は減少する傾向にあるが，Nin と Nout どちらを減らすのが総計算量削減に有効
かは，B と C のバランスにも依存するので一概に言えない．したがって，Nin やNout は各々単独では総計算量の指標にはならない．
そこで，実装やマシンの性能に依存する数値ではあるが，相対比較として，CPU 時間で計算量を示す．CPU 時間はIntel-Corei7(3.2GHz)のパソコン 1 コアでの数値である．
低マッハ数非定常流れとしてカルマン渦列が発生するレイノルズ数=100 の円柱周りの流れを例に最適内部反復回数を調べた．
M=0.001 の場合の，瞬時の圧力分布の例を Fig.1 に示す．

(7.5)
SLAU の質量流束は次で与えられる．

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



2c

(7.6)

p

q  q   q (7.7)

g  g   g   0, 1;   g  min max M ,0 ,1g    max min M  , 0 ,1


(7.8)





  (1  M ) 2

(7.9)

2
2
2
2


1 u   v   u   v  M  min1.0,


2
c



(7.10)CPU time v.s. Inner Iteration # M=0.001 and
M=0.01

上記の関数 g は強い膨張波への対策なので低マッハ数では，簡単に g=0 としてもよい．また，平均圧力はマッハ数を M として下記で与えられる．







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

M 

MFGS M=0.001(7.11)CPU TIME(LOG)p   p     ~p(p  p )22
p  p 1         1
2

100.00pMFGS M=0.001

10.00

TC‐PGS M=0.001
MFGS M=0.01
pMFGS M=0.01

1.00

1

10

20

Nopt

TC‐PGS M=0.01(7.12)
0.10

Fig.2(a) M=0.001 and 0.01

Vnu  xn  v  yn

c
c

5

CPU time v.s. Inner Iteration #M=0.1 and M=0.2(7.13)
1.20

CPU TIME

1.00

MFGS M=0.1

0.80

pMFGS M=0.1

0.60

TC‐PGS M=0.1

0.40

MFGS M=0.2

0.20

pMFGS M=0.2
TC‐PGS M=0.2

0.00

1

5

10

20

nopt

Fig.2(b) M=0.1 and M=0.2
Fig.2 Comparisons of CPU time due to number of inner iterations.
Note that the horizontal axis is not a linear scale and Nopt≦10.
格子形状は O 型，セル数は 100(周方向)x150（半径方向），半径方向の最小格子間隔は 10-2R，計算領域半径は外部境界の影響を避けるため 10000R と大きくとった．
一様流マッハ数は，0.001, 0.01, 0.1,
0.2 とし，
時間刻みは移流速度に基づくクーラン数が約１となるように定めた．その結果，音速ベースのクーラン数は各々のマッハ数で，35,70,700,7000 程度となっている．各時間ステップでは，運動量の残差が初期に比較して 10-4 となるまで外部反復を繰り返し，

Fig.1 Instantaneous pressure distribution around a cylinder in flow
at M=0.001 and Re=100.

5

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C7-3
例して増大する．一方，PC を用いると，マッハ数低下に伴う計算量の増加は一定値にとどまる．
(c) 内部反復回数を適切に選択することで LU-SGS 相当の内部
反復 1 回より 1.4 倍から 4.7 倍程度高速化できる．これは計算量の少ない内部反復を増やして複雑な計算を含む外部反復の収束を早めることが有効なためと考えられる．
(d) MFGS と pMFGS の比較から，PC 有無での違いは M<0.1 で顕著になる．
したがって，マッハ数にかかわらず，適切な内部反復を用いた
TC-PGS が最も効率的である．上記(d)にも述べたように，PC の効果が現れるのは M∞<0.1 以下であり，それ以上での TC-PGS での高速性は内部反復での計算量の削減によるものである．

5 時間ステップでの平均でデータを取った．また，計算開始から
2100 ステップ，
移流速度と直径ベースの無次元時間で 73.5 まで計算を進め準定常状態に達してから，計測している．時間精度は３点後退差分による２次精度，空間精度は（１次元で３次精度再構
築を用いた）２次精度である．
Fig.2 に MFGS, pMFGS, TC-PGS の内部反復回数と CPU 時間の関係を示す．なお，横軸は反復回数に比例しておらず，Nopt≦10（後述）であることに注意．Fig.2(a)は M=0.001 と 0.01 のケース，
Fig.2(b)は M=0.1 と 0.2 のケースであり，Fig.2(a)は縦軸を対数メモ
リとしている．また，Inner Iteration は往復反復の回数なので，GS反復としては，この 2 倍となる．また，Nopt は下記で表わされる
反復回数である．

Nopt  max(1, min(10, N1 / 10 ))(8.2)

Mach # v.s. CPU Time

ただし，N1/10 は内部反復の相対 L2 残差が 1/10 となる内部反復回数であり，したがって，反復回数は 10 回以下となる．なお内部反復の収束判定に，TC-PGS では，エントロピー残差，その他では
密度残差を用いた．この結果から，Nopt は，幅広い条件で，計算時間最小の最適値あるいは，それに近いことが分かる．また最適値近傍では，CPU 時間の変化は内部反復回数に鈍感であることもわかる．これは，内部反復回数増による増加と外部反復回数減による減少が拮抗することによる．
Nopt は低マッハ数非定常流れだけではなく，他の条件でも最適に近い．Fig.3 に NACA0012 翼型，M=0.6，Re=106，Spalart-Allmaras
乱流モデル，局所時間法，MFGS 陰解法（時間一次精度）を用いた定常計算の運動量 L2-Norm の収束履歴を示す．この場合でも，Nopt が計算時間最小となっている．よって，今後の計算では内部
反復回数を Nopt に固定する．

100

pMFGS‐1

Residual of momentum10

CPU(LOG)

MFGS‐Nopt
0.01

0.1

pMFGS‐Nopt1

TC‐PGS‐Nopt

0.1

∝1/M

Fig.4 CPU time (Log scale) v.s. Mach number for one time step ofseveral preconditioned and non-preconditioned schemes.
Mach # v.s. Outer Iteration #10000
MFGS‐11000

Iteration(LOG)

100

TC‐PGS‐1

1

Mach(LOG)

Nopt
LU-SGS
Inner iteration=5
Inner iteration=10
Inner iteration=20

1000

10

0.001

Convergence histry of viscous flow around NACA0012
10000

MFGS‐1

pMFGS‐1
TC‐PGS‐1

100

MFGS‐Nopt10

pMFGS‐Nopt
TC‐PGS‐Nopt

1

1

0.001

0.1
0.01

0.01

Mach(LOG)

0.1

1

∝1/M

Fig.5 Number of outer iteration (Log scale) v.s. Mach number for
one time step of several preconditioned and non-preconditioned
schemes.

0.001
0.0001
1e-05

Mcutoff v.s. Outer Iteration #

1e-06
50

100

150

200

Averaged Outer Iteration

0

CPU Time(sec)Fig.3 Convergence history of steady viscous flow around
NACA0012 at M=0.6
8.2 各手法の計算時間の比較
前項と同じ円柱周りの流れを対象に，MFGS，pMFGS，TC-PGSについて LU-SGS 相当の内部反復が往復 1 回と，前項で定義したNopt 回のケースについて計算時間等の比較を示す．前項と同様に，運動量の相対残差が 10-4 以下となるまでの，１ステップあたりの計算時間及び外部反復の回数を Fig.4,5 に図示する．この結果から次のことがいえる．
(a) マッハ数にかかわらず，最適内部反復回数を用いた TC-PGSの計算時間がもっとも短い．また，M∞=0.2 と M∞=0.001 での計算時間比も 3.6 倍でもっとも小さい．
(b) PC がない場合，CPU 時間も外部反復回数も，ほぼ 1/M に比

250
200

M∞=0.01150

M∞=0.001

100
50
0
1

10

100

Mcutoﬀ/M∞

Fig.6 Effect of cutoff Mach number to number of outer iterations-4
until 10 convergence
8.3 カットオフマッハ数の影響
Fig.6 に M∞=0.01,0.001 のケースについて，Mcutoff/ M∞を 1.0 から 100.0 まで変えた場合の，平均外部反復回数の変化を示す．その他の条件は，前項とおなじである．なお，計算時間もほぼ同様
6

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C7-3
度合いを示す．前項と同じ M=0.01 の円柱周り流れについて，Table.1 に示されるケースについて，揚力変動(Fig.7(a,b))と円柱中
心の上方 100R 位置での圧力変動(Fig8(a,b))の時間履歴を比較する．
計算初期のインパルシブスタートが初期擾乱となるので初期 300
ステップは Case2 のものを共通に用いた．
Fig.7(a)に Case1,2,3,5 の揚力変動を示す．Case5 が少しずれている以外は，ほぼ完全に一致している．Fig.7(b)に示されるように，
Case4,6 の結果は Case1 とは大きく異なっている．Fig.8 の遠方のCp 履歴においても状況は同様である．
もっとも収束判定の厳しい Case1 と大きく結果が変わらない範囲で，計算時間が最も短いのは TC-PGS の反復回数を 2 回とした
Case3 である．収束判定を 10-2 とした Case2 の反復回数は 4 回程度であり，Case3 では残差が 10-1 程度となっている． また，MFGSを用いた収束判定 10-1 である Case5 の精度がそう悪くないことか
ら判断しても，流れ場の時間精度確保には，10-1 程度の収束で許容範囲と考えられる．また，物体から十分離れた 100R 付近の圧力変動を元に音の伝播解析を分離解法で行うならば，この程度の収束で，十分な精度が得られるといえる．

の傾向を示す．前述のように，SLAU 等の全速度スキームを用いる場合には，カットオフマッハ数は，空間離散化への影響はない
ので，収束性のみを考慮して選択できる．図に示されるように，
Mcutoff が 3.0 M∞~10.0 M∞付近で反復回数最小になっている．
Mcutoff= M∞は，最適値からほど遠く，Mcutoff =100 M∞の場合と大差ないのは注目すべき特徴であり，
Mcutoff の選択は過小である
よりは過大である方が良い選択であることを示している．
これは，M∞≧0.1 では Mcutoff=1，即ち PC なし，が最適に近いことも意味しており，8.2 項の検討で，M∞≧0.1 では，PC の効果が見られないことの説明にもなっている．
Table.1 Table of cases for effect of outer iteration convergencecriteria to time accuracy
Case
MethodCriteria
Outer Iteration-4
Case1
TC-PGS10
34.46(Ave.)
Case2
TC-PGS
10-2
4.38(Ave.)
Case3
TC-PGS
2
Case4
TC-PGS
1
Case5
MFGS
10-1
9.86(Ave.)
Case6
MFGS
4

Cp history at (0,100R) M=0.01
0.001

0.0005

Cl history M=0.01
1

Case1
Case2
Case3
Case5

Cl

0
Cp

0.5

Case1
Case2
Case3
Case5

-0.0005-0.0010-0.0015
-0.5

20

30

40

50

60

70

80

90

100

90

100

Non-dimensional Time

Fig.8(a) Case1,2,3,5

-1
20

30

40

50

60

70

80

90

100

Cp history at (0,100R) M=0.01

Non-dimensional Time

0.001

Case1
Case4
Case6

Fig.7(a) Case1,2,3,5
0.0005

Cl history M=0.01
1

Case1
Case4
Case6

Cp

0

Cl

0.5

-0.0005-0.0010-0.0015
-0.5

20

30

40

50

60

70

80

Non-dimensional Time

Fig.8(b) Case1,4,6
Fig.8 Time history of the lift coefficient of cylinder at M=0.01

-1
20

30

40

50

60

70

80

90

100

Non-dimensional Time8.5 陰解法でのクーラン数の音の伝播への影響
ここでは，音の伝播の伝播への陰解法の適用について，2.項の簡易推定の結果を確認する．計算手法は，TC-PGS，3 点後退差分，
外部反復の収束判定は 10-4，カットオフマッハ数は 1.0（PC を用いていないのと同等）であり，その他は前項までと同じである．
Fig.9(a-c)に l/Δx（波長あたりの点数）が 20,40,400 の例を示す．

Fig.7(b) Case1,4,6
Fig.7 Time history of the lift coefficient of cylinder at M=0.01
8.4 外部反復収束の時間精度への影響
ここでは流れ場の時間精度確保のために必要な外部反復の収束
7

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C7-3
8.6 カットオフおよび反復回数の音の伝播への影響
Fig.10 に l/Δx=40，T/Δt=40，クーラン数=1.0，外部反復回数 4回(Fig.10(a))と 20 回(Fig.10(b))の場合について，
カットオフを 0.03，
0.1，1.0 と変化させた場合の音の伝播を示す．PC を用いない
Mcutoff=1.0 では反復 4 回で十分であるが，Mcutoff=0.1 では反復20 回を要している．また Mcutoff=0.03 では反復 20 回でも十分な収束が得られないため，振幅誤差，位相誤差とも大きくなっている．Fig.11 に示す，l/dx=400，，T/Δt=40，クーラン数=10.0，の場
合では，PC を用いない場合には反復４回で十分であるが，Mcutoff≦0.1 では反復 20 回でも十分ではない．
したがって，PC を用いても，外部反復を十分収束させれば音波の計算は可能であるが，クーラン数の増大，Mcutoff の減少とともに収束に必要な反復数は増大し，Mcutoff<0.1 の利用は実用的ではないといえる．8.3 項での検討により， Mcutoff>0.1 で十分な計算効率があるのは M∞>0.01 の場合なので，本手法で流れと音波の同時解析が実用的なのは，M∞>0.01 の場合といえる．

T/Δt>80
音波の一周期を T とすると，
l/Δx やクーラン数によらず，で時間幅の格子収束が得られ，T/Δt>40 程度が減衰や位相誤差の小さい解析に必要であることが分かる．
また，2 次精度では l/Δx =20 は格子の粗いほうの限界であり，より格子の細かい例とは違って，時間幅を小さくとっても厳密解には一致しない．したがって時間空間 2 次精度の陰解法では，T/Δt<40，l/Δx<40 程度が必要なことが分かる．
Sound propagation L/dx=200.0015

Exact Solution
CFL=0.125 T/dt=160
CFL=0.25 T/dt=80
CFL=0.5 T/dt=40
CFL=1.0 T/dt=20

Pressure fluctuation

0.001

0.00050Sound propagation L/dx=40 4 Iteration-0.0005

0.0015-0.001
0

0.2

0.4

0.6

0.8

0.001

1
Pressure fluctuation

X

Fig.9(a) Cases of l/Δx=20
Sound propagation L/dx=400.0015

Exact Solution
CFL=0.25 T/dt=160
CFL=0.5 T/dt=80
CFL=1.0 T/dt=40
CFL=2.0 T/dt=20

Pressure fluctuation

0.001

0.00050-0.0005-0.0010.0005

0

0.5

1

1.5

2

X

Fig.10(a) Computation using 4 outer iterations

0

Sound propagation L/dx=40 20 Iteration-0.0005

0.0015-0.001
0

0.5

1

1.5

Pressure fluctuation

Fig.9(b) Cases of l/Δx=40
Sound propagation L/dx=400
0.0015

Exact Solution
CFL=2.5 T/dt=160
CFL=5.0 T/dt=80
CFL=10.0 T/dt=40
CFL=20.0 T/dt=20

0.001

Exact Solution
Cutoff=1.0
Cutoff=0.1
Cutoff=0.03

0.001

2

X

Pressure fluctuation

Exact Solution
Cutoff=1.0
Cutoff=0.1
Cutoff=0.03

0.00050-0.0005-0.0010.0005

0

0.5

1

1.5

2

X

Fig.10(b) Computation using 20 outer iterations
Fig.10 Effect of cutoff Mach number on accuracy of one
dimensional sound wave propagation of short wave length (l/Δ
x=40 and Courant number=1.0) computed by TC-PGS implicit
scheme for 2 cases of outer iteration numbers.

0-0.0005-0.001
0

5

10

15

20

X

Fig.9(c) Cases of l/Δx=400
Fig.9 Effect of Courant number on accuracy of one dimensionalsound wave propagation computed by TC-PGS implicit scheme
for several mesh densities.
8

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C7-3
なことが分かる．TC-PGS は物理的時間発展を追うことができる
ので，収束解が初期値や時間発展に依存するような複雑な流れ場に対して有効であると考えられる．

Sound propagation L/dx=400 4 Iteration
0.0015

Exact Solution
Cutoff=1.0
Cutoff=0.1
Cutoff=0.03

CONVERGENCE OF FLOW AROUND CYLINDER RE=40
10.00050-0.0005-0.001
0

5

10

15

20

X

0.0001
1e-06
1e-08
1e-10
1e-12

Fig.11(a) Computation using 4 outer iterations

1e-14
0

Sound propagation L/dx=400 20 Iteration
0.0015Pressure fluctuation

200 400 600 800 100012001400160018002000
STEP

Exact Solution
Cutoff=1.0
Cutoff=0.1
Cutoff=0.03

0.001

TC-PGS
TC-off
MFGS

0.01

RESIDUAL OF MOMENTUM(LOG)

Pressure fluctuation

0.001

Fig.12 Convergence of residual of flow around a cylinder at
M=0.01 and Re=40.
CONVERGENCE OF FLOW AROUND CYLINDER RE=40

0.0005
10

TC-PGS
TC-offMFGS

0
8-0.0005CD

6-0.001
0

5

10

15

20

4

X

Fig.11(b) Computation using 20 outer iterations
Fig.11 Effect of cutoff Mach number on accuracy of one
dimensional sound wave propagation of long wave length(l/dx=400 and Courant number=10.0) computed by TC-PGS
implicit scheme for 2 cases of outer iteration numbers.

2

0
0

150

200

Fig.13 Convergence of drag coefficient of flow around a cylinder at
M=0.01 and Re=40.
9. まとめと結論
全速度スキームを用いた低マッハ数の空力音響問題の CFD 解析を目指して，PC(前処理)付きの時間積分スキームの構築と，流れ場，音波の解析に関する基礎的検討を行った．
(1) 時空間の離散化サイズが音波に比べて十分小さい（2 次精度
では 1/40 以下）場合には，音波の解析精度はクーラン数によらず，空力音響問題に関しては大クーラン数を用いた陰解法
でも十分解析可能である．
(2) PC を容易に導入可能であり，また PC を用いても時間無矛盾
で計算量の少ない TC-PGS 陰解法を示し，その有効性を示した．また，反復回数等の調整可能なパラメータに関し，最適値を示した．
(3) M∞≧0.001 の範囲で計算例を示し，流れ場の計算に関しては固有速度の大きな違いによるスティフネスが解消できることを示した．流れの解析だけであれば，外部反復は数回で十分であり，計算量は少ない．
(4) 音の伝搬に関しては，カットオフマッハ数(Mcutoff)が小さい
場合には多数回の外部反復を要し，Mcutoff<0.1 を用いた PC
は音波解析には不適当であることがわかった．

(8.1)

一方，
PC を用いて定常解を求める場合には次のような物理時間発展を無視した方程式が用いられることが多い．

1P~~
Qi   (Ei , j  R i , j ) si , j  0
t
Vi j

100
STEP

8.7 定常問題への適用
TC-PGS は形式的に次の方程式を解いているので，PC がないのと数学的には同じで，物理的時間発展を追うことができる．

1

1
~
~
P  Qi   (Ei , j  R i , j ) si , j   0
Vi j
 t


50

(8.2)

式(8.1,2)の差異は，式(3.4)でΔτ→∞とするかΔt→∞とするかの違いともいえる．定常解の計算例として，8.1,2 項で用いたのと同じ格子を用いた，M=0.01，Re=40 の円柱周りの流れの計算を示す．
ここでは外部反復は 1 回，時間精度は 1 次としている．Fig.12 に残差履歴，Fig.13 に抵抗係数の収束履歴を示す．図中の”TC-off”
は式(8.2)の計算手法を示す．図からわかるように，PC を用いた二つの手法は，収束に至る過程は異なるものの，必要なステップ数
は，ほとんど変わらない．また PC を用いない MFGS より，はるかに収束が早い．したがって，TC-PGS は，PC を用いた非定常計
算のために導入されたものではあるが，定常計算においても有効
9Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C7-3[13] 嶋英志，「構造/非構造 CFD のための簡単な陰解法」，第29 回流体力学講演会講演集 (1997)(5) 本研究の主目的ではないが，TC-PGS は定常計算の場合でも計算効率が高く，収束解が初期値や時間発展に依存するような複雑な流れ場に有効であると考えられる．
この結果から，全速度スキームを用いた空力音響の解析手法に関して次のように結論できる．
(1) M∞≧0.1 では前処理なしの TC-PGS を用いればよい．その程度のマッハ数では PC の効果は小さく，PC なしであれば陰解法でも音の伝播解析に特に支障はないからである．
(2) 0.1＞M∞≧0.01 では Mcutoff=10 M∞とした PC を 20 回程度の十分な外部反復と共に用いるか，音の伝播をクーラン数<1となるような格子で計算することによって対応可能である．
大気中での M=0.01 は約 3.4m/s であることを考慮すると，空力騒音の問題は全速度スキームとTC-PGSの組み合わせで解析可能であるといえる．
(3) M∞<0.01 では，PC を用いてもスティフネスの問題は残されていると言わざるを得ない．音の伝搬を解く場所と流れ場を解く場所で時間積分法を変える，Mcutoff を音の伝搬に適切な 0.1 以上の値に人為的に変えるなどの工夫が必要である．
今後，現実的なケースで有効性を示してゆきたい．
参考文献[1] E.Shima and K.Kitamura, “On New Simple Low-Dissipation
Scheme of AUSM-Family for All Speeds”, AIAA-Paper
2009-136, 2009[2] 嶋英志，「AUSM 族全速度スキームと前処理付き陰解法」，
第 22 回数値流体力学シンポジウム講演論文集(2008)[3] E.Shima and K.Kitamura, ”On AUSM-Family Scheme for All
Speeds with Shock Detection for Carbuncle-Fix”,
AIAA-Paper 2009-3544,2009[4] 嶋英志，北村圭一，「全速度数値流束スキームの改良について」，流体力学講演会講演論文集(2010)[5] E.Turkel, “Preconditioned methods for solving theincompressible and low speed compressible equations”, J.
Comp.Phys., Vol 72, Issue 2, (1987)[6] 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)[7] D. Unrau and D.W. Zingg, “Viscous Airfoil Computations using
Local Preconditioning， AIAA Paper 90-2027 (1997)[8] E. Turkel， Preconditioning Technique in Computational Fluid
Dynamics Annu. Rev. Fluid Mech, Vol 31(1999)[9] E. Turkel， Robust Low Speed Preconditioning for ViscousHigh Lift Flows， AIAA paper 2002-0962 (2002)[10] Y. -H. Choi, C. L. Merkle,” The Application of Preconditioning
in Viscous Flows” J. Comp.Phys., Vol105,(1993)[11] L.D. Dailey and R.H. PLetcher, “Evaluation of Mulgrid
Acceleration for Preconditioned Time-Accurate Navier-Stokes
Algorithms”, Computers&Fluids, Vol. 25, No. 8, (1996)[12] 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)
10

Copyright © 2010 by JSFM