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

AUSM 族スキームの新しい圧力流束：
極超音速空力加熱計算のための衝撃波安定，全速度数値流束の開発に向けて
A New Pressure Flux for AUSM-Family Schemes:
Toward Shock-Stable, All-Speed Flux Functions for Hypersonic Heating○ 北村圭一，JAXA/JEDI，〒252-5210 相模原市中央区由野台 3-1-1, kitamura.keiichi@jaxa.jp
嶋英志，JAXA/JEDI，〒252-5210 相模原市中央区由野台 3-1-1, shima.eiji@jaxa.jp
Keiichi Kitamura, JAXA/JEDI, 3-1-1 Yoshinodai, Chuuou, Sagamihara, Kanagawa 252-5210, Japan
Eiji Shima, JAXA/JEDI, 3-1-1 Yoshinodai, Chuuou, Sagamihara, Kanagawa 252-5210, Japan
In hypersonic flow computations, it is still challenging to design a numerical flux function which satisfies both A)
robustness against shock-induced anomalies, and B) accurate prediction of aerodynamic heating. In this paper, wedevelop such schemes by combining a) a new pressure flux appeared first in SLAU2 scheme having excellent
performance in A), and b) a mass flux of another all-speed AUSM-family scheme (AUSM+-up or LDFSS2001). The
new schemes, named “AUSM+-up2” and “LDFSS2001-2,” are tested through numerical experiments and their desired
performances are demonstrated for a wide spectrum of Mach numbers and aerodynamic problems.
１．はじめに
極超音速流れの数値計算では A)衝撃波異常への堅牢性，およびB)空力加熱の正確な予測が課題となる．A)については“カーバンクル現象(1)-(3)”に代表される衝撃波を数値的に捉える際に生じる異常解（図 1）として，その原因や対策が 20 年以上に亘り研究されているが，いまだに決定的な答えは見つかっていない．こうした衝撃波異常が発生すると，衝撃波の形状そのものが変化したり，数値的に振動したりするため，その背後の流れを大きく乱してしまい，重大な問題となる．複雑形状周りの流体計算では，「そもそも異常解が発生しているかどうか」を検知困難であり，生じる誤差量はもはや予測不能である．
B)については，A)に関連が深く，衝撃波異常が起こればその下流の空力加熱（壁面熱流束）は大きく影響を受けるため，A)と同時に議論される事が多かった(4)-(6)．その中で著者らは以前，空力
加熱を正確に予測するための三つの性質を挙げ，種々の流束関数についてこれらの性質を調べた(4)．その三つとは，
A) 衝撃波異常に対し堅牢B-1) 全エンタルピ保存
B-2) 境界層解像であった．これら全てを満足する手法は今のところ存在せず，特にA)についてはVan Leer のFVS (Flux Vector Splitting)(7)ですら完全には満足しない．しかし，少なくとも限られた条件下では A)は
これらはB-1),
AUSM+(8)などのAUSM族スキームにより満足され，
B-2)も満たすため，条件によっては良い空力加熱分布が得られた．
その後，
著者らは A)を Van Leer の FVS と同等に満足し，
かつ B-1),B-2)の性質も保有する AUSM 族スキームとして SLAU2 を発表した(9)．しかしながら SLAU2 では，後述するように A), B-1), B-2)が満たされた場合でも B)の空力加熱分布に空間的振動が生じる
場合がある事が分かった．よってよどみ点のみにおける空力加熱に注目するのであれば A), B-1), B-2)で十分であるが，それ以外の位置での壁面上の分布も正確に求めるにはもう一つ，
B-3)滑らかな熱流束分布の再現という条件を加える必要がある．
高速流れには近年，A)の性質の強い AUSM 族スキームが用いられる事が多い．AUSM 族スキームはその定式化上，a)圧力流束と b)質量流束に分割できる．それぞれの寄与は基本的に，文字通
り質量および圧力のセル界面における差に起因して計算される．

衝撃波の捕捉については，その前後で質量は保存されるが，圧力差は不連続であるため，a)圧力流束の影響が大きいと考えられる．
一方で，壁面熱流束に直接関わる境界層の解像については，一般に境界層垂直方向には圧力差が存在しない事から，b)質量流束による寄与が支配的である．
a)b)Fig. 1 Carbuncle (top) and “correct” (bottom) solutions.
そこで本稿では，
A)に優れた SLAU2 の a)圧力流束と，
B-1), B-2)を満たし，もし衝撃波異常が発生していなければ（つまり，A が満たされた場合には）B-3)滑らかな熱流束分布を出力する別のAUSM 族スキーム AUSM+-up(10)および LDFSS2001(11)の b)質量流束を組み合わせる事を考え，A), B)双方に優れたスキームを構築する．
なお，SLAU(12), SLAU2, AUSM+-up および LDFSS2001 はいずれも全速度スキーム(All-Speed Schemes)と呼ばれ，低速流れにおけるC) 低速流れにおける低散逸という性質を有する．これらの組み合わせである新しいスキーム
も同様の性質を持つ事を，数値計算例において示す．
1

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
B12-3
響問題等に利用されている AUSM 族スキームである SLAU(12)について示す．セル境界における数値流束を F1/2 と表記すると，

２．支配方程式
支配方程式は圧縮性 Euler もしくは N-S 方程式であり，前処理
行列Γ(13)を含む形で下記のように書ける（3 次元では下添字 k, l, m,
n に 1, 2, 3 が代入される）．なお前処理を行わない場合は，Γ=Ι
と読み替えれば良い．

∂Q ∂Fk ∂Fv kΓ+=∂t ∂xk
∂x k
ρu k
⎡
⎤ ⎫
⎪
F k = ⎢⎢ ρ u l u k + p δ lk ⎥⎥ , ⎪ρu k H
⎣⎢
⎦⎥ ⎪
⎪
⎪
⎡
⎤
⎬
⎢
⎥
0
⎪
⎢
⎥
⎪
Fv k = ⎢
τ lk
⎥
⎪
∂T ⎥
⎢u τ
⎪
κ
+
⎢ m mk∂ x k ⎥⎦⎪⎭⎣

F1/ 2 =(3a)

