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

複数粒子が高マッハ数・低レイノルズ数で移動する流れ解析コードの構築と評価
Development and Validation of a Flow Solver for Moving Multiple Particles at High Mach number
and Low Reynolds number○ 水野裕介，東海大，神奈川県平塚市北金目 4-1-1，1beu2104@tokai-u.jp
高橋俊 ，東海大，神奈川県平塚市北金目 4-1-1，takahasi@tokai-u.jp
野々村拓，宇宙研，神奈川県相模原市中央区由野台 3-1-1，nonomura@flab.isas.jaxa.jp
福田紘大，東海大，神奈川県平塚市北金目 4-1-1，fukuda@tokai-u.jp
Yusuke MIZUNO, Tokai Univ., 4-1-1 Kitakaname, Hiratsuka, Kanagawa, JAPANShun TAKAHASHI, Tokai Univ., 4-1-1 Kitakaname, Hiratsuka, Kanagawa, JAPANTaku NONOMURA, ISAS/JAXA, 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa, JAPANKota FUKUDA, Tokai Univ., 4-1-1 Kitakaname, Hiratsuka, Kanagawa, JAPANThe aim of this work is development and validation of a flow solver that the flow around the multiple particles at highMach number and low Reynolds number using direct numerical analysis of the three dimensions compressible
Navier-Stokes equation for the purpose of high accurate prediction of the acoustic field around a rocket. These particles
are imitation of alumina particles exhausting from the rocket nozzle or water drops surrounding the rocket. The
condition of flow field is the Mach and Reynolds number are set to be 0.3-2.0 and 50-300, respectively. The static
particle is compared the drag coefficient with result obtained from a high-order boundary-fitted-coordinate (BFC)
solver. We investigated the flow around moving particles of number of up to 16.

が数 μm であり，さらに流れがロケット後流を模擬している高速流れである．そのため高 Mach 数，低 Reynolds 数流れを対象とす
る．そのために，検証問題として高 Mach 数，低 Reynolds 数の静止球体まわりの流れ場を対象とし，開発した計算コードの精度を
曲線座標格子(BFC)の結果と比較することで検証する．検証後，複数の移動する球体における流れ場の検証をする．

はじめにロケットノズルから噴出される超音速ジェットは非常に強い音響波を発生させ，フェアリング内の人工衛星や宇宙機のペイロードに損傷を与える可能性がある．これまでは欧米や日本では経験則を用いた予測手法を用いてきた(1)．この手法は米国の莫大な打ち上げ実績や試験データを元に作成されたものである．しかしながら，新ロケットや射場設備の形状が変化すると必ずしも高い精度とは言えず，予測値と実機による測定データとの間に大きな差異が生じる．ロケットの新規開発段階において，音響環境の高い予測精度と定量的な評価ができる手法を必要とする．これに向けてロケット後流内部の物理現象を詳しく理解する必要がある．そのため数値解析による研究が行われてきており，様々な仮定を入れているが定性的な予測が可能になってきた(2-3)．しかしながらより高精度の音響波の予測のためには，現状で考慮できていない物
理現象の切り分け／精査が必要とされる．
実際のロケット打ち上げ時はロケット後流に散水することで音響波を抑制している．福田らは音圧が減衰する原因はロケット後
流内の微粒子と乱流の干渉による音源強さの抑制であるという可能性を示唆した(4)．ここでの微粒子とはロケットノズルから噴出するアルミナ粒子や散水時の液滴である．しかしながら，その現象のメカニズムはいまだに十分にわかっていない．そのためその微粒子周りの流れ場を解析し，基となる物理現象から理解する必要がある．そのため，高速乱流場中における固気二相流の振る舞
いを理解するための直接数値解析の実施が必要となる．
この直接数値解析では複数の運動する微粒子を扱う．この解析には，直交格子を用いた Euler-Euler 型の数値解析が活用でき，この場合微粒子の形状は埋め込み境界法(IBM)を用いた直交格子上で表現される．高橋らは上述の解析を行うための準備として，2次元流れ場を対象としているが高解像度の埋め込み境界法を開発した(5)．
本研究では高橋らが開発した計算コードを 3 次元に拡張し，その性能を検証する．また，本解析の対象である球体周りの流れは
アルミナ粒子まわりや水滴まわりの流れを模擬しており，球直径
1.

