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

新しい堅牢な Weighted Compact Nonlinear Scheme

A new and robust formulation of weighted compact nonlinear scheme○ 野々村拓, 宇宙機構宇宙研, 神奈川県相模原市中央区由野台 3-1-1, nonomura@flab.isas.jaxa.jp
藤井孝藏, 宇宙機構宇宙研, 神奈川県相模原市中央区由野台 3-1-1, fujii@flab.isas.jaxa.jp
Taku Nonomura, ISAS/JAXA, 3-1-1Yoshinodai, Sagamihara, Kanagawa, Japan
Kozo Fujii, ISAS/JAXA, 3-1-1Yoshinodai, Sagamihara, Kanagawa, Japan
In this study, a new weighted compact nonlinear scheme is presented. In WCNS procedure, linear difference scheme ismodified to be using flux on the computational nodes together with flux on mid points. This modification makes thescheme more robust, but more dissipative. The truncation error analysis is conducted and the reasons why thismodification makes scheme robust and dissipative are discussed. The standard shock problems are solved with original
and modified WCNS, and the results show that discontinuity thickness of modified WCNS becomes thicker andover/undershoots at the discontinuities are suppressed better. Also, the stiff shock tube problems, which the original
WCNS can not solve because of negative pressure, can be solved using modified WCNS without blow up. Series ofnumerical tests show the robustness of modified WCNS.
１．緒言
衝撃波を含む流れ場を解析する手法の一つとして Weighted
Compact Nonlinear Scheme(WCNS)1-3 がある． WCNS はコンパク
トスキーム 4 と WENO5 を組み合わせた差分法であり．オリジナ
ルの有限差分 WENO5 に比べ，1)解像度が高い 1,3，2)様々な流束評
価法が使える 1，3)曲線座標系で一様流保持ができる 6，といった利点があり，
これらの多くは WCNS が差分法でありながら変数補間が可能であることに起因する．(有限差分 WENO5 は流束を補間せねばならないことに注意されたい．)
この WCNS のアルゴリズムは，i)非線形補間，ii)近似リーマン解法による流束評価，iii)線形差分となっており，従来の圧縮性の流体解析で良く用いられてきた有限体積の MUSCL7 法と非常に似通っており，構造格子上で書かれた MUSCL 法からの拡張が容易である．
前述のように WCNS は，その WENO に対する様々な利点から
これまでに圧縮性乱流や空力音響問題に適用されてきた．しかしながら，これまでの我々の経験では計算が破綻しやすいという欠点が存在する．我々の経験上，アルゴリズム中で計算が破綻する
部分は，
1) 非線形補間をした際に圧力もしくは密度が負となって計算が破綻する．
2) 線形差分をした後，保存量を更新すると更新した保存量に対応する圧力もしくは密度が負になる．
という 2 点である．
1)に関しては，WCNS 以外の高次精度スキームでも観察される問題であり，これは，補間した圧力もしくは密度が負の場合，1次精度補間に落とすことで計算を安定にすることができる 7．一方 2)の問題は例えば，以下を初期条件とする 1 次元衝撃波問題でも生ずる．
⎧UU(x,t = 0 ) = ⎨ L⎩U R(x ≤ 0.5 )その部分だけ，2 次精度中心差分を使うなどの手法が考えられる．
この方法では線形差分をした後，時間更新をして圧力が負かどうかを判定しなければならず計算コストが大幅に増加する．特に陰解法との相性は悪いように思われる．
また WCNS は線形差分を用いることで保存量の総量を保存しているが，負圧になる部分を 2
次精度などの低次精度の差分に落とす非線形差分を行うと，保存量の総量が保存できなくなる．この問題は衝撃波の近傍で生ずるため不用意にこの処理を行うと衝撃波の捕獲の際，正しい解に収束しないなどの悪影響を及ぼす可能性がある．
以上より，2)の問題に対して，線形差分の枠内で改良を加えて計算が破綻しにくいスキームを構築することでこれを解決することが望まれる．本論文では，アルゴリズムの最後の部分の線形差
分の形を一部変形することでより堅牢なスキームが構成できることを示す．
２．計算手法
本節ではオリジナルの WCNS と改良した WCNS を説明する．
スキームの説明に以下の 1 次元 Euler 方程式を考える．
∂Q ∂E
(2)
+
=0
∂t ∂x
ここで，Q = ( ρ , ρ u , ρ e)T ，E = ( ρ u , ρ u 2 + p, ( ρ e + p )u )T である．
また状態方程式より

ρe =

1

γ −1

p+

1 2
ρu
2

(3)

である．
２．１ オリジナルの WCNS
まずはオリジナルの WCNS のアルゴリズムを説明する．(2r-1)
次精度の非線形補間および 2r 次精度の線形陽的差分を用いた(2r-1)次精度 WCNS を説明する．

(1)(x > 0.5 )２．１．１ 非線形補間
非線形補間では，Q Lj+1/ 2 ，Q Rj+1/ 2 を計算する．ここでは，Q Lj+1/ 2 を計算する手順を示す．まず， Q Lj+1/ 2 を計算するためのステンシル
Sj 内で保存変数から特性量に補間する．
ここで(2r-1)次精度の場合，

U L = ( ρ L , u L , pL ) = (1,0,10)
U R = ( ρ R , uR , pR ) = (1,0,0.1)
0 ≤ x ≤1

