Thermal Science & Engineering Vol.14 No.3 (2006)

粒子法とグリット法の結合によるマルチスケール気液界面解析
石井 英二†

石川 亨‡∗

田辺 好之‡

Simulation of Multi-Scale Free Surfaces within Gas-Liquid flows
using Combined Particle and Grid Methods
Eiji ISHII,† Toru ISHIKAWA‡ andYoshiyuki TANABE‡

Abstract
Gas-liquid flows usually include multi-scale free surfaces, for example, fuel sprays used for automobile engines
become liquid films at the outlet of fuel injectors, and then the liquid films break up into droplets. Fluid flow with freesurfaces, like the liquid films, are mainly simulated using grid methods, where the front of the free surface is directly
captured on regular, fixed grids that cover both liquid and gas domains. In these methods, however, free surfaces are
sometimes lost due to numerical diffusion in case the scale of the free surfaces becomes smaller than the grid size.
Particle methods can avoid the occurrence of such the unreal loss of free surfaces. The particle methods treat the freesurfaces as groups of particles that move in a pattern based on a Lagrangian description. To simulate the multi-scale freesurfaces, we have developed a hybrid particle/grid method where the free surfaces within a sub-grid region are
simulated with particles located near the liquid interfaces. Velocities determined using the particle method are combinedwith those obtained from the grid method. However, the interaction between gas and liquid still depends on the grid size
because the gas regions within the sub-grid regions were calculated by the grid method. Accordingly, the present
method adopts a two-particle model as the particle method; that is, the free surfaces within the sub-grid region are
simulated with two types of particles, one for the gas and the other for the liquid. The sub-grid regions are determinedby using the volume fraction of liquid calculated using the grid method, and then the two-particle model is applied tothe sub-grid regions. We will verify the new method by applying it to prediction of large deformation in free surfaces(the Rayleigh-Taylor instability) and to the fragmentation of a water column. The predicted deformed shape of the watercolumn shows good agreement with measurements reported by Koshizuka et al.
Key Words: Computational fluid dynamics, Multi-phase flow, Multi-scale, Hybrid method

記
C
d
g
n
p
r
R
t
u

w
W

[m/s ]

μ
θ
ρ

[Pa]
[m]
[m]
[s]
[m/s]

τ

:
:
:
:
:
:

添字
CIP
MPS

:
:

号
:
:
:
:
:
:
:
:
:

音速
空間の次元数
重力加速度
粒子数密度
圧力
粒子間距離
補間領域の半径
時間
速度

[m/s]
2

MPS 法の重み関数補間で用いる重み関数
粘性係数
液相の体積分率
密度
粘性応力
CIP 法
MPS 法

[Pa⋅s][kg/m3][N/m2]∗ 受理日: 2006 年 3 月 23 日, 担当エディター: 黒田 明慈†
株式会社 日立製作所 機械研究所 (〒312-0034 茨城県ひたちなか市堀口 832 番地 2)‡
株式会社 日立製作所 オートモーティブシステムグループ (〒312-8503 茨城県ひたちなか市高場 2520 番地)- 51 -© 2006 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.14 No.3 (2006)

１ 緒 言

我々は以前に，マルチスケールの気液界面挙動を解析するために，界面追跡法の CIP 法と粒子法の工業的に利用される気液二相流の多くは，ミクロ

MPS 法を結合した解析手法を開発した[8]．ここではとマクロの気液界面が混在するマルチスケールの流界面近傍のみに粒子法を適用し，他の領域を CIP 法

動現象である．例えば自動車エンジン内に燃料を噴で解析することで，計算負荷の小さい効率的な解析

霧する自動車用インジェクタでは，噴孔出口に液膜手法を構築した．また CIP 法が工業製品の複雑形状

状のマクロ界面が形成され，それが界面不安定性やに対応できるように，非構造格子に対応できる形式

周囲ガス流動による巻き込みにより液滴状のミクロに改良を行った[9]．これにより自動車用スワールイ

