第 28 回数値流体力学シンポジウム
講演番号 C07-2

多成分圧縮性流体解析で速度・圧力・温度平衡を保つ差分 WENO 法
Finite Difference WENO Scheme Maintaining Velocity, Pressure and Temperature Equilibrium in
Multicomponent Compressible Fluid Analysis○野々村拓, 宇宙機構宇宙研, 神奈川県相模原市中央区由野台 3-1-1, nonomura@flab.isas.jaxa.jp
寺門大毅, 東大院, 神奈川県相模原市中央区由野台 3-1-1, terakado@flab.isas.jaxa.jp
藤井孝藏, 宇宙機構宇宙研, 神奈川県相模原市中央区由野台 3-1-1, fujii@flab.isas.jaxa.jp
Taku Nonomura, ISAS/JAXA, 3-1-1Yoshinodai, Sagamihara, Kanagawa, Japan
Daiki Terakado, University of Tokyo, 3-1-1Yoshinodai, Sagamihara, Kanagawa, Japan
Kozo Fujii, ISAS/JAXA, 3-1-1Yoshinodai, Sagamihara, Kanagawa, Japan
In this presentation, the implementation of finite difference WENO scheme maintaining velocity, pressure andtemperature equilibrium in multicomponent compressible fluid analysis is discussed and the results computed by thisfinite difference WENO scheme are demonstrated.

１．はじめに様々な工学的問題で, 異なる気体の混合は重要な物理現象であり，衝撃波を含む高速流体場での予測も重要である．我々はロケット噴流からの音響波の解析を行っている 1,2 が，ロケット噴流と雰囲気気体は異なる気体種であり，これを考慮することで，どのようにせん断層の混合が変化し，結果音響波にどのような影響をあたえるかを数値解析で評価したいと考えている．
これまでに，多成分気体の解析には様々な方法が提案されているが，ここでは最も単純な Diffused-interface 法の内，Single-Fluid
Modeling に着目する．このようなフレームワークでは，Abgrall3が，質量分率を導入した Diffused-interface 法を比熱比の異なる物質の移送問題に適用し，完全保存系では数値振動が発生すると報告している．この問題は，界面で異なる気体種が混合した際にエネルギー式と比熱比の評価によって圧力・速度平衡が崩れることによる．Abgrall3 はこの問題に対し，質量分率のみ(もしくはエネ
ルギー式のみ 4)非保存系で解く準保存系を用いることで振動を消せることを示している．
上記の方法論は Johnsen と Colonius5 によって，高次精度有限体
積 WENO 法 6 へ拡張され，準保存系でも保存変数を利用すると振動が生ずることが示されている．彼らは，有限体積 WENO において基本変数による補間および準保存系の方程式を利用すること
で，数値振動を抑えられることを示した．保存変数を利用した場合の振動は有限差分 WENO7 でも生ずると報告されている．
Johnsen と Colonius5 による実装は振動を抑える観点で非常に優れているが，一方で有限体積ベースであるために有限差分法に比べて計算コストが 10 倍程度かかるために効率的とはいえない．
そこで，Nonomura らは差分法にもかかわらず変数補間が利用できる
WCNS 法 8,9 を利用することで，有限差分法の計算コストで数値振動を抑えた高次精度スキームを提案している 10．彼らは，エラーの発生原因を整理し不適切な数値粘性がより深刻な問題となることを指摘している．さらに Johnsen11 らは，過推定となるように変数を選んだ上で Johnsen と Colonius の方法を拡張することで，速度－圧力－温度平衡を保つ計算方法を提案した．この時，過推定となる他に，基本変数に対する WENO スキームの重みを密度と質量分率で合わせることで温度平衡も保つことができることを報告している．さらに，スプリッティングフォームの中心差分と有限体積 WENO 法のハイブリッド法を提案し，圧縮性多成分乱流向けの計算手法を提案している 12．一方，Terashima らは，Kawaiと Terashima13 が利用した多成分気体向けの LAD 法をベースに，