S j = { j − r , j − r + 1," , j + r}

この問題では，オリジナルの WCNS では 1 ステップ目のルンゲ
クッタ 1 段階が終了する時点で圧力が負となり，計算が破綻してしまう．
この問題に対して，負の圧力に対する同様の処理として，
高次精度差分を行った際に圧力が負になる場所を探し出してきて，(4)となる．保存変数から，特性変数の変換は以下のようになる．

q j +l:m = l j:mQ j +l1

(l = −3, −2," ,3)(5)
Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-1
ここで得られた数値流束は，真の流束の近似となっており，

ここで，lj:m は j 点に於ける m 番目の左固有ベクトルである．
次に，重み付き平均を行う．ここで重み付き平均を行うために定義されるk番目(=1...r)のサブステンシルSj:kは以下のようになる.
S j:k = { j − r + k , j − r + k + 1," , j + k − 1}

2 r −1E WCNS)j +1/ 2 = E j +1/ 2 + O ( Δxが成り立つ．前述の通り WCNS の利点として，流束評価に様々な手法を用いることができる点があり，流束評価には FDS 以外の
FVS(flux vector splitting) 系 9 ， AUSM(advection upwind splitting
method)系 10 や HLL(Harten-Lax-Leer)系 11 を使うことができる．一
般的な(2r-1)次精度有限差分 WENO4 の場合，数値流束 E WENO
は，
j +1/ 2(6)

重み付き平均を行う前に，それぞれのサブステンシル Sj:k から最大精度のU のj の位置でのn(=1...r-1)階微分の値U (j:nk) (k＝1,…,r)を以下の差分の形で求める．

∂E ′j
1
WENO
+ O ( Δx 2 r −1 )( EWENOj +1/ 2 − E j −1/ 2 ) =
Δx
∂x

n

⎛ 1 ⎞ r −1
q (j:nk) = ⎜ ⎟ ∑ akn:l q j − r + k +l
⎝ Δx ⎠ l =0(7)

上式での ank:l を補遺に示す．
それぞれのステンシルで評価された 1~(r-1)階微分を用いると，
k 番目のステンシルを用いたセルノードへの物理量の補間 Uj:k は
r −1⎛ Δx ⎞ ⎛ 1 ⎞q L 1 = q j + ∑ ⎜ ⎟ ⎜ ⎟ q (j:nk)
j + :kn =1 ⎝ 2 ⎠ ⎝ n ! ⎠
2
n

∑

k

k =1

２．１．３ 線形差分スキーム
最後に，線形差分スキームにより微分値を計算する．この線形
差分スキームでは，コンパクトスキームも使えるが，解像度があまり変わらないこと 12 が分かっている．そのため本論文では以下の陽的差分を用いた．

(8)

E j′ =

j+

1 1225 WCNS
( E j +1/ 2 − E WCNSj −1/ 2 )Δx 1024
1 245 WCNS
( E j +3/ 2 − EWCNS−
j −3/ 2 )Δx 3072
1 49
WCNS( E WCNS+
j + 5/ 2 − E j −5/ 2 )Δx 51201 5WCNS
7( E WCNS−
j + 7 / 2 − E j − 7 / 2 ) + O ( Δx )Δx 7168(19)となる．この式は，中点 (Mid point)からノード点(Node point)に差分を行う式となっている．そこで，これ以降 DM (Difference
scheme using Mid points)として記す．

(10)

j + :k
2

k =1

2

(18)

E j′ =

ここで用いられる(2r-1)次精度の Ck は補遺に示す．
しかしながら、実際にはいくつかのステンシルが滑らかでない場合がある．この時に振動を抑えるために最適な重み Ck でなくス
ムースインジケーターISj:k を使って滑らかなステンシルのみの結合となるような非線形な重み wk を計算し、これを用いて重み付き平均を求める．実際にはU Lj+1/ 2 は以下の式で表される．
r1 r −1
WCNS∑ ak ( EWCNS
j + k +1/ 2 − E j − k −1/ 2 )Δx k = 0

係数を Table 1 に示す．例えば，8 次精度差分は

1
j + :k :m2

q L 1 = ∑ wk q L 1(17)と定義されており， E WCNSの定義とは異なっていることに注意
j +1/ 2
されたい．

と表わされる．
すべてのステンシルで物理量の分布が滑らかな場合、それぞれのステンシルから作られた補間の線形結合として(2r-1)次精度の補間が構成できる．この(2r-1)次精度を達成する最適な重みを Ckとすると(2r-1)次精度のU Lj+1/ 2 は以下のように示される．
r(9)qL= C qL1
j + :m2(16)

ここで重み wk は
wk =

αk

(11)

Table 1 Coefficients of the cell-center to cell-node difference scheme.

(12)

coefficients
fourth order
explicit
sixth order
explicit
eighth order
explicit
tenth order
explicit

r

∑α
i =1

i

であり、αk は以下の式で表される．

αk =

Ck

( ISk + ε )p

ここで ISk は
r −1

IS j:k = ∑ ( q (j:nk) )

(13)

2

a0
9
8

a1
1
-24

a2

a3

a4

0

0

0

75
64

25
-384

3
640

0

0

1225
1024

245
-3072

49
5120

5
-7168

0

19845
16384

735
-8192

567
40960

405
-22937635294912n =1でεは 10-6、p は 2 とした．このように設計することで、すべて滑らかなステンシルであった場合 wk は Ck に近づき(2r-1)次精度の補間となることが出来る．
最後に，右固有ベクトル rj を用いて元の変数に変換される．
(14)
QL = r qL(l = −3, −2," , 3)
j +1/ 2∑j

２．２．よりロバストな WCNS
前節でのオリジナルの WCNS に対して，
よりロバストな WCNSを導入する．
よりロバストな WCNS は最後の線形差分の部分のみ異なるだけで他はオリジナルと全く同じである．
よりロバストな WCNS では線形差分を以下のように書く．

j +1/ 2m

ここまでが WCNS の非線形補間の手順である．

E j′ =

２．１．２ 流束評価法
次に，流束評価法に移る．本論文では，特に記述していない場
合 Roe の流束差分法(flux difference splitting：FDS)8 を用いて解析を行った．数値流束 E WCNS
は以下のように計算される．
j +1/ 2LR
E WCNSj +1/ 2 = f ROE (Q j +1/ 2 , Q j +1/ 2 )1 r∑ bk ( E j + k / 2 − E j −k / 2 )Δx k =1

WCNS⎪⎧ E
E j + k / 2 = ⎨ j + k / 2⎪⎩ E j + k / 2(20)( k = 2n − 1)
(k = 2n)

ここで n は整数である．このスキームの係数を Table 2 に示す．例
えば，8 次精度差分は(15)2Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-1
1 8 WCNS( E j +1/ 2 − EWCNS
j −1/ 2 )Δx 5
1 2
( E j +1 − E j −1 )
−
Δx 5
1 8
WCNS
( E WCNS+
j + 3/ 2 − E j − 3/ 2 )Δx 105
1 1( E j + 2 − E j − 2 ) + O( Δx 7 )−Δx 140

∂u ∂f
(22)
+
=0
∂t ∂x
f =u
となる．この式に対し，WCNS を適用する．線形移流方程式での WCNS は，Euler 方程式での実装において，前節での説明に対し Q=q=u，fj+1/2=uj+1/2L，E＝f，l=r=1 と置けばよい．
流れ場が十分に滑らかであることから非線形重みは最適重みになる．最適重みを用いた場合，7 次精度の流束が中点で得られ，この流束がもつ誤差は以下のようになる．
5(23)
f WCNS = f +f (7) Δx 7 + O(Δx8

E j′ =(21)となる．この式は，中点およびノード点(Mid and node point)からノード点(Node point)に差分を行う式となっている．そこで，これ以降 DMN(Difference scheme using Mid and Node points)として記す．
この式は，WCNS の数値流束が WENO のそれとは異なり，真の流束の近似となっているために成立する．本論文後半で示すよ
うに，この差分式はオリジナルの差分式よりロバストとなる． 5
次精度非線形補間とこの表式での 4 次精度差分の組み合わせは既に Point-Flux WENO (PFWENO)13 として発表されている．本論文
では，この手法を WCNS の変形と考えて議論していく．
また，これ以降，スキームは WCNS-IaD bM(N)と記述する．ここで，a は非線形補間の精度，b は差分の精度とする．それぞれのスキームの記述を Table 3 にまとめる．

j+

1
2

j+

1
2

2048

j

)

これを以下の 8 次精度の差分式
1 1225 WCNS
( f j +1/ 2 − f jWCNS−1/ 2 )
Δx 1024
1 245 WCNS
( f j +3/ 2 − f jWCNS−− 3/ 2 )Δx 3072
1 49
WCNS( f jWCNS++ 5 / 2 − f j − 5/ 2 )Δx 5120
1 5WCNS
7( f jWCNS−+ 7 / 2 − f j − 7 / 2 ) + O ( Δx )Δx 7168

f j′ =(24)に代入すると
5Δx 1225  (7)  (7)fj′ = f j′ +( f j − f j −1 )2048 1024
6
5Δx 245  (7)  (7)−( f j +1 − f j − 2 )2048 3072
6

Table 2 Coefficients of the cell-center to cell-node difference scheme.
coefficients
fourth order
explicit

b0
4
3

b1
1
-6

b2

b3

b4

0

0

0

sixth order explicit

3
2

3
-10

1
30

0

0

eighth order
explicit

8
5

2
-5

8
105

1
-140

0

tenth order explicit

5
3

10
-21

5
42

5
-252

1
630

5Δx 6 49  (7)  (7)( f j + 2 − f j −3 )2048 5120
5Δx 6 57 (7)( fj(7)−+ 3 − f j − 4 ) + O ( Δx )2048 7168+

5Δx  (8)(25)
fj′ = f j′ +f j + O( Δx 7 )2048が得られる．
式(22)の右辺第 2 項，4 項(下線部)は負の数値粘性が入っている
ことがわかる．本解析では，最適重みで考えたが，衝撃波前後でこの負の数値粘性が強く働くことが考えられ，スキームの安定性を考える上では好ましくない．これが原因でスキームが脆弱であると予想される．一方，式(23)から，負の数値粘性は正の数値粘性よりも弱く全体としては正の数値粘性となることがわかる．
一方，式(20)を，新しい WCNS の差分式
7Table 3 Schemes tested in this paper
Scheme
WCNS-I3-D4M
WCNS-I5-D6M
WCNS-I7-D8M
WCNS-I9-D10M
WCNS-I3-D4M
WCNS-I5-D6M
WCNS-I7-D8M
WCNS-I9-D10M

interpolation
3rd order
5th order
7th order
9th order
3rd order
5th order
7th order
9th order(22)

difference
Eq. 18 with 4th order coefs.
Eq. 18 with 6th order coefs.
Eq. 18 with 8th order coefs.
Eq. 18 with 10h order coefs.
Eq. 20 with 4th order coefs.
Eq. 20 with 6th order coefs.
Eq. 20 with 8th order coefs.
Eq. 20 with 10h order coefs.

( )1 8 WCNS( f j +1/ 2 − f jWCNS−1/ 2 )
Δx 5
1 2−
( f j +1 − f j −1 )
Δx 5
1 8
WCNS
+
( f jWCNS+ 3/ 2 − f j − 3/ 2 )Δx 105
1 1−( f j + 2 − f j − 2 ) + O ( Δx 7 )Δx 140f j′ =

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

(26)に代入すると

３．評価
３．１ 誤差解析
本節では WCNS-DM と改良した WCNS-DMN の誤差解析を行い，WCNS-DMN がよりロバストになる理由を簡単に説明する．
ここでは，7 次精度の非線形補間および 8 次精度の陽的線形差分を用いた全体で 7 次精度のスキームを考える．
線形移流方程式かつ流れ場が十分に滑らかな場合の誤差を解析する．
線形移流方程式は

5Δx 6 8  (7)  (7) 5Δx 6 8  (7)  (7)f j′ = f j′ +( f j − f j −1 ) +( f j +1 − f j − 2 ) + O (Δx 7 )2048 52048 105(27)

5Δx ⎛ 64  (8) ⎞fj′ = f j′ +f j ⎟ + O(Δx 7 )⎜2048 ⎝ 35
⎠

(28)

7

が得られる．
前述のオリジナルの WCNS での負の係数に相当する場所(式(24)，下線部)は，真の流束を用いており，このため数値粘性項が入っておらず，結果的に負の数値粘性とならない．これによりス
キームが安定になると考えられる．一方で，式(26)より数値粘性の合計はオリジナルの WCNS よりも数値粘性が大きいことがわ

3Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-1
かる．これによりスキームはより散逸的になると考えられる．

3.0

３．２ フーリエ解析
本節では WCNS-DM と改良した WCNS-DMN のフーリエ解析を行い，WCNS-DMN では数値粘性過多となり，解像度が落ちる
ことを議論する．
本章では高次精度 WCNS の解像度を解析的な評価を示す．
すべて滑らかな領域と仮定し、理想的な重みを用いたと仮定して解析を行った．解析の方法は以下のようである．

Re(ω'), Im(ω')

1.5

(29)

Uj=exp (i w xj)

をΔx=1 の格子上でスキームを用いて数値的に微分した時に、数値的な一階微分 Uj’は、

0.0

Uj’=i w* exp (i w xj)

(30)

0.0

と表される．この w*を w に対して解析的に計算した．理想的には w*=w となることが望ましい．結果を図 1 に示す．横軸が計
算前の波の波数で縦軸が微分スキームの評価した波数である．こ
こで、Ideal で表される直線に近いほうが、解像度が高いことを示している．厳密解や中心差分スキームでは w*の虚数部分は表れな
いが、風上差分となることによって w*の虚数部分も現れる．
Figure 1 に WCNS-DM と WCNS-DMN，同じ精度の WENO，低い精度の WCNS-DM の結果を示す．3，5，7 および 9 次精度において，WCNS-DMN のほうが WCNS-DM よりも解像度が低くなる
ことがわかった．また同じ精度の WENO と比較すると，ほぼ同
等か若干低い解像度であること，また低い精度の WCNS-DM より
は解像度が高いことが伺える．前述の誤差解析から得られたように WCNS-DMN がより散逸的であることがわかった．

3.0

WCNS-I5-D6M
WCNS-I5-D6MN
WCNS-I3-D4M
WENO5
Ideal

-1.5

1.5

3.0

ω
(b)5th order WCNS

3.0

WCNS-I3-D4M
WCNS-I3-D4MN
WENO3
Ideal

WCNS-I7-D8M
WCNS-I7-D8MN
WCNS-I5-D6M
WENO7
Ideal

Re(ω'), Im(ω')

1.5

Re(ω'), Im(ω')

1.5

0.0
0.0

1.5

3.0

0.0
0.0

1.5

3.0

-1.5

ω
(c)7th order WCNS
Fig. 1 Continued

-1.5

ω

(a)3rd order WCNS
Fig. 1 Comparison of resolution using modified wave number

4

Copyright © 2010 by JSFM

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

1.2
3.0

WCNS-I9-D10M
WCNS-I9-D10MN
WCNS-I7-D8M
WENO9
Ideal

WCNS-I7D8M
WCNS-I7D8MN
Exact

density

1
0.8
0.6
0.4

Re(ω'), Im(ω')

1.5

0.2
0
0

1

0.5
x

0.0
0.0

1.5

3.0

Fig . 2 Density distribution of Sod problem computed by 7th order
WCNS

0.45

-1.5

ω

WCNS-I7D8M
WCNS-I7D8MN
Exact

0.4

th

density(c)9 order WCNS
Fig. 1 Continued
３．３ 標準的な衝撃波問題
本 節 で は 標 準 的 な 衝 撃 波 問 題 に WCNS-DM お よ び
WCNS-DMN を適用して，衝撃波厚さ，アンダー/オーバーシュー
トの大きさを比較する．
まず，以下の初期条件を持つ Sod 問題 15 を解いた．
⎧U
U(x,t = 0 ) = ⎨ L⎩U R(x ≤ 0.5 )

0.35
0.3
0.25
0.6

0.65

(31)

0.7

0.75

0.8

x

(x > 0.5 )(a)Contact surface

U L = ( ρ L , u L , pL ) = (1,0,1)

WCNS-I7D8M
WCNS-I7D8MN
Exact

0.3

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

0.25
density

格子点数は 101 点とし，CFL 条件を 0.6 として t=0.2 まで解析した．WCNS-I7-D8M および WCNS-I7-D8MN の結果を Fig．2 に示す．また接触面および衝撃波面の拡大図を Fig. 3 にしめす．
WCNS-I7-D8MN のほうが衝撃波の厚さが厚いことがわかる．

0.2
0.15
0.1
0.75

0.8

0.85

0.9

0.95

x
(b)Shock wave
Fig. 3 Enlarged figure of Fig.2
より定量的に，不連続での特性を議論するため t=0.1 から t=0.2までの間において，不連続厚さ，アンダー/オーバーシュートを接触面および衝撃波において平均して比較する．ここで不連続厚さ
は，

5Copyright © 2010 by JSFM

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

δ=

Δρ

(32)

0.70%

と定義される．ここで， Δρ は厳密解の不連続での値の変化値である．
結果を Fig. 4 に示す．
まず不連続厚さは WCNS-DM に比べ，
WCNS-DMN のほうが厚くなっていることがわかる．これは，不連続に於いてより大きな数値粘性が加えられたためだと考えられる．次に，アンダー/オーバーシュー卜を見ると， WCNS-DMN
はオーバーシュートは大きいものの，アンダーシュー卜を抑えていることがわかる．計算が破綻しやすいのは強い衝撃波のアンダーシュートであることから，WCNS-DMN はよりロバストであることが期待される．

0.60%

)

Discontinuity thickness

5

0.30%

0.10%
0.00%
1

3

5
7
Order of Accuracy

911

次に以下の初期条件を持つ Lax 問題 16 を解いた．
⎧U
U(x,t = 0 ) = ⎨ L
⎩U R

2
1

(x ≤ 0.5 )(33)(x > 0.5 )U L = ( ρ L , u L , pL ) = (0.445,0.698,3.528)

1

3

5
7
Order of Accuracy

9

U R = ( ρ R , u R , pR ) = (0.5,0,0.571)

110 ≤ x ≤1(a)Discontinuity thickness

格子点数は 100 点とし，CFL 条件を 0.6 として t=0.13 まで解析した．WCNS-I7-D8M および WCNS-I7-D8MN の結果を Fig．5 に示す．また接触面および衝撃波面の拡大図を Fig. 6 にしめす．Sod
問題同様 WCNS-I7-D8MN のほうが衝撃波の厚さが厚いことがわ
かる．

0.80%
WCNS-DM, contact
WCNS-DMN, contact
WCNS-DM, shock
WCNS-DMN, shock

0.70%
0.60%

1.4

0.40%

1.2

0.30%

1
density

0.50%

0.20%
0.10%
0.00%

WCNS-I7D8M
WCNS-I7D8MN
Exact

0.8
0.6
0.4

1

Fig. 43(c)OvershootFig. 4 Continued

40Undershoot

0.40%

0.20%

WCNS-DM, contact
WCNS-DMN, contact
WCNS-DM, shock
WCNS-DMN, shock

6

WCNS-DM, contact
WCNS-DMN, contact
WCNS-DM, shock
WCNS-DMN, shock

0.50%
Overshoot(max ρ j +1 − ρ j

3

5
7
9
Order of Accuracy

11

0.2
0

(b)Undershoot
Discontinuty thickness, under/overshoot in Sod problem.

0

0.5

1

x
Fig . 5 Density distribution of Lax problem computed by 7th order
WCNSs.

6

Copyright © 2010 by JSFM

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

6Discontinuity thickness

1.3

density

1.1
0.9

WCNS-I7D8M
WCNS-I7D8MN
Exact

0.7
0.5
0.3
0.6

0.65

0.7

0.75

WCNS-DM, contact
WCNS-DMN, contact
WCNS-DM, shock
WCNS-DMN, shock

5
4
3
2
1
0

0.8

1

x

3

WCNS-I7D8M
WCNS-I7D8MN
Exact

11

WCNS-DM, contact
WCNS-DMN, contact
WCNS-DM, shock
WCNS-DMN, shock

0.035%
0.030%
Undershoot

1.1
density9(a) Discontinuity thickness(a)contact surface

1.3

5
7
Order of Accuracy

0.9
0.7

0.025%
0.020%
0.015%
0.010%

0.5

0.005%
0.3
0.75

0.8

0.85

0.9

0.000%

0.95

1

x

3

5
7
9
Order of Accuracy

11(b)Undershoot(b)Shock wave
Fig. 6 Enlarged view of Fig. 50.25%より定量的に，不連続での特性を議論するため t=0.065 から
t=0.13 までの間において，不連続厚さ，アンダー/オーバーシュー
トを接触面および衝撃波において平均して議論する．
結果を Fig. 7に示す．まず不連続厚さは Sod 問題と同様 WCNS-DM に比べ，
WCNS-DMN のほうが厚くなっていることがわかる．これは，不連続に於いてより大きな数値粘性が加えられたためだと考えられる．次に，アンダー/オーバーシュー卜を見ると，接触面のオーバ
ーシュートは高次精度の時により大きくなるが，それ以外でのアンダーオーバーシュートはほぼ同じか，小さくなることがわかる．
特に，本問題で生ずる強い衝撃波において，計算が破綻しやすい
アンダーシュートが抑えられており，これがロバスト性向上につながると考えられる．

WCNS-DM, contact
WCNS-DMN, contact
WCNS-DM, shock
WCNS-DMN, shock

Overshoot

0.20%
0.15%
0.10%
0.05%
0.00%
1

Fig. 773579Order of Accuracy

11(c)Overshoot
Discontinuty thickness, under/overshoot in Lax problem.

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-1
３．４ 高周波の解像度
本節では Shu-Osher 問題 5 に於いてスキームの高周波の波の解
Shu-Osher 問題は以下の初期条件で定義される．
像度を議論する．
(x ≤ 4 )(x > 4 )(34)

density⎧U
U(x,t = 0 ) = ⎨ L⎩U R

U L = ( ρ L , uL , pL ) = (3.857143, 2.629369,10.3333)
U R = ( ρ R , uR , pR ) = (1 + 0.2sin(0.5π x),0,1)
0 ≤ x ≤ 10

CFL が 0.6 となるように時間刻みを決定し，
t=1.8 まで計算した．
WCNS-I7-D8M および WCNS-I7-D8MN の結果を Fig．
8 に示す．
どちらも大きな数値振動を起こさずに計算ができていることがわ
か る ． Figure 9 に 拡 大 図 を 示 す ． 拡 大 図 か ら ， 新 し い
WCNS(WCNS-DMN)はオリジナルの WCNS よりも(WCNS-DM)よりも解像度が低くなることがわかる．これは，3.1 で説明したように，数値粘性がより大きく加えられるためである．

densitydensity

WCNS-I7D8M
WCNS-I7D8MN
Exact

0
10

5

7

8

7

8

7

8

(b)5th orderWCNS

3

0

6
x

4

1

WCNS-I5D6M
WCNS-I5D6MN
Exact

5

5

2

5
4.8
4.6
4.4
4.2
4
3.8
3.6
3.4
3.2
3

WCNS-I7D8M
WCNS-I7D8MN
Exact

5
4.8
4.6
4.4
4.2
4
3.8
3.6
3.4
3.2
3
5

6

x

x
th

5
4.8
4.6
4.4
4.2
4
3.8
3.6
3.4
3.2
3

(c)7th orderWCNS

WCNS-I3D4M
WCNS-I3D4MN
Exact
densitydensityFig. 8 Density distribution of Shu-Osher problem computed by 7 order
WCNS

5

6

7

8

WCNS-I9D10M
WCNS-I9D10MN
Exact

5
4.8
4.6
4.4
4.2
4
3.8
3.6
3.4
3.2
3
5

x

6
x

(a)3rd orderWCNS
Fig.. 9 Enlarged figures of density distribution of Shu-Osher problem
computed by WCNSs(d)9th order WCNS
Fig.. 9 Continued

8

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-1
３．５ より厳しい衝撃波問題
以下の初期条件のより厳しい衝撃波問題を解析した．
(x ≤ 0.5 )

10

(35)

(x > 0.5 )WCNS-I7D8M
WCNS-I7D8MN

8pressure⎧U
U(x,t = 0 ) = ⎨ L
⎩U R

12

U L = ( ρ L , uL , pL ) = (1,0,0.1PR)
U R = ( ρ R , uR , pR ) = (1,0,0.1)

64
2

0 ≤ x ≤ 100
ここで，圧力比 PR＝10, 100, 1000, 10000, 100000 のケースを扱った. CFL が 0.6 となるように時間刻みを決定した. 計算が破綻するケースに関しては，より小さな時間刻みでも計算したが，破綻する時刻はおおよそ変化しなかった．
様々な精度の WCNS-DM と WCNS-DMN がこの問題を解ける
るかどうかを Table 3 にまとめた．緒言でも示したようにオリジナ
ルの WCNS-DM では，PR=100 で計算が 1 ステップでも破綻してしまう．一方で，WCNS-DMN では全てのケースで計算が安定にできることがわかった．このようなシビアな条件では，新しいWCNS の表式を用いることでより安定に計算ができることがわ
かった．
Figure 10 に PR=100 の時の WCNS-I7-D8M の 1 ステップ目
Runge-Kutta の 1 ステージ後の圧力分布を示す．不連続で圧力が負になっており，これが計算の破綻につながる．一方
WCNS-I7-D8MN ではこのような負の圧力は観察されておらず，
Fig．11 のように最後まで，破綻することなく計算できた．

-20Fig. 10 Stiff shock tube problem with PR=100. Pressure distribution
at 1st stage of 1st step.

6

Exact
density4
32that problem is solved without blow-up and a cross denotes that the0

computation is blow-up.

WCNS-I3-D4M
WCNS-I5-D6M
WCNS-I7-D8M
WCNS-I9-D10M
WCNS-I3-D4M
WCNS-I5-D6M
WCNS-I7-D8M
WCNS-I9-D10M

10
○
○
○
○
○
○
○
○

100
×
×
×
×
○
○
○
○

WCNS-I7D8MN

5

1

Scheme0
10,000
×
×
×
×
○
○
○
○

1

x

Table 3 Computation of stiff shock tube problems. A circle denotes

PR
1,000
×
×
×
×
○
○
○
○

0.5

0.2

0.4

0.6

0.8

1

x

100,000
×
×
×
×
○
○
○
○

Fig. 11 Stiff shock tube problem with PR=100. Pressure distribution
solved by WCNS-DMN at t=0.1.
この問題では，初期の数ステップで計算が破綻してしまう．初期の数ステップでは，不連続の厚さが非常に薄いためであると考えられる．そこで，次のように，初期解から不連続がある程度の厚さを持った状態で解析を始めた．ここで，不連続での分布はKawai と Terashima の研究 17 を参考にした．
U(x,t = 0 ) = (1 − f smooth ) U L + f smoothU R
f smooth = 0.5 + 0.5* erf (

(36)

s
)
Cε Δx

ここでも同様に，PR＝10, 100, 1000, 10000, 100000 のケースを解析した．様々な精度の WCNS-DM と WCNS-DMN がこの問題をとけるかどうかを Table 4 にまとめた．オリジナルの WCNS-DM
は PR=100 まで解けるようになったが，PR=1000 以上の圧力比でこれらの問題では，は計算が 30 ステップ程度までの間に破綻した．
初期に不連続を滑らかにしていても，衝撃波が非常に強いため衝撃波付近で不連続が薄くなっていき，その後負圧が観測されるようである．当然ではあるが，シャープな初期条件と同様に
WCNS-DMN では破綻することなく計算が終了した．
Figure 12 に PR=100 の時の WCNS-I7-D8M の計算が破綻する直前の圧力分布を示す．衝撃波の前方で圧力が負になっており，これが計算の破綻につながる．一方 WCNS-I7-D8MN ではこのよう

9

Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
D10-1
な負の圧力は観察されておらず，Fig. 13 のように最後まで，破綻することなく計算できた．

7
6

WCNS-I7D8MN

5

Exact

Table 4 Computation of stiff shock tube problems with smootheddensity

initial condition. A circle denotes that problem is solved without blow-up
and a cross denotes that the computation is blow-up.
Scheme

10
○
○
○
○
○
○
○
○

WCNS-I3-D4M
WCNS-I5-D6M
WCNS-I7-D8M
WCNS-I9-D10M
WCNS-I3-D4M
WCNS-I5-D6M
WCNS-I7-D8M
WCNS-I9-D10M

PR
1,000
×
×
×
×
○
○
○
○

100
○
○
○
○
○
○
○
○

10,000
×
×
×
×
○
○
○
○

100,000
×
×
×
×
○
○
○
○

4
3
2
1
0
0

0.2

0.4

0.6

0.8

1

x
Fig. 13 Stiff shock tube problem with PR=1000 and smoothed initialcondition. Pressure distribution solved by WCNS-DMN at t=0.036.

120
100

WCNS-I7D8M
WCNS-I7D8MN

pressure80

３．６ ２次元衝撃波問題
以下の Woodward と Collela の問題 18 を解析した．
π
1
⎧

60

20
0

U L = ( ρ L , uL , vL pL ) = (8.0,1.516, −0.8751,116.5)

-20

(37)

U R = ( ρ R , u R , vR , pR ) = (1.4,0,0.1)

0

0.5

1
0 ≤ x ≤ 4, 0 ≤ y ≤1x

ここで y=0, x>1/6 において壁面となる．ここでは滑り有り壁を模擬するため，仮想格子点を用いて反射条件を課した．また，上側の境界条件では衝撃波の動きに合わせて厳密解を与えた．本問
題では、Roe の FDS では衝撃波不安定が生じたため流束評価に
SLAU(Simple Low-dissipative AUSM)19 を用いた．時間刻みは CFL
が 0.6 となるように決定した．
この問題では，緒言で議論した QL,, QR,を求めた際にも負圧が観測されて計算が破綻することがあったため，このような場合１次精度に落とすことにしている．
この計算を行うと，オリジナルの WCNS では１ステップ目で”更新した後の流体量が負圧となるため”に計算が破綻する．
これに対して，WCNS-DMN では計算が破綻することなく終了
した．結果を Fig. 14 に示す．

(a)Entire view

1

WCNS-I7D8M
WCNS-I7D8MN

0.8
pressure( x − y * tan( ) ≤ )36π1( x − y * tan( ) > )
3
6

⎪⎪U L
U(x,t = 0 ) = ⎨
⎪U
⎪⎩ R

40

0.6
0.4
0.2
0
-0.2
0

0.2

0.4

0.6

0.8

1

x
(b)Enlarged view for negative pressure computed by WCNS-I7-D8M.
Fig. 12 Stiff shock tube problem with PR=1000 and smoothed initialcondition. Pressure distribution before blow up of WCNS-DM.

Fig. 14 Double Mach reflection problem solved by
WCNS-I7-D8MN.

10Copyright © 2010 by JSFM

以上のように，WCNS-DMN とすることでより堅牢な高次精度スキームが構成できた．
k=1
k=2

４．議論
本論では，ノード点の流束も活用する WCNS-DMN を導入し，その利点であるロバスト性と欠点である散逸的になる効果を議論した．
この新しいスキームを導入する利点欠点をここで議論する．
WCNS は WENO に比べ，1)解像度が高い 1,3，2)様々な流束評価法が使える 1，3)曲線座標系で一様流保持ができる 4，という利点があったが，WCNS-DMN にすることによって，1)の解像度が高いという利点は無くなる．一方で，2)様々な流束評価法が使える
1，3)曲線座標系で一様流保持ができる 4，といった利点は残ったままであり，実用計算での使用する際には WENO に比べて非常に良い特性であるといえる．
また計算コストに関しては，ノード点での流束の計算は特に近似リーマン解法等を必要としないため，
オリジナルの WCNS と比べほぼ変わらない．
さらに，緒言では述べなかったが，オリジナルの WCNS では，使用するステンシルが非常に大きく，多くの点で境界スキームを用いなければならず，実装が手間であったが，新しい WCNS-DMN
はよりコンパクトなステンシルとなり，境界条件の実装の手間や，並列化した場合に重ねあわせる点数を減らすことができると考えられる．
一方で，本論で示したようにスキームの解像度は落ちるため，注意が必要となる．
５．結論
新しい差分式を導入し，より堅牢な WCNS の構築に成功した．
これに新しい WCNS は衝撃波のアンダーシュートを抑えており，
より堅牢になると考えられる． 特にオリジナルの WCNS では解けなかった Stiff な衝撃波問題でも安定して解けることがわかった．
一方でこの表式を用いるとスキームの解像度は落ちてしまう．
ただし，新しい(2r-1)次精度の WCNS は従来の(2r-3)次精度のWCNS よりも高い解像度であること，および新しい(2r-1)次精度のWCNS が(2r-1)次精度の WENO とほぼ同様の解像度であることを確認した．
最後に実問題に近い例として， 7 次精度で Double Mach
Reflection 問題を解いた．補間時の負圧の補正もいれることで，破綻することなしに解けることを確認した．
以上のように，本論文で導入した新しい WCNS はオリジナルと比べて堅牢なスキームとなった．新しい WCNS は WENO に比べた際の解像度が高いという利点は失ったものの，2)様々な流束評
価法が使える 1，3)曲線座標系で一様流保持ができるといった利点は変わらず持っており，曲線座標かつ強い衝撃波が存在する問題などでの適用が考えられる．