Ψ = (1, u , v, w, H ) , N = (0, n x , n y , n z , 0 )
T

(1a)

T

(3b)

と書け，a) 圧力流束 ~
pは

⎡ ρ ⎤
Q = ⎢⎢ ρ u l ⎥⎥ ,⎣⎢ ρ E ⎦⎥⎛ ∂u∂u ⎞ 2 ∂uτ lk = μ ⎜⎜ l + k ⎟⎟ − μ n δ lkx∂
⎝ k ∂xl ⎠ 3 ∂x n

m& + m& + m& − m& − ~Ψ +Ψ + pN22

p + pR β + α = 0 − β − α = 0~( pL − pR )
p= L
+
2
2
p + pR+ (1 − χ ) β + α = 0 + β − α = 0 − 1 L
2

(

(1b)

)

)

(1c)

ここで，ρ は密度，ui は速度成分, E は全エネルギ，p は圧力，H は全エンタルピ(H = E + (p/ρ)), そして T は温度である．
，プラントル数は Pr=0.72
気体は完全気体の空気（比熱比γ =1.4）
である．分子粘性μ は一定とし，熱伝導係数κ とはκ=cpμ/Pr の関．更に，乱流粘性μt と乱流プラントル係がある（cp は定圧比熱）数(Prt =0.89)からκt=cpμt/Prt の関係がある．
これはデルタ形式で次のようになる．

(3c)χ = (1 − M )2

(3d)

2
2
2
2
2
2
⎛
)
1 u L + v L + wL + u R + v R + wR ⎞⎟M = min⎜1.0,
⎜
⎟
c
2
⎝
⎠

(3e)

⎧1
if M ≥ 1⎪ (1 ± sign(M )),β± α =0 = ⎨ 2
1⎪ (M ± 1)2 (2 m M ), otherwise
⎩4

(3f)

V n un x + vn y + wn z
=
c
c

(3g)

cL + cR
2

(3h)

M =
c=

&は
であり，b) 質量流束 m

Vi
ΔQ i + Γ i−1 ∑ (Fi , j − Fv i , j ) S i , j = 0
Δt i
j

(2)

m& =

ここで，Vi はセル i の体積，Δti は局所時間刻み，ΔQi は保存量の時間変化，Fi,j および Fvi,j はセル i とその隣接セル j との界面
Si,j を通るそれぞれ非粘性（数値）流束および粘性流束である（図2 を参照）
．

) (

(

)

+
−
1⎧
χ ⎫
⎨ ρ L VnL + Vn + ρ R VnR − Vn − Δp ⎬2⎩c⎭| Vn | + = (1 − g ) | Vn | + g | VnL |,(3j)| Vn | − = (1 − g ) | Vn | + g | VnR |
Vn =

(3i)

ρ L VnL + ρ R VnR
ρL + ρR(3k)

g = − max [min (M L , 0 ),− 1](3l)⋅ min [max (M R ,0 ),1] ∈ [0, 1]

である．
SLAU の発表後，著者らはまずセル境界音速 c の与え方に任意性が有る事に注目し，この衝撃波異常への影響について調べた(9)．
そして SLAU の a)圧力流束の散逸量を制御する式(3c)の最終項に注目し，式(3e)では流速がどんなに大きくても M>1 では散逸量が一定であった事から，散逸量をマッハ数に比例する形へと変更を行った．更に理想気体以外への拡張性を考慮し，式(3c)をFig. 2 Schematic of cell geometric properties.
３．数値流束ここで，本稿で扱う数値流束について説明する．本稿では，Liouらによって発表された AUSM 族スキーム(9)-(12)を扱う．AUSM 族
スキームは Van Leer の FVS に端を発し，
これまでに数多くのバージョンが発表されており，主要なものを図 3 にまとめた．

p + pR β + α = 0 − β − α = 0~( pL − pR )
p= L
+
2
2
2

2

2

2

2

u L + vL + wL + u R + vR + wR2
⋅ β + α = 0 + β − α = 0 − 1 ρ c1 2+

３－１．SLAU および SLAU2
まず，Shima らによって発表され，以降，様々な空力問題，音

(

2

2

(4a)

)

Copyright © 2010 by JSFM

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

ρL + ρR

ρ=

(4b)

2

+

(

c L R = c *2 max c * , Vn

c

±

)

(4c)

2(γ − 1)H=
(γ + 1)fとした．この新しい圧力流束を用いた SLAU を SLAU2 と呼ぶ．
SLAU2 は衝撃波異常に対し非常に堅牢であり，対応する数値実験の結果は後に示す．なお，これらの手法では一様流マッハ数 M∞，
参照速度 Vref などの参照パラメタは不要である．

m& + m& + m& − m& − ~Ψ +Ψ + pN
2
2

ρ=

T

T

(5b)

と表現するとき，a) 圧力流束 ~pは~
p = β+

αpL + β −

α

(5l)

−

p R + pu(5c)⎧1
⎪ (1 ± sign(M )),β± α = ⎨221⎪ (M ± 1)2 (2 m M ) ± αM M 2 − 1 ,
⎩4

(

)

⎧ 1
(M ± M ),if M ≥ 1
⎪
=⎨ 2
2
1
1
⎪± (M ± 1)2 ± M 2 − 1 , otherwise
8
⎩ 4

(

Kp
fa

(

max 1 − σM 2 , 0)) p ρc− p ,
R

L

2
12

(5o)

ρL + ρR
2~
F = m& + Ψ + + m& − Ψ − + ~
pN

(5d)

(6a)

Ψ = (1, u , v, w, H ) , N = (0, n x , n y , n z , 0 )T(

pu = − K u β + β − (ρ L + ρ R )( f a c1 2 ) Vn − Vn
−

+

)

(
(

R

c *2 =α=)np + + p − β + α =0 − β − α =0 +~p=+
p − p−
22+ ρ V ref2 ,1 2 β + α = 0 + β − α = 0 − 1

(

(

(5f)

)

2(γ − 1)H
(γ + 1)

(5g)

(

(5h)

3
− 4 + 5 f a2
16

)

(

(

M o = min 1, max M 2 , M ∞2

M2 =

Vn

+2

+ Vn
2c12 2

))

)

)