2. 数値計算法
2. 1 支配方程式
本研究では支配方程式を式(1)に示す 3 次元圧縮性 Navier-Stokes
方程式とする．

Q E F G Ev Fv Gv






t x y zxyz u  v   w  vu   u 2 u  p   wu  
2Q   v , E   uv , F   v  p , G  wv ,



 

 
2
 uw  vw   w  p  w( e  p)u ( e  p)v  ( e  p) w e 




 

(1)







0
0
0






 yx
 xx
 zx






, G  
, F  Ev  
 xy


yy
zy
 v 

 v 







 xz
 yz
 zz






 yx u   yy v   yz w  Ty  zxu   zy v   zz w  Tz  xx u   xy v   xz w  Tx 

ここで E, F, G, Ev, Fv, Gv はそれぞれ x, y, z 軸方向の非粘性流束，粘性流束であり，Q は保存変数ベクトルである．式中の圧力 p と単位質量当たりのエネルギ e は，密度 ρ，各軸方向の速度 u, v, w 比
熱比 γ を用いて式(2)の状態方程式で関係づける．

1

Copyright © 2014 by JSFM

e 

p
1
  (u 2  v 2  w 2 ) 1 22
32
 yy   (2v y  wz  u x )32
 zz   (2wz  u x  v y )3 xy   yx   (u y  v x )第 28 回数値流体力学シンポジウム
講演番号 C04-3
スの圧力・応力のみで空気力を評価する方法を利用した．本手法
は，ゴーストセルに値が入っていれば，体積力評価時に改めて物体を表現するレベルセット関数やポリゴンを使った評価を行う必要がない．合わせて，その精度もポリゴン等を使って正確に求めた空気力とほぼ変わらないことがわかっている．(8)(2)

 xx   (2u x  v y  wz )2. 4 スポンジ境界
本手法では完全等間隔直交格子並びに Skew-symmetric 型のス
イッチングスキームを導入しているため，亜音速流れを解析する
場合に音波の反射によって解析の不安定化を招く場合がある．それを防ぐために，外部境界近傍に Mani によって開発されたスポンジ境界を設置した．(9-10)本計算手法では以下の式(6)によってス
ポンジ層の強さを流れ場に応じて設定する．スポンジの強さは減衰率を表す ηtarget の式(7)を元に，スポンジの層数を表す Lsp，スポンジの強さを制御するパラメータ σsp，β を境界ごとに設定し，最終的に式(8)によって減衰率の概算値を求める．式(8)の ηtarget はデシベル単位の減衰率であり，Mani の論文では 3 桁程度の減衰率を設定することで流れ場への影響を非常に小さく抑制できることが示されている．本研究では等間隔直交格子による計算負荷の増大を抑制すべく，このスポンジ層を活用して計算領域を適切な大きさに設定している．

(3)

 yz   zy   (v z  w y ) zx   xz   ( wx  u z )C pT 

e  p 1 2
 (u  v 2  w 2 )

2

(4)

また，
粘性応力テンソル τ と温度T は式(3)と式(4)から求められる．
ここで，μ, Cp, κ はそれぞれ粘性係数，定圧比熱，熱伝導係数である．
上記の式の全ての変数は解析に当たって一様流の密度，
音速，
代表長さによって無次元化される．またこれらの方程式はセル中心の等間隔直交格子上で離散化される．式(1)の非粘性流束は圧力波と乱流の干渉場を解くために 3 次精度 MUSCL スキーム(6)と 2
次精度の Skew-symmetric 型のスキーム(7)から構築されるハイブリ
ットスキームによって計算される．

target  20

 x x x    sp  sp
 x 
 sp 
xsp  x  Lsp