界面へと形態を変える．このような燃料噴霧の分裂

ンジェクタの燃料噴霧の液膜がリング状に分裂する

挙動はエンジン内の燃焼に影響を及ぼすことがわか

挙動を定性的に再現することが可能となった[10]．

っている．そして近年の自動車エンジンの低排気しかし本手法は MPS 法の計算において，（気相と液化・低燃費化の流れに伴い，燃料噴霧の分裂挙動解

相の両方ではなく）液相のみに粒子法を適用したた

析の必要性が増加してきている．
そこで本研究では，め，気液間の相互作用を十分に計算できない場合が

燃料噴霧の分裂挙動解析技術構築の第一ステップとあった．例えば Rayleigh-Taylor 不安定問題では液相して，マルチスケール界面挙動の解析技術についての沈降挙動を正確に再現できないケースがあった．

検討を行った．

そこで本研究では界面近傍に気相用と液相用の二種

従来，気液界面挙動の計算では界面追跡法のよう類の粒子を配置することで，気液間の相互作用を高に固定された計算格子により界面移動を計算する手精度に計算することを試みた．そして解析手法の検法が主に用いられてきた．界面追跡法としては証の第一ステップとして，本解析手法（これよりVolume of fluid (VOF)法[1]，Level set 法[2]，Cubic

Hybrid 法と呼ぶ）を基礎的な流れ場（Rayleigh-Taylor

interpolated propagation (CIP)法[3]等があり多くの成不安定問題とダム崩壊問題）に適用した結果につい果を挙げている．しかしこれらの手法では，界面サて述べる．

イズが計算格子サイズに近くなると，数値拡散によ
２ 解析手法

り界面が消失してしまう問題がある．一方，数値拡散による界面消失のない解析手法として，Moving
particle semi-implicit (MPS)法[4] や Smooth particle

図 1 に Hybrid 法のコンセプトを示す．計算格子よhydrodynamic 法[5]のような粒子法がある．これらのりも大きな界面挙動は主に CIP 法で計算し，計算格手法では，計算格子の代わりに空間に配置した粒子

子よりも小さな界面挙動は主に MPS 法で計算を行を Lagrange 的に移動するので，数値拡散による界面う．

消失が無いのが特徴である．しかし自動車用インジ
ェクタの内部から大気領域までを粒子法で計算する
ためには，粒子数が膨大となり計算負荷が大きくな

Grid

る問題がある．

Gas

CIP method

Micro surface

MPS method

上記の界面追跡法と粒子法の問題に対して，界面追跡法と粒子法を結合し，計算の効率化を図った手法がある．Liu ら[6]は MPS 法に差分格子を併用した．
また Enright ら[7]の Hybrid particle level set 法では，気液界面に散布したマーカー粒子を Lagrange 的に用いることで，数値拡散による界面の消失が少ない
高精度な解析を実現している．これらの解析手法は気相と液相（もしくは固相）のどちらかを粒子法でMacro surface

CIP method

置き換えるか，または粒子自体はマーカーとしての

LiquidCIP method

み使用するため，粒子サイズ（粒子間距離）は計算
Fig. 1 Assignment of CIP and MPS methods in

格子サイズと同程度であり，マルチスケールの界面

gas-liquid flowsに対応するために十分とはいえない．
- 52 -© 2006 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.14 No.3 (2006)

MPS 法の粒子は界面近傍のみに配置される．界面近傍以外の領域は計算格子を利用して CIP 法で計算する．以下では MPS 法と CIP 法について簡単に説明した後，MPS 法と CIP 法の結合方法について説明すGGG
uin +1 = ui* + δ ui

(9)

G
G
G
ri n +1 = ri* + Δ t δ ui

(10)

る．

ただし，

2.1 Moving Particle Semi-implicit (MPS)法[4](∇p )i = d ∑ ⎢⎢⎡n0 j ≠ ir −r⎣⎢ j i
ここで記号*は粒子移動後の値を示す．