⎧1
if M ≥ 1⎪ (1 ± sign(M )),β± α =0 = ⎨ 2
1⎪ (M ± 1)2 (2 m M ), otherwise
⎩4

(6b)

(6c)

(6d)

であり，セル境界音速や参照速度等は下記で与えられる．

M L R = Vn± c~1 2 , c~1 2 = c~L ⋅ c~R2

f a (M o ) = M o (2 − M o )2Tと書け，a) 圧力流束 ~
pは

(5e)

であり，

c1 2 = min (c~L , c~R ),+c~L = c *2 max c * , V n ,−c~ = c *2 max c * , − V(5n)３－３．LDFSS2001 および LDFSS2001-2
最後に，Edwards らによって提案された LDFSS2001(11)について説明する．LDFSS2001 は，

if M ≥ 1otherwise(5m)である．ここで，Ku=0.75, Kp=0.25, σ=1.0 とする．また，一様流マ
ユーザが指定する必要がある．
ッハ数M∞は参照パラメタであり，
AUSM+-upでは圧力流束の散逸はSLAUや後述するLDFSS2001と異なり付加的な項 pu（式 5e）によって与えられる．この項は，M>1 ではゼロとなるため，あくまでも低速流れにおける散逸量を与えるためのものと解釈できる．SLAU2 と同様に超音速流れ（M>1）においてマッハ数に比例する散逸量を与えるためには，
pu 以外の項に修正を加える必要がある．これは圧力流束そのものを大幅に変更する事になり，様々な方法が考えられるが，ここでは SLAU2 と同じ式(4a)を採用する．そしてこのときの AUSM+-up
を AUSM+-up2 と呼ぶ事にする．

(5a)Ψ = (1, u , v, w, H ) , N = (0, n x , n y , n z , 0 )±M

Mp = −

３－２．AUSM+-up および AUSM+-up2
次に，Liou によって提案された AUSM 族スキームの最新版である AUSM+-up(10)について述べる．SLAU 等と同じく数値流束を

F1/ 2 =

if M 1 2 > 0otherwiseM1 2 = fM + fM + M p

c1 2 = min (c L , c R ),*2⎧ρ
m& = M 1 2 c1 2 ⎨ L
⎩ρ R

(

)

2

(5i)

2
Vn± 1 − M ref+ 4Vref2
c~L R =
2
1 + M ref

(5j)

2
M ref
= V ref2 c1 2 ,

(6e)

2

(

(

v2
2
V ref2 = min c1 2 , max V , V ∞2

−2

(5k)

v2
V = u 2 + v 2 + w2

)

(6f)

&は
である．また，b) 質量流束 m
3Copyright © 2010 by JSFM

c1 2 = min (c L , c R ),(c L R = c *2 max c * , Vn

c

*2

±

)

第 24 回数値流体力学シンポジウム
B12-3で与えられる．なお，参照速度 Vref はユーザ指定パラメタである．
ここで，圧力流束（式 6c）における SLAU2（式 3c）との違い
は第3項のみである．
この第3項はここでも散逸量を制御するが，
Vref の与え方から考えると，このままの形では SLAU と同様，M>1
では散逸量が一定となる．そこで LDFSS2001 に対しても M>1 で散逸量がマッハ数に比例するよう，第 3 項を SLAU2 のものに置き換える事を考える（表記上は，圧力流束をそのまま SLAU2 のもの（式 4a）に変える事と同義である）．この数値流束を，
LDFSS2001-2 と呼ぶ事にする．
このように，SLAU, AUSM+-up, LDFSS2001 の圧力流束に式(4a)を用いる事で，それぞれ SLAU2, AUSM+-up2, LDFSS2001-2 を構築した（図 3 参照）．以降，数値計算例によりこれらの有効性を示
す．

(6g)

2(γ − 1)=H
(γ + 1)

また b) 質量流束 m& + , m& − は

(

)

(

m& + = ρ L c~1 2 f M+ − f 1+2 , m& − = ρ R c~1 2 f M− + f 1−2

f

±
M

)

⎧1
⎪ (M ± M ), if M ≥ 1=⎨2
1⎪± (M ± 1)2 , otherwise
⎩ 4

(6h)

(6i)

４．極超音速流れの計算における衝撃波安定性および空力加熱分
布予測の両立
ここでは新しい圧力流束を用いた流束関数SLAU2, AUSM+-up2,
および LDFSS2001-2 の性質として，冒頭で述べた A)衝撃波異常への堅牢性，および B)空力加熱の正確な予測を調べるため，A)
著者らによる極超音速流れの 1.5 次元問題および 2 次元鈍頭物体
周りの流れの計算，そして B) 2 次元鈍頭物体壁面における空力加熱計算を行う．

(6j)

４－１．1.5 次元問題（非粘性，空間 1 次精度）
これは 1.5 次元問題と呼ばれ，2 次元空間に 1 次元衝撃波を配置し，選択した計算手法がそれを安定に捉えられるかを調べるた
めに著者らによって行われた数値実験である(3)．この問題は 2 次
元鈍頭物体よどみ点前方の流れを抽出したものと見なせるため，
ここでの結果は 2 次元鈍頭物体周りの計算の安定性，堅牢性にほぼ近い．ここでは問題の概要と結果のみ述べ，詳細は文献(3)に譲る．
図 4 に示すように，計算格子は擾乱が全く入っていない 50×25⎛ [Δp + δ Δp ] ⎞⎟,f 1+2 = f 1 2 ⎜1 −
2
⎜
⎟
2
ρ
V
L
ref
,
1
2
⎝
⎠
⎛ [Δp − δ Δp ] ⎞⎟f 1−2 = f 1 2 ⎜1 +
2
⎜
⎟
2
ρ
V
R ref ,1 2 ⎠⎝Δp = p L − p R , δ ∈ [0, 4]
2⎧1 ⎡ 1⎤⎪ ⎢M L2 + M R2 − 1⎥ , if M L < 1, M R < 1⎪⎦
f1 2 = ⎨ 4 ⎣ 2⎪⎪ 0,otherwise⎩()

