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

衝撃波・乱流干渉流れ解析のための埋め込み境界法を用いた数値解析法
High Resolution Flow Simulation based on Immersed Boundary Method for Compressible
Turbulent Flows with Shocks○ 高橋俊 , 東海大, 神奈川県平塚市北金目 4-1-1, takahasi@tokai-u.jp
野々村拓, 宇宙研, 神奈川県相模原市中央区由野台 3-1-1, nonomura@flab.isas.jaxa.jp
福田紘大, 東海大, 神奈川県平塚市北金目 4-1-1, fukuda@tokai-u.jp
Shun 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, JAPANA computational code adopting immersed boundary methods for compressible gas-particle multiphase turbulent flows
is developed and validated through two-dimensional numerical experiments. The turbulent flow region is modeled by asecond-order pseudo skew-symmetric form with minimum dissipation, while the monotone upstream-centered scheme
for conservation laws (MUSCL) scheme is employed in the shock region. The present scheme is applied to the flow
around a two-dimensional cylinder under various freestream Mach numbers. Compared with the original MUSCL
scheme, the minimum dissipation enabled by the pseudo skew-symmetric form significantly improves the resolution of
the vortex generated in the wake while retaining the shock capturing ability. In addition, the resulting aerodynamic
force is significantly improved. Also, the present scheme is successfully applied to moving two-cylinder problems.

１．はじめにロケット後流から生じる音波は非常に強く，ロケットのフェア
リング内に搭載されている衛星に損傷を与える可能性がある．これまでの研究では音圧レベルは経験的に予測されていたが(1)，新型のロケット発射場に対するこの手法の予測精度は必ずしも高くない．これを向上させるにはロケット後流内部の物理現象に基づいた革新的な近似モデルが必要である．そのために数値解析を応用した研究がこれまで行われてきた(2-5)．ロケット後流で生じる音
波は高温かつ高速の固気液混相流から構成されているため，正確に予測することが困難である．実際のロケット発射場ではロケッ
ト後流に散水することでこの音波を抑制している．福田ら(6)はロケット後流内にある液滴やその他微粒子が音波の減衰にはほとん
ど影響しないことを示すことで，音圧が減衰する原因は微粒子と乱流の干渉による音源強さの抑制であるという可能性を示唆した．
しかしながら，その現象はいまだよくわかっておらず，この音源
強さの抑制現象を高精度に見積もるためにはより基礎的な解析が必須である．そこで液滴を含むマルチスケールな高速二相流を解析することで，この物理現象理解を狙う．本解析ではロケット後流が解析対象であるため，流れ場中では多数の衝撃波を同時に扱う必要がある．一連の解析では，液滴と乱流が干渉するマイクロスケール流れを DNS により，液滴の振る舞いをモデル化した中間スケールを LES により，そしてマクロスケールの乱流場をRANS によって解析する．本研究はこの中のマイクロスケール流れ解析に相当し，高速乱流場中における固気二相流の振る舞いを理解するための直接数値計算を用いた研究である．
この直接数値計算では多数の運動する微粒子を扱う必要がある．この解析には微粒子の運動にともない格子の再生成等が必要となるため，
構造格子や非構造格子を用いた解析は容易ではない．
このような場合には直交格子を用いた Euler-Euler 型の数値解析が活用でき，この場合微粒子の形状は埋め込み境界法を用いて直交格子上で表現される．本研究では直交格子に単純かつ高解像度の埋め込み境界法を用いて，微粒子を含む固気二相高速乱流の解析手法を開発する．
埋め込み境界法は元来非圧縮性流体解析(7)から発展されてきており，その後圧縮性流体解析に適用されてきた(8-9)．そしてこれまでの多くの埋め込み境界法では流れ場を安定化させるために風上近似が適用されてきた．しかし，風上近似はその過大な数値粘性から我々の主目的の一つである乱流場の高精度な解像に対して好ましくなく，乱流場を高精度に解像するために数値粘性を最小化する必要がある．また本研究では流れ場中に衝撃波も含まれるため，衝撃波由来の数値振動を防ぐために風上化による数値粘性は不可欠である．これらの乱流と衝撃波を同時に解析するために，
数値粘性のスイッチングスキームを導入する．本研究では 2 次元の検証計算によってスイッチングスキームの性能を検討する．
２．数値計算法
２．１ 支配方程式
本研究では支配方程式を2 次元圧縮性 Navier-Stokes 方程式とする．直接数値計算を行うため平均化操作やフィルタ操作は行わな
い．