第 24 回数値流体力学シンポジウム
D10-1
Table A.1 a1k:l of third order (r=2)WCNS.
a1k:1
a1k:2
-1
1
-1
1

k=1
k=2
k=3

Table A.2 a1k:l of fifth order (r=3)WCNS.
a1k:1
a1k:2
1/2
-2
-1/2
0
-3/2
2

a1k:3
3/2
1/2
-1/2

k=1
k=2
k=3

Table A.3 a2k:l of fifth order (r=3)WCNS.
a2k:1
a2k:2
1
-2
1
-2
1
-2

a2k:3
1
1
1

Table A.4 a1k:l of seventh order (r=4)WCNS.
a3k:1
a3k:2
a3k:3

a3k:4

k=1

-1/3

3/2

-3

11/6

k=2

1/6

-1

1/2

1/3

k=3

-1/3

-1/2

1

-1/6

k=4

-11/6

3

-3/2

1/3

Table A.5 a2k:l of seventh order (r=4)WCNS.
a3k:1
a3k:2
a3k:3
k=1

-1

4

-5

2

k=2

0

1

-2

1

k=3

1

-2

1

0

k=4

2

-5

4

-1

Table A.6 a3k:l of seventh order (r=4)WCNS
a3k:1
a3k:2
a3k:3