Fig. 3 Van Leer’s FVS and AUSM-family schemes.
4

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
B12-3•‘0’は計算が不安定(unstable)となった事を示し，このときは残差がある程度の大きさからは下がらなかった．この場合に得られた解は，衝撃波が±2 セルを超える振動もしくは非対称性を示し，多くの場合，衝撃波が崩壊した（カーバン
クル解）．
今回得られた結果については，表 1 にまとめた．表の最右列の数字が大きい程，その手法が衝撃波に対し堅牢という事になる．
この結果より，
SLAU2, AUSM+-up2, LDFSS2001-2 (δ=4)はいずれもオリジナル手法よりも衝撃波に対し堅牢であり，Van Leer の FVSと同等の安定性を持つ事が分かる（以降，特に断らない限り．
LDFSS2001-2 では散逸パラメタδ=4 とする）

個のセルから成る等間隔格子である．衝撃波の左(L: i≤12)と右(R:i≥14)の状態の初期条件は垂直衝撃波に対するランキン・ユゴニオ(Rankine-Hugoniot)条件を満たすように与えられる．
ishock=12+ε

M∞=6
j

i

４－２．2 次元鈍頭物体周りの流体計算（非粘性，空間 1 次精度）
1.5 次元問題において示された新手法の A)衝撃波異常への堅牢性を確認するため， 2 次元の鈍頭物体周りの流体計算を行う．気流条件は名古屋大学衝撃風洞を模擬し（表 2），
円柱半径は r=20mmとする．計算格子は，160（壁に沿う方向）×80（壁垂直方向）セルから成り，各方向において格子間隔は変化させていない．計算格子および座標系を図 6a, b にそれぞれ示す．

Fig. 4 Computational grid and conditions for 1.5D problem.
ここで衝撃波位置パラメタεを導入し，ε =0.0, 0.1, ..., 0.9 とする．
例えばε =0.0 のとき，衝撃波はちょうどセル i=12 と i=13 の間のセル境界に初期配置される．一方でε =0.5 の場合には，i=13 のセル
中心上に衝撃波が存在する事になる．一様流マッハ数は M∞=6.0とする．図 4 の j 方向境界には周期条件を課している．
数値計算結果の例を図 5 に示す．
a) Good (2)b) Fair (1)Table 2 2D blunt-body flow conditionsc) Poor (0)Fig. 5 Typical solutions for 1.5D test, a) Good (2), b) Fair (1), and c)

計算は空間，時間共に 1 次精度で行い，陽解法で CFL=0.5 にて時間発展させた．いずれのケースも 30,000～40,000 ステップ程度で収束解が得られた．流れ場の可視化結果を図 6c-i に示す．新しい圧力流束を用いた 3 つの手法 SLAU2, AUSM+-up2, およびLDFSS2001-2，そして SLAU ではいずれも衝撃波の形状が滑らかに捉えられている（図 6c-f）．一方で図 6g の AUSM+-up では計算
格子線の影響を受けて衝撃波形状がやや振動しており，図 6h のLDFSS2001 では流れ場が非対象となっている．よってこれらの手

Poor (0).
••.1
1.31×105
370.6
63.73

M∞
Rer
P∞ [Pa]T∞ [K]‘2’は安定(stable)かつ上下対称な解が得られた事を示し，この場合には残差が 3 桁以上低下し，指数関数的に収束に向かった．
‘1’は衝撃波が±2 セルの限られた範囲において振動(oscillatory)もしくは上下非対称性(asymmetry)を示した場合を指す．

Table 1 1.5D test results for various numerical flux functions (ε: shock location parameter (3); the initial shock is imposed exactly on a cell-interface whenε=0.0, and at the cell-center when ε=0.5.)
Numerical Flux Functions

ε=0.0

0.1

0.

0 3

0.4

0.5

0.6

0.7

0.8

0.9

Total

SLAU2

2

2

2

2

2

2

2

2

2

2

20

AUSM+-up2

2

2

2

2

2

2

2

2

2

2

20

2

2

2

2

2

2

2

1

2

2

19

2

2

2

2

2

2

2

2

2

2

20

SLAU

1

2

2

2

1

1

1

1

1

1

13

AUSM+-up

2

2

2

2

2

1

1

1

1

2

16

2

2

2

2

2

2

1

1

1

2

17

2

2

2

2

2

2

1

1

2

2

18

Van Leer FVS

2

2

2

2

2

2

2

2

2

2

20

Roe

0

2

2

2

0

0

0

0

0

0

6

LDFSS2001-2

LDFSS2001

δ=1
δ=4

δ=1
δ=4

5

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
B12-3b) Coordinatesa) Gridφc) SLAU2d) AUSM+-up2h) LDFSS2001i) Roe (E-fix)e) LDFSS2001-2yxM∞=8.1ijg) AUSM+-upf) SLAU

Fig. 6 Hypersonic, inviscid, 2D blunt-body problem a) grid (160×80; every other grid lines are shown), b) coordinates, and results (Cp contours ofconverged solutions) of c) SLAU2, d) AUSM+-up2, e) LDFSS2001-2, f) SLAU, g) AUSM+-up, h) LDFSS2001, i) Roe (E-fix).
b) Coordinatesa) GridφM∞=8.1d) AUSM+-up2c) SLAU2e) LDFSS2001-2

y
x

i

j

Fig. 7 Hypersonic, viscous, 2D blunt-body problem a) grid (160×160; every other grid lines are shown), b) coordinates, and results (Mach number contours
at 100,000 steps) of c) SLAU2, d) AUSM+-up2, e) LDFSS2001-2.
6

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
B12-3ても，この性質が損なわれていない事を確認する．
計算対象としては NACA0012 翼を扱う．この周りの M=0.01 の
非粘性流れを解き，抵抗が十分小さい事を確かめる．本来，抵抗はゼロとなるべきであり，数値計算で現れる抵抗値は誤差の大きさとみなせる(19)．
数値流束にはSLAU2, AUSM+-up2, LDFSS2001-2, そして比較のため SLAU, Roe を用い，時間積分には LU-SGS もしくは WeissSmith の前処理付き LU-SGS（pLU-SGS）(13)を用いた．CFL=20 とし，2,000 ステップ計算を行ったところ，密度残差は約 5 桁低下した．
このときの可視化結果を図 9 に，抵抗値を表 3 に示す．
Roe/LU-SGS の結果では非物理的な圧力分布が現れてしまい，大きな抵抗値（誤差）が生じてしまっているが，新スキームではいずれも誤差が Roe/LU-SGS の 5%未満と小さく，SLAU/pLU-SGSと同等であり，良好な解が得られている（ちなみに Roe/pLU-SGSの組み合わせでは計算が発散する）．従って新手法は，オリジナルの全速度スキームと同じく低速流れへの適用が可能である．