Q E F Ev Fvt x yxy u  u  u  u 2  p , F   vu ,Q   , E  2 v 
 uv  v  p 
 




 e ( e  p )u ( e  p )v 
0
0








 xx
 yx
, F  Ev  
v




 xy
 yy




 xx u   xy v  Tx  yx u   yy v  Ty (1)

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

e 
1

p
1
  (u 2  v 2 ) 1 2(2)
Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
C02-2
式(9)中の Ʌ は Ʌ=LAR を満たす対角行列である．式(8)中のコンシステント部を式(7)中の skew-symmetric 型に置き換えることで，式の精度を落とすことなく，式(10)のように変形できる．

式(1)中の粘性応力項は以下のように表される．

2
3

2
3

 xx   (2u x  v y ),  xy   yx   (u y  vx ),  yy   (2v y  u x )(3)

Ei 1 2, shock  Ei 1 2, cent 

また，式(1)中の熱流束項は，以下の式(4)と式(5)を用いて最終的に式(6)の形で計算される．ここで Cp は定圧比熱，T は温度，κ は熱伝導係数である．

C pT  e 
Pr 

p



1
 (u 2  v 2 )
2

Tx 

Ei1 2  (1   i1 2 ) Ei1 2, turbulen   i1 2 Ei1 2, shoc
 Ei1 2, cent  (1   i1 2 ) i1 2   i1 2(5) (10)

式(7)と式(10)を組み合わせて，本研究で用いるハイブリッドス
キームを式(11)のように構築する．

(4)

Cp



1
Ai 1 2, Roe Qi 1 2, R  Qi 1 2, L 
2





1
Ai1 2,Roe Qi1 2,R  Qi1 2,L 
2

(11)



  p   p  , Ty  Pr   1    x
Pr   1    y(6)

このようにして，衝撃波もしくはポテンシャル領域では数値粘性を付加すると同時に，その他の領域では数値粘性を最小化するスキームが構築される．ここで式中の βi+1/2 は Ducros 型の衝撃波センサー(15)で，次の式(12)で定義される．

この手法の詳細は文献(10)に示されている．上記の式の全ての変数は解析に当たって一様流の密度，音速，代表長さによって無次元化される．またこれらの方程式はセル中心の等間隔直交格子上で離散化される．
衝撃波近傍とポテンシャル流れに近い領域以外での過大な数
値粘性を抑制するため，非粘性流束項は skew-symmetric 型の差分スキーム(11)と MUSCL-Roe スキーム(12-13)のハイブリッドスキーム
によって計算される．乱流領域において用いられる
skew-symmetric 型のスキームは以下の式(7)で示される．ただしこ
こでは x 方向の非粘性流束のみを示す．

 i1 2  min(1, i  i1 )

0

 
1





(7)

 u

(12)

2

 u  u  
2



2



(13)

