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

Weighted Compact Nonlinear Scheme の多成分圧縮性流体への適用
Implementation of Weighted Compact Nonlinear Scheme toCompressible Multicomponent Flows○ 野々村拓, 宇宙機構宇宙研, 神奈川県相模原市中央区由野台 3-1-1, nonomura@flab.isas.jaxa.jp
寺島洋史, 宇宙機構情報・計算工学センタ, 茨城県つくば市千現 2-1-1, ishihara.hiroshi@jaxa.jp
藤井孝藏, 宇宙機構宇宙研, 神奈川県相模原市中央区由野台 3-1-1, fujii@flab.isas.jaxa.jp
Taku Nonomura, ISAS/JAXA, 3-1-1Yoshinodai, Sagamihara, Kanagawa, Japan
Hiroshi Terashima, JEDI/JAXA, 2-1-1 Sengen, Tsukuba, Ibaraki, Japan
Kozo Fujii, ISAS/JAXA, 3-1-1Yoshinodai, Sagamihara, Kanagawa, JapanA weighted compact nonlinear scheme (WCNS) was applied to simulations of compressible mulitcomponent flows,
and four different implementations (fully or quasi-conservative forms and conservative or primitive variables
interpolations) are examined in order to investigate numerical errors generated due to each implementation. The results
show that the different types of error in pressure are generated, when fully-conservative form or conservative variables
interpolation is selected, while quasi-conservative form generally shows poor mass conservation property. The error
due to fully-conservative form is generated only in initial flow fields, while the error due to conservative variables
interpolation exists during the computations, which leads to significant spurious oscillations near interfaces. The error
due to fully-conservative form can be reduced simply by smooth initial interface. The primitive variables interpolationis therefore considered to be the best choice for compressible mulitcomponent flows in the framework of WCNS,
either with fully or quasi-conservative forms. The choice of fully or quasi-conservative forms depends on problems
since there are advantages or disadvantages for each problem, such as good mass conservation property or spurious
oscillations.
よる．Abgrall17 はこの問題に対し，質量分率のみ(もしくはエネル
ギー式のみ 18)非保存系で解く準保存系を用いることで振動を消せることを示している．さらに Johnsen と Colonius19 が，高次精度へ拡張し，有限体積 WENO を用いた場合，準保存系でも保存変数を特性量補間すると振動が生ずることを示している．彼らは，有限体積 WENO において基本変数の特性量補間を行うことで，数値振動を抑えられることを示した．このような振動は有限差分
WENO でも生ずるようである 19．しかし，準保存系のため，質量保存性が悪いこともわかっている 22．Kawai と Terashima21 はコン
パクトスキームと局所人工粘性(LAD ：Localized Artificial
Diffusivity) を用い，完全保存系ではあるが初期界面を滑らかにして計算を始めることで，Abgrall17 が指摘した振動を最小限に抑えて解析している．彼らの手法は完全保存系であることから質量保存が完全に満足される．
Johnsen と Colonius19 らの手法は急激な界面を持っていてもWENO ベースで数値的な振動を抑えながら解析できることから幅広い条件での多成分気体の計算に適していると考えられるが，以下の欠点がある：

１．緒言
異なる気体の混合は様々な工学的問題で取り扱われ，衝撃波を含む高速流体場での予測も重要である．我々はロケット噴流からの音響波の解析を行っている 1-2 が，ロケット噴流と周りの気体は異なる気体種であり，これを考慮することで，どのように音響波の特性が変わるかを数値解析で評価したいと考えている．
そこで，
数値解析で評価する前段階として，本研究では，我々がこれまでに改良してきた高次精度不連続捕獲差分法である Weighted
Compact Nonlinear Scheme (WCNS)3-5 を圧縮性多成分気体へ適用した際の特性を明らかにすることを目的とする．
WCNSはWeighted Essentially Nonoscillatory Scheme (WENO)系6-7のスキームであり，オリジナルの有限差分 WENO に比べ，1)解像度が高い 3,5，2)様々な流束評価法が使える 1，3)曲線座標系で一様流保持ができる 8，といった利点があり，これらの多くは WCNSが差分法でありながら変数補間が可能であることに起因する．(有限差分 WENO は流束を補間しなければならない．)
これまでの多成分圧縮性流体の解法を大別すると，
Sharp-interface 法および Diffused-interface 法に分けられる．
Sharp-interface 法では界面で混ざらない流体を仮定し，Level-set9-11法や Front-tracking12-14 法などで界面を追跡して，各流体の方程式をそれぞれ解くものである．界面での境界条件の与え方に，
Ghost-fluid method13,15 やRiemann 問題を解く方法 14,16 などが挙げられる．これらの解法では界面をシャープに扱える一方で，界面のトポロジーが複雑になる場合には適用が難しい．
一方 Diffused-interface17-21 法は流体同士が混ざり合うことを許容し，質量分率などを導入し，界面を特に意識しない．我々の最終的な目標であるロケット噴流からの音響波の評価では流体同士が混ざり合い，複雑な界面になることが予想されるため，本研究では質量分率を導入した Diffused-interface 法を WCNS で実現する．
Abgrall17 は，質量分率を導入した Diffused-interface 法を比熱比の異なる物質の移送問題に適用し，完全保存系では数値振動が生じることを示した．この問題は，界面で異なる気体種が混合した際にエネルギー式と比熱比の評価によって圧力平衡が崩れることに

1) 有限体積法であるため，差分法に比べ 2-3 次元問題で 4-10 倍程度計算コストが高い．
2) 準保存系であるため，それぞれの気体種の質量の総和が保存されない．
1)に対して，本研究で導入する WCNS は差分法でありながら変数補間(補間に用いる変数は自由である)が可能であるため，基本変数を特性補間することで，Johnsen と Colonius19 が指摘した振動を抑えつつ，計算の低コスト化が期待できる．また 2)の欠点に関しては，化学反応などが伴う場合など好ましくない場合がある．
この欠点は保存系/準保存系の選択に起因するものであり，Abgrall17 が指摘した振動とトレードオフの関係にあるが，このトレードオフを含めた上での高次精度スキームに対する議論は十分になされていない．また，多成分気体に対する高次精度スキーム
1

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-2
また，準保存系の実装の場合，式(1)に付け加えて以下の非保存系の質量分率の式を解く．

の取り組みそのものが非常に少ない．
本研究では，WCNS を多成分気体に適用する際に，完全保存系/準保存の選択及び基本/保存変数補間の選択など，どのような実装を選択すべきか明らかにすることを目的とする． Kawai とTerashima21 が行ったように初期条件の影響も調べることで，それぞれの実装で現れる数値振動や誤差がどのように変化するかということにも言及する．

∂u∂Y1 ∂Y1ui+− Y1 i = 0
∂t
∂xi
∂xi

非保存系の質量分率の式は本来，∂Y1∂Y+ ui 1 = 0∂t∂xi

２．計算手法
本研究では，WCNS を前述のように完全保存系/準保存の 2 通り
および基本/保存変数補間の 2 通りの計 2×2=4 通りを実装した．
完全保存系および保存変数補間の実装を Fcnsv (Fully conservative
form with characteristic interpolation of conservative variables)，完全保存系および保存変数補間の実装を Fprim (Fully conservative form
with characteristic interpolation of primitive variables)，準保存系および基本変数補間の実装を Qcnsv (Quasi conservative form withcharacteristic interpolation of conservative variables)，準保存系および基本変数補間の実装を Qprim (Quasi conservative form withcharacteristic interpolation of primitive variables)として以後議論していく．ここで補間にはすべて特性補間(Characteristic interpolation)を用いた．界面の特性は成分毎補間(Componentwise interpolation)
でも同様と考えられるが，成分毎補間では衝撃波で数値振動が生ずるためである．これらの実装を Table 1 にまとめる．