2. 2 埋め込み境界法
本研究では物体境界はレベルセット法とゴーストセル法を利用
して埋め込み境界で取り扱う．レベルセット関数は各セル中心の座標から物体境界までの垂直距離に符号がついた符号付距離の値である．本研究では，このレベルセット法を複数物体の解析に適用するために，重合格子法の最小距離アルゴリズムを拡張した．
レベルセット関数の値からセルを 3 種類(流体セル，物体セル，ゴーストセル)に分類する．ゴーストセルは流体領域と物体領域の間で境界条件を与えるガードセールとして機能し，本研究では 2 層分定義するように式(5)に基づいてセルを分類した．下付き文字のFC，GC，OC はそれぞれ流体セル，ゴーストセル，物体セルを指
している．

(6)



(7)



2 log e xsp  xsp  x  target  20
 sp dx
1  M 2 0
 xsp  20

(8)

2 log e  sp xsp
1 M 2  1

3. 静止球体の検証
3. 1 計算条件
本研究では，三次元球体周りの流れ場において，BFC の結果と比較して検証する．比較する計算ケースを表 1 で示す．格子解像度は球直径を D として，0.1D，0.05D の 2 種類，Mach 数は 0.3，
0.8，1.2，2.0 の 4 種類，Reynolds 数は 100，200，300 の 3 種類とし，合計 21 ケースの比較を行う．
ここで，BFC には Nonomura らによる WENOCU6 法を曲線座標に拡張したソルバーを利用した．
格子点数は 107*48*177 である．

d PC  0
d GC  0 and d GC  3x

2 log e xsp
 x dx
1  M 2 0

(5)

d OC  3x
本研究ではイメージポイントを用いるゴーストセル法を採用する．
イメージポイントはゴーストセルから物体壁面に対して垂直方向に伸びたプローブの先端にある点である．このプローブ長さは本
研究では格子幅の 1.75 倍(1.75Δx) で固定である．
この 1.75Δx という長さは再帰参照をさけるためセルの斜辺の長さ√3Δx より大きくなるように設定した．このイメージポイントを囲む流体セル群からイメージポイントに流れの諸量を内挿する．本研究ではこの
内挿には Tri-linear 内挿を用いた．このイメージポイントの値からゴーストセルにおける値を決定する．

(11-12)

計算領域を x-y 平面で図 1 に示す．z 軸方向においても y 軸方向と同じ計算領域の量である．流入境界では全ての諸量にDirichlet 条件，
上下境界と流出境界の諸量に Neumann 条件を課す．
ただし，亜音速流れである Mach0.3 の場合の流出境界は密度のみ一様流の値で固定する．スポンジ層は計算領域の外部境界から 10層分に配置する．ただし，Mach1.2 と Mach2.0 の場合はスポンジ層を配置していない．計算格子を図 2 に示す．黒線がセル中心をつないだ格子線，赤線が物体境界，青色，灰色，橙色の領域がそれぞれ流体セル，ゴーストセルセル，物体セルである．

2. 3 空気力計算
空気力の評価には流体セルとゴーストセルの界面のセルフェイ
2

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C04-3
Table 1 Test cases.
Re

Uniform flow Mesh size Mesh number
0.3
0100D
100.100.100
0100D
100.100.100
0.8
0050D
100.100.100
100
0100D
100.100.100
1.2
0050D
100.100.100
0100D
100.100.100
2.0
0050D
100.100.100
0.3
0100D
100.100.100
0100D
100.100.100
0.8
0050D
100.100.100
200
0100D
100.100.100
1.2
0050D
100.100.100
0100D
100.100.100
2.0
0050D
100.100.100
0.3
0100D
100.100.100
0100D
100.100.100
0.8
0050D
100.100.100
300
0100D
100.100.100
1.2
0050D
100.100.100
0100D
100.100.100
2.0
0050D
100.100.100