また y 方向の偏微分項∂F/∂y 等も同様にして求める．
粘性項は二次精度有限差分的に求め，時間進行には三段階のTVD Runge-Kutta スキーム(18)を用いる．
２．２ 埋め込み境界法
本研究では物体境界はレベルセット法とゴーストセル法による埋め込み境界法で扱う(19-20)．レベルセット関数は各セル中心の座標から物体境界までの垂直距離に符号がついた符号付距離の値である．レベルセット関数の模式図を図 1 に示す．レベルセット法の利点はレベルセット関数値の勾配から物体の法線方向が容易に求められることである．本研究では，このレベルセット法を多数物体の解析に適用するために，重合格子法の最小距離アルゴリ
ズム(21)を拡張した．次にレベルセット関数の値から図 1 のようにセルを三種類（流体セル，物体セル，ゴーストセル）に分類する．
ゴーストセルは流体領域と物体領域の間で境界条件を与えるガー
ドセルとして機能し，本研究では 2 層分定義するように式’(14)に



ここで A は非粘性流束のヤコビアン，下添え字の Roe は左右の状態量の Roe 平均された値を指す．ここで，|A|は左右の固有値行
列 L と R を用いて式(9)のように書ける．

A  RL u  u  2E Ei 1 2  Ei 1 2x
x

ここで下添え字の i は i 番目の格子点を指し，L と R はそれぞれ左側と右側の状態量を表し，三次精度の MUSCL 法(13)と vanAlbada のリミッター(14)を用いた内挿により求められている．
また，衝撃波とポテンシャル流れの領域では，非粘性流束は次の式(8)で表わされる二次精度のMUSCL-Roe スキームから求めら
れる．
1
(8)
Ei 1 2, shock  Ei  Ei 1  Ai 1 2, Roe Qi 1 2, R  Qi 1 2, L
2

2

ここで，式(12)中の ε は ε=10-15 程度のゼロ割を防ぐための小さな値であり，スイッチングのしきい値となる θ には本研究では
θ=0.6 の値を用いる．式(12)内の速度の発散と回転を求める際は，各格子点において 2 次精度有限差分的に求める．本研究ではDucros 型のセンサーのみで衝撃波とポテンシャル流れの領域を
センシングしているが，過去の研究例には衝撃波領域のセンシン
グに Jameson 型センサーと組み合わせた研究例が存在する(15-17)．
また，後述する本研究では，スイッチングパラメータである β の値を物体近傍で固定する手法も提案する．本研究における非粘性
流束の偏微分項は，最終的に以下の式(13)のように求められる．

Ei 1 2, turbulent  Ei 1 2, cent   i 1 2, cent1( i  i 1 )(ui  ui 1 )


4
1

1
2
 ( i  i 1 )(ui  ui 1 )  ( pi  pi 1 )
2

Ei 1 2, cent   8
 1 (    )(u  u )(v  v ) 
i 1
i
i 1
i
i 1
 8 i

1

 ( i ei  i 1ei 1  pi  pi 1 )(ui  ui 1 ) 
4

0


u


u
1
i 1 2 , R
i 1 2 , L  i 1 2, cent   ( i  i 1 )(ui  ui 1 )  vi 1 2, R  vi 1 2, L 4 ei 1 2, R  ei 1 2, L 

 u
2

(9)
2

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
C02-2
基づいてセルを分類した．下付き文字の FC，GC，OC はそれぞれ流体セル，ゴーストセル，物体セルを指している．

Δx

VGC  VIP 

Fluid Cell

Ghost Cell

dOC
y

dGC
Object Cell
Boundaryx(15)

２．３ 空力係数算出法
本研究では，多数の物体の流体力を同時に見積もる必要がある
ため，適用したゴーストセル法を拡張して用いる．本手法では物体境界はゴーストセルと流体セルのどちらにも含まれているため，それぞれのセルで流体力を見積もる必要がある．ゴーストセルでは通常の流体計算と同様にイメージポイントを用いて流体力を見積もる．一方，流体セルで流体力を見積もる場合，セル中心と物体境界が非常に近くなることで粘性力が過大評価される場合がある．そこで図 3 のように流体セルからもプローブを伸ばしてイメ
ージポイントを定義し，粘性力を見積もる．ただし，この場合のプローブの長さは√2/2Δx とする．このようにすることで，粘性力は式(16)から求められる．ただし式(16)の速度成分 VIP’と VIB は物
体壁面に平行な方向の速度成分である．