A. 補遺
非線形補間における係数を示す．これらは Nonomura et al.6 が示したものであるが，参考文献内 6 に係数の間違いが合ったため，
ここで改めて示す．

11

a3k:4

k=1

-1

3

-3

1

k=2

-1

3

-3

1

k=3

-1

3

-3

1

k=4

-1

3

-3

1

k=1
k=2
k=3
k=4
k=5

謝辞
JAXA 石向桂一氏には研究に対するコメントを頂いた．ここに記して感謝する．

a3k:4

Table A.7
a3k:1
1/4
-1/12
1/12
-1/4
-25/12

a1k:l of ninth order (r=5)WCNS.
a3k:2
a3k:3
a3k:4
-4/3
3
-4
1/2
-3/2
5/6
-2/3
0
2/3
-5/6
3/2
-1/2
4
-3
4/3

a3k:3
25/12
1/4
-1/12
1/12
-1/4

Copyright © 2010 by JSFM

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

k=1
k=2
k=3
k=4
k=5

Table A.8 a2k:l of ninth order (r=5)WCNS.
a3k:1
a3k:2
a3k:3
a3k:4
11/12
-14/3
19/2
-26/3
-1/12
1/3
1/2
-5/3
-1/12
4/3
-5/2
4/3
11/12
-5/3
1/2
1/3
35/12
-26/3
19/2
-14/3

