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

弧状衝撃波不安定性の三次元数値シミュレーション
Three-dimensional numerical simulation of bow-shock instability
佐藤陽介，東北大院, 仙台市青葉区荒巻字青葉 6-6-01, E-mail：yosato@rhd.mech.tohoku.ac.jp
保江かな子，JAXA, 東京都調布市深大寺東町 7-44-1, E-mail：yasue@chofu.jaxa
大西直文，東北大院, 仙台市青葉区荒巻字青葉 6-6-01, E-mail：ohnishi@rhd.mech.tohoku.ac.jp

Yosuke Sato, Tohoku University, 6-6-01 Aramaki-Aza-Aoba, Aoba-ku, Sendai 980-8579, Japan
Kanako Yasue, Japan Aerospace Exploration Agency, 7-44-1 Jindaiji-Higashi, Chofu, Tokyo 182-8522, Japan
Naofumi Ohnishi, Tohoku University, 6-6-01 Aramaki-Aza-Aoba, Aoba-ku, Sendai 980-8579, Japan
Three-dimensional numerical simulations were conducted using discontinuous Galerkin (DG) ﬁnite-element
method to clarify mechanism of bow-shock instability which has been experimentally observed for a blunt body ina low-γ gas. Shock-surface formation was computed with a blunt body in Mach 3.9 ﬂow. A well-known numericalinstability, carbuncle phenomenon, often occurs simultaneously in such a ﬂow condition. Suﬃcient resolution of
stream-ward grids is required to obtain carbuncle-free solution. In a bow shock ahead of a cone, an instability
was observed, and it may be physical instability induced by the edge eﬀect.

1.

緒論

超音速飛行する鈍頭物体前方に形成される弧状衝撃波について, これまでに多くの実験や数値計算が行われてきた. 一般に, 一様流中に形成される弧状衝撃波は安定であることが知られている. しかしながら, 弾道法における実験において, 比熱比の小さな気体中では, 鈍頭物体前方に形成される弧状衝撃波が不安定であることが確認されている(1, 2) . 様々な条件による実験から, この不安定性はマッハ数や気体の圧力, 物体の大きさ, 縁の丸み, 前面の曲率に依存して発生することが分かっている. Fig. 1 にBaryshnikov らによる実験で得られた不安定な弧状衝撃波を示す. この現象の要因として, 比熱比が小さいことによる衝撃波後方の力学的非平衡と, 急激な温度上昇による化学的非平衡が考えられているが, 実験では衝撃波後方の流れ場の構造を詳しく解析することが困難であるため, 原因は解明されていない.
現在, 超音速で飛行する物体の解析には数値流体力学(CFD) がよく用いられる. CFD を用いることによって,
衝撃波後方の流れ場の構造を詳しく解析することができ,
不安定性の原因を解明することが可能になると考えられる. 一方で, 強い衝撃波が形成される数値計算において,衝撃波が計算格子に平行に近い場合に, 衝撃波捕獲スキー
ムは数値的に不安定になることが知られている. この不安定性はカーバンクル現象と呼ばれる(3) . カーバンクル
現象は主流条件, 計算格子, 計算手法等の様々な要因によって生じるため, 適切な計算条件を整えることが重要である.
本研究では, 鈍頭物体前方に形成される弧状衝撃波不安定性のメカニズムの解明を目指し, 非構造格子を用いた
Discontinuous Galerkin (DG) 法(4) による三次元数値流体解析を行う. カーバンクル現象の存在により, 発生した
不安定性が物理的なものか数値的なものかの判別が難しくなる. そこで, カーバンクル現象を生じにくいと考えられる数値流束関数を用いて計算する(5) . また, 計算格子には, カーバンクル現象を生じにくくするために主流方
向の解像度が高い計算格子を用いて計算を行う(6, 7) . 本
稿では, 物体形状が不安定性に与える影響に焦点を当て,物体の曲率による影響や角を含む物体について数値計算
結果に基づく考察を行う.

2.

(a) large-amplitude
perturbation(b) complete destruction