dFC

Δx

d IP  dGC(VIP  VIB )d IP

Figure 1 Definition of level set function and cell classification
Δx
Δx

d FC  0

Fluid Cell

(14)

dGC  0 and dGC  2 2x
dOC  2 2x

Ghost Cell

本研究ではイメージポイントを使うゴーストセル法を採用するが，この手法の長所はアルゴリズムの簡便さと計算の安定性にある．図 2 にイメージポイントの模式図を示す．イメージポイン
トとはゴーストセルから物体壁面に対して垂直方向に伸びたプローブの先端にある点である．
このプローブの長さ dIP も重要なパラメータとなりうるが，本研究では格子幅の 1.75 倍（1.75Δx）で固定した．この 1.75Δx という長さは三次元への拡張を考えた場合のセルの斜辺の長さ√3Δx に基づいている．このイメージポイントを囲む流体セル群（図 2 における橙色で囲まれたセル群）からイメ
ージポイントに流れの諸量を内挿する．本研究ではこの内挿にはbilinear 内挿を用いた．このイメージポイントの値から（必要に応じてプローブ内での分布を考慮して）ゴーストセルにおける値を決定する．本研究では，速度はプローブに沿って線形に分布していると仮定して，式(15)に基づいて値を決定する．

Image Point’Δxxf IB  

VIP

(16)

２．４ 外部境界条件
超音速流れを解析する場合は外部境界条件に Neumann 条件を適用すればよいが，それだけでは亜音速流れを解析する場合には不安定化を招く．それを防ぐために本研究ではスポンジ境界(23)を適用し，音波の反射を抑制した．同時に亜音速の外部境界条件では密度を一定として，その他の変数に Neumann 条件を課した．これにより，外部境界に渦が干渉する場合も安定に解析が行える．

VIB

dGC

y
Object Cellx

VIP  VIBSd IP '  d FC

また，レベルセット関数の情報無しに表面力をより簡易な手法
で求める計算法も存在する(22)が，ここでは扱わない．

dIP
VGC

Boundary

Figure 3 Image point for estimation of surface force

Image Point

Ghost CelldIP’

Object Cell

Δx

Fluid CellΔSy２．５ フローチャート
図 4 に本解析のフローチャートを示す．時間進行のループ内に
三段階 Runge-Kutta 法の内部ループが存在する．流束計算の直前に内部境界条件として作用するゴーストセルの値を求めている．
移動物体を含む場合には，Runge-Kutta ループの前にレベルセット

Boundary

Figure 2 Schematic of image point

3

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
C02-2
関数やセル属性，イメージポイントの再定義などを毎回行う．流
体力を求める前にも一度，粘性力を正しく求めるためにゴーストセルの値を更新している．

top boundary
inflow boundary

初期条件入力
レベルセット関数・セル属性定義
イメージポイント定義

Runge-Kutta法

時間進行のループ

非粘性・粘性流束計算

O

イメージポイント定義

18

D=1

28

x

bottom boundary
Figure 5 Computational domain

レベルセット関数・セル属性定義

変数の更新

8

y

外部境界条件ゴーストセルの更新

10

1
outflow boundary

sponge layer

物体移動

ゴーストセルの更新

移動物体の場合
流体力の計算

計算終了

Figure 4 Flowchart

３．二次元円柱周りの検証計算
３．１ 計算条件
本研究で提案したスキームを，二次元円柱周りの流れ場によって MUSCL-Roe スキームと比較して検証する．比較する計算スキームは表1 の3 種類，
格子解像度は円柱直径をD として0.2D，
0.1D，
0.05D，0.025D の 4 種類，Mach 数は 0.3，1.2，2.0 の 3 種類とし，合計 36 ケースで比較を行う．Reynolds 数は全ケースにおいて 300で固定する．
Figure 6 Computational grid
Table 1 Numerical scheme for comparison(top-left: 0.200D top-right: 0.100D
bottom-left: 0.100D bottom-right: 0.025D)