非圧縮性と運動量の式は次式で表される．
∂ρ
=0
∂t
G
G
Du
ρ
= −∇p + μ∇ 2 u
Dt(1)示す．表面張力は無視した．式(2)左辺は Lagrange 的
ターは次式で計算する．
j0 j ≠i

i

j

i

(12)

0
⎞
⎛
⎟
⎜
⎛F⎞
⎜ ⎟
− ρ u j, j
⎟
⎜
ρ
⎜ ⎟
⎟
1
(13)
f =⎜ ⎟ g =⎜ 1
⎜ − p ,i + τ ij , j ⎟
ui
⎜ ⎟
ρ
ρ
⎟
⎜
⎜ p⎟
⎝ ⎠
⎜ − ρ C s2 u j , j ⎟⎠⎝τ は粘性応力，Cs は音速を表す．表面張力は無視しな粒子の扱いを示す．式(2)のラプラシアンのオペレ

i⎦⎥

流れの支配方程式は次式で表される．
∂f+ u j f, j = g∂t

ここで u は速度，p は圧力，ρ は密度，μ は粘度を

2

(11)

2.2 Cubic Interpolated Propagation (CIP)法[3](2)(∇ φ ) = λ2nd ∑ (φ − φ )w( r − r )⎤( r j − ri )w( rj − ri )⎥⎥2

p j − pi(3)

ここで φi は物理変数，ri は粒子座標，d は次元数，n0 は基準となる粒子数密度を示す．w は重み関数でた．F は液相の体積分率 θCIP を用いて次式で定義さあり次式で定義される．

れる．

(rij ≤ re )(rij > re )

⎧ re
⎪ −1
w(rij ) = ⎨ rij⎪0⎩θ CIP =
(4)

音速は次式で求める．

示す．粒子数密度 ni と変数 λ は次式で定義する．

(

j ≠i

)

(5)

∑ r j − ri w( r j − ri )
2

λ=

j ≠i

∑ w( r j − ri )(6)

μ = μlθCIP + μ g (1 − θCIP )(16)

Cs = CslθCIP + Csg (1 − θ CIP )(17)式に分割して計算される．
∂f+ u j f, j = 0∂t

的に移動させる．次に式(1)の非圧縮条件を満足する
ように次の圧力 Poisson 方程式を解く．

i

(15)

す．式(12)は次式のように移流方程式と非移流方程まず式(2)の粘性力と体積力を求め，粒子を Lagrange(∇ p ) =ρ = ρlθ CIP + ρ g (1 − θCIP )ここで添え字の l と g は液相と気相での物性値を示

j ≠i

ρ ni − n 0

2

(14)

θCIP は気相で 0，液相で 1 の値をとる．密度，粘度，

ここで rij は粒子間距離，re は粒子間の影響範囲を
ni = ∑ w r j − riarctan(F )
+ 0.5
0.99π

Δt

2

n0

∂f
=g
∂t

(7)

(18)
(19)

CIP 法では式(18)を解く際，3 次スプライン補間を用式(7)で得られた圧力を用いて，速度と粒子座標は次いることで，3 次精度で移流方程式を解くことが可式で修正される．
GΔt
δ ui = − ∇p

能である．式(18)を計算した後に，式(19)を有限体積ρ法で離散化し，fractional-step 法で時間積分した．

(8)- 53 -© 2006 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.14 No.3 (2006)

2.3 MPS/CIP hybrid 法

体 2 の界面からのメッシュ数（距離）で定義し，矢MPS 法と CIP 法の結合方法を以下に示す．

印方向を正とした．また基準位置は CIP 法で求める密度関数 θ （流体 1 で 0，流体 2 で 1）が 0.5 の位① MPS 法と CIP 法で別々に移流項を計算して移流置で，δ (i) をゼロとした．粒子の移流速度は，界面速度を求める．
② MPS 法で使用する粒子の粒子濃度と，CIP 法で近傍（δ が 0～δmin ）では粒子法で計算し，他の領求めた気液界面からの距離を変数とする重み関域（δ が δmax 以上）では CIP 法で計算した．ただし，数で，各々の移流速度を修正する．

