第26 回数値流体力学シンポジウム
講演番号 E05-3

非常に低いマッハ数流れのための全速度スキームの改良について
On the improvement of the all speed flux scheme for very low Mach number flows○嶋英志，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. However, it
is known that, they sometimes produce wiggle at very low Mach number (M<0.01). The new numerical flux scheme, named
WS(Wiggle-Sensing) SLAU is formurated in this paper, aiming to deliminate wiggles with keeping other favorable feature of
SLAU, such as capability to compute sound propagation. The effectiveness of the new scheme is shown by numerical examples,
including 1D sound propagation, very weak shock tube and flow around a airfoil.
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]を提案した．また，これらの手法は時間微分前処理との併用により低 Mach 数での定常解に関してはスティフネスも解決可能である．更に，非定常計
算にも容易に適用できる陰解法 TC-PGS1(Time ConsistentPre-conditioned Gauss-Seidel)および，それを行列前処理に用いた
FGMRES(Frexible Generalized Minimum Residual Method)を提案し，後者は低マッハ数流れと音波を同時に効率良く解ける可能性があること，よりシンプルな前者でも流れ場のみの非定常計算には十分であることを示した．[8,9]一方，未解決の問題点としては，マッハ数の非常に低い（M<0.01）の条件では，解に数値的な Wiggle が発生することがあることが報告されている．[10]例えば，翼周り定常計算に関する例としては
Fig.7 参照．円柱周り O 型格子のような，直交等間隔格子に近い，
きれいな格子で，問題が発生しないケースも考え合わせると，
SLAU 等は低マッハ数では，数値散逸が減少するため，格子不均一によるノイズが消えないという，無散逸の中心差分と同様の現象が生じていると考えられる．2.2 節で，この状況を簡単に示す．
また，Hosangadi[11]らは，等間隔格子を用いた，1 次元衝撃波間問題において，接触不連続の移流速度が M≒10--4 程度となる，小さな圧力差の場合に，様々な全速度流束スキームが Wiggle を発生することを示した．この場合も，計算初期の相対的に大きな擾乱から発生するノイズを消すのに十分な数値散逸が含まれていないことが原因と考えられる．
本論の目的は，低マッハ数の非定常流れと音波を同時に計算可能であり，かつ，上記のノイズを消すことのできる数値流束関数を構築することである．基本的には，数値散逸を増加すればよい
ので，独立な数値粘性項の追加によっても，実現は可能であるが，数値流束関数に含めることで既存コード上での改良が容易になることを狙っている．
以下，2 節では基本的な手法と理論的な分析，更に新手法の定式化を行い，3 節では数値例により有効性を示す．
2. 理論的な検討
2.1 基礎方程式および基本的な解法
保存系の圧縮性 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)

j はｉ番目のセルのｊ番目の境界面(辺)もしくは，
ここで添え字i，
その面で接する隣のセル，n は時間ステップ，V，s はセルの体積~ ~，R は非粘性項と粘及び境界の面積，Δt は時間刻み幅を示す．E
性項の数値流束である．なお，構造格子の場合は，単に，4 角形
あるいは 6 面体の検査体積に適用すれば良いので，非構造格子との区別は必要ない．低マッハ数においても，物理粘性項の取り扱
いに変化は無いので，本論では非粘性数値流束のみに注目する．
AUSM 族スキームによる数値流束は次のように書ける．
（詳細は Appendix 参照）~ m  m  m  m  ~E   pN22(2.4)

Copyright © 2012 by JSFM   u , v, w, h T

第26 回数値流体力学シンポジウム
講演番号 E05-3
考に，次のように表わせる．

(2.5)

   xn , y n , z n ,0(2.6)

h  (e  p ) / 

(2.7)

T

Ei 1/ 2 

ここで次の性質がある．
「行列（D±A）が正定値であれば，上記のスキームは，適切な時間積分との併用でモノトーンになりうる」
例えば，Roe 法も，その一例である．

ただし，(xn,yn,zn)は左から右方向の法線ベクトル，u,v,w は x,y,z 方向の速度，ρ,p, e は密度，圧力，全エネルギーである．また上付
き±は，左(+)右(-)のセルから正負方向に伝わる成分を意味していれる．本論では質量流束の修正に焦点を当てるが，境界圧力 ~pにp  u 