Fig. 1: 不安定な弧状衝撃波(1, 2) .
件と計算手法を用い, 計算格子のみを変更してその影響を調べる必要があるため, 計算格子には物体の形状適合性が高く, 格子の配置を比較的自由に変更できる非構造格子を用いる. 計算手法としては, 有限要素法の一種で,高次精度を達成できる DG 法を用いて数値計算を行う.
2.1 Discontinuous Galerkin 法
基礎方程式
三次元 Euler 方程式は次式で与えられる.
∂F
∂G
∂Q ∂E
+
+
+
=0
(1)
∂t
∂x
∂y
∂z
ここで, Q は保存関数, E, F , G はそれぞれ x, y, z 方向の流束関数であり, 以下のように表される.
 
ρ
 
 ρu  Q= ρv 
 
ρw







e




ρu
ρv
ρw
 2





 ρu + p  ρuv  ρuw 







 2



E=
 ρuv  , F =  ρv + p  , G =  ρvw 




 2

 ρuw  ρvw  ρw + p (e + p)u(e + p)v(e + p)w

数値計算法

弧状衝撃波の物理的不安定性を調べるために, 上で述
べたカーバンクル現象を防ぐ必要がある. カーバンクル
現象を検証するには, 同形状の物体について, 同じ主流条(2)

1

Copyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
B12-5ρ は密度, u, v, w はデカルト座標系における速度成分, eは単位体積あたりの全エネルギーである. また, 圧力 p が次式で表される理想気体における状態方程式に従うことを仮定する.
{}
1 2
2
2
p = (γ − 1) e − (u + v + w )(3)2

直交基底関数
前述のように, 基底関数に直交性があれば, 自由度ごとに独立な時間発展の方程式が得られる. 非構造格子を用いて計算を行う場合には, 各セルの物理空間における座標 (x, y, z) を基準セルの計算空間の座標 (ξ, η, ζ) およびテンソルセルの写像空間座標 (r, s, t) に座標変換を行うことで直交性が実現される. 三次元の計算セルにおけるテンソルセル座標の直交基底関数 φlmn (r, s, t) および
基準セル座標の直交基底関数 φlmn (ξ, η, ζ) は Jacobi 多
項式 Pγα,β (z) を用いて記述される.(8) 三角柱セルの場合
は次式のようになる.

ただし, γ は比熱比である.
DG 法は有限要素法の一種であるので, 次式のような弱形式にした積分方程式を解くことで解を得る.
()∫∂Q ∂E ∂F ∂G
w(x, y, z)+++dΩ = 0 (4)
∂t
∂x
∂y
∂z
Ω

0,0
φlmn (r, s, t) =Pl0,0 (r)Pm(s)(1 − t)l Pn2l+1,0 (t)
(
)
(11)
1+ξ
0,0
0,0
φlmn (ξ, η, ζ) =Pl2− 1 Pm(η)(1 − ζ)l1−ζΩ は計算空間での領域, w(x, y, z) はテスト関数を表す.
式 (4) に対して Gauss の発散定理を適用し, 展開すると次式を得る.
∫Ω∂wQ∂t

ここで, Jacobi 多項式 Pγα,β (z) は次式で与えられる.

∫Pγα,β (z) =dΩ +w (Enx + F ny + Gnz ) d∂Ω
∂Ω
)
∫ (
∂w
∂w
∂w
−
E
+F
+G
dΩ = 0 (5)
∂x
∂y
∂z
Ω

∫
K2i =

Ω

∂t

dΩ

E(Qh )Ω∂vi∂x+ F (Qh )∂viQ̂ = J −1 Q

(7)

∂y

+ G(Qh )

∂vi
∂z

)

Ê = J −1 (ξx E + ξy E + ξz E)
F̂ = J −1 (ηx F + ηy F + ηz F )J −1 は座標変換におけるヤコビアンであり次式で表さ
れる.

dΩ

j

dt

vi vj dΩ =Ω∂(x, y, z)

J −1 =

さらに式 (6) を K1i に代入すると∑ dQj ∫

∑
j

Iij

dQj
dt

(14)

Ĝ = J −1 (ζx G + ζy G + ζz G)

(8)

(9)

K1i =

(13)

ここで, Q̂ は保存関数, Ê, F̂ , Ĝ はそれぞれ ξ, η, ζ 方向の流束関数である.