法は，新しい圧力流束により堅牢となり，衝撃波をより滑らかに捉えられるようになったと言える．なお，図 6i の Roe (E-fix) (14),(15)ではカーバンクル現象が発生し明らかに異常な解が現れた．
４－３．2 次元鈍頭物体の空力加熱（粘性，空間 2 次精度）
次に，極超音速流中に置かれた円柱壁面の熱流束を考える．気流条件は非粘性の場合と同じく表 2 の通りとし，円柱半径は（等温壁条件），r=20mm である．
今回，
物体壁温度はTw=300K とし壁近傍に格子を集め，格子レイノルズ数（最小格子幅に基づくレイノルズ数）は Recell=1.3 とした．計算格子および座標系を図 7a, bにそれぞれ示す（160×160 セル）．
計算には流束関数として SLAU2, AUSM+-up2, LDFSS2001-2 を使用した．MUSCL(16)により空間 2 次精度とし（Van Albada の制限，時間発展に関数(17)を使用；粘性項も中心差分により 2 次精度）
は LU-SGS を用いて CFL=200 で 100,000 ステップ行った．このとき，残差（密度の L2 ノルム）は初期のほぼ 3 桁低下した．可視化結果を図7c-eに示す．
手法によって解にやや違いがあるものの，いずれも基本的な流れ場は再現できていると言える．
対応する壁面熱流束分布を図 8 に示す．横軸に当たる角度φは
図 7b で定義されており，φ=0˚はよどみ点である．このとき，熱流束値は，よどみ点における Fay-Riddell(18)の値（qFR =17.5 W/cm2）で規格化している．この図より，冒頭でも述べたが，SLAU2 はA) 衝撃波に対し堅牢であり，かつ，B-1)全エンタルピ保存，B-2)
境界層解像の性質を有しているにも関わらず壁面熱流束分布は波打ってしまっている．一方で，A) の性質を持つ SLAU2 の圧力流
束を B-1)全エンタルピ保存，B-2)境界層解像，そして B-3) 滑らかな熱流束分布を出力，の性質を併せ持つ他の質量流束と組み合わせた AUSM+-up2, LDFSS2001-2 では，実際に今回の問題で良好な熱流束分布が得られている．
新しい圧力流束（式
従ってAUSM+-up2 およびLDFSS2001-2 は，(4a)）を用いる事でオリジナルの手法の衝撃波への堅牢性を向上させ，他の性質は残した事で，A), B-1), B-2), B-3)の全ての性質を満足する事が確かめられた．

a) SLAU2/ pLU-SGSb) AUSM+-up2/ pLU-SGS

M∞=0.01c) LDFSS2001-2/ pLU-SGSd) Roe/ LU- SGS(19)1.4

Fig. 9 Cp contours (0 < Cp < 390) for low speed flow over
NACA0012 airfoil (M∞=0.01), a) SLAU2/ pLU-SGS, b) AUSM+-up2/
pLU-SGS, c) LDFSS2001-2/ pLU-SGS, d) Roe/ LU- SGS(19).

1.2

q/q FR

1.0
0.8

SLAU2

0.6

Table 3 Computed drag coefficients CD for NACA0012 airfoil.

AUSM+-up2
0.4

Numerical Flux/
Time Integration Method
SLAU2/ pLU-SGS
AUSM+-up2/ pLU-SGSLDFSS2001-2/ pLU-SGSSLAU/ pLU-SGS
Roe/ LU- SGS

LDFSS2001-2

0.2
0.0
-75

-60

-45

-30

-15

0

15

30

45

60

75

φ [deg.]Fig. 8 Aerodynamic heating profiles over blunt-body at Mach 8.1.
５．数値計算例
新しいスキームを用いた数値計算例として，低速流れ，3 次元における鈍頭物体空力加熱問題，そして実問題への適用として複雑形状から成るロケット周りの数値解析例を示す．

CD
0.00320.00370.00480.0037(19)0.0720(19)５－２．3 次元鈍頭物体の空力加熱（粘性，空間 2 次精度）
ここでは NASA LaRC（ラングレー研究所）の Gnoffo 博士らが
LAURA コードなどで用いている M∞=17 の検証問題(20),(21)を扱う．
計算条件は V∞ = 5000 m/s, Re = 376,930 /m, ρ∞= 0.001 kg/m3, T∞ =200 K, Tw = 500 K である．計算格子は 30×64×10 セル，壁近傍の格子レイノルズ数は Recell = 1.00 であり，Gnoffo 博士らにご提供いただいた（私信，2009）．CFL=200 で 100,000 ステップ計算したと
ころ，AUSM+-up2 では密度残差が 4 桁低下，LDFSS2001-2 および Roe (E-fix)ではマシンゼロまで低下した．

５－１．2 次元翼周り低速流れの計算（非粘性，空間 2 次精度）
ここでは新しい手法の低速流れにおける性質について調べる．
圧力流束を変更する前の流束関数はいずれも全速度スキーム(All-Speed Schemes)として低速流れの計算に対応できるように作られていた．このため，圧力流束変更後の新しいスキームについ7

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
B12-3a) AUSM+-up2b) LDFSS2001-2c) Roe (E-fix)(21)

a)
1.2
1.0

q/q FR

0.8
0.6

AUSM+-up2, 3D

0.4

AUSM+-up2, 2D

0.2

Ref. (LAURA), 2D

0.0
-90

-75

-60

-45

-30

-15

0

15

30

45

60

75

90

45

60

75

90

φ [deg.]

b)
1.2
1.0

M∞=17

q/q FR