含まれる数値散逸も AUSM 族スキームの安定化に大きく寄与しており，SLAU 等では，これを制御することで低マッハ数での低散逸を実現している．
SLAU 系の質量流束 m は次のように書ける．

p p 
1   (Vn )  (Vn )  Vn       1 g  f p
m

c 
2


Vn 



 Vn   Vn  




(2.9)

g

(2.10)
 max(min( M ,0),1)  min(max( M ,0),1)
ここで，関数 g は高マッハ数の強い膨張でのみ値を持つので，低マッハ数の場合には，単に 0 としてもよい．fp は圧力差分項の効果を制御する関数で，SLAU の場合次のようなり，マッハ数が小さい場合には，fp は，ほぼ１となる．


f pSLAU   ( M )  (1  M ) 2(2.11)

M  min 1 .0, M (2.12)M  (u 2  v 2  v 2 ) / c 2(2.13)

Kp
1 (Mo )max(1M 2 ,0)

Mo  min1, maxM 2 , Mc2 2 u  u  c

(2.19)

S  u(2.20)

 pTCPGS1 | Vn | c /(1   ( M ' ))(2.21)

 uTCPGS 1   sTCPGS 1 | Vn |  (1   ( M ' ))c(2.22)



(2.23)

M<<1 では，(1-χ)≈2M であり，Δu，Δs は移流速度のオーダーとなり，数値散逸が適切なスケールになっていることが分かる．
一方，対応する SLAU の数値散逸は，近似的に次のようにあらわされる．1

 pSLAU | Vn | c

(2.24) uSLAU   sSLAU | Vn |  (1   ( M ))c(2.25)

エントロピー変数での対角形式の数値散逸の範囲で，式(2.18-20)が単調性のための必要十分であることから，SLAU の質
量流束の圧力差分項に関する数値粘性は，単調性を満足するには，不足していることが分かる．このため，低マッハ数では中心差分に近づくため，ノイズが消えない，と考えられる．
したがって，SLAU の数値散逸を式式(2.18-20)の形式に近づけることで，改良をはかる．[10]圧力差分項に関する数値散逸の増強には，次のような問題点が挙げられる．
(1) ゼロ割を防ぐための TC-PGS1 におけるのと同様に，カット
オフマッハ数の調整が必要である．
(2) 低マッハ数では圧力に関する数値散逸が大きなものとなり，圧力伝搬(=音波)の減衰が極めて早い．
（3.1 節）
このうち，前者は，一様流が存在する問題では，一様流マッハ数に固定することで，大きな問題は生じない．
（2）の問題点を回避し，ノイズを除去する質量流束を次節で提案する．

(2.14)(2.15)

ここで，Kp,σ は O(1)の定数，カットオフマッハ数 Mc は一様流
マッハ数程度である．M<<1 では，(1-χ)≈2M であり，低マッハ数
で，圧力差分項の寄与が極めて大きくなるのが SLAU との違いである．次節では，この項の分析を示し，新スキームを導く．
2.2 適切な数値散逸の検討
数値散逸項の適切なスケールの考察のために，ここでは，一次元非粘性流を考える．様々な変数 U で，方程式は次のように書ける．

U t  AU x  0(2.18)

M '  min 1.0, M 2  Mc 2

一方，Liou によって提案された全速度数値流束である
AUSM+-up[12]では，完全に対応するものではないが，fp は，近似的に下記のように表わされる．
（詳細は Appendix 参照）

f pAUSMUP c

ここで，α は任意の正数である．即ち，数値散逸の選択には α の分の自由度があるが，圧力と速度に対するものは独立では無い．
全速度スキームに対する陰的時間積分法 TC-PGS1[8,9]は，この変数変換と対角形式の数値散逸を応用したものである．TC-PGS1 では，低マッハ数での速度場の低散逸と，淀み点での無散逸回避を狙って，次のような修正数値散逸を用いている．

(2.8)


T

D を対角行列 D=diag(Δp,Δu,ΔS)であると仮定すると（半）正定値性の必要十分条件は次のように表わせる．

あり，これらの違いにより様々な AUSM 族数値流束関数が定義さ

次に，エントロピー変数，U  p , u , p  c 2  ) をとり，る．また，m ，~p はセル境界で定義された質量流束および圧力で1
A i1/ 2 U i  U i1   Di1/ 2 U i1  U i  (2.17)
2(2.16)

1

FVM 形式で，数値散逸 D を加えた数値流束は Roe 法の表記を参2

定義式からの類推であり，正確な固有値解析の結果では無い．
ただし，TC-PGS1 との組み合わせで安定に計算できることから，
この値を超えることは無いことは数値的には云える．
Copyright © 2012 by JSFM

第26 回数値流体力学シンポジウム
講演番号 E05-3
うな，各座標方向の最高波数一次元的 Wiggle では，2.3 新しい数値散逸項の提案
モノトーンスキームとするため，AUSM+-up と同様に，fp を次のように設定することを考える．

f pMT 1
 f pSLAU
1   (M ' )(  p i / p i )  1
を満足し，目的を達成することが分かる．

(2.26)

2.4 運動量項の修正
前述のように，単調性の保証のためには，TC-PGS1 陰解法の数
値散逸相当の式(2.21-23)を用いればよい．式(2.22)と式(2.25)の比較
より，速度に関する数値散逸も下記のような修正が必要なことが分かる．

しかし，これは音波を強く減衰してしまう．
（3.1 節参照）そこ
で，平滑な部分では数値散逸を抑制し，Wiggle の発生する部分だけで増強するために，次のような形式を用いる．

 (M ' )   f pSLAU
f pNEW  1  f ws
1   (M ' )   1 at wiggle
f ws  
 0 otherwise

 uSLAU   sSLAU

そこで，SLAU の運動量項中での平均圧力を下記のように変更することでこの効果を実現する．

max  p L ,  p R 
f ws  min 1, C ws


max

p

,

p

,



L
R



   
p p
~
 
p
( p  p )22  2
 1   ( M ' )       1
c
2






(2.36)SLAU では不要であった，カットオフマッハ数の定義が必要になることが新手法の難点であるが，質量流束で必要であるものと同じであり，更なる新パラメータの導入には当たらない．前節で述
べた，質量流束の利用と合わせて，新スキームを WS(Wiggle
Sensing)-SLAU と呼ぶことにする．なお，修正後も，LHS の数値
散逸が RHS のそれより大きいから，TC-PGS1 陰解法は，何ら変更なく利用できる．

(2.28)

ここで，ε は，ゼロ割と丸め誤差のノイズを防ぐための小さな正数定数，Cws は O(1)の正数定数，下付き L,R はセル境界の左右のセル中心を示す．境界値は， pi ,  pi 各々を０次外挿とした．
定数は，以下を標準とするが，比較のために，Cws について Table.1のケースを検討する．

C ws  1.0,   0.1* q ref(2.35)


| Vn |  (1   ( M ))c | Vn | (1   ( M ' )) c

(2.27)

高マッハ数では，χ→0 なので高速での性質は変わらない．
Wiggle は，隣り合う格子点が互い違いに変動するような，高波数の波である．したがって，振幅に関わらず，高波数の波だけを捉えるような関数(Wiggle Sensor)を設計する必要がある．ここでは，圧力を用いた次のような関数を利用する．



(2.34)

Table 1 Cws setting in each case
Note
Cws
Case0
0.0
Original SLAU
Case1
1.0
Standard fws
fws =1, i.e. monotone
CaseInf
Very Large
scheme of Eq.(2.27)(2.29)

ここで，qref は参照動圧である．また，  は，Laplacian を模擬した，セル中心で評価される，次のような差分オペレータである．

pi 

Nd
 pi  p 
2

1D Sound propagation Cutoff=0.01(2.30)

0.006
0.005

Wiggle の多次元格子内でのスケール調整のために導入されてい0.004
Velocity fluctuation

ここで， p は，周囲のセルでの平均，Nd は空間次元で，一次元的る．   は同じオペレータを 2 回を適用することで計算される．
次に，fws の性質を格子間隔 Δx の一次元の例で説明する．Taylor
展開により，

Cws=0
Cws=1
Cws=inf

0.003
0.002
0.001

x 2
x 4
p i  
p xx , pi 
p xxxx4
16(2.31)

0-0.001

したがって，0 pi / pi   ( x / 4) p xxxx / p xx
2

2

0.4

0.6

0.8

1

1.2

1.4

X

(2.32)

Fig.1 Velocity distribution of 1D sound propagation. Each results
is shifted.

である．p にフーリエ成分を代入すると，波数を k として， pi / pi  (1 / 4)kx 

0.2(2.33)となり， (  pi / pi ) は波数センサーとして働くことが分かる．
離散式で解析できる波では，|kΔx|<<1 であるから，Wiggle のない
所では，十分小さい値である．一方， 0,1 が互い違いに表れるよ3

3. 数値例
本節では，3種類の計算で提案手法WS-SLAUの有効性を示す．
定数は，式(2.29)を用いるが，比較のために，極端に小さい場合と大きい場合も検討する．Table.1 で Case0 は fws =0.0 であり，オリジナルの SLAU と同等である．また，C aseInf は至る所で fws =1.0
Copyright © 2012 by JSFM

第26 回数値流体力学シンポジウム
講演番号 E05-3
であり，式(2.26)相当のモノトーンスキームとなる．

Verocity distribution of weak shock tube

3.1 一次元音波の伝搬
x=0 で正弦波の条件を与え，正方向に伝播する，音波の解析を行った．計算条件は下記である．
・空間離散化：3 次精度 MUSCL
・制限関数：無し
・時間積分法：2 次精度 2 段階陽解法
・CFL 数=0.01・比熱比γ=1.4・一波長あたり格子点：40(ppw)
・ステップ数：25000
・速度振幅：1x10-3
・カットオフマッハ数：0.01
・qref=5x10-7
なお，40ppw は，減衰なしに音波を捉える，ほぼ，最低限の格子密度である．[6,7]Fig.1 に各ケースでの速度分布を示す．Case0，Case1 では，音波は減衰せずに伝搬しているのに対し，CaseInf では減衰が大き
い．Wiggle 検知が想定通りに働き，滑らかな部分での数値散逸増加を防いでいることが分かる．

Cws=0
Cws=1
Cws=inf

0.000140.00012

Verocity

0.0001
8e-05
6e-05
4e-05
2e-05
0
0

20 40 60 80 100 120 140 160 180 200
Grid number

Fig.2 Velocity distribution in weak shock at 5000 step.
Density distribution of weak shock tube
1.0004

Cws=0
Cws=1
Cws=inf

1.00035
1.0003

3.2 弱い一次元衝撃波
理想気体の一次元衝撃波管で圧力差が小さく，発生する接触不連続のマッハ数が非常に低い場合を考える．界面左右の初期条件は次で与えられる．

  L  1.0002804    R  1.0    uL   c  
 L 

01.0    ,  u R    0   c  1.0   R  

Density

1.00025
1.00021.00015
1.0001(3.1)

1.00005
1

接触不連続のマッハ数は約 1.0x10-4 であり，qref=0.5x10-8 とした．
その他計算条件は，下記とした．
・空間離散化：2 次精度 MUSCL
・制限関数：van Leer・時間積分法：2 次精度 2 段階陽解法
・CFL 数=0.01・比熱比γ=1.4
Fig.2,3 に 5000step 後の速度，密度の分布を示す．Case0 では接触不連続を挟んで速度場に，Wiggle が発生しているが，Case1,
CaseInf では消去されている．密度分布には顕著な Wiggle は，どのケースでも発生していないが，Case0 ではショック・膨張波でオーバーシュートが見られる．省略するが，圧力も密度と同様である．Fig.4 に示されるように，5000step 目では，ショックおよび膨張波の伝播した部分では fws は 1 に近い大きさになっている．一方，Fig.5 に示される 50000step 目では，伝播の先端でやや大きな値を取る他は，全体に小さくなっており，初期の擾乱を減衰させた後は，数値散逸の増強なしでも Wiggle の無い解が得られることを示している．つまり，大部分では数値散逸は十分小さいので，音波の減衰も小さいことが期待される．
このように，このケースでは，WS-SLAU は必要な場所で，数値散逸を強化し，Wiggle 抑制に有効に作用していることが分かる．

020 40 60 80 100 120 140 160 180 200
Grid number

Fig.3 Density distribution in weak shock at 5000 step.
fws distribution of weak shock tube at 5000 step
Cws=1
1

fws

0.8

0.6

0.4

0.2

0
0

20

40

60

80

100 120 140 160 180 200

Grid number

Fig.4 Distribution of fws(Wiggle sensor) of Case1 in weak shock at
5000 step.

4

Copyright © 2012 by JSFM

第26 回数値流体力学シンポジウム
講演番号 E05-3
fws distribution of weak shock tube at 50000 step
0.06

Cws=1

0.05

fws

0.04
0.03
0.02
0.01
0
0

200 400 600 800 1000 1200 1400 1600 1800 2000
Grid number

Fig.5 Distribution of fws(Wiggle sensor) of Case1 in weak shock at
50000 step.
3.3 2 次元翼型の定常計算
NACA0012，M=0.001，迎角=0 度の非粘性流れの例を示す．勾配制限関数無し，3 次精度 MUSCL，TC-PGS1 陰解法による定常解の比較である．後縁付近±0.1C 範囲の格子図(Fig.6)に示される
ように，後縁周りの O 型格子のため，格子の歪が生じている．これに起因する Wiggle を，圧力係数のカラーマップで Fig.7-9 に示
す．
Fig.7 に示される，Case0(Cws =0)の結果では，後流部分に激しい
Wiggle が発生している．なお，絶対圧力での変動は僅かなので計算の安定性には影響しない．本論の新手法の結果を，Fig.8(Case0)
および Fig.9(CaseInf)に示す．Fig.5 のケースでは，至る所で数値散逸が増強され，式(2.27)に相当する．Fig.8,9 では Wiggle は，ほぼ抑制されている．
このように，このケースでは，WS-SLAU は必要な場所で，数値散逸を強化し，Wiggle 抑制に有効に作用していることが分かる．
また，流れ場のみの解析で，音波の解像が必要なければ，式(2.26)の利用で十分であり，あえて fws を導入する必要は無いが，WS-SLAU の利用も悪影響は無いともいえる．

Fig.7 Color map of the pressure coefficient(-0.1<Cp<0.1).
Case0(Cws=0.0), which corresponds to the original SLAU.

Fig.8 Color map of the pressure coefficient(-0.1<Cp<0.1).
Case1(Cws=1.0)

Fig.6 Enlarged view new trailing edge of the mesh around
NACA0012 airfoil. x 1.0C 0.1C, y
0.1C
5

Copyright © 2012 by JSFM

第26 回数値流体力学シンポジウム
講演番号 E05-3[5] E.Shima and K.Kitamura, ”On AUSM-Family Scheme for All
Speeds with Shock Detection for Carbuncle-Fix”,
AIAA-Paper 2009-3544, 2009[6] 嶋英志，北村圭一，「全速度スキームと前処理付陰解法を用いた空力音解析手法について」，第 24 回数値流体力学
シンポジウム講演論文集, 2010[7] E.Shima and K.Kitamura, ”CFD Method for Aero-acoustics
using All-speed Numerical Flux and Preconditioned Implicit
Time Integration”, AIAA-2011-3045, 2011[8] 嶋英志，「圧縮性 CFD による低マッハ数流れ計算のための新しい陰的時間積分法」，第 24 回数値流体力学シンポ
ジウム講演論文集, 2011[9] E.Shima and K.Kitamura, “New approaches of efficient andtime accurate compressible CFD for unsteady low Mach number
flows”, AIAA-P2012-0571, 2012[10] A. Hosangadi, Personal Communication, 2011[11] A. Hosangadi, J. Sachdev and V. Sankaran, “Improved Flux
Formulations for Unsteady Low Mach Number Flows”,
ICCFD7-2202, Proceedings of 7th International Conference on
Computational Fluid Dynamics , 2012

Fig.9 Color map of the pressure coefficient (-0.1<Cp<0.1).
20
CaseInf(Cws=10 ), which corresponds to the monotone scheme.
4. 結論
全速度数値流束関数 SLAU に対して，非常に低いマッハ数(M<0.01)で発生することのある，Wiggle を除去し，解を改良するために，改良を行った．
・固有値分析から導出され，既に AUSM+-up 等で使われているのと類似の，質量流束項での圧力差分起因の数値散逸項強化は，特に定常解では Wiggle 除去に有効である．
・流れ場のみであれば問題はないが，この数値散逸の強化は，音波への減衰効果が大きく，低マッハ数流れと音波の同時解析には不適当である．
・音波の計算との両立するために，Wiggle 発生箇所を検知する関数を用いた，Wiggle でのみ数値散逸を増強する，新手法（WS-SLAU）の提案を行い，一次元音波，弱い衝撃波，翼型周りの流れ場で有効性を確認した．
・一次元音波では，オリジナル SLAU と同様に，減衰のない計算が実施できるとともに，衝撃波，翼型周りの計算で発生していた
Wiggle を除去することに成功した．

[12] Liou, M., “A Sequel to AUSM, Part II, AUSM+up for All
Speeds,” Journal of Computational Physics, Vol. 214, No. 1, pp.
127-170, 2006.

謝辞
CRAFTTech 社の Ashvin Hosangadi 博士，名古屋大学の北村圭一
博士との貴重な議論に感謝します．
参考文献[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] 嶋英志，北村圭一，「全速度数値流束スキームの改良について」，流体力学講演会講演論文集(2010)
6

Copyright © 2012 by JSFM

第26 回数値流体力学シンポジウム
講演番号 E05-3
Appendix A. WS-SLAU 定式化まとめ
数値流束関数 WS-SLAU は AUSM 族の共通の形式で次のように書ける．

~ m  m  m  m  ~FΦ Φ  pN22(A.1)

Φ  1, u,v , w, h 

(A.2)

T

Ν  0, xn , y n , z n ,0(A.3)

h  (e  p ) / 

(A.4)

T

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



Vn 

 Vn   Vn  




(A.6)


f pSLAU    (1  M )2

(A.7)


M  min 1 .0, M 

(A.9)



(A.10)(A.8)

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

f pSD  SLAU  max  p L ,  p R 
f ws  min 1, C ws


max

p

,

p

,



L
R



(A.12)

C ws  1.0,   0.1* q ref(A.13)

(B.3)

2

(B.4)

(B.5)

CSD1  0.1
一方，AUSM+-up[12]の質量流束は下記で与えられる．

1  cM1/ 2        c M1/ 2      m
2
M1/ 2  M(4)  M(4) 

は参照動圧である．また，  は，左右のセル中心で評価される，

Kpfamax(1M 2 ,0)p   p
c 2

(B.6)

(B.7)

fa  Mo(2  Mo) [0,1]次のような差分オペレータである．

Nd pi  p 2

M 1  M 1  2 M

  C p / p  C  2 
SD1    min1,  SD 2  p max / p  CSD1  
CSD 2  10

ここで M はマッハ数，Vn は垂直外（右）向き速度，c は音速，qref

pi 

(B.2)

M  Vn / c(A.11)

(B.1)

SD-SLAU' 
f pWS  SLAU  1  f ws
 '1  ' '  (1  M ' ) 2

SLAU 系各スキームの質量流束の差は fp の違いに帰着する．
SLAU：

 max(min( M  ,0),1)  min(max( M  ,0),1)


M '  min 1.0, M 2  Mc 2

(A.17)1p 
  (Vn )  (Vn )  Vn  1 g  f p 
2
c



g

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


m

1
  (Vn )  (Vn )  Vn  1 g  f pWSSLAU 2c(A.16)

Appendix B. 質量流束項の比較
SLAU，SD-SLAU の質量流束は SLAU 系共通の形式で，与えられる．

p  (A.5)







M 

面の外向き法線ベクトル，e,p, ρ,h は密度，体積当たりの全エネル
ギー，圧力，密度，エンタルピー，である．また上付き±は，左(+)右(-)のセルから正負方向に伝わる成分を意味している．
WS-SLAU の質量流束は次で与えられる．





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

ここで，  u,v, w  はデカルト座標の速度， x n , y n , z n  はセル表

m



(A.15)(B.8)

Mo2  min(1, max(M 2 , Mc2 )(A.14)



M 2  u2  u2 / 2c 2
ここで，Kp,σ は O(1)の定数 M(4) ,M(4) は左(+)右(-)のマッハ数を  は同じオペレータを 2 回を適用することで計算される．
境界圧力 ~
p は次で定義される．

用いた四次関数である．近似的に SLAU 系と同等の形式で fp は，
次のように書ける．

f pAUSMUP 
7

Kp
1 (Mo )max(1M 2 ,0)(B.9)

Copyright © 2012 by JSFM