35/12
11/12
-1/12
-1/12
11/12

k=1
k=2
k=3
k=4
k=5

Table A.9 a3k:l of ninth order (r=5)WCNS.
a3k:1
a3k:2
a3k:3
a3k:4
3/2
-7
12
-9
1/2
-3
6
-5
-1/2
1
0
-1
-3/2
5
-6
3
-5/2
9
-12
7

a3k:3
5/2
3/2
1/2
-1/2
-3/2

k=1
k=2
k=3
k=4
k=5

Table A.10
a3k:1
1
1
1
1
1

a4k:l of ninth order (r=5)WCNS.
a3k:2
a3k:3
a3k:4
-4
6
-4
-4
6
-4
-4
6
-4
-4
6
-4
-4
6
-4

a3k:3
1
1
1
1
1

(1996), pp. 200-212(6) 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.
(7) Xu, Z.; Liu, Y. and Shu, C.-W., “Hierarchical Reconstruction for
Discontinuous Galerkin Methods on Unstructured Grids with a
WENO-type Linear Reconstruction and Partial Neighboring Cells,”
Journal, Vol.228, (2009), pp. 2194-2212.
(8) Roe, P. L., “Approximate Riemann Solvers, Parameter vectors, andDifference Schemes,” Journal of Computational Physics, 43, (1982),
pp357-372.
(9) van Leer, B., “Flux Vector Splitting for the Euler Equations,”
ICASE Report 82-293 (1982).
(10) Liou, M-S. and Steffen, C.J., “A New Flux Splitting Scheme,”
Journal of Computational Physics, 107, (1993), pp357.
(11) Harten, A., Lax, P.D. and van Leer, B., “On upstream differencing
and Godunov-type schemes for hyperbolic conservation laws,”
SIAM Review. 25 (1983), pp 35-61.
(12) Nonomura, T., and Fujii, K., “Effects of Difference Scheme Type in
High-order Weighted Compact Nonlinear Schemes,” Journal ofComputational Physics 228, (2009), pp. 3533-3539.
(13) Kudryavtsev, A. N. and Khotyanovsky, D.V., “Application of
WENO Schemes for Numerical Simulations of High Speed
Flows,” Book of Abstracts of International Conference on
Computational Fluid Dynamics, (2006), pp. 303.
(14) Gottlieb, S. and Shu, C.-W., “Total Variation Diminishing
Runge–Kutta Schemes, Mathematics of Computation 67 (1998), pp.
73–85(15) Sod, G.A., “A Survey of Several Finite Difference Methods for
Systems of Nonlinear Hyperbolic Conservation Laws,” Journal ofComputational Physics 27, (1978), pp. 1–31.
(16) Lax, P. D., “Weak Solutions of Nonlinear Hyperbolic Equations andTheir Numerical Computation. Communications on Pure andApplied Mathematics 7, (1954), 159–193.
(17) Kawai, S. and Terashima, H., “A High-resolution Scheme forCompressible Multicomponent Flows with Shock Waves,”
International Journal of Numerical Methods in Fluids, 2010.
(18) Woodward P. and Collela P., “The Numerical Simulation of Two
Dimensional Fluid with Strong Shock,” Journal of Computational
Physics, 54, (1984), pp.115-173.
(19) Shima, E. & Kitamura, K., “On AUSM-Family Scheme forAll
Speeds with Shock Detection for Carbuncle-Fix,” AIAA Paper
2009-3544.

a3k:3

Table A.11 Ck of third order (r=2)WCNS.
C1
C2
1/4
3/4

C1
1/16

Table A.12 Ck of fifth order (r=3)WCNS.
C2
C3
10/16
5/16

C1
1/64

Table A.13 Ck of seventh order (r=4)WCNS.
C2
C3
C4
21/64
35/64
7/64

C1
1/256

Table A.14 Ck of ninth order (r=5)WCNS.
C2
C3
C436/256
126/256
84/256

C5
9/256

参考文献
(1) Deng, X. G. and Zhang, H., “Developing High-order Weighted
Compact Nonlinear Schemes,” Journal of Computational Physics，
165, (2000), pp.22-44.
(2) Nonomura, T., Iizuka, N. and Fujii, K., “Increasing Order ofAccuracy of Weighted Compact Nonlinear Schemes,” AIAA Paper
2007-893, (2007).
(3) 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.
(4) Lele, S. K., “Compact Finite Difference Schemes with
Spectral-Like Resolution,” Journal of Computational Physics,
103(1992), pp. 16-42(5) Jiang, G.-S. and Shu, C.-W., “ Efficient Implementation of
Weighted ENO Schemes”, Journal of Computational Physics，126,

12

Copyright © 2010 by JSFM