0.8
0.6
LDFSS2001-2 (d=4), 3D

0.4

LDFSS2001-2 (d=4), 2D

0.2

Ref. (LAURA), 2D

0.0
-90

-75

-60

-45

-30

-15

0

15

30

φ [deg.]c)
2.0
Roe (E-Fix), 3D
Roe (E-Fix), 2D

1.5

q/q FR

Ref. (LAURA), 2D
1.0

0.5

Fig. 10 Pressure (top; 0 < P/P∞ < 390) and Mach number (bottom; 0< M∞ < 17) contours for 3D circular-cylinder (second-order in space;freestream Mach number M∞=17), a) AUSM+-up2, b) LDFSS2001-2,c) Roe (E-fix)(21).

0.0
-90

-75

-60

-45

-30

-15

0

15

30

45

60

75

90

φ [deg.]Fig. 11 Aerodynamic heating profiles over blunt-body at Mach 17, a)
AUSM+-up2, b) LDFSS2001-2, c) Roe (E-fix)(21).

図 10，図 11 に計算結果を示す．まず図 10 より，この問題に
AUSM+-up2 または LDFSS2001-2 を用いれば，対称かつ滑らかな衝撃波形状が得られるが，Roe (E-fix)を用いると，スパン方向に空間的振動が生じてしまう事が分かる．
空力加熱分布についても同様である．図 11 に今回の計算結果
を，比較のため対応する 2 次元計算を行った結果と，LAURA コードによる計算結果 (20) と共に示した．AUSM+-up2 およびLDFSS2001-2では対称，
滑らかかつ正確な空力加熱分布が得られ，2 次元の解，LAURA コードの解やよどみ点の Fay- Riddell 解析値とほぼ一致している．図中では全壁面要素における熱流束を示しているが，スパン方向への値の変動は見られない．一方で Roe(E-fix)ではスパン方向の（3 次元的な）振動が現れており，他の解とは大きな誤差が生じてしまっている．

子を作製し（図 12，約 6,400 万セル）
，これを同じく JEDI 内製の非構造格子ソルバーLS-FLOW(24),(25)で計算を行った．計算条件は対応する風洞試験（以下，風試）に合わせて M∞=1.5, Re=1.2×107，
迎角 5 度，ロール角 60 度である．
流束関数に SLAU2, AUSM+-up2 または LDFSS2001-2 を用い，Green-Gauss 法により空間 2 次精度とし，Venkatakrishnan の制限関数(26)を利用した．粘性項は Wang の方法(27)で 2 次精度で算出し，乱流モデルには S-A モデル(28)を用い，LU-SGS と局所時間刻み（CFL=1）を利用して定常計算した．計算には JAXA スーパーコ
ンピュータ JSS（M システム）を利用し，64 並列にて実行した．
以降，風試結果(29)や従来法(SLAU, AUSM+-up)による計算結果(30)との比較を行う．
図 13 に計算結果の可視化結果を示す．流れ場の詳細な議論は既に文献(30)で行われているので割愛するが，いずれの結果も各突起部から発生する衝撃波や境界層剥離の様子，渦の形成など流れ場の詳細を捉えている．機体重心付近の突起（B2 ケーブルダク

５－３．複雑形状周り超音速流れの計算（粘性，空間 2 次精度）
最後に，現在 JAXA で開発中のイプシロンロケットを対象とした計算例を示す．イプシロンロケットは複雑形状を有するため，
JAXA/JEDI にて開発中の LS-GRID(22),(23)を用いて周囲に非構造格
8Copyright © 2010 by JSFM

(B2) Cable Duct

第 24 回数値流体力学シンポジウム
B12-3ト）を乗り越えた流れが機体に衝突する位置で機体表面に局所的な圧力上昇が見られ，更にその流れが下流の SMSJ と呼ばれる突起で馬蹄渦を形成し分岐して流れ去る様子が流線で描かれている．
これらの様子や衝撃波，境界層剥離の強さ，位置について，計算手法による目立った差異は見られない．また，機体先端部におけるカーバンクルも見られない．
図 14 では，機体のうち手法の影響を受け易いと考えられる下流部の突起（SMSJ）付近の表面流線パターンを示している．ここでも，手法による本質的な流れ場の違いは見られない．また，いずれの手法を用いた場合も，SMSJ 周囲で流れがその外側へ向かう様子など，オイルフロー可視化試験結果(29)と良く一致する．
更に，各ケースで得られた全機ローリングモーメント Cmx の値を表 4 に示した．全機ローリングモーメントは SMSJ やケーブルダクトなどの各突起によってお互い逆方向に作られる力が相殺し合った結果生まれ，風試で計測できない程小さい値となるため，計算方法への感度が非常に高いと考えられる．それでも表 4 に

SMSJ

Fig. 12 Computational grid (Cartesian/body-fitted hybrid,
unstructured grid) generated by using LS-GRID for Epsilon Launch
Vehicle(30).
a) SLAU2

M∞=1.5b) AUSM+-up2

Fig. 13 Computational results for rocket configurations at Mach 1.5, a) SLAU2, b) AUSM+-up2, c) LDFSS2001-2, d) SLAU(30). Body: Cp contours(-0.5<Cp<0.5) with body-constraint streamlines; Flow: absolute density gradient with streamlines (continued).
9

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
B12-3c) LDFSS2001-2d) SLAU(30)Fig. 13 Computational results for rocket configurations at Mach 1.5, a) SLAU2, b) AUSM+-up2, c) LDFSS2001-2, d) SLAU(30). Body: Cp contours(-0.5<Cp<0.5) with body-constraint streamlines; Flow: absolute density gradient with streamlines (concluded).
よれば，Cmx はいずれの手法を用いた場合でも，SLAU の場合(1.57×10-4)とほぼ一致する事が分かり，少なくとも異常な解には至っていないと考えられる．
以上より，本稿で紹介した新しい圧力流束を用いた手法は，複雑形状の実問題にも適用可能であると言える．