Case
100-03-10-100
100-08-10-100
100-08-20-100
100-12-10-100
100-12-20-100
100-20-10-100
100-20-20-100
200-03-10-100
200-08-10-100
200-08-20-100
200-12-10-100
200-12-20-100
200-20-10-100
200-20-20-100
300-03-10-100
300-08-10-100
300-08-20-100
300-12-10-100
300-12-20-100
300-20-10-100
300-20-20-100(a) 0.1D(a) D=1(b) 0.05D
Figure 2 Computational grid around a particle.
3. 2 Mach 数 0.3 の球体周りの比較
Mach0.3 の場合の全体の抗力係数，圧力抗力係数，摩擦抗力係
数を Reynolds 数ごとに図 3 に示す．赤線は球直径に対して格子を10 分割して得られた値であり，黒線が BFC によって得られた参考値である．抗力係数は近い値の結果である．しかし，BFC の結果と比較すると圧力抗力係数は高く，摩擦抗力係数は低い値である．これは直交格子で球体を再現しているため，凹凸ができてしまい解像度が足りていないためである．また，圧力分布を図 4 に示す．球体の後流は非定常流れになっている．その結果，球体周りの分布は似ているが後流は異なっている．なお(a)において後方に圧力分布が生じているがこれは外部境界に配置するスポンジ層の影響である．

(b) D=2
Figure 1 Computational domain.
3

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C04-3(a) Average drag coefficient(a) 300-03-10-100(b) Average pressure drag coefficient(b) BFC-03
Figure 4 Pressure contours.
3. 3 Mach 数 0.8 の球体周りの比較
Mach0.8 の場合の全体の抗力係数，圧力抗力係数，摩擦抗力係
数を Reynolds 数ごとに図 5 に示す．橙色は球直径に対して格子を20 分割して得られた値である．10 分割に比べると 20 分割した結果が BFC の結果により近い値である．また，抗力係数の値が近い
20 分割の結果の圧力分布を図 6 に示す．球体の後流は非定常流れとなっている．(a)は球体が外部境界に接近しているため，流れが球体の前方から後方に押されて分布が変化している．また，これはスポンジ層の影響もある．

(c) Average friction drag coefficient
Figure 3 Aerodynamic coefficients at Mach 0.3.

4

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C04-3(a) Average drag coefficient(a) 300-08-20-100(b) Average pressure drag coefficient(b) BFC-08
Figure 6 Pressure contours.
3. 4 Mach 数 1.2 の球体周りの比較
Mach1.2 の場合の全体の抗力係数，圧力抗力係数，摩擦抗力係
数を Reynolds 数ごとに図 7 に示す．M0.8 同様，20 分割の結果が近い値である．ただし，圧力抗力係数が近い値の結果を示しているのに対して，摩擦抗力係数が低い値を示しており，これにより全体の抗力係数が低い値となる．これは BFC に比べ，境界層の解像度が足りないためである．20 分割の結果の圧力分布を図 8 に示
す．M0.8 同様，外部境界に接近しているため，流れが球体の前方から後方に押されて分布が変化している．しかし，この解析では(超音速であるため)スポンジ層を与えていないため，後方の分布には大きく影響を与えていない．

(c) Average friction drag coefficient
Figure 5 Aerodynamic coefficients at Mach 0.8.

5

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C04-3(a) Average drag coefficient(a) 300-12-20-100(b) Average pressure drag coefficient(b) BFC-12
Figure 8 Pressure contours.
3. 5 Mach 数 2.0 の球体周りの比較
Mach2.0 の場合の全体の抗力係数，圧力抗力係数，摩擦抗力係
数を Reynolds 数ごとに図 9 に示す．Re300 の 10 分割したケース
は計算の途中で不安定化する．M0.8，M1.2 の結果と同様 20 分割の結果が近い値である．異なる Reynolds 数では 10 分割の結果を可視化すると球前方に正しく離脱衝撃波が形成されない．そのため，M2.0 では 10 分割のケースは正しく流れ場を解けていないことがわかった．一方で，20 分割したケースは解けており，物体と離脱衝撃波の間の格子点数が十分にあることが正しい流れ場を得
るうえで必要と考えられる．また，20 分割の結果の圧力分布を図10 に示す．衝撃波の形状並びに後流の流れもほぼ同じ分布を示し
た．

(c) Average friction drag coefficient
Figure 7 Aerodynamic coefficients at Mach 1.2.