δ が δmin ～ δmax の領域は MPS 法と CIP 法の結合領域とした．また粒子の散布範囲は，δ が 0～3（界面から 3 メッシュの範囲）とし，δmin は 0.4，δmax は③ 最後に圧力ポアソン式を用いて（MPS 法と CIP
法で別々に）圧力を計算する．

MPS 法と CIP 法の通常の解析手順に，上記②の移流

0.5 として界面近傍の約 1 メッシュ分を MPS 法の計

速度の修正を加えたのみのシンプルな解析手法であ算に割当てた．CIP 法での界面消失を MPS 法で補うる．移流速度の修正は次式で行う．
G*
G
G
u MPS
= (1 − λ MPS ) u MPS + λ MPS u (CIP )という目的から δmax は 0.5 から 1.0 程度が望ましく，

(20)

G*
G
G
u CIP
= (1 − λCIP ) u ( MPS ) + λCIP u CIP(21)

これより大きいと界面挙動が平均化されて鈍った傾向を示す．また粒子濃度 β は次式で定義した．

β = ∑ W ( Rij ) N 0

(28)

i≠ j

ここで uMPS は MPS 法で計算した（粒子座標上の）

ここで N0 は基準の粒子数密度である．

移流速度，uCIP は CIP 法で計算した（グリッド座標上の）移流速度，u(CIP) は uCIP の粒子座標上への補間値，u(MPS) は uMPS のグリッド座標上への補間値を(a) Particles → Grids．
表す．それぞれの補間は次式で行う（図 2 参照）
G∑ u MPS W ( Rij )Gi≠ ju ( MPS ) =(22)∑ W ( Rij )(b) Grids → Particles

Re

Re

i≠ j

G

∑ u CIP W ( Rij )Gi≠ ju (CIP ) =(23)∑ W ( Rij )Fig. 2 Interpolation between particle and gridi≠ jcoordinates

⎧ Re
−1
⎪
ただし， W ( Rij ) = ⎨ Rij⎪0⎩(Rij ≤ Re )(Rij > Re )(24)

Fluid 1

ここで Re は有効半径（1.0～1.5 格子分）
，Rij は粒子

(1)
δ max

座標とメッシュ座標の距離．一方，重み関数は次式(i )(δ (i ) < δ min)

ただし， β = θ CIP βθ = 0.5

Fluid 2

(25)
(26)
(1)

+ (1 − θ CIP ) β( 2)

δ (1) = 0
δ ( 2) = 0

(i )
(i )
(δ min≤ δ (i ) ≤ δ max
)
(i )
(i )
(δ > δ max )λCIP = λ MPS (1 − β )

CIP

(1)
δ min