数値粘性の形を適切に取ることで，速度－圧力－温度平衡を保つ
数値粘性のベクトルフォームを提案している 14．具体的には，セクション 2 の WENO スキームの導入前に簡単にその形に触れた
い．
差分法ベースのスキームを選び，(Johnsen らは差分法ベースで局所的に 2 次性度有限体積法を導入しているが，その妥当性は明らかでない．)計算コストを抑えた上で非線形重み付けスキームで多成分圧縮性流体を解析する方法の一つとして，前述の WCNS法があるが，様々な亜種があるスタンダードな差分 WENO 法 15による実装も同様に価値があると考えられる．実用的な観点で，Johnsen らのスキームとほぼ同様の結果になる可能性もあるが，ここで適切な数値粘性の形を議論することは，今後より複雑な問題への適切なアプローチを指し示すことにもなろう．
そこで本研究では，差分 WENO 法により，速度－圧力－温度
平衡に関して考察し，その実装を示す．今回の実装は，それまで難しかった差分 WENO 法での一様流保持を実現した方法 16 と類似の方法であり，一様流保持の時と同様，ある種の平衡条件を満たすために差分 WENO 法の数値粘性を工夫した形である．
２．数値解析手法
本研究では 2 つの実装を検討した．
保存系 FC(Fully Conservative)，
数値粘性を修正した過推定準保存系 OQCmod(OverestimatedQuasi-conservative-modified)である.
２．１．支配方程式
まず，Single-Fluid Modeling における保存系(FC)の支配方程式を以下に示す．
  u0tx ui  (  uu  p )

0
t
x
 e  ((  e  p )u )

0
t
x

(1)

ここで，は密度，u は速度，p は圧力，e は質量あたりのエネル
ギー，ij はクロネッカーのデルタを示している．完全保存の実装の場合，式(1)に付け加えて以下の質量分率を保存系で解く．
Y1 Y1u
(2)

0
t

x

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

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C07-2
Y0  1  Y1d1  1 u

(3)と書ける．このとき，比熱比は


1
1
Y M
Y1 M
 0
 1  1 M1  0  1 M 0

(4)

Y
Y1
 0
M1 M 0

(6)

x 

x 

 ui  (  uu  p )     u

t
x
x 
x 

t

p 0 1

(10)

1 
u
1
1 1
1 
b1  b4  b2  b3     b1u  
c
2
c 2
2
 l1   2 l   1  b  b
b
u

b
b

2
2
3
1
1
4

L


l3   1 
u
1
1 1
1 
b

b


b
u

b

b
    2 3
1
4

 1

2
c
2
c 2
2 
l4 


Y1
0
0
1 


２．２．Terashima らの実装
Terashima らは，Localized Artifical Diffucsivity を利用して準保存系に対し圧力速度温度平衡を満たすために，以下のように実装した．
  u    

 

x

d3   0 0(9)

この時，支配方程式は以下のように書ける．
Q F
(11)

0
t x
ここで流束ヤコビアンの左固有ベクトルと右固有ベクトルを導入する．