∂Ω∫ ([(1 − z)α+γ (1 + z)β+γ ]∂ Q̂ ∂ Ê
∂ F̂
∂ Ĝ
+
+
+
=0
∂t
∂ξ
∂η
∂ζ

vi (E(Qh )nx + F (Qh )ny + G(Qh )nz ) d∂Ω

K3i =

dz γ

離散化
式 (1) の物理空間 (x, y, z) における支配方程式を計算
空間 (ξ, η, ζ) へ変換すると, 次式のように記述できる.

DG 法 で は, テ ス ト 関 数 w(x, y, z) を 基 底 関 数
vj (x, y, z) とし, 式 (5) において厳密解 Q(t) を式 (6) で表される近似関数 Qh (x, y, z, t) に置き換える. そして, w
が vi のときの各積分項を以下のように定義する.
∂vi Qhdγ

三次元二次精度の場合, 自由度は 4 個必要となる. また,
添字 l, m, n は多項式の次数であり, (l, m, n) の組み合
わせは定数項では (0, 0, 0) となり, 1 次項では (1, 0, 0),(0, 1, 0), (0, 0, 1) となる.

j

∫

2γ γ!

(1 − z)−α (1 + z)−β(12)

ここで, ∂Ω は領域 Ω の境界, nx , ny , nz はそれぞれ境界における x, y, z 方向の法線ベクトルを表す. 有限要素法では, 次式のような自由度 Qj (t) と基底関数 vj (x, y, z)の積の線形和で表される近似関数 Qh (x, y, z, t) を用いてセル内の物理量の分布を記述する.
∑
Qh (x, y, z, t) =Qj (t) vj (x, y, z)
(6)

K1i =

(−1)γ

(15)

∂(ξ, η, ζ)

式 (13), (14) に対して DG 法を適用する. 弱形式にした
積分方程式は次式のようになる.

(10)

(

∫

となる
. ここで, 基底関数の積のモーメントを Iij =∫vvdΩ と表す. モーメント Iij は基底関数に直交性がijΩある場合には i = j 以外では 0 となる. その場合, 自由度の数だけ時間積分を独立に実行することが可能となり,自由度ごとに解を得る.

w(ξ, η, ζ)
Ω̂

∂ Q̂
∂t

+

∂ Ê
∂ξ

+

∂ F̂
∂η

+

∂ Ĝ
∂ζ

)
dΩ̂ = 0(16)

Ω̂ は計算空間での領域, w(ξ, η, ζ) はテスト関数を表す.
式 (16) に対して Gauss の発散定理を適用し, 展開すると
2

Copyright c 2011 by JSFM

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

次式を得る.

∫

∂wQ̂
Ω̂

∂t

を数値流束関数で置き換える必要がある. その際の流束計算には Riemann 解法を適用する(.
式 (20) において, 流束関数 Ê(Q̂h )nξ +F̂ (Q̂h )nη)+Ĝ(Q̂h )nζ はセル内部の物理量 Q−
h と隣のセルの物∫()
w Ênξ + F̂ nη + Ĝnζ d∂ Ω̂
∂ Ω̂
)
∫ (
∂w
∂w
∂w
Ê
dΩ̂ = 0
−
+ F̂
+ Ĝ
∂ξ
∂η
∂ζ
Ω̂

dΩ̂ +

−
+
理量 Q+
h を用いて表される流束関数 ĥ(Qh , Qh ) に置き換えられる. よって, 式 (20) は次式のように記述できる.
∫
+
K2i =
(26)
φi ĥ(Q−h , Qh )d∂ Ω̂

(17)

∂ Ω̂

∂ Ω̂ は計算空間領域 Ω̂ の境界, nξ , nη , nζ はそれぞれ境界における ξ 方向, η 方向, ζ 方向の法線ベクトルを表す.
ここで, 次式のような自由度 Q̂j (t) と前述した直交基底関数 φj (ξ, η, ζ) の積の線形和で表される近似関数
Q̂h (ξ, η, ζ, t) を用いてセル内の分布を記述する.
∑Q̂h (ξ, η, ζ, t) =Q̂j (t)·φj (ξ, η, ζ)(18)+
ここで, 本研究では, 流束関数 ĥ(Q−h , Qh ) に数値流束(9)SLAU (Simple Low-dissipation AUSM) を用いる.
K2i の積分の評価は, 各境界上に精度に見合った Gauss
点をとり, 二次元の Gauss 求積法を用いて面積分を行い
K2i の近似を行う. また, K3i の評価は, セル内に精度に見合った Gauss 点をとり, 三次元の求積法を用いて体積積
分により近似を行う. それにより式 (26) および式 (21) はそれぞれ次式のように置き換えられる.

j

式 (17) において厳密解 Q̂(t) を式 (18) で表される近似関数に置き換え, テスト関数を直交基底関数とする. i 番目の基底関数を乗じた各積分項を以下のように定義する.

∫∂φi Qˆh

K1i =

∫
K2i =

∂t

Ω̂

(

dΩ̂

Kˆ2i =

∑∑(
e

(19)

Kˆ3i =)∑(
Ê(Qˆh )l+
φi ĥl (Q−h , Qh )σl

l

∂φi
∂ξ

+ F̂ (Qˆh )

∂φi
∂η

)
(27)
e

+ Ĝ(Qˆh )

∂φi

(20)

ここで, e はセル境界面, l はセル境界面の Gauss 点, σl は各 Gauss 点における重みを表す. 二次精度の場合, 三角
形および四角形のセル境界面の Gauss 点の数は 4 個, セル内の体積ガウス点は 8 個となる.

Ê(Qh )Ω̂∂φi∂ξ+ F̂ (Qh )∂φi∂η+ Ĝ(Qh )

∂φi
∂ζ

)

dΩ̂

2.2 勾配制限法
本研究で定義した自由度 Qj は, 二次精度の場合には
Q1 がセル内の物理量の平均値, Q2 , Q3 , Q4 がセル内の物理量の一階微分を表す. よって, 勾配制限は Q2 , Q3 , Q4に対して以下のように行う.
φQ
2,3,4 , if pl max > pm max or pl min > pm min
Q̃2,3,4 =Q, otherwise(21)
式 (19)–(21) より, 式 (17) は次式のように記述できる.

K1i + K2i − K3i = 0(22)

また, 式 (18) を式 (19) に代入すると∑ dJ −1 Qj ∫
j

l

(28)

∫ (

K1i =

σl

φi Ê(Qh )nξ + F̂ (Qh )nη + Ĝ(Qh )nζ d∂ Ω̂
∂ Ω̂

K3i =

∂ζ

)

dt

φi φj dΩ̂ =

Ω̂

∑

Iij

j

dJ −1 Qj
dt

2,3,4

φ = min(φmax , φmin )pm max − pn
φmax =
pl max − pnpm min − pnφmin =
pl min − pn(23)となる. ここで, 基底関数の積のモーメントを Iij =∫φ φ dΩ̂ と表す. モーメント Iij は基底関数の直交性にΩ̂ i jより i = j 以外で 0 となるので, 以下のように表せる.

K1i = Iii

dJ −1 Qidt

ここで, pn はセル内の圧力, pl max , pl min はそれぞれセル境界の Gauss 点における圧力の最大値および最小値,
pm max , pm min はそれぞれ各面に隣接するセルの圧力の最大値および最小値を表す.

(24)

式 (22) と式 (24) より, 最終的に次式に対して時間積分を行うことで解を得る.

dQi
dt

=−

1
J −1 Iii(K2i − K3i )2.3 TVD Runge-Kutta 法
本計算で時間積分に用いた二段階 TVD Runge-Kutta
法(11) は以下に示す通りである.

(25)

Q(1) = Qn + ∆tR(Qn )
(
)
1
n
(1)
(1)
n+1
Q + Q + ∆tR(Q )Q=
2

K2i の評価において, DG 法では各セル毎に独立した近似関数を用いてセル内部の物理量を記述しているため, セル境界における流束関数の値は不連続となり, 直接求めることは出来ない. そのため, セル境界における流束関数

(29)
(30)

ここで, R は式 (25) の右辺を表す.

3

Copyright c 2011 by JSFM

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

3.

計算結果および考察

弧状衝撃波の力学的不安定性を調べるために, カーバ
ンクル現象を防ぐ必要がある. 以前の結果から, カーバン
クル現象を防ぐには, 主流方向の解像度を均一かつ高くすることが重要であることが分かった(6, 7) . 本研究では,そのような計算格子を用いて計算を行うことにより, 物
理的な不安定性を調べる.

3.1

物体の曲率による衝撃波面への影響

計算条件
物体の曲率が不安定性に与える影響を調べるために, 曲
率の異なる 2 つの物体について計算を行った. 計算対象として, 短径と長径の比がそれぞれ 1:10, 1:20 の楕円体を用いた. 計算格子はそれぞれセル数を 199,200, 749,500とし, 計算にはプリズムセルを用いた. カーバンクル現象を防ぐために, 上に述べたような主流方向の解像度が均一かつ高い格子を用いた. それぞれの計算格子を Fig. 2に示す. 主流条件は主流マッハ数を 3.9, 迎え角を 0, 比
熱比を 1.01 とした. また, CFL 数は 0.5 とした.

(a) 1:10得られた衝撃波面と亜音速領域のマッハ数分布
Fig. 3 にある断面における圧力分布を示す. 色は圧力の大きさを表す. 赤と青の境界が衝撃波面であり, どちらの曲率の物体を用いて計算を行ったときにも衝撃波面は定常になった.
弧状衝撃波不安定性において, 亜音速領域における乱れ
は衝撃波面に影響を与える可能性があると考えられる. 2つの物体における亜音速領域のマッハ数分布を Fig. 4 に示す. 短径と長径の比が 1:10 の場合に比べて 1:20 の場
合には亜音速領域に擾乱が残っていることが分かる. また, 1:10 の場合にも, わずかに非対称性が見られる. しかし, この亜音速領域の乱れは衝撃波面の歪みを招くまでには至っていない.
この 2 つの形状で不安定性は見られなかったが, より大きな短長径比をもつ楕円体を用いて計算を行った場合には, 亜音速領域のマッハ数分布の乱れはより大きくなることが予想される. しかし, 1:10 の場合と 1:20 の場合のマッハラインの位置の比較から, 短長径比が大きくなってもマッハラインの位置はほとんど変わらず, 衝撃波面に大きな変化が起きる可能性は低いと考えられる.

(b) 1:20
Fig. 3: 異なる短長径比をもった楕円体前方の圧力分布.

(a) 1:10(a) 1:10(b) 1:20
Fig. 4: 亜音速領域のマッハ数分布.
3.2 物体の角による衝撃波面への影響
前節で曲率の異なる 2 つの楕円体について計算を行ったが, 衝撃波面に不安定性は見られなかった. しかし, 短
長径比が 1:10, 1:20 の場合の比較から, より大きな曲率をもつ場合に亜音速領域の乱れが大きくなることが分かった. 次により大きな曲率をもつ極限として, 前面が平坦な物体を考える. そのような, 物体が角をもつ場合には, 物
体の角付近で生じた擾乱が音波として亜音速領域中を伝わり, 衝撃波面に達することで不安定性を引き起こす可能性がある. そこで, 角を含む物体についての計算結果を示す.

(b) 1:20
Fig. 2: 各計算格子の物体中心を通る断面図.

4

Copyright c 2011 by JSFM

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

計算条件
計算対象として円柱, 円錐を用いた. 計算格子は, それ
ぞれセル数が 896,340, 375,436 とし, 計算にはプリズムセルを用いた. カーバンクル現象を防ぐために主流方向の解像度が均一かつ高い格子を用いた. それぞれの計算
格子を Fig. 5 に示す. 主流条件は主流マッハ数を 3.9, 迎
え角を 0, 比熱比を 1.01 とした. また, CFL 数は 0.5 とした.

(a) 円柱(a) 円柱(b) 円錐Fig. 6: 角を含む物体の圧力分布.
(b) 円錐Fig. 5: 各計算格子の物体中心を通る断面図.

得られた衝撃波面と亜音速領域のマッハ数分布
Fig. 6 に流れ場がある程度定常になった後の圧力分布を示す. 色は圧力の大きさを表す. どちらの物体形状で計算を行ったときにも衝撃波面は定常になったが, 円錐を用
いて計算を行った場合は衝撃波面に歪みが見られる. 2 つの物体における亜音速領域のマッハ数分布を Fig. 7 に示
す. どちらの場合にも衝撃波内のマッハ数分布には非対称性が残った. また, どちらの物体を用いた場合にもマッ
ハラインが物体の角にかかり, 角を含まない物体に比べて分布の乱れが大きくなっていることが分かる. さらに,円錐においてはマッハラインに歪みが見られる. 前節において楕円体を用いて計算した結果や円柱を用いた結果との比較から, この不安定性には円錐の角付近が影響していると考えられる. この衝撃波面の歪みは, Baryshnikov
らによる実験で得られた不安定性 (Fig. 1 (a)) と同じ不安定性である可能性がある.

(a) 円柱結論
本研究では, 弧状衝撃波不安定性の力学的メカニズムを調べるために, DG 法を用いて鈍頭物体周りの衝撃波を含む流れ場を調べた.
弧状衝撃波の物理的不安定性を調べるために, カーバ
ンクル現象を生じにくいと考えられる数値流束, 計算格子を用いて, 不安定性の物体形状依存性を調べた. 曲率が異なる 2 つの楕円体形状物体を用いて計算を行った結果,
どちらの場合にも衝撃波面は定常になったが, 曲率が大きくなるにつれて亜音速領域の乱れは大きくなった. しかし, 衝撃波面の歪みは確認できなかった. 次に, 物体の角が不安定性に与える影響を調べるために, 円柱と円錐

4.

(b) 円錐Fig. 7: 角を含む物体の亜音速領域のマッハ数分布.

5Copyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
B12-5を用いて計算を行った. その結果, 円錐を用いて計算した
場合に衝撃波面に歪みが見られた. したがって, 衝撃波
面の歪みには円錐の角付近が影響する可能性がある. また, この歪みは, Baryshnikov らによって得られた弧状衝
撃波不安定性と同じ不安定性である可能性がある.
今後, 得られた不安定性が物理的なものか確かめるために, 更に計算格子依存性を調べ, スペクトル解析を通してメカニズムの解明を目指す.
参考文献(1) A. S. Baryshnikov, A. P. Bedin, V. G. Maslennikov
and G. I. Mishin, “Stability of a bow shock,” Sov.
Tech. Phys. Lett., 5 (1979), pp. 113–114.
(2) A. S Baryshnikov, “Baroclinic turbulence of dissociating gas as a reason for instability of bow shock
wave,” Phys. Scr., T132 (2008), 014007.
(3) J. Quirk, “A contribution to the great Riemannsolver debate,” International Journal for Numerical
Methods in Fluids, 18 (1994), pp. 555–574.
(4) B. Cockburn and C.-W. Shu, “The Runge-Kutta
discontinuous Galerkin method for conservation
laws V: Multimensional systems,” Journal of Computational Physics, 141 (1998), pp. 199–224.
(5) 鈴木佑一郎, 保江かな子, 大西直文, 第 42 回流体力学講演会, (2010), 2E10.
(6) Y. Sato, Y. Suzuki, K. Yasue and N. Ohnishi,“Three-dimensional simulation of bow-shock instability using discontinuous Galerkin method,”
Proceedings of 28th International Symposium on
Shock Waves, (2011), 2713.
(7) 佐藤陽介, 鈴木佑一郎, 保江かな子, 大西直文, 第 43
回流体力学講演会, (2011), 1E09.
(8) J. Sherwin, “Hierarchical hp finite elements in hybrid domains,” Finite Elements in Analysis andDesign, 27 (1997), pp. 109–119.
(9) E. Shima and K. Kitamura, “On New Simple
Low-Dissipation Scheme of AUSM-Family for All
Speeds,” AIAA Paper 2009-136 (2009).
(10) Y. Wada and M. S. Liou, “A Flux Splitting Schemewith High-Resolution and Robustness for Discontinuities,” AIAA Paper 94-0083 (1994).
(11) S. Gottlieb and C.-W. Shu, “Total Variation Diminishing Runge-Kutta Schemes,” Mathematics ofComputation, 67 (1998), pp. 73–85.

6

Copyright c 2011 by JSFM