となる．
0
⎧
⎪ (i )
(i )
(i )
(i )
)
λMPS = ⎨(δ − δ min ) /(δ max− δ min⎪1⎩

Grids δ (1)

( 2)
δ min
( 2)
δ max

δ ( 2)

MPS

CIP

(27)
Fig. 3 Computational model

図 3 は解析モデルであり，δ (i) (i =1, 2)は流体 1 と流

- 54 -© 2006 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.14 No.3 (2006)

３ 計算結果3.1

Rayleigh-Taylor 不安定問題

解析手法の定性的妥当性を確認するため，Hybrid

Heavy fluid

法 を Rayleigh-Taylor 不 安 定 問 題 に 適 用 し た ．

Rayleigh-Taylor 不安定とは，軽い流体の上に重い流体が存在する場合に起こる界面の不安定現象である．

1
Particles for heavy fluid

計算の初期条件を図 4 に示す．軽い流体を領域の下側 2/3 に，重い流体を上側 1/3 に配置し，界面の初Particles for light fluid

期振幅を 0.1 とした．軽い流体と重い流体の密度比は 1：4，動粘度は共に 0.01，重力加速度は 10 とした．計算領域の幅は 1，高さを 3 とし，外周は壁面とした．計算格子は 60×180（格子間隔は 0.017），

Light fluid2

初期粒子間距離は 0.005 として界面近傍に MPS 法用の粒子を発生させた．重い流体用の粒子は，CIP 法で計算した液相率が 0.5 の位置から重い流体側に 3
メッシュの範囲で発生させた．この 3 メッシュより外の重い流体領域に進入した重い流体用の粒子は削1除した．軽い流体用の粒子は軽い流体側 3 メッシュの範囲で発生させた．また軽い流体中で重い流体用の粒子の濃度が上昇した場合には，その周囲に軽いFig. 4 Initial arrangement of light and heavy fluids

流体用の粒子を発生させた．

(a) Hybrid method(b) CIP method only(c) CIP method only

60×180 grids

t = 0.8

t = 1.1

t = 1.5

t = 0.8

120×360 grids

t = 1.1

t = 1.5

t = 1.44

Fig. 5 Distribution of light and heavy fluid: (a) Hybrid method: (b) CIP method (60×180 grids): (c) CIP method (120×360 grids)- 55 -© 2006 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.14 No.3 (2006)

図 5(a) に Hybrid 法の計算結果を表す．また図(a)5(b)(c)は CIP 法を二種類の格子数で適用した結果で

(b)

(c)

あり， Hybrid 法の計算結果の定性的検証，およびHybrid 法の粒子挙動により補われる界面について検討した．図 5 (c)（t = 1.44）は液相沈降位置が一致する時刻（(b)の t = 1.5）での結果を示す．図 5 よりHybrid 法では重い流体が軽い流体に沈降して拡散する過程を定性的に再現できていることがわかる．また Hybrid 法の計算結果を，CIP 法単独での計算結果と比較すると，数値拡散に伴う界面の消失が，Hybrid
法の粒子挙動により補われていることがわかる．
図 6 は計算結果の粒子数依存性を示したものである（重い流体用粒子の表示は省略している）．初期粒
子間距離はそれぞれ，(a) 0.01 (b) 0.0067 (c) 0.005 である．粒子間距離が短い（粒子数が多い）と，分布の対称性が改善される傾向がある．この理由は粒子Fig. 6 Particle/number-of-density sensitivity study.

間距離が短いと，粒子発生時に界面形状（CIP 法かInitialら求まる液相率 0.5 の位置）に沿って粒子を配置しchanged: (a) 0.01, (b) 0.0067 and (c) 0.005

易くなり計算が安定化するためと考えられる．

(presently adopted)

distance

betweenparticles

were

3.2 ダム崩壊問題
複雑な気液界面挙動の予測精度を検証するため，

Hybrid 法をダム崩壊問題で検証した．図 7 に解析モaデルを示す．計算条件は越塚ら[4]の実験条件と同一

Water column，高さ 29.2 cm
として，横 14.6 cm（図 7 中の記号 a）

8a（2a）の水柱が幅 58.4 cm（4a）の容器内で崩壊する過程を計算した．計算では容器の高さは 116.8 cm

2a（8a）とした．空気密度は 1.2 kg/m3，水密度は 998

Wall

kg/m3，空気粘度は 1.8×10-5 Pas，水粘度は 1.0×10-3
Pas，重力加速度は 9.8 m/s とし，水の表面張力は無4a視した．また外周は壁面境界とした．計算格子は 96
×192（格子間隔は 0.608 cm），初期粒子間距離はFig. 7 Initial geometry for the simulation of collapse of

0.162 cm とした．粒子の発生範囲は Rayleigh-Taylor

water column

不安定問題の場合と同様の条件で設定した．
図 8 は水柱崩壊過程の計算結果を示したものである．水柱は，時刻 0.1 秒から 0.2 秒では重力により台

図 9 は格子解像度の異なる二種類の CIP 法の計算

形状に変形しながら右側へ進行する．その後，時刻結果を比較して，Hybrid 解析における MPS 法の利

0.3 秒付近で右壁面に衝突し，水膜が壁面を上昇し始用効果を評価した．図 9 (a)は越塚と岡[4]による実験

める．水膜は 0.5 秒付近を境に重力により下降を始

結果，(b)は 96×192 の計算格子での CIP 法の計算結める．そして時刻 0.7 から 0.9 秒では左側へ移動する

果，(c)は 192×384 の計算格子での CIP 法の計算結ようになる．これらの計算結果より，Hybrid 法を使果，(d)は Hybrid 法による計算結果（96×192 の計算用したことにより，液相の気相への飛散および気相格子を使用）である．図 9(b)と(c)の CIP 法の結果での液相への巻き込みが，界面付近に配置された粒子

は，格子解像度を 2 倍にすることで，0.6 秒と 0.8 秒の挙動により表現されていることがわかる．

では水柱形状が改善される傾向が見られる．一方，
図 9(d)の Hybrid 法の結果は(b)と同じ 96×192 の計算- 56 -© 2006 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.14 No.3 (2006)

格子使用しているにも関わらず，2 倍の格子解像度

また Hybrid 法に関しては，初期粒子間距離が 0.162の結果に近い計算結果が得られていることがわかる．

（図中の Hybrid (1)）と 0.122 cm（図中の Hybrid (2)）
cm

また計算時間に関しては，Hybrid 法は同じ格子解像の 2 ケースについても結果が示されている．図は無の CIP 法の約 1.5 倍の計算時間であり，度（96×192）次元数で表示し，t は時間，g は重力加速度を示す．

計算負荷が大幅に増加しないことがわかった．

Hybrid 法の結果は検討した二種類の粒
図 10(b)より，

図 10 は Hybrid 解析の計算結果を定量的に評価す子間隔では同一の結果であった．また VOF 法の結るために，水柱右端位置の時間変化を Martin と果と比較するとほぼ一致することがわかった．ただMoyce [11]の実験結果と比較したものである．(a)はし無時限時間で 2.0 付近から僅かに差が発生する傾時刻 0.24 秒での水柱形状を，(b)は水柱右端の座標を向がみられた．この原因は，図 10(a)に示すように，示したものである．
図 10(b)では Martin と Moyce [11]Hybrid 法では図中の丸印の右側に薄く染み出すようの実験結果に加えて，VOF 法による計算結果[1]も比な部分が発生するためである．

較のために示されている．

Particles for gas
Particles for liquid
Water column

0.1 sec

0.2 sec

0.3 sec

0.5 sec

0.7 sec

0.9 sec

Fig. 8 Simulation of collapse of a water column using 96×192 grids with grid size of 0.608 cm. Initial distance
between particles is 0.162 cm- 57 -© 2006 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.14 No.3 (2006)(a) Koshizuka and Oka [4](b) CIP only (96×192 grids)(c) CIP only (192×384 grids)(d) Hybrid (96×192 grids)

0.4 sec

0.6 sec

0.8 sec

Fig. 9 Grid/number-of-density sensitivity study through the estimation of recovered scale with MPS by comparing
cases using CIP: (a) measurements of Koshizuka and Oka [4], (b) CIP only (96×192 grids), (c) CIP only (192×384 grids) and (d) Hybrid (96×192 grids)

図 10(a)において丸印が付いているところの座標が解析手法検証の第一ステップとして，

VOF 法の計算結果と一致する座標である．よってこRayleigh-Taylor 不安定問題とダム崩壊問題に適用しの薄い染み出しを取り除くことが今後の課題といえて以下の結論を得た．

る．また VOF 法の計算結果が Martin と Moyce [11]の実験結果から外れる点に関しては，容器底面の濡・ Rayleigh-Taylor 不安定問題に適用して定性的にれ性などにより実験では水柱右端先端が丸まるため

妥当な結果が得られた．また CIP 法単独での計と考えられる．

算結果と比較して数値拡散による界面の消失が低減されることがわかった．

４ まとめ・ ダム崩壊問題での複雑な界面挙動を解析し，定マルチスケールの気液界面挙動を解析するため，性的に妥当な結果が得られた．また Hybrid 化し界面近傍を気相用と液相用の二種類の粒子を配置しても計算負荷が大幅に増加しない（1.5 倍程度でた粒子法（MPS 法）により計算し，他の領域は計算ある）ことがわかった．

格子を用いる界面追跡法（CIP 法）で計算する解析手法を検討した．

- 58 -© 2006 The Heat Transfer Society of Japan

Thermal Science & Engineering Vol.14 No.3 (2006)

201-225.

(a) Collapse of water column using hybrid method.

[2] Sussman, M. Smereka, P., and Osher, S., “A levelset

approachforincompressible

computingtwo-phase

solutions

to

Journalofflow”,

Computational Physics, 114 (1994), 146-159.

[3] Yabe, T. and Aoki, T., 1996, “A dream to solve
dynamics of all materials together”, Proceedings ofInternational Conference on High-performance

Leading edge

Computing(b) Leading edge position.

in

Automotive

Design,(1996),

2105-2108.

4.0[4] Koshizuka, S. and Oka, Y., “Moving-particle
Martin and Moyce
VOF
Hybrid (1)Hybrid (2)3.5

semi-implicit

methodforfragmentationofincompressible fluid”, Nucl. Sci. Eng., 123 (1996),
421-434.
[5] Gingold, R. A. and Monaghan, J. J., “Kernel

Leadingedge position
X* = x /aestimates as a basis for general particle methods in
3.0hydrodynamics”, Journal of Computational Physics,
46 (1982), 429-453.

[6] Liu, L., Koshizuka, S. and Oka, Y., “Computation

2.5of multiphase flow by coupling the MPS methodwith mesh method”, 日本原子力学会 2003 年春の年会 (2003), 479.

2.0y[7] Enright, D., Fedkiw, R., Ferzger, J., and Mitchell, I.,“A hybrid particle level set method for improved

1.5

interface capturing”, Journal of Computational

Leadingedge

2aPhysics, 183 (2002), 83-116.

x[8] Ishii, E., Ishikawa, T and Tanabe, Y., “Particle/CIP

a

1.0
0.0

1.0

2.0

hybrid method for predicting motions of micro andmacro free surfaces”, Proceedings of ASME

3.0

HT-FED 2004, Paper No. HT-FED2004-56142

Nondimensional time(2004).

1/2

T = t (2g /a )[9] Ishii, E., Kameya, H., Watanabe, A. and Urashin,M., “Calculation of two-phase flows in an oil

Fig. 10 Movement of a leading edge caused by theseparator by using a cubic interpolated propagation

collapse of a water column

code on unstructured grid spaces”, Proceedings of
５ 謝辞

ASME

FEDSM’02, Paper

31229(2002).

No. FEDSM2002-[10] Ishii, E., Ishikawa, T and Tanabe, Y., “Simulation of
ダム崩壊問題の実験写真を東京大学越塚教授よりliquid jet breakup using a combination of particle

ご提供いただいたことに心から謝意を表する．

and

grid

methods”,

Proceedings

of

ASME

FEDSM’05, Paper No. FEDSM2005-77001 (2005).
参考文献[11] Martin J. C. and Moyce W. J., “An experimental[1] Hirt, C. W. and Nichols, B. D., “Volume of fluid

study of the collapse of liquid columns on a rigid(VOF) method for the dynamics of free boundaries”,

horizontal plane”, Phil. Trans. Royal Soc. of

Journal of Computational Physics, 39 (1981),

London, 244 (1952), 312-324.

- 59 -© 2006 The Heat Transfer Society of Japan