ここでの比熱比 p は圧力をエネルギーから決定する際にのみ使う新たな変数である．このような定式化にして，合わせていくらかのテクニックを利用することで，圧力平衡が保てることがわかっている． 式(2)，(4)，(5）の代わりに式(6)を利用することを準保存系というが，式(2)，式(6)を同時に解く過推定システムで温度平衡が実現できることが報告されている．

t

t


 


 u  u 2  p Q   , F   e 
u (  e  p )   Y   uY 

次に，準保存系の支配方程式を示す．準保存系の場合，エネル
ギーと圧力を結びつける比熱比(を変形したもの)を移流方程式で解いている．

 1 
  1   u 0
t   p  1 x   p  1 

d2   0 0 0 1 0t

２．３．差分 WENO 法
２．３．１ 保存系(FC)での実装
まず，基本となる保存系で標準的な WENO 法の構築方法を示
す．1 次元問題を考え，解ベクトルと流束を以下の用に選ぶ．

(5)

1


0 0


である．本論文では Terashima らが利用した数値粘性の形をベク
トル型の数値粘性と呼ぶこととする．

となる．ここで M は混合気体の分子量であり，気体種 1 の分子量M1 と気体種 0 の分子量 M0 を用いて以下のように書ける．
M

1 2
u
2

R   r1 r2

r3

(7)

 1
 u c
r4    H  uc

 Y1

1
u

1
uc

1 2
u
2

H  uc

Y1

Y1

0 
0 
1 e 
 Y1 
1 

(12)

(13)

 e  ((  e  p )u )   2      1  u    p 
t
x
x 
x  x 
x   p  1  である．式中の c, H , b1, b2, b3, b4, (∂e/∂Y)は以下のように書ける：

Y1 Y1u   Y1 

 

t
x
x  x c

 1 
  1      1 
  u    t   p  1 x   p  1  x  x   p  1  H

ここで は密度の高階微分から求めた数値粘性の強さであり，は比熱比の高階微分から求めた数値粘性の強さである．この式をベクトル表記すると，


  


1

 u  u 2  p 






 1

 e 
u (  e  p ) , F  Q1, M   Y   uY 
1


 1 


1

u 




   1  1 

p


b1 

b2 

(8)



p

 1 



 1

u2
2

(15)

(16)

c2

u 2b1
2

b3  b1Y1

b4  b1

Y1   QF       1 M  d1    d2    d 2 tx x x  x 
x  x 
x   p  1  

(14)

(17)
1 e

 Y1

1 e

 Y1

(18)

(19)

となり，ここで，2

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C07-2

 1 M
e
1 M 
 p


Y1


1
M

1
0 1 M 2 
 1
 1
1 



 Y
Y
1
1   M1 M 0 
 p 1
 0

2
  1  1 M1  0  1 M 2    Y
Y   M M  
0  
 1

F
1
WENO
(29

 FjWENO
1/2  Fj 1/2 x x
本論文では，理解を深め，かつ準保存形や後述のテクニックを導入する際に見通しをよくするため，中心差分と数値粘性部分をわけた分離型 WENO スキームを導入する．上記方式のまま，準保存型や後述のテクニックを導入することも検討したが，その実装は必ずしも自明ではない．以下の分離型 WENO スキームは，
オリジナルの WENO スキームと同一のスキームである．
F
1
1
CENT
DISP
(30)

 FjCENT FjDISP
1/2  Fj 1/2  
1/2  Fj 1/2 
x x
x
ここで，
1
(31)
FjCENT Fj 2  8Fj 1  37Fj  37Fj 1  8Fj 2  Fj 3 
1/2 
60(20)

以下では，時間方向には離散化しない半離散化の状態で空間差分に対して議論を進める．差分 WENO 法ではまず，特性変数に変換した低散逸 Lax 風上流束を構築する．
1(21)
f j ,m  lm  Fj  mQ j 
2
これらの風上流束から非線形補間で，特性変数上での数値流束を求める．一般的には，以下のように書かれるが，

 , DISP
 , DISPFjDISP1/2  Fj 1/2  Fj 1/2



r
DISP
 ,m
 ,m
 ,m
Fj,1/2
  m  20w1  1 f j3,1/2,1
 10  w1  w2   5 f j3,1/2,2
 f j3,1/2,3
m 60

WENO
WENO , m
f j,1/2
  wk f j,1/2,k

1  ,m 7  ,m 11  ,m
f j  2  f j 1  f j
3
6
6
1 ,m 5 ,m 1 ,m
 , m ,WENO
f j 1/2,2   f j 1  f j  f j 1
6
6
3
1  ,m 5  ,m 1  ,m
 , m ,WENO
f j 1/2,3  f j  f j 1  f j  2
3
6
6

m ,WENOf j,1/2,1


(22)



r
DISP
,m
,m
,m
Fj,1/2
  m  20w1  1 f j3,1/2,1
 10  w1  w2   5 f j3,1/2,2
 f j3,1/2,3
60
m
 ,m




f j3,1/2,
k   f j 3 k  3 f j  2  k  3 f j 1 k  f j  k
,m




f j3,1/2,
k   f j 1 k  3 f j  2  k  3 f j 3 k  f j  4  k

このように書くと，右固有ベクトルは，数値粘性部分にのみ働いており，Terashima らが利用した数値粘性ベクトルフォームを表していると理解できる．また，このように考えると，左固有ベク
トルは非線形重みの決定のためのスカラーを作っているとも考えられ，流れ場の中の特性波を示すセンサとして利用されているというように見える．
上述で一般的な差分WENO 法と分離型WENO 法の手順を説明した．ただし，これらの方法では，準保存型では無いため，圧力平衡が満たされないことを指摘しておく．平衡が満たされない理由の一つは，過去の研究で指摘されたように，完全保全形では圧力平衡が保たれないことと，もう一つは次節で説明するように適切で無い数値流束が混入するためである．

これらの 3 つの線形結合で、滑らかな部分で 5 次精度になるよ
うに重みが求められ、不連続がある場合にはそのステンシルを含まないように重みが変化すればよい．5 次精度になる理想的な重み Ck はC1 

1
63, C2  , C3 
10
10
10

(23)

となる．これらから衝撃波付近で振動しないように以下のように重みを求めるk  3 k
 l

(24)

l 1

ここで、

k 

Ck

  ISk 2

２．３．２ 数値粘性を修正した準保存系(OQCmod)での実装
ここで提案する平衡を保つ実装が満たすべき要請として，

保存変数を保存する
スムースインジケータで，音響波がなく比熱比が変化する
場所で，音響波の変化があると判定されない

Terashima らが指摘したような適切な数値粘性になっていることが重要と考えており，前者に関しては差分を中心差分部分と数値粘性部分に分けて，数値粘性部分にはスムースインジケータが誤作動しないような人工的な流束を利用することで対応した．
後者に関しては，特に新たな提案はないが，特性波がつくる数値
粘性の形がそのままの形で Terashima らが指摘した「適切な」形となっていることを確認した．
まず本実装では，過去の研究と同様に，過推定準保存系の以下の式を考える．

(25)

であり ISk はスムースインジケータで、

13  , m
1( f j  2  2 f j,1m  f j , m ) 2  ( f j,2m  4 f j,1m  3 f j , m ) 2
12
4
13  , m
1
IS2  ( f j 1  2 f j , m  f j,1m ) 2  ( f j,1m  f j,1m ) 2
(26)
12
4
13
1
IS3  ( f j , m  2 f j,1m  f j,2m )2  (3 f j , m  4 f j,1m  f j,2m ) 2
124である．次に，特性変数に右固有ベクトル rm をかけることで，実の流束に戻す．
ここで風下流束は，
風上流束と対称に評価する．
 ,WENO ,WENO , m
(27)
F
 r f
IS1 

j 1/2



m


  


1

 u  u 2  p 
 1







 e 
u (  e  p ) , F  Q1, M   Y   uY 
1 

 1 


1

u 





   1  1 

j 1/2m

WENO
WENO , mFj,1/2
  rm f j,1/2
m

最後にこれらを利用して差分を取る．
 ,WENOWENO
FjWENO Fj,1/2
1/2  Fj 1/2



(28)
3

(32)

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C07-2
比熱比式が uj+1/2 を含んだ形である

エネルギー式の比熱比が固定されている．
1 点目は，エネルギー式と比熱比式の数値粘性が「適切」になり，かつエネルギー式を保存系で書くためである．また，2 点目は，比熱比のみが変化した場合に，1,3 番目の音波の特性波の変動を無くすためである．
後者に関しては，
例をあげ以下で説明する．
速度圧力一定で，密度比熱比のみ不連続をまたいで変化する場合を考える．一定の値は L/R の上付き添字を省略する．また，左固有ベクトル，右固有ベクトルを作成する j+1/2 の u および p の値
も，一定値を利用する．この場合，通常の流束 F で不連続の左側と右側で 1 番目の特性量の値を確認すると，

Q
F
M
0
t
x
ただし，本式を直接特性空間で実装された WENO 法で離散化する方法は自明でない．そこで，この式を変形して
Q
Q
(33)
B
0
t
x
とし，B に対する左固有ベクトルと右固有ベクトルを利用することを考えた．このとき B の左固有行列と右固有ベクトルは以下の
用になる．

1 
u
1
1 1
b1
 b2     b1u  
c
2
c 2
 l1   2 
l   1  b
b1u
b1
2
 2 L  l3   1 
u
1
1 1
b1
    b2     b1u  
l
2
c
2
c 2



 4 
l5   Y1
0
0

0
0
0

R   r1 r2

r3

r4

 1
 u c
r5    H  uc

 Y1
 0


1
u
1 2
u
2
Y1
0


0 0

0 0

0 0


1 0
0 1 

0 0
u  c 0 0 H  uc 0 p 

Y1
1 0
0
0 1 

f1L 

(34)




Mj

x

CENT
DISP
DISP
 M j  FjCENT
1/2  Fj 1/2    Fj 1/2  Fj 1/2 

1
1 
1  R Ru 2 R 2pu  b1, j 1/ 2u j 1/ 2    u  p   b1  R2 
c j 1/ 2 
2   1
2 


1
1 
1  Rpu  b1, j 1/ 2u  p  b1  R
2 
c j 1/ 2 
2   1

(38)となり，特性量の値が変化する．第一不変量は音波の特性波であり，音圧変化がない場合は本来現れるべきでないが，それが現れた．この値を利用して，WENO スキームを使用すると，不適切な
数値粘性が入れられる．一方で，提案する比熱比を固定した FDISPを使うと，以下のように不連続の左右で 1 番目の特性波は値が変化しない．
f1L , DISP 



(36)



(37)

1
1 
1   j 1/ 2
 Lu 2 L 2
pu  b1, j 1/ 2u j 1/ 2    u  p   b1 2
c j 1/ 2 
2   j 1/ 2  1
2 

u j 1/ 2  R1 b2, j 1/ 2 
 u
2 
c j 1/ 2 

数値粘性項にのみ使う風上化する流束を

u j 1/ 2  L1 b2, j 1/ 2 
 u
2 
c j 1/ 2 


1
1 
1   j 1/ 2pu  b1, j 1/ 2u  p  b1 2 
c j 1/ 2 
2   j 1/ 2  1


f1R 

ここで， FjCENT， FjDISPの定義は式とほぼ同じであるが，j+1/2 の
1/2
1/2
1
f j ,m  lm  Fj  mQ j  for m  1, 2,3, 4
2
 f j ,m  0
u j 1/2  0
  ,mDISP
u j 1/2  0
 f j  lm Fj

u j 1/ 2  R1 b2, j 1/ 2 
 u
2 
c j 1/ 2 


(35)


1
1 
1  L Lu 2 L 2pu  b1, j 1/ 2u j 1/ 2    u  p   b1  L2 
c j 1/ 2 
2   1
2 


1
1 
1  Lpu  b1, j 1/ 2u  p  b1  L
2
c j 1/ 2 
2   1


f1R 

1

である．前述したように，R の形は，数値粘性のベクトルフォー
ムを示している．
第2，
第3 および第5 ベクトルをみると，
Terashima
らが主張している適切な数値粘性の形に非常に似通っていることがわかる．(単相流であれば，完全一致である．) この形を利用することで，平衡を保つスキームを構築することを考える．
ただし，上記の特性量を利用した場合，保存系の式に対しては一般的な数値流束の構築方法があきらかであるが，準保存系ではその構築方法は確立されていない．そのため以下の式変形は試行
錯誤を重ね得たもので，分離型のWENO スキームをベースとし，
数値粘性側に対していくらか任意の変形を加える事で保存変数の保存性を保ちかつ，平衡を保つスキームの基本型として提案する．
Fju j 1/ 2  L1 b2, j 1/ 2 
 u
2 
c j 1/ 2 

1
1 
1   j 1/ 2
 Ru 2 R 2pu  b1, j 1/ 2u j 1/ 2    u  p   b1 2 
c j 1/ 2 
2   j 1/ 2  1
2 


1
1 
1   j 1/ 2pu  b1, j 1/ 2u  p  b1 2 
c j 1/ 2 
2   j 1/ 2  1(39)
このため，不必要な数値粘性が導入されない．
これらの式変形は，数値粘性パートのみに行っている．数値粘性パートは計算が安定であり，かつ細い格子を使った場合に数値粘性が適切に(べき乗則にしたがって)減っていけばそれで良い．
そのため，本論文で提案するような変更を行っても，数値テストの結果振動するなどの問題が置きなければ，精度上も実計算上でも十分な性能といえる．また，上記記述に関して，M が数値粘性項の外側に着くような実装も考えたが，数値粘性項が適切にならない，保存系にならないなどの問題があり，ここでは採用しなかった．
上記のスキームの構築により，速度，圧力平衡が保たれるスキームが構築できる．さらに温度平衡のために，第 4 特性変数(Y1)の解析時に第 2 特性変数(エントロピー)と同じ重み関数を選ぶと

 jk


  u 2pj k j k
j k




p jk
 p j k )
u j  k (
DISP


1
Fj

j 1/2 j k u j kYj ku j 1/2






1
j

k


である．この FjDISP は，以下の 2 点で通常の流束と異なっている．
4

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C07-2
実質的に Johnsen らが主張するように密度と各気体種の密度を同じ重みで計算したことに一致し．これが実現される．(Johnsen ら
は基本変数で実装しており彼らの場合には，密度の重み関数と各気体種の密度の重み関数を合わせることに一致)．本論文では，第2 特性変数で得た重み関数を第 4 特性変数に利用している．ただし，WENOCU6 スキームを利用した場合は振動が起こったため逆に，第 4 特性変数で得た重み関数を第 2 特性変数に利用した．
３．数値実験
本章では，3 つのテスト問題に適用した結果を示す．テスト中の時間積分には CFL0.4 で 3 段 TVD Runge-Kutta 法を利用した．
３．１ 移動物質境界問題
以下の Terashima らと同様の問題を考える：



u

p 




 7.0 1.0 1.0 /1.4 1.4 28 0.25  x  0.25
M
othrwise1.0 1.0 1.0 /1.4 1.66 4 (b)velocity

計算領域を-0.5<x≤0.5，格子点を等間隔 100 点とし，周期境界条件を利用した．格子幅はx=0.01 である．1 点ジャンプで不連続を表現した．
結 果 を 図 1 に 示 す ． FC-WENO5 ， OQCmod-WENO5 ，
OQC-WENOCU6 のすべてのスキームで，ほぼ圧力，速度，温度でエラーがないように見える．しかしながら，よりレンジを絞った図 2 では，FC では 10-4 のオーダーのエラーを持っていることがわかる．一方で OQC の実装では，このレンジでもエラーがなく，
より小さなレンジでも調べた上で平衡が取れていることを確認した．ここで，FC に関して，従来研究よりもエラーが少ない理由を考える．本論文で指摘したように，このケースはエントロピー波と Y1 の特性量に関して，同じ数値粘性が入れられれば，温度平衡が保てる．この場合にかぎり，Y1 の方程式が，1/(p-1)の式と同値になり，準保存系のメカニズムを利用して圧力平衡も保てるよ
うになる．(すべての点でのモル数一定の条件下の場合に，生ずる．本論文での FC は温度平衡を保つようには設計していないが，ほぼ同様の密度と質量分率でほぼ同様の不連続があり，これらに対して，同様の数値粘性が入ったために，このような結果になったと思われる．最後に，WENOCU6 の効果に関して議論する．
本解析から，WENOCU6 でも OQC の実装は可能であり，テスト
問題では，密度分布から不連続面をよりシャープに捉える．すな
わち，WENOCU6 により平衡を保つ高解像度スキームが実現できるといえる．
(c)pressure

(d)temperature
Fig. 1 Continued.
(a)density
Fig. 1 Results of the material interface problem

5

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C07-2で生ずる大きなエラーを OQCmod では，適切に取り除けること，
またこの方法は WENOCU6 方へ拡張でき，その際に高解像度が実現できるといえる．

(a)velocity(1)density(b)pressure(b)velocity(c)temperature
Fig. 2 Enlarged results of the material interface problem(c)pressure
Fig. 3 Results of the contact surface problem

３．２ 接触面境界問題
ここでも Terashima らのテストケースを参考にして，以下の問
題を考える：



u

p 


14.0 1.0 1.0 /1.4 1.4 28 0.25  x  0.25
M
othrwise 1.0 1.0 1.0 /1.4 1.66 4 

計算領域，格子点，不連続の与え方は上記と同様である．
結果を図 3 に示す． FC では，平衡が取れないことによるエラーが確認できる．一方で OQCmod-WENO5，OQC-WENOCU6 では完全に圧力，速度，温度平衡が取れている．このテストで，FC
6

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C07-2(c)temperature
Fig. 3 Continued.

(c)pressure

３．３ 衝撃波問題
以下の Terashima らの衝撃波問題を考える：



u

p 

x0
 1 0 1.0 /1.4 1.4 28
M
0.1
0
0.125
/1.4
1.66
4
othrwise

計算領域を-0.5≤x≤0.5，格子点は 101 点とした． 結果を図 4 に示
す． 今回の問題では，すべての実装で特に大きなエラー無く解析ができることがわかった．WENUCU6 に関しても前節と同様の特性を持っていると言える．

(d)specific heat ratio
Fig. 4 Continued.
４．考察－Terashima らの実装との比較－
最後に Terashima らの実装との比較を考える．２で説明したよ
うに今回 WENO スキームで出てきた，特性変数をベースにした
数値粘性は，左固有ベクトルと流束もしくは変数をかけたものが「センサ」になっており，そのセンサから得られた数値粘性の強さに右固有ベクトルを書けた形で数値粘性のべクトルフォームが書かれていると理解できる．この時に，Terashima らが主張している数値粘性のベクトルフォームと右固有ベクトルから得られるベ
クトルフォームが一致していることが興味深い．寺島らは，密度および比熱比のみ変化する特定のシチュエーションで平衡を保つことを意図してこのフォームを提案しているが，結局のところ，
このようなシチュエーションを考えることは，ある特定の特性波が存在する単純波の問題を考えることに一致しており，このような単純波の問題での適切な数値粘性が右固有ベクトルで書かれる可能性が高いことを示唆している．また寺島らは，物理現象の考察からエントロピー波の数値粘性の決定に密度を選んでいるが，密度は音波でも変化するため，本論文で提案するスキームが利用する左固有ベクトルを流束／変数にかけたものを「センサ」として利用するほうが良いことを示唆している．
すなわち，Terashima らの実装も含め，特性空間での定式化はそのままでも「適切」な数値粘性の形に繋がる可能性が高く，任意のスキームで利用することや参考にしてその数値粘性の形を決めることが，より複雑なシステムでの，適切な数値粘性の導入に役立つと考えられる．

(a)density(b)velocity
Fig. 4 Results of the shock problem
7

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C07-2
Diffusion Terms for Simulating Compressible Multicomponent
Flows,” Computers & Fluids , Vol.88, (2013), pp. 484 – 495.
(15) Jiang, G.-S. and Shu, C.-W., “ Efficient Implementation of
Weighted ENO Schemes”, Journal of Computational Physics，126,(1996), pp. 200-212(16) Nonomura, T., Terakado, D., Abe, Y., and Fujii, K., “A New
Technique for Freestream Preservation of Finite-difference WENO
on Curvilinear Grid,” Computers & Fluids , to appear (2015).

５．まとめ
本論文では,圧力-速度-温度平衡を保つ差分 WENO 法の実装を提案した．提案手法により，１次元問題ではあるが，圧力-速度温度平衡が保てることがわかった．また，近年の高解像度スキー
ムである，WENOCU6 法にも利用できることがわかった．
また一連の解析から，特性量変換する左固有ベクトルは流場のセンサとして，右固有ベクトルは適切な数値粘性の形として有望であることがわかった．今後，WENO スキームだけでなく，様々な支配方程式，数値解法で上記のアイデアを利用して，適切な数
値粘性の入れ方を検討することが有用であるといえる．
参考文献(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) 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.
(4) Abgrall, R. and Karni, S., “Computations of Compressible
Multifluids,” Journal of Computational Physics, 169, (2001), pp.
594-623.
(5) Johnsen, E. and Colonius, T., “Implementation of WENO Schemesin Compressible Multicomponent Flow Problems” Journal ofComputational Physics, 219, (2006), pp. 715-732.
(6) Titarev, V. A. and Toro, E. F., “Finite-volume WENO Schemes for
Three-dimensional Conservation Laws,” Journal of Computational
Physics, Vol. 201, (2004), No. 1, pp 238-260.
(7) Johnsen, E., “Spurious oscillations and conservation errors ininterface-capturing schemes,” Annual Research Briefs 2008, (2008),
Center for Turbulence Research.
(8) Deng, X. G. and Zhang, H., “Developing high-order weighted
compact nonlinear schemes,” Journal of Computational Physics，
165, (2000), pp.22-44.
(9) 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.
(10) Nonomura, T., Morizawa, S., Terashima, H., Obayashi, S., and Fujii,
K., “Numerical (error) issues on Compressible Multicomponent
Flows Using a High-order Differencing Scheme: Weighted
Compact Nonlinear Scheme,” Journal of Computational Physics,
Vol. 231, (2012), pp. 3181-3210.
(11) Johnsen, E. and Ham, F. E. “Preventing Numerical Errors
Generated by Interface-capturing Schemes in Compressible
Multi-material Flows,” Journal of Computational Physics, Vol. 231,
(2012), pp. 5705-5717.
(12) Movahed, P. and Johnsen, E., “A Solution-adaptive Method for
Efficient Compressible Multifluid Simulations, with Application to
the Richtmyer-Meshkov Instability,” Journal of Computational
Physics, Vol. 239, (2013), pp. 166-186.
(13) Kawai, S. and Terashima, H., “A High-resolution Scheme forCompressible Multicomponent Flows with Shock Waves,”
International Journal for Numerical Methods in Fluids, Vol. 66,(2011), No. 10, pp. 1207-1225.
(14) Terashima, H., Kawai, S. and Koshi M., “Consistent Numerical8

Copyright © 2014 by JSFM