6

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C04-3(a) Average drag coefficient(a) 300-20-20-100(b) Average pressure drag coefficient(b) BFC-20
Figure 10 Pressure contours.
3. 6 本解析コードの評価
空気力係数は M0.3 では直径に対して 10 点，M0.8，M1.2，M2.0
では直径に対して 20 点用意することで BFC に近い値である．流
れ場は M0.3，M0.8 では球後流の非定常流れ，M1.2，M2.0 では球
前方に離脱衝撃波を形成した．各ケースにおける評価を表 2 にまとめた．BFC の結果と比べて，抗力係数の誤差が 5%より小さいケースを「◎」，抗力係数の誤差が 5%以上のケースを「○」，計算が不可能なケースを「×」と示す．これらの結果を踏まえると，
本解析コードは高 Mach 数，低 Reynolds 数の静止球体まわりの流れ場を対象とした解析において適応性が十分にあることがわかる．

(c) Average friction drag coefficient
Figure 9 Aerodynamic coefficients at Mach 2.0.

7

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C04-3
4. 1 球体 4 個の計算結果
各軸方向の力を移動速度から求められる基準動圧と球体直径で正規化した軸力係数の無次元時間 t の変化を図 12 に示す．t=0 で球体が動きだし，t=7 付近で物体が計算領域外に出るため，値がゼロになっている．また，t=3 付近で物体同士がすれ違うために，
各軸方向において変動がある．また，t=1.2，2.4，3.6，4.8 の時の圧力分布を 3 次元場で図 13 に示す．
各時刻において球前方で衝撃波が生じ，t=1.2 で物体が移動することにより圧力分布に変動が生
じ，t=2.4 で球体の前方の流れ並びに衝撃波の干渉，t=3.6 で衝撃波の干渉，t=4.8 で球体の後方の流れの干渉がわかる．

Table 2 Evaluation of test cases.
Re

Uniform flow Mesh size Mesh number
0.3
0100D
100.100.100
0100D
100.100.100
0.8
0050D
100.100.100
100
0100D
100.100.100
1.2
0050D
100.100.100
0100D
100.100.100
2.0
0050D
100.100.100
0.3
0100D
100.100.100
0100D
100.100.100
0.8
0050D
100.100.100
200
0100D
100.100.100
1.2
0050D
100.100.100
0100D
100.100.100
2.0
0050D
100.100.100
0.3
0100D
100.100.100
0100D
100.100.100
0.8
0050D
100.100.100
300
0100D
100.100.100
1.2
0050D
100.100.100
0100D
100.100.100
2.0
0050D
100.100.100

Case
Evaluation
100-03-10-100○
100-08-10-100○
100-08-20-100◎
100-12-10-100○
100-12-20-100○
100-20-10-100
×100-20-20-100◎
200-03-10-100◎
200-08-10-100◎
200-08-20-100◎
200-12-10-100○
200-12-20-100○
200-20-10-100
×200-20-20-100◎
300-03-10-100◎
300-08-10-100○
300-08-20-100◎
300-12-10-100○
300-12-20-100○
300-20-10-100
×300-20-20-100◎

4. 複数の移動する球体の検証
計算ケースを表 3 に示す．また，球体を 4 個(O1~O4))配置した
場合の計算領域と計算条件を図 11 に示す．z 軸方向においても y軸方向と同じ計算領域の量である．4 個の球体は同時刻に動きだ
し，同じ速度 U1=Mach1.2 で移動する．球体を 16 個配置した場合
は速度 M-0.8～0.8 の間の速度を各軸方向に与え移動をする．初期
場で流れを与えず，四方の境界を Neumann 条件と設定する．

(a) x-axis

Table 3 Test Cases.

Particle number Mesh size Mesh number
40100D
100.100.100
160050D
200.200.200

Case
04-10-100
16-20-200(b) y-axis

Figure 11 Computational domain.

(c) z-axis
Figure 12 Temporal variation of force coefficient.
8

Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C04-3(a) Time = 1.2(d) Time = 4.8
Figure 13 Pressure contours.
4. 2 球体 16 個の計算結果
t=1.2，2.4，3.6 の時の圧力分布を 3 次元場で図 14 に示す．時間が経過すると同時に球体が移動し始める．その影響を受けて，圧力場が変化していることがわかる．