A

全領域を本手法[式(11)]で解く

B

物体の隣接セルをMUSCL-Roe[式(8)]
その他の領域を式(11)で解く

C

全領域をMUSCL-Roe[式(8)]で解く

３．２ マッハ数 0.3 の円柱周りの比較
スキーム間における流れ場の違いを見るために，図 7 に格子解
像度が 0.200D の同時刻の密度分布を示す．図から C のMUSCL-Roe スキームの渦が拡散しているのに対し，A の提案のスキームでは渦が高精細に捉えられているのが確認できる．一方で図 8 から物体近傍での圧力場を見ると，A では渦の保存特性が良好な一方で中心差分の影響から物体付近で若干の圧力振動が確認される．対して C では滑らかな圧力分布が見られるが同時に渦も散逸的な様相を示している．B ではこれら両者の特性が同時に確認でき，物体近傍では圧力振動を抑制しつつも，後流の渦の解像度を担保している．

計算領域と計算格子を図 5 と図 6 に示す．流入境界では全ての諸量に Dirichlet 条件，上下境界と流出境界の諸量に Neumann
条件を課す．ただし，亜音速流れである Mach 0.3 の場合の流出境
界の密度のみは一様流の値で固定する．スポンジ層は計算領域の外部境界近傍に配置する．計算格子の図において，黒線がセル中
心をつないだ格子線，赤太線が物体境界，青色，白抜き，赤色の領域がそれぞれ流体セル，ゴーストセル，物体セルである．また，計算結果の妥当性を検討するため，同様の計算条件で境界適合格
子の計算(24-25)から得られた結果と比較する．その計算格子は(imax,jmax)=(208, 177)の計算格子で，およそ 2 倍の解像度の計算格子を用いて格子収束を確認済みである．

4

Copyright © 2013 by JSFM

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

scheme-A Δx=0.200D

scheme-A Δx=0.050D

scheme-B Δx=0.200D

scheme-B Δx=0.050D

scheme-C Δx=0.200D

scheme-C Δx=0.050D

Figure 7 Instantaneous density distribution at Δx=0.200D

Figure 8 Instantaneous pressure distribution at Δx=0.050D(top: scheme-A middle: scheme-B bottom: scheme-C)(top: scheme-A middle: scheme-B bottom: scheme-C)

A と B の違いを明確にするために，式(12)で示されるスキームのスイッチングパラメータ β を図 9 に可視化する．図における黒と白の領域がそれぞれβ=0 と1 を表し，
式(11)より黒の領域が式(7)の skew-symmetric 型スキーム，白の領域が式(9)の MUSCL 型スキームが適用されていることを示す．
格子解像度が粗い 0.100D の場
合は B の方に広い MUSCL 領域が見られるが，格子解像度が細かな 0.025D の場合には同様の分布を呈している．
提案手法は衝撃波とポテンシャル流の領域のみに MUSCL スキームを適用することを目的としており，格子解像度が細かい場合にはこの目的が満足されている．一方で格子解像度が粗い場合には，物体付近の圧力振動から生じる数値誤差の影響により式(12)のセンサがうまく機能せず，MUSCL 領域が狭くなったと考えられる．この時に B の方が正確にポテンシャル領域を探知できているのは，物体境界付
近に MUSCL-Roe スキームを用いることで圧力振動を抑制し，数値誤差の伝搬を抑えたことによると考えられる．

Mach 0.3 の場合の抗力係数，圧力抗力係数，摩擦抗力係数を図10 に示す．赤線が境界適合格子によって得られた参考値である．
圧力抗力には格子収束が見られつつあるが摩擦抗力はまだ収束していない．Re 数 300 における境界層厚さは 1/√Re からおよそ
0.058D 程度と見積もられ，最も細かい格子解像度の 0.025D であっても 2 点が境界層に入るかどうかで摩擦抗力の格子収束には解像度不足であると思われる．また，A における圧力の過大評価が著しいが B においてかなりその傾向が緩和されており，C の値とも比較して B のスキームが良好な性能を示していると言える．