Fully/Quasi

Interplolation

Fully

Conservative

Conservative

variables.

Fprim

Fully

Primitive

Conservative

variables.

Qcnsv

Quasi

Conservative

Conservative

variables.

Quasi

Primitive

Conservative

variables.

Qprim∂ρ ui ∂ ( ρ ui u j + pδ ij )
+
=0
∂t
∂x j

ρe =

∂xi

p

1
+ ρ ui ui(7)

γ −1 2Governing Eqs.

２．２．差分スキーム
本研究では，7 次精度の非線形補間，8 次精度陽的差分からな
る WCNS を用いた．流束評価法に HLLC 法 23 を用いた．Euler 方
程式に対する WCNS の手順を以下に示す．WCNS は非線形補間,,
流束評価,線形差分の 3 つの部分からなる．

Eqs (1), (2), (3), (6), (7)Eqs (1), (2), (3), (6), (7)Eqs (1), (3), (4), (6), (7)２．２．１非線形補間
ここでは，j 方向の Q Lj+1/ 2 を求める手順を示す．まず用いるス
テンシル内を中心となる格子点の左固有ベクトルを用いて，

Eqs (1), (3), (4), (6), (7)q j +l ,m = l j ,mQ j + l(l = −3, −2," , 3)(8)として，特性量に変換する．本研究では，保存変数補間と基本変数補間を扱うが，保存変数補間の場合 Q は保存変数で l は流束ヤ
コビアン A の左固有ベクトル，基本変数補間の場合 Q は基本変数
で l は基本変数 Q に対して基礎方程式を ∂Q / ∂t + M ∂Q / ∂x = 0と書いた場合の係数行列 M の左固有ベクトルとなる．
その後、m 番目の特性量 q Lj +, k1/ 2, m (k = 1,…,4)が k 番目の 4 点のサブステンシルから補間される．
5
21
35
35
q Lj+,11/ 2 = − q j −3 + q j − 2 − q j −1 + q j
16
16
16
16
1
5
15
5
q Lj+,21/ 2 = q j − 2 − q j −1 + q j + q j +1
16
16
16
16
1
9
9
1
(9)
q Lj+,31/ 2 = − q j −1 + q j + q j +1 + q j + 2
16
16
16
16
5
15
5
1
q Lj+,41/ 2 = q j + q j +1 − q j + 2 + q j +3
16
16
16
16
次に， q Lが q Lj +, k1/ 2, m の重み付き平均として計算される．

(1)

=0

ここで Y1 は気体種 1 の質量分率である．本研究では 2 成分気体を対象にするため，気体種 0 の質量分率 Y0 は
Y0 = 1 − Y1

Y0

Table 1 にそれぞれの実装がどの支配方程式を用いたかを併せて示す．

ここで，ρは密度，ui は速度ベクトル，p は圧力，e は質量あたりのエネルギー，δij はクロネッカーのデルタを示している．完全保存の実装の場合，式(1)に付け加えて以下の質量分率を保存系で解く．
(2)∂ρY1 ∂ρY1ui++

ここで，γ1は気体種1の比熱比，γ0は気体種0の比熱比である．この
比熱比を用いて，以下の理想気体の状態方程式より圧力とエネル
ギーの関係が決まる．

∂ρ e ∂ (( ρ e + p)u j )
+
=0
∂t
∂x j

∂t

Y1

(6)

+1

γ1 −1 γ 0 −1

２．１． 基礎方程式
本研究では，保存系の Euler 方程式および質量分率の式を基礎方程式とした．保存系の Euler 方程式は以下のようになる．
∂ρ ∂ρ ui
+
=0
∂t
∂xi

1

γ=

Conservative
Fcnsv

(5)

のように書けるが，Johnsen と Colonius19 にならって後述の流束評
価法において，移流速度を合わせるために式(4)の形とした．2種
気体から決まる比熱比は，界面での温度平衡を仮定すると，以下のように定義される．

Table 1 Summary of implementation of WCNS
Name

(4)

j +1/ 2, m

q Lj +1/ 2, m = w1, m q Lj +,11/ 2, m + w2, m q Lj +,21/ 2, m + w3, m q Lj +,31/ 2, m + w4, m q Lj +,41/ 2, m

(3)

(10)
ここで wk , m は非線形重みであり、

と書ける．
2

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-2
wk , m =α k ,m∑ α l ,m

(11)

Ck

(12)

1 1225 WCNS{
( E j +1/ 2 − E WCNSj −1/ 2 )Δx 1024
245 WCNS−( E j + 3/ 2 − E WCNS
j − 3/ 2 )
3072
49
WCNS
+
( E WCNS
j +5 / 2 − E j −5 / 2 )
5120
5
WCNS
−
( E WCNS
j + 7 / 2 − E j − 7 / 2 )}7168

E j′ =lα k ,m =

( IS

k ,m

+ε )

2

として求まる．ここで(C1 , C2 , C3 , C4 ) = (1/ 64, 21/ 64,35 / 64, 7 / 64)(12)は最適重みであり，ISk ,m は以下のように定義されるスムースイン
ジケータである．
4⎛ 4⎞
IS k , m = ∑ ⎜ ∑ cn , k ,l q j + k + l −5,m ⎟n =1 ⎝ l =1⎠

２．３．時間積分スキーム
時間積分には 3 次精度 Total-Variation-Diminishing Runge-Kutta(Strong Stability Preserving Runge-Kutta)26 を用いて解析を行った．

2

(13)

Deng と Zhang3，Nonomura ら 4 は式(9)および(13)を計算する際に，それぞれのステンシルで計算される n 階微分の計算値を用いて計算を効率化している．
最後に，右固有ベクトル rj を用いて元の変数に変換される．
(14)
QL = r qL(l = −3, −2," , 3)
j +1/ 2∑j

２．４．準保存系および基本変数補間による圧力平衡の維持
最後に，圧力，速度一定，密度および質量分率が変化する場合にQprim の実装を用いることでJohnsen とColonius19 と同じように
完全に圧力平衡が保持されることを示す．
1 次元で速度 u>0 一定，圧力 p 一定の流れ場を考える．基本変数補間を用いると，uL=u，pL=p が成り立つ．(保存変数補間ではこれが満たされない．) HLLC での流速評価を考えて WCNS 差分
および 1 次精度 Euler 陽解法を仮定すると，密度，運動量，エネ
ルギーの変化は以下のように示される．

j +1/ 2m

ここまでが WCNS の非線形補間の手順である．
２．２．２ 流束評価法
次に，流束評価法に移る．本研究では，Johnsen と Colonius にならって，HLLC23 を用いた．ここでは簡単のため 1 次元 HLLCを説明する．
1 + sign( s* ) L − *L
⎡ E + s ( Q − Q L )⎤⎣⎦2*
1 − sign( s ) R + *R
⎡ E + s (Q − Q R )⎤+⎣⎦2

EWCNS =

(15)(

( ρu ) jn +1

u nj +1 =である．ここでs − = min ( 0, s L )s = max ( 0, s ),u − c )(s = max ( ( u + c ) , u + c )
R

24

ROE

)

ρ nj

=u

(21)

L

L

ROE

R

R

2
⎧
⎡( ρ u )n +1 ⎤ ⎫
j ⎦ ⎪
n +1
n +1
n +1
⎪
n +1
p nj +1 = (γ n +1 − 1) ⎨( ρ e ) j − ⎣⎬ = ( γ − 1) ⎡⎣( ρ e ) j − u ( ρ u ) j ⎤⎦ = p
2 ρ nj +1 ⎪
⎪
⎩
⎭

(17)

(22)
が得られる．よって，これらの式から，Johnsen と Colonius の有限体積法と同様に，圧力平衡および速度平衡が保てる高次精度差分法が構築できると期待される．

である．ここで，c は音速であり，ROE の上添字は Roe 平均値を示す．また接触面の速度は以下のように定義される．

s* =

(20)

n+1 ステップでも速度を維持することがわかる．この結果を使い，圧力に関しても若干の計算を行うと，

R

s L = min ( u − c ))ここで bk，s は差分スキームによって決まる変数である．ここでは簡単のため陽的スキームを用いた．この式より，n+1 ステップでの速度が求まり，

t

(16)

+

)