(b) Time = 2.4(a) Time = 1.2(c) Time = 3.6
9Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 C04-3
参考文献(1) K. M. Eldred, “Acoustic Loads Generated by Propulsion System”,
NASA SP-8072, 1971
(2) T. Nonomura, K. Fujii, S. Tsutsumi, “Toward accurate simulation
and analysis of strong acoustic wave phenomenon – A review fromthe experience of our study on rocket problems”, International
Journal for Numerical Methods in Fluids, 1412-1432, 2010
(3) K. Fukuda, S. Tsuysumi, K. Fuji, K. Ui, T. Ishii, H. Oinuma，J.
Kazawa and K. Minesugi, “Acoustic Measurement and Predictionof Solid Rockets in Static Firing Tests”, The 15th AIAA//CEAS
Aeroacoustics Conference, AIAA paper 2009-3368, 2011.
(4) K. Fukuda, S. Tsutsumi, T. Shimizu, R. Takaki, and K. Ui,“Examination of Sound Suppression by Water Injection at Lift-off
of Launch Vehicles”, The 17th AIAA/CEAS Aeroacoustics
Conference, AIAA paper 2011-2814, 2011.
(5) S. Takahashi, T. Nonomura and K. Fukuda, “A Numerical Scheme
Based on an Immersed Boundary Method for compressible
Turbulent Flows with Shock: Application Two-dimensional Flowsaround Cylinders”, Journal of Applied Mathematics, vol. 2014
(6) B. van Leer, “Towards the Ultimate Conservative Difference
Scheme. IV. A New Approach to Numerical Convection.”, Journalof Computational Physics, Vol. 43, no. 2, pp. 357-372, 1981
(7) K. Georgesa, G. Winckelmans, P. Geuzainea, “Improving
shock-free compressible RANS solver for LES on unstructured
meshes” Journal of Computational and Applied Mathematics, Vl.
215, pp. 419-428, 2008
(8) T. Nonomura, J. Onishi, K. Fujii, “A simple Hydrodynamic force
evaluation method for immersed boundary methods”, 第 27 回数
値シンポジウム講演集, 2013
(9) A. Mani, “Analysis and optimization of numerical sponge layers as
a nonreflective boundary treatment”, JCP, Vol. 231, pp. 704-716,
2012(10) J. B. Freund, “Proposed Inflow/Outflow Boundary Condition for
Direct Computation of Aerodynamic Sound”, AIAA Journal, Vol.
35, No. 4, 1997(11) T. Nonomura, D. Terakado, Y. Abe, K. Fujii, “A new Technique forFreestream Preservation of Finite-difference WENO on Curvilinear
Grid”, Computers & Fluids Journal, 2014(12) X.Y. Hu, Q. Wang, N. A. Adams, “An adaptive central-upwind
weighted essentially non-oscillatory scheme”, Journal ofComputational Physics, Vol. 229, pp. 8952-8965, 2010(b) Time = 2.4(c) Time = 3.6
Figure 14 Pressure contours.

5. 結言
本研究は 3 次元圧縮性 Navier-Stokes 方程式に対して，高 Mach
数，低 Reynolds 数の流れ場における球体周りの流れを解析することによって，性能を検証し，その適応性を調査した．静止球体周
りで Mach 数を 0.3，0.8，1.2，2.0，Re 数を 100，200，300 として，
2 つの異なる格子解析度によって検証を行った．空気力係数は参考値と比較することで直径に対して 10 点，
ケースによっては直径に対して 20 点用意することで良好な結果を得た．また，複数の移動する球体を解析することによって，複雑問題の適用性を確認した．今後は，球体同士の衝突判定並びに運動方程式を導入し，物体移動の計算の性能を検証すると共に流れ場の解析を行う．

謝辞
本研究は科研費(24656522)の助成を受けたものである．ここに感謝の意を示す．
10

Copyright © 2014 by JSFM