5

Copyright © 2013 by JSFM

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

scheme-A Δx=0.100D

A

B

C

BFC

average drag coefficient

2.0
1.8

1.6
1.4
1.2
1.0
0.01

0.10

1.00

mesh size

scheme-B Δx=0.100D
A

B

C

BFC

average pressure drag coefficient

2.0
1.8
1.6

1.4
1.2
1.0
0.01

scheme-A Δx=0.025D

0.10

1.00

mesh size

A

B

C

BFC

average friction drag coefficient

0.3

scheme-B Δx=0.025D

0.2

0.1

0.0
0.01

0.10
mesh size

1.00

Figure 10 Drag coefficient in Mach 0.3(top: total middle: pressure bottom: friction)

Figure 9 Instantaneous distribution of β of Eq. (11)6

Copyright © 2013 by JSFM

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

scheme-A Δx=0.100D

scheme-A Δx=0.100D

scheme-B Δx=0.100D

scheme-B Δx=0.100D

scheme-C Δx=0.100D

Figure 12 Instantaneous distribution of β of Eq. (11)

A

B

C

BFC

drag coefficient

1.8

1.6

1.4

1.2
0.01

Figure 11 Converged density distribution at Mach 1.2

0.10

1.00

mesh size(top: scheme-A middle: scheme-B bottom: scheme-C)A

３．３ マッハ数 1.2 の円柱周りの比較
スキーム間における流れ場の違いを見るために，
図 11 に格子解
像度が 0.100D の物体付近の密度分布を示す．
Mach 0.3 の場合と同様に，A，B，C の順で散逸的な流れ場になっていることが確認できる．またスイッチングパラメータ β を図 12 のように Mach 0.3の場合と同様に可視化すると，超音速のため MUSCL の領域が広がっている．だが A と B の分布において Mach 0.3 ほどの違いは見られない．物体の近傍にのみ若干ながら異なる分布が見られる
が，ここでは紙面の都合上から図は割愛する．その結果として図13 のように圧力抗力に差が生じて抗力係数に差が生じているものの，その差は Mach 0.3 ほど大きくはない．圧力抗力の収束値が境界適合格子の結果と若干ずれているが，この原因は調査中であ
る．

B

C

BFC

pressure drag coefficient

1.8

1.6

1.4

1.2

0.01

0.10
mesh size

1.00

Figure 13 Drag coefficient in Mach 1.2(top: total bottom: pressure)

7

Copyright © 2013 by JSFM

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

scheme-B Δx=0.100D

A

B

C

BFC

drag coefficient

1.6

1.4

1.2

0.01

0.10

1.00

mesh size

scheme-C Δx=0.100D

A

B

C

BFC

pressure drag coefficient

1.6

1.4

1.2
0.01

0.10

1.00

mesh size

Figure 14 Converged density distribution at Mach 2.0
A(top: scheme-B bottom: scheme-C)

B

C

BFC

friction drag coefficient

0.3

scheme-B Δx=0.100D

0.2

0.1

0.0
0.01

0.10
mesh size

1.00

Figure 16 Drag coefficient in Mach 2.0(top: total middle: pressure bottom: friction)

３．４ マッハ数 2.0 の円柱周りの比較
スキーム間における流れ場の違いを見るために，
図 14 に格子解
像度が 0.100D の物体付近の密度分布を示す．
Mach 2.0 の本ケースでは A の計算が途中で発散したため，B と C のみで比較を行う．
他の場合と同様にBよりもCの方が散逸的な流れ場になっている．
しかしながら抗力係数で見ると図 16 のように B と C の間に大きな差は見られない．一方でスイッチングパラメータ β を図 15 のように異なる格子解像度間で比較すると，格子解像度の向上に伴って明瞭にスイッチングが行われており，今後多数の物体が移動する問題を解析する場合には大きく影響することが予想される．