(

sL/ R − uL/ R
s L / R − s*
⎛ L/ R L/ R * L/ R
pL/ R ⎞ ⎞*L/ R ⎛ L/ R *⎜ ρ , ρ s , E + ( s − u )⎜ ρ s + L/ R⎟⎟
s − uL/ R ⎠ ⎠
⎝
⎝

(

Δx
s −1
u ∑ k =0 bk ( ρiL+1/ 2+ k − ρiL−1/ 2− k )
Δt
s −1
Δx
n +1
n
( ρ u ) j = ( ρ u ) j − u 2 ∑ k =0 bk ( ρiL+1/ 2+ k − ρiL−1/ 2−k )Δt⎡ u 2 s −1⎤
LL⎢ ∑ k =0 bk ( ρi +1/ 2+ k − ρi −1/ 2− k )⎥Δx ⎢ 2n +1n⎥( ρe) j = ( ρe) j − u ⎢
⎛
⎞
⎥
s −1
1
1
Δt
− L
⎟⎟ ⎥
⎢ + p ∑ k =0 bk ⎜⎜ L⎢⎣
⎝ γ j +1/ 2 + k − 1 γ j −1/ 2− k − 1 ⎠ ⎥⎦
s
Δ
x
−
1
Y jn +1 = Y jn +u ∑ k =0 bk (Yi +L1/ 2+ k − Yi −L1/ 2 − k )Δtρ nj +1 = ρ nj +ここで，EL,ER はセル境界での左側,右側の物理量から計算された流束であり, QL,QR は左側,右側の保存量, q*はq* L / R =(19)

p R − p L + ρ Lu L ( s L − u L ) − ρ Ru R ( s R − u R )ρ L ( sL − uL ) − ρ R ( sR − uR )(18)
３．１次元テスト問題
３．１ 精度検証
まず4つの実装が全て7次精度を保っていることを確認する．不連
続がある場合，スキームが 1 次精度に落ちるため，精度検証は滑らかな流れ場に対して行った．初期条件として，一定速度を持ち，非常に滑らかに気体種の割合が変わる問題を解析した．初期条件を以下のように設定した．

２．２．３ 線形差分スキーム
最後に，線形差分スキームにより微分値を計算する．この線形
差分スキームでは，コンパクトスキームも使えるが，解像度があまり変わらないことが分かっている．そのため並列化，ベクトル化のため以下の陽的 8 次精度差分 25 を用いた．
3

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-2
U = ( ρ , u , p, Y1 ) = (5.5 + 0.9sin(2π x ),0.5,1/1.4,0.5 + 0.1sin(2π x))

γ 0 = 1.4, γ 1 = 1.6

L2 error(23)−0.5 ≤ x ≤ 0.5

density

初期条件を Fig. 1 に示す．境界条件は周期境界である．時間刻み
Δt はΔt＝(Δx)3 として t＝2.0 まで計算した．この問題において t=2.0 の解析解は初期条件と同じであるため，初期条件からのずれをみることで誤差を調べた．格子点数は 8,16,32,64 点とした．Figure 2 に誤差を示す．全てのスキームで 5.5～6.5 次精度が達成された．7 次精度にならないのは臨界点で最適重みにならないことなどが原因 27 である．これらの問題を解決する Mapped technique27 を用いて，行った同様のテストの誤差は Fig. 3 に示すようにおおよそ 7 次精度となった．
これらの結果より 4 つ全ての実装で精度が適切に出せることがわかった．以降のテストでは，通常の重みの評価方法を用いている．

Fcnsv
Fprim
Qcnsv
Qprim
7th order
0.01

0.1

1

x
Fig. 3 Error of WCNS implementation with Mapping technique for the
problem examining order of accuracy

３．２ 物質移送問題
次に，Johnsen と Colonius19 や Kawai と Terashima21 らと同様に 1 次
元物質移送問題を解くことで 4 つ実装の特徴を明らかにする．初期条件は以下のように設定した．

7
6
5
4
3
2
1
0

⎧U
U(x,t = 0 ) = ⎨ material⎩U ambient(0.25 ≤ x ≤ 0.75)
otherwise(24)

U material = ( ρ m , um , pm , Y1m ) = (10,0.5,1/1.4,1.0)
-0.5

-0.3

-0.1

0.1

0.3

U ambient = ( ρ a , ua , pa , Y1a ) = (1,0.5,1/1.4,1.0)0.5x

γ 0 = 1.4, γ 1 = 1.6
0 ≤ x ≤1(a)density

ただし，誤差を低減させるために用いる滑らかな初期条件を作る
際には，Kawai と Terashima21 に習い，Error Function を用いて，以下のように定義した．

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

U(x,t = 0 ) = (1 − f smooth ) U ambient + f smoothU material

Velocity
Pressure
Y1
Y0
-0.5

-0.3

-0.1

0.1

0.3

f smooth = 0.5 + 0.5* erf ((b)Velocity, pressure, mass fraction
Fig. 1 Initial condition of the problem examining order of accuracy
1.00E+00
1.00E-01
1.00E-02
1.00E-03
1.00E-04
1.00E-05
1.00E-06
1.00E-07
1.00E-08
1.00E-09
1.00E-10

Fcnsv
Fprim
Qcnsv
Qprim
7th order
0.01

0.1

s
)
Cε Δx

(25)

ここで，s は符号付きの界面からの距離であり，Δx は格子幅，Cεは滑らかさを変化させるパラメタである．Cεを変化させた場合の密度分布を Fig. 4 に示す．
境界条件は周期境界である．時間刻みは CFL 条件が 0.6 となるよ
うに選び，t＝2.0 まで計算した．この問題において解析解は圧力および速度が一定値を保つ．格子点数は 100 点とした．
まず，シャープな界面を初期条件にもつ場合の t=2.0 における密度，圧力，気体 1 の質量分率の分布を Fig. 5 に示す．Qprim は完全に圧力平衡を保っており，Qcnsv が小さな数値振動，Fcnsv, Fprim が大きな数値振動を示している．シャープな界面を初期条件にした場
合は，完全保存系にすることによる誤差が大きいことを示している．
次に，Fig. 6 に気体種 1 の質量の総量，圧力，速度の L2 誤差の時間
履歴を示す．当然ながら Fcnsv，Fprim は気体種を完全に保持している．一方で，圧力，速度の誤差はどの時刻でも，Fig．5 の考察同様，
Qprim が完全に圧力平衡(速度平衡)を保ち，Qcnsv が小さな数値振動，Fcnsv, Fprim が大きな数値振動をもっている．すなわちシャープな界面を初期条件にした場合は，完全保存系にすることによる誤差
がどの時刻でも大きい．

0.5

xL2 error

1.00E+00
1.00E-01
1.00E-02
1.00E-03
1.00E-04
1.00E-05
1.00E-06
1.00E-07
1.00E-08
1.00E-09
1.00E-10

1xFig. 2 Error of WCNS implementation for the problem examining orderof accuracy
4

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-2
12

density108
64|(M1-M1|t=0)/M1|t=0|

sharp
1
2
3
4
5

2
0
-0.5

-0.3

-0.1

0.1

0.3

1.00E-01
1.00E-03
1.00E-05

Fcnsv

1.00E-07

Fprim

1.00E-09

Qcnsv

1.00E-11

Qprim

1.00E-13
1.00E-15
1.00E-17

0.5

0.00

x

Fcnsv

10

Fprim

1.00E-01

Qcnsv

1.00E-03

L2 norm of velocity

density

12

Qprim

6

Exact

4
2
0
-0.5

-0.3

-0.1

0.1

0.3

1.00E-05

Fprim

1.00E-09

Qcnsv

1.00E-11

Qprim

1.00E-13
1.00E-17
0.00

Qcnsv
Qprim
Exact

0.1

0.3

0.5

1.00E-05

Fcnsv

1.00E-07

Fprim

1.00E-09

Qcnsv

1.00E-11

Qprim

1.00E-13
1.00E-15

Fprim

1

Qcnsv

0.8

Qprim

0.6
0.2
0
-0.2
0.3

1.00

1.50

2.00

次に，滑らかな界面を初期条件に持つ場合の t=2.0 における圧力場を Fig. 7 に示す．ここでの結果は Cε =3 の結果を示す．Qprim は同様に完全に圧力平衡を保っているが，Fprim も同様に圧力が変動している様子は見られない．一方で Fcnsv と Qcnsv は振動が残っていることがわかる．すなわち，滑らかな初期条件で始めた場合，保存変数で補間することの誤差が相対的に大きくなることがわかる．Figure 8に気体種 1 の質量の総量，圧力，速度の L2 誤差の時間履歴を示す．
Fcnsv，Fprim は気体種を完全に保持している．一方で，圧力，速度の誤差はどの時刻でも，Fig．7 の考察同様，Qprim が完全に圧力平
衡(速度平衡)を保ち，Fprim が小さな数値振動，Fcnsv, Qcnsv が大きな数値振動をもっている．すなわち滑らかな界面を初期条件にした場合は，保存変数を補間することによる誤差が相対的に目立つよう
になる．

Exact

0.4

0.1

0.50

t
(c)L2 error in pressure field
Fig. 6 Time history of error in material convection problem with sharp
initial interface (t=2.0)

Fcnsv

1.2

mass fraction

2.00

1.00E-03

0.00(b)pressure

-0.1

1.50

1.00E-17

x

-0.3

1.00

1.00E-01

Fprim

L2 norm of pressurepressure

Fcnsv

-0.5

0.50

t
(b)L2 error in velocity field

0.755
0.75
0.745
0.74
0.735
0.73
0.725
0.72
0.715
0.71
0.705
-0.1

2.00

Fcnsv

1.00E-07(a)density

-0.3

1.50

1.00E-15

0.5

x

-0.5

1.00

t
(a) Mass conservation error of gas 1.

Fig. 4 Effects of Cε on the initial profile of material convection problem

8

0.50

0.5

x
(c) Mass fraction of gas 1.
Fig. 5 Results of material convection problem with sharp initialinterface (t=2.0)

5

Copyright © 2010 by JSFM

1.00E+01

12

Fcnsv

10

Fprim

L2 norm of velocity

density

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

Qcnsv

8

Qprim

6

Exact

4
2
0
-0.5

-0.3

-0.1

0.1

0.3

L2 norm of pressurepressure

Fprim
Qcnsv
Qprim
Exact

-0.3

-0.1

0.1

0.3

Qprim

1.00E-13
1.00E-15
0.50

1.00

t

1.50

2.00

Qprim

0.6

Exact

0.4
0.2
0
0.1

0.3

0.5

x
(c) Mass fraction of gas 1.
Fig. 7 Results of material convection problem with smoothed initialinterface (t=2.0)
1.00E+01
1.00E-01
1.00E-03
1.00E-05

Fcnsv

1.00E-07
1.00E-09

Fprim
Qcnsv
Qprim

1.00E-13
1.00E-15
1.00E-17
0.50

1.00

t

1.00E-07
1.00E-09

Fprim

1.00E-11

Qprim

1.50

Qcnsv

1.00E-13
1.00E-15
0.50

1.00

t

1.50

2.00

最後に，Fig. 9 に示す，シャープな界面から界面を滑らかにしていった場合の t=2.0 での誤差の変化をみる． 基本的に全ての誤差は界面を滑らかにすることで減ることがわかる．また，界面の滑らかさに関わらず，完全保存系の Fcnsv と Fprim は質量の総和の誤差がマシ
ンゼロ程度，Qprim は速度，圧力の誤差がマシンゼロとなっている．
Figure 9 の Qcnsv と Qprim の質量の総和に対する誤差は，界面を滑らかにすることで 2 桁程度軽減できることがわかる．一方で速度，圧力に関する誤差は，Fprim は界面を滑らかにすることでその誤差を 2桁程度減らせるのに対し Fcnsv および Qcnsv はあまり誤差が減らな
いことがわかる．これは以下のように考察できる．Fprim は完全保存
系にすることによる誤差のみを持つと考えられ，この誤差は界面を滑らかにすることで急激に減る．一方 Qcnsv は保存変数補間をすることによる誤差のみを持つと考えられるが，この誤差は界面を滑らかにしてもあまり減らない．Fcnsv は両方の誤差を持つが，Fig. 9 より，Fcnsv
が Qcnsv もしくは Fprim の大きい方とほぼ同じ程度の誤差を持つことから，これらの誤差はそれぞれ独立なものとして考えて良いといえる．
以上の知見をまとめる．
z
JohnsenとColoniusと同様に準保存系および基本変数補間を用いることで，高次精度差分法のフレームワークの中で圧力・速度平衡を保つスキームが開発できた．
z
質量の総和，圧力平衡，速度平衡に対する誤差は，初期界面を滑らかにすることによって緩和できることがわかった．
z
完全保存系にすることによる誤差はシャープな界面で顕著であるため，初期条件および境界条件で滑らかな分布を与えることでかなり緩和できる．
z
保存変数補間による誤差は界面を滑らかにしてもあまり減らないため，基本変数補間が望ましいと考えられる．

Qcnsv

0.8

1.00E-11

Fcnsv(c)L2 error in pressure field (Qcnsv is identical to Fcnsv)
Fig. 8 Continued

Fprim

1

-0.1

1.00E-05

1.00E-17

0.5

Fcnsv

-0.3

1.00E-01
1.00E-03

0.00

1.2

mass fraction

Qcnsv

1.00E-11

1.00E+01
Fcnsv(b)pressure|(M1-M1|t=0)/M1|t=0|

Fprim(b)L2 error in velocity field (Qcnsv is identical to Fcnsv)

x

0.00

Fcnsv

1.00E-07
1.00E-09

0.00(a)density

-0.2
-0.5

1.00E-05

1.00E-17

0.5

x

0.755
0.75
0.745
0.74
0.735
0.73
0.725
0.72
0.715
0.71
0.705
-0.5

1.00E-01
1.00E-03

2.00(a) Mass conservation error of gas 1.
Fig. 8 Time history of error in material convection problem with sharp
initial interface (t=2.0)
6

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-2
波および界面近くの圧力分布を拡大して観察すると，Fcnsv は Qcnsvに比べて，接触面，膨張波が始まる部分で大きな振動をしている．同様に Fprim は Qprim に比べて膨張波が始まる部分で大きな振動をしている．Fcnsv でみられた接触面の振動は Fprim ではみられておらず，これは保存変数補間による誤差だと考えられる．一方，Qcnsv とQprim を比べると，Qcnsv は接触面での誤差が見られるが，Qprim では逆に衝撃波では弱い振動をしていることがわかる．すなわちこのテストで，基本変数補間は界面捕獲には良いが，衝撃波捕獲には悪い影響を及ぼすことがわかった．この問題でも，滑らかな界面を初期条件として解析を始めたが，この場合，完全保存系であるために生ずる膨張波が始まる場所での誤差が大きく改善され，その他の誤差は変わらなかった．

1.00E+01|(M1-M1|t=0)/M1|t=0|

1.00E-01
1.00E-03
1.00E-05

Fcnsv

1.00E-07

Fprim

1.00E-09

Qcnsv

1.00E-11
1.00E-13

Qprim

1.00E-15
1.00E-17
0.00

1.00

2.00

3.00

Cε

4.00

5.00

(a) Mass conservation error of gas 1.
1.2

1.00E-01

1

1.00E-03

0.8

Fcnsv

0.6

Fprim

1.00E-05
1.00E-07
1.00E-09
1.00E-11
1.00E-13

density

L2 norm of velocity

1.00E+01

Fcnsv
Fprim

Qcnsv

0.4

Qcnsv

0.2

Qprim

0

1.00E-15

Qprim
Exact
0

0.2

0.4

0.00

0.6

0.8

1

x

1.00E-17
1.00

2.00

3.00

Cε

4.00

5.00

Fig. 10 Density distribution of multi-component Sod problem(t=0.2)(b)L2 error in velocity field
0.315

1.00E+01
1.00E-03
1.00E-05
1.00E-07
1.00E-09
1.00E-11
1.00E-13

pressure

L2 norm of pressure

1.00E-01

Fcnsv
Fprim

Fcnsv
Qcnsv

Qcnsv

Exact

Qprim

0.307

1.00E-15

0

0.2

0.4

1.00E-17
0.00

1.00

2.00

Cε

3.00

4.00

0.8

1

0.8

1

5.00

(a)Fcnsv and Qcnsv(c)L2 error in pressure field
Fig. 9 Effects of smoothness of interface on errors in material convection
problem.

pressure

0.315

３．３ 衝撃波問題
以下の 2 成分 Sod 問題 18,19.21 を解析した．
⎧U
U(x,t = 0 ) = ⎨ L
⎩U R

0.6

x

(x ≤ 0.5 )

Fprim
Qprim
Exact

(26)

0.307

(x > 0.5 )

0

0.2

0.4

0.6

x

U L = ( ρ L , u L , pL , Y1L ) = (1,0,1,0)(b)Fprim and Qprim
Fig. 11 Enlarged figure of Fig. 10

U R = ( ρ R , u R , pR , Y1R ) = (0.125,0,1,1)

γ 0 = 1.4, γ 1 = 1.6

次にかの研究と同様 21,28 な以下のより厳しい条件の改良 Sod 問題を解析した．

0 ≤ x ≤1

本問題で，シャープな界面を初期条件とし，格子点数を 201 点，CFL
条件は 0.6 とした．ｔ＝0.2 の密度分布を Fig. 10 に示す．図より，数値
振動はあまり顕著でなくほぼ同じ程度である．Figure 11 に示す衝撃⎧U
U(x,t = 0 ) = ⎨ L⎩U R

7(x ≤ 0.5 )(27)(x > 0.5 )Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-2
した，2 つの誤差を切り分けるのに適しているためである．まず全体の概要を示すため，Qprim で 5001 点を用いた解析の密度の x-t ダイ
アグラムを Fig 13 に示す．

U L = ( ρ L , uL , pL , Y1L ) = (10,0,10,0)
U R = ( ρ R , u R , pR , Y1R ) = (0.125,0,1,1)

γ 0 = 1.4, γ 1 = 1.6
0 ≤ x ≤1

本問題では，シャープな界面を初期条件とした場合計算が破綻した
ため，Cε =1 の滑らかな初期条件で解析した，格子点数を 201 点，
CFL 条件は 0.6 とした．ｔ=0.16 の密度分布を Fig. 12 に示す．本問題の場合，実装による差は，解析解からのずれから比べて十分に小さくなったため，実装によらずおおよそ同程度の誤差を持った解が得られると考えられる．シャープな界面を初期条件とした場合計算が破綻する問題に関しては，別途，本シンポジウムで発表するより堅牢なWCNS を用いることで解決できるが，解像度の若干の低下が同時に生ずるため本論文では議論しない．

t12

density

10
8

Fcnsv

6

Fprim

x

Qcnsv

4

Qprim

2

Fig. 13 x-t diagram of density field of the one-dimensional
shock-bubble interaction problem．

Exact

0
0

0.2

0.4

0.6

0.8

1

x

Figure 13 より，この問題では，衝撃波にヘリウムバブルが入っていき，そこで膨張波，衝撃波が発生する．その後，界面付近から弱い圧
縮波の発生が見られるが，これらはバブル界面に於ける衝撃波の反射によるものである．
この問題をシャープな界面及び 1001 点，CFL 条件を 0.6 として解析した．本問題は 201 点の格子点でも解析しているが，誤差の切り分
けが難しかったため，1001 点を選んだ．t=4.0 における結果を Fig. 14に示す．
Fprim と Qprim を比較した Fig．14(b)より，完全保存系にすることによる誤差は図内丸で囲んだ部分でみられる．これらは，シャープな初期条件で急に界面が動き始めたときに発生したものであり，音速でこの位置まで伝播してきている．一方，Qcnsv と Qprim を比較した
Fig．14(c)より，保存変数補間による誤差は，同じく丸で囲んだ部分に存在することがわかる．この部分では，高周波の振動が存在することがわかる．この誤差は，初期だけではなく，シミュレーションを通して生成しており，その特性は好ましくない．Fcnsv は Fprim と Qcnsv のもつ誤差が両方現れている(Fig．14(a))．最後に，ヘリウムの質量の総和の時間変化を Fig. 15 に示した．準保存系では，計算終了時の質量の総和の誤差は 2.5％程度見られた．

Fig. 12 Density distribution of multi-component severe Sod
problem(t=0.16)以上の知見をまとめる．
z
2成分Riemann問題では実装の差はあまりあらわれない．特に，条件の厳しい改良 Sod 問題では大きな差はなかった．一方でSod 問題では以下の若干の誤差の出方の変化が見られた．
z
基本変数補間を用いることで，界面での若干の誤差は減らせるが，衝撃波での誤差が生ずる．
z
完全保存系を用いると，膨張波が始まる場所で誤差が現れる．
ただし滑らかな界面の初期条件ではこの誤差は観察できな
い．
３．４ バブル衝撃波干渉問題
次に，以下の初期条件をもつバブル衝撃波干渉問題の結果を示
す．
⎧ U bubble⎪
U(x,t = 0 ) = ⎨ U preshock⎪U⎩ postshock( x + 4) 2 < 0.5
x < −3otherwise(28)

U bubble = ( ρ B , uB , pB , Y1B ) = (0.1819,1.22,1/1.4,1)

U preshock = ( ρ pre , u pre , p pre , Ypre ) = (1,1.22,1/1.4,0)
U postshock = ( ρ post , u post , p post , Ypost ) = (1.3764,0.8864,1.5698 /1.4,0)

γ 0 = 1.4, γ 1 = 1.648−5 ≤ x ≤ 5

本問題は，一般的に行われているテスト 18,21 とは観測する座標系
を変えて空気中での衝撃波が静止するようにした．これは，3.2 で示
8Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-2|(M1-M1|t=0)/M1|t=0|

Error due to fully conservative form

1.2

pressure

1.1
1Error due tointerpolation ofconservative
variables

0.9
0.8
0.7

Fcnsv
Qprim

1.00E-03
1.00E-05

Fcnsv

1.00E-07

Fprim

1.00E-09

Qcnsv

1.00E-11

-5

-3

-1

1

3

1.00E-15

50.00(a) Fcnsv and Qprim

Error due to fully conservative form
1.2

1.0
0.9
Fprim

0.7

Qprim

0.6
-5

-3

-1

1

3

5

x

pressure

1
0.9
0.8

1.1

0.7

1

0.6

Fcnsv
Qprim
-5

0.9
0.8

-3

-1

0.7

Qprim
1

3

3

5

(a) Fcnsv and Qprim

0.6
-1

1

x

Qcnsv

-3

4.00

1.1

Error due to interpolation ofconservative variables

-5

3.00

1.2

(b) Fprim and Qprim

1.2

2.00

次に，Cε =3 で滑らかな界面を初期条件にして解析を行った．
Figure 16 に結果を示す．シャープな界面でみられた Fprim の持っていた誤差は非常に小さくなり，分布からは判断ができない．一方で
Qcnsv が持っていた誤差はほぼ同程度で残っている．この誤差はシ
ャープな界面を初期条件にもつ場合とほぼ同じで，計算初期だけでなくシミュレーションを通して生成しており，その特性は好ましくない
ことがわかる．Fcnsv でみられる誤差は，Qcnsv でみられる誤差と同じである．最後に，ヘリウムの質量の総和の時間変化を Fig. 17 に示した．準保存系では，質量の総和の誤差は 0.5％程度見られた．

1.1

0.8

1.00

t
Fig. 15 Time history of error in mass conservation of the
one-dimensional shock-bubble interaction problem with sharp initialinterface．

x

pressure

Qprim

1.00E-13
1.00E-17

0.6

pressure

1.00E-01

1.2

5

1.1

pressurex(c) Qcnsv and Qprim
Fig. 14 Pressure distribution of the one-dimensional shock-bubble
interaction problem with sharp initial interface (t=4.0)．

1
0.9
0.8

Fprim

0.7

Qprim

0.6
-5

-3

-1

1

3

5

x
(b) Fprim and Qprim
Fig. 16 Pressure distribution of the one-dimensional shock-bubble
interaction problem with smooth initial interface (t=4.0)．

9

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-2⎧ U bubble⎪
U(x, y , t = 0 ) = ⎨ U preshock⎪U⎩ postshock

Error due to interpolation ofconservative variables

1.2

pressure

1.1

x<2otherwise(29)

U bubble = ( ρ B , uB , vB , pB , Y1B ) = (0.1819,1.22,0,1/1.4,1)

1

U preshock = ( ρ pre , u pre , v pre , p pre , Ypre ) = (1,1.22,0,1/1.4,0)

0.9
0.8

Qcnsv

0.7

U postshock = ( ρ post , u post , v post , p post , Ypost )Qprim= (1.3764,0.8864,0,1.5698 /1.4,0)

0.6
-5

-3

-1

1

3

5

γ 0 = 1.4, γ 1 = 1.648x

0 ≤ x ≤ 16 , −0.89 ≤ y ≤ 0.89(c)Qcnsv and Qprim
Fig. 16 Continued．

初期条件の模式図を Fig. 18 に示す

1.00E+01|(M1-M1|t=0)/M1|t=0|( x − 1) 2 + y 2 < 0.5

Stationary Shock Wave

1.00E-01
1.00E-03
1.00E-05

Fcnsv

1.00E-07
1.00E-09

Fprim

1.00E-11

Qprim

Qcnsv

1.00E-13
1.00E-15

UpreshockFig. 18 Schematic of initial condition of the two-dimensional
shock-bubble interaction problem.

1.00E-17
0.00

Upostshock

1.00

2.00

t

3.00

4.00

本研究では，誤差の切り分けを容易にするために，一般的なバブル衝撃波干渉問題 19,21 と異なる座標系で解析した．格子点数は1601×179 であり，バブル直径に 100 点が入るようにした．上下の境界は滑り壁面とし，左右の境界条件は初期条件と同じである．
シャープな界面と Cε=3 の滑らかな界面の初期条件に対して
CFL0.6 として解析した．
Qprim でシャープな界面の初期条件で解析した結果をFig. 19 に示す．衝撃波の中にバブルが入って行き，その後，バブルの形状が変形していく様子がわかる．全般的なバブルの形状変化は実験結果や他の数値解析結果とよく一致している．
次に，この問題で生ずる誤差を議論する．t=1.086 のときの速度発散を Fig. 19 に示す．速度発散を用いることで生じた数値的な誤差が観察できる．Fprim と Qprim を比較することで完全保存系にすることによる誤差を議論する．Fprim と Qprim の差は Fprim のみに見られる外側に存在する音波である．これは，計算の初期に発生したものである．それ以外は特に大きな差異は見当たらず，
完全保存系にすることによって生ずる誤差は計算初期以外ではほとんど見えない．さらにこの誤差は滑らかな界面の初期条件にすると，大幅に低減され，Fig. 19(b)の滑らかな界面の初期条件の結果の可視化では確認できなくなる．次に，Qcnsv と Qprim を比較することで保存変数補間を行うことによる誤差を議論する．
Qcnsv と Qprim の差はバブル内部の振動である．これは，計算の初期のみでなく，計算を通して発生している．これらの誤差は滑らかな初期条件で解析した場合も存在しているため，保存変数補間よりも基本変数補間が良いと考えられる．最後に，Fcnsv の誤差は，Fprim と Qprim に存在する誤差を足し合わせたものとなっており，これらの誤差はほぼ無関係であるといえる．

Fig. 17 Time history of error in mass conservation of the
one-dimensional shock-bubble interaction problem with smooth initialinterface．
バブル衝撃波干渉問題では以下の知見が得られた．
z
完全保存系にすることによる誤差はシャープな界面から発生し，計算を進めると誤差の発生はほとんど無くなる．また界面を滑らかにしておけば，この誤差はほとんど見られなくなる．
z
保存量補間による誤差は界面の滑らかさに関わらずどの時点でも通じて発生しており，望ましくない．
ここまでの 1 次元テストより，保存変数補間よりは基本変数補間が望ましいと考えられる．また完全保存系にすることによる誤差は初期のシャープな界面から強く発生するのみであることがわかった．これらの誤差は準保存系および基本変数補間を用いた実装で，完全に消すことができる．これは Johnsen と Colonius の有限体積 WENO と同様の結果であるが，WCNS は差分法であるため，多次元問題で高効率になることが期待できる．ただし準保存系にすることによる気体
種の質量の総和の誤差(今回の問題では最大 2.5%程度)が許容できないシミュレーションを行う際は，滑らかな初期条件を与えて，完全保存系および基本変数補間で実装することで，気体種の質量の総和を完全に維持しつつ，大きな誤差が発生しない状態でシミュレーショ
ンできることがわかった．
４．2 次元テスト問題
４．１．バブル衝撃波干渉問題
次に，2 次元バブル衝撃波干渉問題を解析した．初期条件は以下の式で表わされる．

10

Copyright © 2010 by JSFM

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

Fcnsv
Error due toconservative form
Fprim

Error due to Qcnsv
interpolation ofconservative variables
Qprim(a)Sharp initial interface

Fcnsv

Error due to Fprim
interpolation ofconservative variables
Qcnsv

Fig. 19 Example of computational results of the two-dimensional
shock-bubble interaction problem. Density distribution. The result is
computed by Qprim. Top to bottom; t=0.410, 0.622, 0.821, 1.086, 2.002,
3.242. 4.880 and 6.929. The figures show the regions at 0.9<x<11.7.

Qprim(a)Smooth initial interface
Fig. 20 Divergence of velocity fields of the two-dimensional
shock-bubble interaction problem. Top to bottom: Fcnsv, Fprim, Qcnsv,
Qprim．The results at t=1.086 are shown. Contour range normalized byinitial bubble diameter is set from -1 to 1. The figures show the regions at
0.9<x<11.7.
次に，バブルの形状に対する影響をみる．t=6.929 における密度分布を Fig. 20 に示す．図より，Fcnsv と Qcnsv は界面に振動を持っており，界面の動力学を議論する際に悪影響となることが予想される．
この傾向は滑らかな界面の初期条件でも変わらないため，
保存変数補間はこのような問題に対して悪影響を及ぼすことがわ
かる．一方で，Fprim の界面は，Qprim のそれとほぼ同じである．
Fprim の誤差は界面がシャープであるときのみに初期にでる音波だけと考えてよく，全体として良い精度を保っていると考えられる．
11

Copyright © 2010 by JSFM

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

Fprim

|(M1-M1|t=0)/M1|t=0|

Fcnsv

Error due tointerpolation ofconservative variables

1.00E-01
1.00E-03
1.00E-05

Fcnsv

1.00E-07

Fprim

1.00E-09

Qcnsv

1.00E-11

Qprim

1.00E-13
1.00E-15
1.00E-17
0.00

Qcnsv

2.00

6.00

8.00

|(M1-M1|t=0)/M1|t=0|t(a)Sharp initial interface

Qprim(a)Sharp initial interface

Fcnsv

1.00E-01
1.00E-03
1.00E-05

Fcnsv

1.00E-07

Fprim

1.00E-09

Qcnsv

1.00E-11

Qprim

1.00E-13
1.00E-15
1.00E-17
0.00

Fprim

4.00

2.00

4.00

6.00

8.00

t
(a)Smooth initial interface
Fig. 22 Time history of error in mass conservation of thetwo-dimensional shock-bubble interaction problem.

Error due tointerpolation ofconservative variables

４．２．Rhcitmayer-Meshkoff 不安定
次に Rhcitmayer-Meshkoff 不安定問題を解析した．初期条件は以下の式で表わされる．

Qcnsv⎧ U heavy⎪U(x, y , t = 0 ) = ⎨ U preshock⎪U⎩ postshock

Qprim

x < 0.4 − 0.1sin(2π ( y + 0.25))
x < 0.7
otherwise(30)

U heavy = ( ρ h , uh , vh , ph , Y1h ) = (5.04,1.24,0,1/1.4,1)
U preshock = ( ρ pre , u pre , v pre , p pre , Ypre ) = (1,1.24,0,1/1.4,0)
U postshock = ( ρ post , u post , v post , p post , Ypost )(a)Smooth initial interface
Fig. 21 Density fields of the two-dimensional shock-bubble interaction
problem. Top to bottom: Fcnsv, Fprim, Qcnsv, Qprim．The results at
t=6.929 are shown. Contour range is set from 0.2 to 1.4. The figures show
the regions at 0.9<x<11.7.

= (1.4112,0.8787,0,1.6272 /1.4,0)

γ 0 = 1.4, γ 1 = 1.093
0 ≤ x ≤ 16 , 0 ≤ y ≤ 1

最後に Fig. 22 に示す質量保存を議論する． 2 次元問題での計算終了時の Qcnsv, Qprim の準保存系の質量保存の誤差はシャープな界面(Fig. 22(a))で 9％の質量の総和の誤差，滑らかな界面(Fig.
22(b))で 6％の質量の総和の誤差となっている．Fcnsv，Fprim の完全保存系では誤差はマシンゼロ程度である．準保存系での質量保存の誤差が許容範囲外の場合，質量保存の誤差がなく初期界面を滑らかにした状態で Fprim を用いることが勧められる．一方これらが許容範囲内であれば，Qprim は圧力平衡を保つため，この実装がより良いと考えられる．

初期条件の模式図を Fig. 23 に示す．
Stationary Shock Wave
Heavy fluid

Light fluid
Upreshock

Upostshock

Fig.23 Schematic of initial condition of the two-dimensional RM
instability problem.
12

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-2
本研究では，誤差の切り分けを容易にするために，バブル衝撃
波干渉問題と同様，一般的な RM 不安定性問題 11,13,21 と異なる座標系で解析した．格子点数は 2049×128 であり，1 波長に 128 点が入るようにした．上下の境界は周期境界を課し，左右の境界条
件は初期条件と同じである．シャープな界面と Cε=3 の滑らかな界面の初期条件に対して CFL0.6 として解析した．
Qprim でシャープな界面の初期条件で解析した結果をFig. 24 に示す．RM 不安定性により，スパイクが生成していく様子が見られる．この形状変化は過去の計算と良い一致を示している．

Fcnsv
Error due toconservative form
Fprim

Qcnsv

Error due to
Qprim
interpolation ofconservative variables(a) Sharp initial interface
Fig. 19 Example of computational results of the two-dimensional
shock-bubble interaction problem. Density distribution. The result is
computed by Qprim. Top to bottom; t=0.275 at 2＜x＜6, t=0.550 at 4＜x＜8, t=0.825 at 6＜x＜10 and t= 11.0at 8＜x＜12．

Fcnsv

この問題でも生ずる誤差を議論する．t=0.275 のときの速度発散
を Fig. 25 に示す．速度発散を用いることで生じた数値的な誤差が観察できる．Fprim と Qprim を比較することで完全保存系にすることによる誤差をみる．Fprim と Qprim の差は Fprim のみに見られる衝撃波前方にある音波である．これは，計算の初期に発生したものである．それ以外は特に大きな差異は見当たらず，完全保存系にすることによって生ずる誤差は初期以外ではほとんど見えない．さらにこの誤差は滑らかな界面の初期条件にすると，大幅に低減され，Fig. 25(b)の可視化では確認できなくなる．
次に，
Qcnsv と Qprim を比較することで保存変数補間を行うことによる誤差をみる．Qcnsv と Qprim の差はスパイク内部で見られる振動である．これは，計算の初期のみでなく，計算を通して発生している．
これらの誤差は滑らかな初期条件で解析した場合も存在しているため，
保存変数補間よりも基本変数補間が良いと考えられる．最後に，Fcnsv の誤差は，Fprim と Qprim に存在する誤差を足し合わせたものとなっており，
これらの誤差はおおよそ無関係であるといえる．

Error due toFpriminterpolation ofconservative variables
Qcnsv

Qprim

(b) Smooth initial interface
Fig. 25 Divergence of velocity fields of the two-dimensional RM
instability problem. Top to bottom: Fcnsv, Fprim, Qcnsv, Qprim．The
results at t=0.275 are shown. Contour range normalized by wave length isset from -0.5 to 0.5. The figures show the regions at 2<x<6.
最後に，界面の形に対する影響をみる．t=0.825 における密度分
布を Fig. 26 に示す．図より，Fcnsv と Qcnsv は界面に振動を持っており，界面の動力学を議論する際に悪影響となることが予想される．この傾向は滑らかな界面の初期条件でも変わらないため，
保存変数補間はこのような問題に対して悪影響を及ぼすことがわ
かる．一方で，Fprim の界面は，Qprim のそれとほぼ同じである．
Fprim の誤差は界面がシャープであるときのみに初期にでる音波だけと考えてよく，全体としては良い精度を保っていると考えられる．

13

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-2間の 2×2 通りを試した．
まず，従来の結果と同様に，完全保存系による誤差と保存変数補間を行うことによる誤差が界面で発生することを確かめた．
完全保存系による誤差は，シャープな計算初期界面を用いた場
合にのみ発生し，滑らかな初期界面からはあまり発生しないことがわかった．一方，保存変数補間にすることによる誤差は界面の滑らかさに関わらず発生し，計算の精度に悪影響を及ぼすことがわかった．
Johnsen と Colonius の有限体積 WENO と同様に，WCNS でも準保存系および基本変数補間を行うことでこれらの誤差は完全に消し去ることができた．WCNS は有限差分であるため，有限体積
WENO よりも多次元問題で 4-10 倍程度効率が良いスキームとなっている．ただし，この方法では，Johnsen と Colonius の WENOと同様に準保存系で気体種の質量保存が悪く，今回のテストでは
気体種の質量の総和の誤差が最大 10％程度であった．燃焼問題のように質量の総和が重要となる問題では，滑らかな初期界面を導入し，完全保存系かつ基本変数補間による WCNS の実装も考慮にいれるべきであろう．
最後に，本論文では詳細に議論しないが，Yang の人工圧縮 29
を，音波以外の特性波に用いた場合では，基本変数および完全保存系でも計算を通して強い誤差が界面から生じた．これは Yangの方法を用いることで，界面が常に急峻になり，本論文中での完全保存形の計算において初期で生じた誤差が出続けるためである．
このようなテクニックを併せて使う場合は，準保存および基本変数補間が必須になると考えられる．

Fcnsv

Fprim

Qcnsv

Error due tointerpolation ofconservative
variables

Qprim(a) Sharp initial interface

Fcnsv

Fprim

Qcnsv

謝辞
JAXA 石向桂一博士には研究に対するコメントを頂いた．東京大学大学院小嶋亮次氏には本論文で用いた計算コード LANS3D への HLLC ソルバーの導入を手伝っていただいた．また計算の一部
は JAXA, JSS システムを用いた．ここに記して感謝する．

Error due tointerpolation ofconservative
variables

参考文献(1) Fujii, K., Nonomura, T. and Tsutsumi, S., “Toward Accurate
Simulation and Analysis of Strong Acoustic Waves,” International
Journal of Numerical Methods in Fluids. (to appear)(2) Nonomura, T. and Fujii, K., “Recent Efforts for Rocket Plume
Acoustics,” Computational Fluid Dynamics Review 2010, (2010)
pp. 421-446, World Scientific Publishing Company.
(3) Deng, X. G. and Zhang, H., “Developing high-order weighted
compact nonlinear schemes,” Journal of Computational Physics，
165, (2000), pp.22-44.
(4) Nonomura, T., Iizuka, N. and Fujii, K., “Increasing Order ofAccuracy of Weighted Compact Nonlinear Schemes,” AIAA Paper
2007-893, (2007).
(5) Zhang, S., Jiang, S. and Shu, C.-W., “Development of Nonlinear
Weighted Compact Schemes with Increasingly Higher Order
Accuracy” Journal of Computational Physics, 227, (2008),
pp.7294-7321.
(6) Liu, X., Osher, S. and Chan, T., “Weighted EssentiallyNon-oscillatory Schemes,” Journal of Computational Physics, 115,
(1994), pp. 200-212.
(7) Jiang, G.-S. and Shu, C.-W., “ Efficient Implementation of
Weighted ENO Schemes”, Journal of Computational Physics，126,(1996), pp. 200-212(8) Nonomura, T., Iizuka, N. and Fujii, K., “Freestream and Vortex
Preservation Properties of High-order WENO and WCNS on
Curvilinear Grids,” Computers & Fluids, 39, (2010), pp. 197-214.

Qprim(b) Smooth initial interface
Fig. 26 Density fields of the two-dimensional RM instability problem.
Top to bottom: Fcnsv, Fprim, Qcnsv, Qprim．The results at t=0.825 areshown. Contour range is set from 1.5 to 8.5. The figures show the regionsat 2<x<6.
ここで得られた結論はバブル衝撃波干渉問題とほぼ同じであった．2 次元問題で得られた知見をまとめると，z
完全保存系にすることによる大きな誤差は計算初期に発生する音波のみである．
z
保存変数補間による誤差は計算を通して発生しており，界面の形状にも悪影響を及ぼすため，基本変数補間が望
ましい．
５．結論
WCNS を多成分気体問題に対して適用した．実装方法として，
完全保存系もしくは準保存系，保存変数補間もしくは基本変数補14

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-2
Runge-Kutta Schemes,” Mathematics of Computation, Vol.67,(1998), pp. 73-85.
(27) Henrick, A. K., Aslam, T. D. and Powers, J. M.,“Mapped Weighted
Essentially Non-oscillatory Schemes: Achieving Optimal Order
Near Critical Points,” Journal of Computational Physics, Vol.207,(2005), pp. 542-567.
(28) Wackers, J., and Koren, B. “A Fully Conservative Model for
Compressible Two-fluid Flow, International Journal of Numer.ical
Methods in Fluids 47 (2005) pp. 1337–1343.
(29) Yang, H., “An Artificial Compression Method for EN0 Schemes:
The Slope Modification Method,” Journal of Computational
Physics, Vol.89, (1990), pp. 125-160.

(9) Osher, S. and Serthian, J.A., “Fronts Propagating with
Curvature-dependent Speed: Algorithms based on Hamilton-Jacobi
Formulations,” Journal of Computational Physics, 79, (1988), pp.
12-49.
(10) Karni, S., “Multicomponent Flow Calculation a Consistent
Primitive Algorithm,” Journal of Computational Physics, 112,
(1994) pp. 31-43.
(11) Nourgaliev, R. R., Dinh, T. N., and Theofanous, T. G. “Adaptive
Characteristics-based Matching for Compressible Multifluid
Dynamics,” Journal of Computational Physics, 213, (2006), pp.
500-529.
(12) Tryggvason, G., Bunner, B., Esmaeeli, A., Juric, D., Al-Rawahi, N.
Tauber, W., Han. J., Nas, S., and Jan, Y.-J., “A Front tracking
methods for the computations of multiphase flow,” Journal ofComputational Physics, 169, (2001), pp. 708-759.
(13) Terashima, H. and Tryggvason, G., “A Front-tracking/Ghost-fluid
Method for Fluid Interfaces in Compressible Flows,” Journal ofComputational Physics, 2009, 228, 4012-4037.
(14) Terashima, H. and Tryggvason, G., “A Front-tracking Method with
Projected Interface Conditions for Compressible Multi-fluid
Flows,” Computers and Fluids, 39, (2010), 1804-1814.
(15) Fedkiw, R. P., Aslam, T., Merriman, B. and Osher, S. “A
Non-oscillatory Eulerian Approach to Interfaces in Multimaterial
Flows (the Ghost Fluid Method), Journal of Computational Physics,
152, (1999), pp. 457-494.
(16) Hu, X. Y. and Khoo, B. C., “An Interface Interaction Method forCompressible Multifluids,” Journal of Computational Physics, 198,
(2004), pp. 35-64.
(17) Abgrall, R. M., “How to Prevent Pressure Oscillations in
Multicomponent Flow Calculations: A Quasi Conservative
Approach,” Journal of Computational Physics, 125, (1996), pp.
150-160.
(18) Abgrall, R. and Karni, S., “Computations of Compressible
Multifluids,” Journal of Computational Physics, 169, (2001), pp.
594-623.
(19) Johnsen, E. and Colonius, T., “Implementation of WENO Schemesin Compressible Multicomponent Flow Problems” Journal ofComputational Physics, 219, (2006), pp. 715-732.
(20) Chang, C. H. and Liou, M. S. “A Robust and Accurate Approach toComputing Compressible Multiphase Flow: Stratified Flow Model
and AUSM+up Scheme,” Journal of Computational Physics, 225,
(2007), pp. 840-873.
(21) Kawai, S. and Terashima, H., “A High-resolution Scheme forCompressible Multicomponent Flows with Shock Waves,”
International Journal of Numerical Methods in Fluids, 2010.
(22) Johnsen, E. “Spurious Oscillations and Conservation Errors inInterface-Capturing Schemes, ”Annual Research Briefs, (2008), pp.
115-126.
(23) Toro, E., F.; Spruce, M. and Speares, W., “Restoration of the Contact
Surface in the HLL-Riemann Solver,” Shock Waves, Vol.4, (1994),
pp. 25-34.
(24) Roe, P. L., “Approximate Riemann Solvers, Parameter Vectors, andDifference Scheme,” Journal of Computational Physics, Vol.43,(1981), pp. 357-372,(25) Nonomura, T. and Fujii, K., “Effects of Difference Scheme Type in
High-order Weighted Compact Nonlinear Schemes,” Journal ofComputational Physics, Vol.228, (2009), pp. 3533-3539.
(26) Gottlieb, S. and Shu, C.-W., “Total Variation Diminishing
15

Copyright © 2010 by JSFM