６．まとめ
AUSM 族スキームの新しい圧力流束を提案し，これを 3 種の
AUSM 族全速度スキーム SLAU, AUSM+-up, LDFSS2001 に適用した．こうして構築された新スキーム SLAU2, AUSM+-up2, およびLDFSS2001-2 は下記の性質を満足する．
高速流れA) 衝撃波異常への堅牢性（Van Leer の FVS と同等）
B-1) 全エンタルピ保存
B-2) 境界層解像
B-3) 滑 ら か な 熱 流 束 分 布 を 出 力 （ AUSM+-up2 お よ び
LDFSS2001-2）
これらを満たす事で，極超音速空力加熱が精度良く計算できる．
また，Table 4 Computed rolling moment coefficients Cmx for Epsilon
Launch Vechicle.
Numerical Flux
SLAU2
AUSM+-up2
LDFSS2001-2
SLAU
AUSM+-up

Cmx
1.65×10-4
1.70×10-4
1.68×10-4
1.57×10-4(30)1.69×10-4
10Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
B12-3a) SLAU2b) AUSM+-up2c) LDFSS2001-2d) Experiment(29)SMSJ

Fig. 14 Computational and reference experimental results for rocket configurations at Mach 1.5 (close-up views for SMSJ), a) SLAU2, b) AUSM+-up2,c) LDFSS2001-2, and d) oilflow experiment(29). Cp contours (-0.5<Cp<0.5) with body-constraint streamlines.
イプシロンロケットの解析は JAXA スーパーコンピュータ JSS（M システム）にて実行した．またこのとき利用した JAXA/JEDI
内製コードLS-GRIDおよびLS-FLOWは同センター所属の藤本圭一郎氏，葛生和人氏と共同で開発し，アイオワ州立大学 Z.J. Wang
教授，㈱計算力学研究センターの青野淳也氏のご協力をいただいた．藤本氏には本稿の内容についてもご意見をいただいた．
また本研究で用いたイプシロンロケット風洞試験データの取得には，JAXA/JEDI 藤本氏，葛生氏，宇宙科学研究所の野中聡助
教，入門朋子氏，大屋友児氏，大山聖助教（現 准教授）
，イプシ
ロンロケット研究チームの福添森康氏，三保和之氏，田村誠氏のご協力をいただいた（風洞試験を実施した 2009 年 5 月当時）．
ここにご協力いただいた全ての皆様へ感謝の意を表する．

低速流れC) 低速流れにおける低散逸の性質も有するため，全速度スキームとしての有用性は残されている．
数値計算例によりこれらの手法が幅広い速度域，そして複雑形状を有する実問題に利用できる事が示された．
ただし，AUSM+-up2 および LDFSS2001-2 には低速流れ用に一様流マッハ数等の参照パラメタが含まれており，内部流れ等でこれらを定義する事は困難である．一方で，SLAU2 にはこれらのパラメタが不要である．
更に，AUSM+-up2 および LDFSS2001-2 は Van Leer の FVS と同等に堅牢だが，衝撃波異常が全く起きないわけではない．近年研究が進められている多次元散逸項(20),(31)-(37)を導入すれば，より堅牢な手法を構築できる可能性がある．
従って今後は，AUSM+-up2 および LDFSS2001-2 より更に衝撃波に対し堅牢であり，かつ SLAU2 のように低速流れ用の参照パラメタを要しないような統合的手法の開発を行う．そして将来的には，多相流への拡張(38)等を行い，統合的手法の適用性を拡げて行く予定である．
謝辞
本研究は NASA LaRC（ラングレー研究所）の Dr. Hiroaki
Nishikawa（附属機関 NIA 所属）および Dr. Jeff White との議論，
そしてJAXA 宇宙科学研究所の野々村拓氏のご指摘をきっかけとして始まったものである．
本研究の中で行った 1.5 次元問題は，第一著者が米国ミシガン
大学 Philip L. Roe 教授の下で実施した研究である．LAURA ベン
チマーク（M∞=17 の円柱）のデータは NASA LaRC の Dr. Peter
Gnoffo よりご提供いただいた．
本研究で用いた計算コードの一部は，名古屋大学中村佳朗教授の研究室で開発された．
11