scheme-B Δx=0.025D

Figure 15 Instantaneous distribution of β of Eq. (11)8

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
C02-2
４．超音速ですれ違う二次元円柱周りの解析
４．１ 計算条件
計算領域と計算条件は図 17 のように設定し，二つの円柱 C1 とC2 は同時にスタートして同じ速度 U1=Mach 1.2 で移動しながらすれ違う．
この速度と円柱直径から求められる Re 数を 300 とする．
初期場では流れが無いとし，四方を周期境界条件と設定する．格子解像度には 0.050D と 0.025D の二通りを用い，計算スキームには表 1 の B を採用する．またこの計算ではセルの属性が物体・流体・ゴーストの三種類に動的に変化し，一般的にフレッシュセル(23)と呼ばれるセルの切り替わりの問題が生じる可能性がある．それを効率的に扱うために様々な対処法が提案されているが，ここでは特別な処理を行わずに計算を行う．

y
O

10
1

C2

3
U1

U1

C1
10

20

Periodic boundary

Periodic boundary

Periodic boundary

20

x Periodic boundary
Figure 17 Computational setup

４．２ 計算結果
無次元時間tが7.8，
12.9，
18.1，
22.0.のときの格子解像度が0.025Dの密度場の可視化図を図 18 に示す．t=7.8 では円柱前方の衝撃波同士の干渉が見られ，t=12.9 では円柱 C1 が円柱 C2 の後流に干渉している．その後 t=18.1 では後流同士が干渉し，t=22.0 では後流に渦が形成されているのが確認できる．t=22.0 のとき β は図 19 のように分布しており，円柱後流部にのみ skew-symmetric 型スキー
ムが適用されていることが確認できる．各軸方向の力を移動速度から求められる基準動圧と円柱直径で正規化した軸力係数の時間変化を図 20 に示している．t=0 で円柱が動き始めてから t=20 付近
で円柱が計算領域外に出て値がゼロになっている．t=10 付近にあるC2 のx方向の軸力係数のピーク値には格子解像度の差によって約 10%程度の差があるが，他の時間における格子解像度による差はさほど大きくない．このようにして，本手法を用いれば多数物体が運動する流れ場において効率的に流体力を求めることが可能である．
本計算において，セル属性がある時急に流体セルに切り替わる
フレッシュセルは間違いなく生じているが，本手法の場合には物体内部の 2 層分をゴーストセルとすることで，物体セルから流体セルに切り替わることなく，必ずゴーストセルから流体セルに切り替わっている．そのためフレッシュセルの影響が緩和され，数値誤差による圧力振動等は確認されずに安定に計算がなされたと考えられる．

Figure 18 Distribution of instantaneous density distribution(t=7.8, 12.9, 18.1, 22.0)

9

Copyright © 2013 by JSFM

第 27 回数値流体力学シンポジウム
C02-2
謝辞
本研究は科研費（24656522）の助成を受けたものである．また，
本計算は東北大学サイバーサイエンスセンターNEC SX-9 にて行った．ここに感謝の意を記す．
参考文献(1) K. M. Eldred, NASA SP-8072, 1972.
(2) K. Fujii, T. Nonomura, S. Tsutsumi, International Journal of
Numerical Methods in Fluids, Vol. 64, No. 10-12, pp. 1412-1432,
2010.
(3) T. Nonomura, K. Fujii, CFD Review 2010, pp. 421-446, World
Scientific Company, 2010.
(4) T. Nonomura, K. Fujii, AIAA Journal, Vol. 49, pp. 2292-2294,
2011.
(5) T. Nonomura, Y. Goto, K. Fujii, International Journal of
Aeroacoustics, Vol. 10, No. 4, pp. 401-426, 2011.
(6) K. Fukuda, S. Tsutsumi, K. Ui, T. Ishii, R. Takaki, K. Fujii,
Transaction of JSASS, Vol. 54, No. 184, pp. 120-129, 2011.
(7) R. Mittal, G. Iaccarino, Annual Review of Fluid Mechanics, Vol. 37,
pp. 239-261, 2005.
(8) S. K. Sambasivan, H. S. Udaykumar, AIAA Journal, Vol. 47, pp.
2923-2937, 2009.
(9) X. Y. Hu, B. C. Khoo, N. A. Adams, F. L. Huang, Journal ofComputational Journal, Vol. 219, No. 2, pp. 553-578, 2006.
(10) 藤井孝蔵：流体力学の数値計算法，東京大学出版会，1994(11) L. Georgesa, G. Winckelmans, P. Geuzainea, Journal ofComputational and Applied Mathematics, Vol. 215, pp. 419-428,
2008.
(12) P. L. Roe, Journal of Computational Physics, Vol. 43, No. 2, pp.
357-372, 1981.
(13) B. van Leer, Journal of Computational Physics, Vol. 23, No. 3, pp.
276-299, 1977.
(14) G. van Albada, B. van Leer, W. Roberts, Astronautics andAstrophysics, Vol. 108, pp. 76–84, 1982.
(15) F. Ducros, V. Ferrand, F. Nicoud, C. Weber, D. Darracq, C.
Gacherieu, T. Poinsot, Journal of Computational Physics, Vol. 152,
pp. 517-549, 1999.
(16) A. Jameson, W. Schmidt, E. Turkel, AIAA Paper, 81-1259, 1981.
(17) S. Teramoto, AIAA Journal, Vol. 43, No, 11, pp. 2354-2363, 2005.
(18) S. Gottlieb, C. W. Shu, Mathematics of Computation, Vol. 67, No.
221, pp. 73-85, 1998.
(19) M. T. McGurn, K. P. Ruggirello, P. E. DesJardin, Journal ofComputational Physics, Vol. 241, pp. 364-387, 2013.
(20) K. Nakahashi, AIAA Paper, 2011-3386, 2011.
(21) S. Takahashi, I. Monjugawa, K. Nakahashi, Transactions of JSASS,
Vol. 51, No. 172, 2008.
(22) J. Onishi, T. Nonomura, “Notes on the Simple Evaluation of theForce on Bodies in Immersed Boundary Methods for Fluid
Computation,” (in preparation.)(23) J. B. Freund, AIAA Journal, Vol. 35, No. 4, 1997.
(24) X. Y. Hu, Q. Wang, N. Adams, Journal of Computational Physics,
Vol. 229, pp. 8952-8965, 2010.
(25) T. Nonomura, D. Terakado, Y. Abe, K. Fujii, AIAA paper,
2013-2569, 2013, (submitted for publication.)(26) J. H. Seo, R. Mittal, Journal of Computational Physics, Vol. 230, pp.
7347-7363, 2011.

Figure 19 Instantaneous distribution of β of Eq. (11) at t=22.0

Figure 20 Temporal variation of force coefficient

５．結言
本研究では圧縮性 Navier-Stokes 式に対して，等間隔直交格子法とレベルセット法，埋め込み境界法を用い，skew-symmetric 型とMUSCL-Roe 型のハイブリッドスキームを適用しての直接数値計算を行ってその性能検討と応用性を調査した．二次元円柱周りでMach 0.3，1.2，2.0，Re 数を 300 と固定して，4 つの異なる格子解像度によって検証を行った．その結果，物体周りを強制的にMUSCL-Roe，その他を提案スキームとする場合が物体境界付近の圧力振動を抑制して適切な解を導くことが確認された．
また，2 つの二次元円柱が相対運動する解析により，本解析の定性的，定量的な有効性が示された．今後はこのスキームを三次元に拡張し，多数の粒子が移動する複雑固気混相流解析を実施する．
10

Copyright © 2013 by JSFM