参考文献(1) Peery, K.M. and Imlay, S.T., “Blunt-Body Flow Simulations,”
AIAA Paper 88-2904 (1988).
(2) Pandolfi, M. and D’Ambrosio, D.: Numerical Instabilities in
Upwind Methods: Analysis and Cures for the “Carbuncle”
Phenomenon, Journal of Computational Physics, Vol. 166, No. 2,
2001, pp.271-301.
(3) Kitamura, K., Roe, P., and Ismail, F.: Evaluation of Euler Fluxes forHypersonic Flow Computations, AIAA Journal, Vol. 47, No.1, 2009,
pp.44-53.
(4) Kitamura, K., Shima, E., Nakamura, Y., and Roe, P.: Evaluation ofEuler Fluxes for Hypersonic Heating Computations, AIAA Journal,
Vol.48, No.4, 2010, pp.763-776.
(5) Henderson, S.J.. and Menart, J.A., “Grid Study on Blunt Bodies
with the Carbuncle Phenomenon,” AIAA Paper 2007-3904 (2007).
(6) Mazaheri, A.R. and Kleb, B., “Exploring Hypersonic,
Unstructured-Grid Issues through Structured Grids,” AIAA Paper
2007-4462 (2007).
Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
B12-3(7) Van Leer, B.: Flux Vector Splitting for the Euler Equations, Lecture
Notes in Phys., Vol. 170, 1982, pp. 507-512.
(8) Liou, M.-S.: A Sequel to AUSM: AUSM+, Journal ofComputational Physics, Vol. 129, 1996, pp.364-382.
(9) Kitamura, K., and Shima, E., “Improvements of Simple
Low-dissipation AUSM against Shock Instabilities in considerationof Interfacial Speed of Sound,” Paper No.1283, Proceedings of
ECCOMAS CFD 2010, Lisbon, Portugal (2010).
(10) Liou, M.S.: A Sequel to AUSM, Part II: AUSM+-up for All Speeds,
Journal of Computational Physics, Vol. 214, 2006, pp. 137-170.
(11) Edwards, J.R., “Towards Unified CFD Simulation of Real Fluid
Flows,” AIAA Paper 2001-2524 (2001).
(12) Shima, E. and Kitamura, K., “On New Simple Low-Dissipation
Scheme of AUSM-Family for All Speeds,” AIAA Paper 2009-136(2009).
(13) Weiss, J.M. and Smith, W.A.: Preconditioning Applied to Variable
and Constant Density Flows, AIAA Journal, Vol. 33, No.11, 1995,
pp. 2050-2057.
(14) Roe, P.L.: Approximate Riemann Solvers, Parameter Vectors, andDifference Schemes, Journal of Computational Physics, Vol. 43,
1981, pp.357-372.
(15) Harten, A.: High Resolution Schemes for Hyperbolic Conservation
Laws, Journal of Computational Physics, Vol. 49, 1983,
pp.357-393.
(16) Van Leer, B.: Towards the Ultimate Conservative Difference
Scheme. V. A Second-Order Sequel to Godunov’s Method, Journalof Computational Physics, Vol. 32, 1979, pp.101-136.
(17) Van Albada, G.D., Van Leer, B. and Roberts, Jr., W.W.: A
Comparative Study of Computational Methods in Cosmic Gas
Dynamics, Astron. Astrophys., Vol. 108, 1982, pp.76-84.
(18) Fay, J.A. and Riddell, F.R.: Theory of Stagnation Point Heat
Transfer in Dissociated Air, J. Aeronautical Sciences, Vol. 25, 1958,
pp.73-85.
(19) Kitamura, K., Shima, E., Fujimoto, K., and Wang, Z.J.:
Performance of Low-Dissipation Euler Fluxes and Preconditioned
LU-SGS at Low Speeds, Communications in Computational
Physics (Accepted for Publication).
(20) Gnoffo, P.A., “Multidimensional, Inviscid Flux Reconstruction for
Simulation of Hypersonic Heating on Tetrahedral Grids,” AIAAPaper 2009-599 (2009).
(21) Kitamura, K., Shima, E., and Roe, P., “Three-Dimensional
Carbuncles and Euler Fluxes,” AIAA Paper 2010-1273 (2010).
(22) Fujimoto, K.: Study on the Automated CFD Analysis Tools forConceptual Design of Space Transportation Vehicles, PhD
Dissertation, Univ. of Tokyo, 2006.
(23) Fujimoto, K., Fujii, K., Wang, Z.J., Kitamura, K., “Advanced
Cartesian Grid-based CFD Approach for High Reynolds number
Turbulent Flows,” IMECE2008-67110, ASME International
Mechanical Engineering Congress & Exposition, Boston, MA(2008).
(24) Kitamura, K., Fujimoto, K., Shima, E., Kuzuu, K., and Wang, Z.J.:
Validation of an Arbitrary Unstructured CFD Code for
Aerodynamic Analyses, Transactions of the Japan Society for
Aeronautical and Space Sciences, (Accepted for publication).
(25) Kuzuu, K., Kitamura, K., Fujimoto, K., and Shima, E., “Numericalflow simulation of a reusable sounding rocket during nose-up
rotation,” MX.00005, American Physical Society 63rd Annual
Division of Fluid Dynamics (DFD) Meeting, Long Beach, CA, Nov.

2010.
(26) Venkatakrishnan, V.: Convergence to Steady State Solutions of theEuler Equations on Unstructured Grids with Limiters, Journal ofComputational Physics, Vol. 118, 1995, pp.120-130.
(27) Wang, Z.J.: A Quadtree-based Adaptive Cartesian/Quad Grid Flow
Solver for Navier-Stokes Equations, Computers and Fluids, Vol. 27,
1998, pp. 529-549.
(28) Spalart, P.R. and Allmaras, S.R., “A One-Equation Turbulence
Model for Aerodynamic Flows,” AIAA Paper 92-0439 (1992).
(29) 北村圭一，藤本圭一郎，野中聡，入門朋子，福添森康，嶋英
志：次期固体ロケット空力特性についての風洞試験，航空宇
宙技術，Vol. 9, 2010, pp.9-14．
(30) 北村圭一，藤本圭一郎，葛生和人，野中聡，入門朋子，福添
森康，嶋英志：次期固体ロケット空力特性についての数値解
析，航空宇宙技術（査読中）；もしくは“イプシロンロケッ
トの空力特性，”第 54 回宇宙科学技術連合講演会，1F05(2010).
(31) Loh, C.Y. and Jorgenson, P.C.E.: Multi-dimensional Dissipation forCure of Pathological Behaviors of Upwind Scheme, Journal ofComputational Physics, Vol. 228, 2009, pp.1343-1346.
(32) Yoon, S.H., Kim, C., and Kim, K.H.: Multi-Dimensional Limiting
Process for Three-Dimensional Flow Physics Analyses, Journal ofComputational Physics, Vol. 227, 2008, pp.6001-6043.
(33) Nishikawa, H. and Kitamura, K.: Very Simple, Carbuncle-Free,
Boundary-Layer Resolving, Rotated-Hybrid Riemann Solvers,
Journal of Computational Physics, Vol. 227, 2008, pp. 2560-2581.
(34) Shima, E. and Kitamura, K., “On AUSM-Family Scheme for All
Speeds with Shock Detection for Carbuncle-Fix,” AIAA Paper
2009-3544 (2009).
(35) Roe, P.L. and Kitamura, K., “Artificial Surface Tension to Stabilize
Captured Shockwaves,” AIAA Paper 2008-3991 (2008).
(36) Kim, K.H., Kim, C., and Rho, O.H.: Cures for the shock instability Development of a shock-stable Roe scheme, Journal ofComputational Physics, Vol. 186, 2003, pp.342-374.
(37) Kim, S.S., Kim, C., Rho, O.H., Hong, S.K.: Methods for the
Accurate Computations of Hypersonic Flows I. AUSMPW+
scheme, Journal of Computational Physics, Vol. 174, 2001,
pp.38-80.
(38) Liou, M.S., Chang, C.H., Nguyen, L., and Theofanous, T.G.: How
to Solve Compressible Multifluid Equations - A Simple, Robust,and Accurate Method, AIAA Journal, Vol. 46, No.9, 2008, pp.
2345-2356.

12

Copyright © 2010 by JSFM