粒子法による流れの数値解析

ながれ ）

〔特集〕流れと粒子

粒子法による流れの数値解析

Numerical Analysis of Flow using Particle Method

∗_{東大院・工学系研究科附属原子力工学研究施設}

_越

_塚

_誠

_一

†

Seiichi KOSHIZUKA

1 はじめに ここで紹介する粒子法は，粒子を用いた連続体 の数値解析手法である．すなわち，連続体の運動 を有限の数の粒子の運動として離散化する．有限 体積法や有限要素法と比較すると，（1）完全ラグ ランジュ法であり，（2）格子を全く用いない，と いう特徴がある．完全ラグランジュ法ということ は，移流項を離散化しないことを意味するから， 移流項に起因する数値拡散が生じない．そのた め，自由液面などの界面を特別なスキームを用い なくとも明瞭に追跡することができる．もちろん 有限体積法でも格子を界面に合わせて移動すれば よいわけだが，界面の変形が大きいと格子のゆが みが著しくなり，計算が破綻してしまう．しかし ながら，粒子法では格子を用いないので界面の変 形が大きくても問題は生じないし，流体の分裂や 合体が生じるようなトポロジーの変化に対して特 別なアルゴリズムを用いる必要も無い． 粒子法の研究の歴史は古く，米国Los Alamos 研究所でのPAF（Particle-and-Force）法が最初の 粒子法であろう．その後にMAC（ Marker-and-Cell）法が提案されることになるが，ここでは粒 子は液面のトレーサーとしての役割しか果たして いない．なお，MAC法における格子での非圧縮 性流れの計算アルゴリズムが今日の有限体積法に 発展している．PIC（Particle-in-Cell）法は，移流 項を粒子で，その他の項を格子で計算する方法で ∗_{〒319-1188 茨城県那珂郡東海村白方白根 2-22} † _{E-mail: koshi@utnl.jp} ある．粒子を用いて完全ラグランジュ法で移流を 計算すれば，その時点では確かに数値拡散は発生 しないが，PIC法では粒子と格子の間で変数のや り取りが必要であり，この過程で大きな数値拡散 が発生する．なお，Los Alamos研究所における 一連の粒子法の研究はHarlowのレビューに詳し い1）_． 粒子のみを用いるという意味での純粋な粒子 法は，PAF法の後にSPH（Smoothed Particle Hy-drodynamics）法が提案された2）．SPH法は圧縮 性流れの計算手法で，これまで主に宇宙物理学の 領域で使われてきた．最近では自由液面流れや固 体力学にも適用範囲を広げている．粒子法では各 粒子がある一定の質量を保持しながら運動するの で，流体の密度が高くなる場所には粒子が集中す る．従って，密度の高い場所の空間解像度が自動 的に高くなるので，宇宙物理学にように真空から 星まで極めて大きな密度変化を扱わなくてはなら ない問題に対して，粒子法は適している． 粒子法に非圧縮性流れの計算アルゴリズムを導 入したのがMPS（Moving Particle Semi-implicit）

法である3–5）_{．非圧縮条件を粒子数密度一定の条} 件とし，ここから圧力のポアッソン方程式を導出 した．自由液面は，粒子数密度の低下によって判 定するので液面形状を描く必要が無く，流体の分 裂や合体も容易に計算できる． 粒子法の利点は，（1）界面の大変形を扱うこと が容易，（2）煩雑な格子生成作業が必要ない，こ とである．これらは，複雑な自由液面流れや混相 流の詳細解析においては革新的である．これまで に，砕波6, 7），沸騰8），蒸気爆発9），などの数値

解析に成功している．さらに，固体力学10）や， 流体-構造連成問題11）にも適用されている． 以降では，まずMPS法による非圧縮性流れの 数値解析手法を説明し，次に計算例を示す．計算 例については，計算結果の概要とその意義を紹介 するにとどめるので，定量的な議論に関してはそ れぞれの引用文献を参照して欲しい． 2 MPS法 2.1 非圧縮性流れの支配方程式 非圧縮性流れの支配方程式として以下のものを 用いる． ∂ρ ∂t =0 （1） Du Dt = − 1 ρ ∇P+ ν∇2u+ f （2） 質量保存則には，通常の有限体積法では速度の発 散を用いるが，MPS法では密度一定の条件を用 いる．これは後述する非圧縮性流れの計算アルゴ リズムと関係している．また，MPS法は完全ラ グランジュ記述の計算手法であり，運動量保存則 の時間微分にはラグランジュ微分を用いればよ く，移流項を陽に表記する必要は無い． 2.2 粒子間相互作用モデル 粒子法の基本的な課題は，格子を用いずにどの ように微分方程式を離散化するかということであ る．MPS法では，勾配，発散，回転，ラプラシ アンといった微分演算子に対してそれぞれ粒子間 相互作用モデルを用意し，これらを用いて微分方 程式を離散化する． 重み関数wを導入し，粒子間相互作用モデル にはこの重み関数を利用する． w (r) =_ re r − 1 0 ≤ r < re 0 re≤ r （3） ここでrは粒子間距離である．従って，式（3）の重 み関数を用いると，粒子間距離がパラメータreよ り短い場合のみ相互作用することになる（図1）． e r 図1 粒子間相互作用 r= 0でwは無限大になるが，これは計算を安定 にする効果がある． 現在のMPS法では，計算精度の観点から重み 関数のパラメータreを粒子間距離の2∼4倍とし ている4）_{．そのため，相互作用する近傍粒子の数} が2次元では12∼44個程度になる．有限体積法 では隣接するセルの数が4∼8個であるので，粒 子法の方が有限体積法よりも計算時間が長くな る．ただし，3次元化については，有限体積法で は展開された項の数が著しく増加するが，MPS 法の粒子間相互作用モデルでは，基本的に粒子間 距離の計算に第3番目の次元を追加するだけで あり，プログラミングが簡単である． さて，粒子iの近傍粒子 jに対して，重み関数 の和を取ったものを粒子数密度と呼ぶ． ni=  j
i wrj− ri （4） 各粒子の保持する質量が一定であるとすると，粒 子数密度は流体の密度と比例する．従って，非圧 縮条件（式（1））より，粒子数密度も一定でなけ ればならない．この一定値をn0と書く． 支配方程式（1）（，2）には微分演算子として勾配 とラプラシアンが含まれており，後述する圧力の ポアッソン方程式でもラプラシアンが現れるの で，これら2種類の微分演算子の粒子間相互作用 モデルが必要である．粒子iの位置における物理 量φの勾配モデルは， _∇φ i= d n0  j
_i   _rφj− φi j− ri 2  rj− ri  wrj− ri   （5）

粒子法による流れの数値解析

i

j

I

_j i

I

re



I

_j



_i

r

j



r

i 2

r

j



r

_i





I

図2 勾配モデル とする．これは，図2に示すように，粒子iとそ の近傍粒子 jとの間で単純な勾配ベクトルを定義 し，重み関数wを掛けて平均を取ったものであ る．ここでdは空間の次元数である． SPH法の場合には，まず物理量の分布をカー ネル関数（MPS法の重み関数に相当）の重ね合 わせで表現し，この分布を微分する．そのため， 空間微分はカーネル関数の微分の重ね合わせにな る．SPH法とMPS法の空間微分を比較すると， 粒子配置が乱雑である場合には，MPS法の方が 微分値がばらつきにくい3）． MPS法の勾配モデルは，有限体積法の考え方 を踏襲している．有限体積法では，セル中心に配 置された変数の1階微分はセル境界に現れ，セル 中心で1階微分が必要なときは，セル境界での1 階微分を平均化する．粒子法では，粒子がセルに 相当すると考えると，1階微分は粒子間に現れる と考えるのが順当であろう．そして，粒子位置で の1階微分は粒子間で得られる微分値を重み関 数で平均化したものとする． ラプラシアンモデルは，  ∇2_φ i= 2d λ n0  j
i  φj− φi  wrj− ri  （6） とする．これは図3に示すように，粒子iの物理 量を近傍粒子 jに重み関数の分布で分配すること を意味している．ここで係数λは，分布の統計的 な分散の増加を解析解と一致させるために導入 する． i j r_e 図3 ラプラシアンモデル λ =  j
i rj− ri 2_w rj− ri  j
i wrj− ri （7） ラプラシアンは物理的には拡散を意味してい る．初期分布をデルタ関数とした場合の非定常拡 散方程式の解析解はガウス分布である．任意の初 期分布をデルタ関数の重ね合わせと考えれば，そ の解析解はガウス分布の重ね合わせになる．従っ て，各粒子が物理量を保持していることを，デル タ関数の重ね合わせと見なせば，ラプラシアンモ デルとしてガウス分布による分配を採用すればよ いことになる．しかしながら，ガウス分布は中心 付近で急峻であるため粒子への離散的な分配に用 いると誤差が大きくなることと，分布が無限遠方 まで達するので計算時間が膨大になるという問題 がある． 中心極限定理によれば，ガウス分布とは異なる 分布であっても，分散の増分を一致させておけ ば，その分布による分配を繰り返していくとガウ ス分布に収束する．そこで，誤差を減らすために あまり急峻でなく，計算時間を節約するため有限 の距離までしか達しない分布で分配することを考 えた．これがMPS法のラプラシンアンモデルで ある． なお，拡散に対する格子を用いない計算手法は 渦法でよく検討されている12）．拡散に対して最 も広く用いられている粒子計算モデルはランダム

ウォークであろう．しかし，この場合は粒子の確 率的な挙動として拡散が表されるので，統計的に 十分多数の粒子が必要であり，差分などの決定論 的なスキームと比較すると計算時間が膨大になっ てしまう．MPS法のラプラシアンモデルは決定 論的である． ラプラシアンは，勾配にさらに発散を作用させ たものである．しかしMPS法では，勾配モデル に発散モデルを作用させたものと，ラプラシアン モデルは一致しない．これも有限体積法と似てい る．1階微分を中心差分で離散化し，これを2回 作用させても，2階微分の中心差分とは一致しな い．ただし，MPS法の勾配モデルとして，粒子 位置における平均化された勾配ではなく，平均化 前の粒子間の位置における勾配に対して発散モデ ルを作用させると，ラプラシアンモデルと一致す る13）． 2.3 計算アルゴリズム 非圧縮性流れの計算アルゴリズムには，有限体 積法におけるSMAC（Simplified MAC）法と同様 に，半陰的アルゴリズムを用いる．まず陽的なス テップでは，運動量保存則（式（2））の圧力勾配 項以外の項を陽的に計算し，粒子の仮の速度u∗ と位置r∗を得る． u∗= un+ ∆tν∇2un+ fn （8） r∗= rn+ ∆t u∗ （9） 非圧縮条件より，粒子数密度は一定値n0_とす る必要がある．しかし，陽的なステップが終了し た段階での粒子数密度n∗はn0_{になっていない．} そこで，次の陰的なステップで粒子数密度をn0 に戻すことを考える． n0 = n∗+ n （10） 速度の修正量uは陰的な圧力勾配項によって生 じるとする． u= −∆t_{ρ ∇}Pn+1 （11） 速度の修正量と粒子数密度の修正量は，連続の式 より，次の関係がある． n n0∆t+ ∇ · u _{= 0 （12）} 式（10）∼（12）より，次の圧力のポアッソン方程 式が得られる． ∇2_Pn_{+1= −} ρ ∆t2 n∗− n0 n0 （13） 式（13）の左辺を各粒子位置においてラプラシ アンモデル（式（6））で離散化すると，連立一次 方程式が得られる．これを解けば，圧力Pn+1_が 求まる．この圧力を式（11）に代入すると速度の 修正量が得られ，新しい時刻における粒子の速度 および位置が確定する． un+1= u∗+ u （14） rn+1= r∗+ ∆t u （15） 粒子法では，各粒子がある一定の質量を保持し ているので，質量保存は特に何の条件も課さずと も成立する．非圧縮条件は密度一定，すなわち一 定の質量が一定の空間を占めることが条件であ るので，粒子法の場合には一定の空間の中に一定 の数の粒子が存在するように圧力場を計算する． 粒子数密度（（4）式）は重み関数の空間的広がり 中に何個の粒子が入っているかを表しており，こ れが一定値n0になることを非圧縮条件とする． 従って，圧力のポアッソン方程式の右辺が粒子数 密度になる．式（13）より，右辺のn∗がn0_と等 しければ圧力場は発生しないことがわかる．一 方，有限体積法では，格子を用いることから，本 質的に保存されているのは空間である．従って， 密度を一定にするためには，一定の空間の中での 質量が保たれる必要がある．そのため，圧力のポ アッソン方程式の右辺が速度の発散になる．この 場合，速度の発散がゼロであればポアッソン方程 式の右辺はゼロになり，圧力場は発生しない． 粒子数密度一定の条件から圧力のポアッソン方 程式を導く手順は，PIC法で試みられており14）_， MPS法ではこれを参考に非圧縮性流れのアルゴ リズムを構築したものである． 既に述べたように，式（3）の重み関数はr= 0 で無限大になるようにしている．こうすると，2

粒子法による流れの数値解析 個の粒子が互いに接近する場合に，距離がゼロに なるまで限りなく粒子数密度は増加する．従っ て，2個の粒子が完全に重なり合う前に必ず粒子 数密度がn0より大きくなり，圧力が上昇し，粒 子間には反発力が働くことになる．これは，粒子 のすれ違いや粒子配置のむらを防ぐのに効果があ り，計算を安定化する4）．一方，PAF法やSPH 法では人工粘性を加えることで同様の安定化をお こなっていると思われる． SPH法では完全に陽的なアルゴリズムを用い る．粒子数密度から流体の密度を求め，状態方程 式より圧力を計算する．そのため，SPH法で非 圧縮性流れを計算する場合は，わずかな密度変化 から大きな圧力変化が生じるように状態方程式を 設定する．これは，擬似圧縮性の計算アルゴリズ ムに相当する． 2.4 自由液面の境界条件 MPS法では自由液面の判定条件にも粒子数密 度を用いる．自由液面の外側には粒子が配置され ていないので，自由液面上の粒子は粒子数密度が 低下する．そこで， n∗< β n0 （16） を満たす粒子が自由液面上に存在すると判定す る．ここで，βは0.97などの1.0未満の値を用い る．自由液面上に存在すると判定された粒子は， 圧力のポアッソン方程式を解く際に，Pn+1 _{= 0} のディリクレ境界条件を与える．βの影響を見る ために，同じ問題に対してβの値を変えた計算を おこなった．計算結果はβの値によってあまり 変化しなかった4）_． 粒子数密度を利用する自由液面の境界条件は単 純であり，液面形状を描く必要が無い．流体が分 裂や合体を生じる場合にも，特別な取り扱いは不 要である．ただし，液体内部で負圧が生じる場合 には，粒子数密度が低下することにより流体内部 に自由液面が発生する．これはキャビテーショ ンを模擬しているとも言えるが，通常キャビテー ションが発生する圧力は自由液面の圧力よりも 低い．その場合，MPS法では液体内部で負圧の 発生を許すような特別なアルゴリズムが必要に なる． 2.5 さまざまな現象の粒子計算モデル 本稿では詳細は省略するが，これまでに，沸騰 と凝縮，凝固と融解，表面張力，の粒子計算モデ ルを開発した．大気圧の空気と水のような，密度 差が大きな2種類の流体を扱う場合には，圧力の ポアッソン方程式をそれぞれの流体で別に解く ことで安定に計算できるアルゴリズムも開発し た9）． 固体力学にもMPS法を適用し，剛体，弾性体， 塑性体，粘塑性体，薄肉弾性体の計算モデルを開 発した．特に，厚肉弾性体では，各粒子に回転の 自由度も与えることによって，ハミルトンの正準 方程式による定式化に成功した15）．これは，動 的問題において仮想仕事の原理を満たすことを 意味する．この離散化式は粒子同士が垂直バネ と剪断バネで接続されていると解釈することが でき，計算アルゴリズムとしてはDEM（Distinct Element Method）とほぼ同じになる． 粒子法は完全ラグランジュ記述の方法であるの で，移流項の計算は必要無いが，粒子の移動を追 跡する必要がある．これは従来の有限体積法には 無い計算であり，むしろ天体の運動の計算や分子 動力学と共通である．弾性体ではハミルトンの正 準方程式として定式化されているので，例えばエ ネルギー保存に優れたシンプレクティックスキー ムを使うことができる．しかし，非圧縮性流れで はこれまでのところハミルトンの正準方程式とし て定式化することに成功しておらず，現状では研 究課題である． 3 MPS法による計算例 3.1 水柱の崩壊4） 水柱の崩壊の実験とMPS法による計算を示す （図4）．水柱は支えの板を上方に引き抜くことで 崩れ始める．水は底面を右方向に進み，水槽の右 壁に衝突して高く跳ね上がる．跳ね上がった水は

図4 水柱の崩壊の実験と計算（間隔0.2秒） やがて落下し，今度は左上方に跳ねる．自由液面 の大変形のみならず，流体の分裂や合体が生じて も計算できることがわかる． 3.2 砕波6） 斜面を進む波の砕波を計算した．砕波は，海岸 工学や船舶工学の分野で重要な現象でありなが ら，これまではシミュレーションが困難であっ た．図5にMPS法による計算結果を示す．この 例では，波の先端部が巻きながら落下する「巻き 砕波」が発生した．計算結果を見ると，前の波が 斜面を引いていくときに次の波がやってくるの で，波の先端が立ち上がり，おおいかぶさるよう に砕波していることがわかる．波の波長や斜面 の勾配を変化させると，「巻き砕波」だけでなく， 前面での回転があまり大きくない「崩れ砕波」も 見られた．どちらの砕波が生じるかは砕波パラ メータによってほぼ整理でき，実験とも一致する 結果が得られた． 図5 巻き砕波の発生メカニズム

粒子法による流れの数値解析 3.3 流体─構造連成振動11） 円筒水槽に水を入れ，底面を左右に強制振動さ せる計算をおこなった．ここで，水槽の固有振動 数と液面の固有振動数を一致させ，さらに強制振 動周波数もこれに合わせた．水も水槽も粒子で表 し，それぞれ非圧縮性流れの方程式と薄肉弾性体 の方程式をMPS法で離散化している．図6に結 果を示す．液面と水槽のどちらも連成しながら大 変形している．また，振動が1次元的であるにも かかわらず，液面は「スワール」と呼ばれる回転 運動になっている．MPS法は微分方程式の解法 として一般性があるので，流体解析，構造解析， さらには，流体─構造連成解析にも適用すること ができる． 3.4 核沸騰における単一気泡の上昇と離脱8） 粒子法にメッシュレスな移流項計算スキー ムを加えると，任意ラグランジュ─オイラー （ALE，Arbitrary Lagrangian-Eulerian）法が可能 になる．メッシュレス移流項計算スキームと し てMAFL（Meshless Advection using Flow-directional Local-grid）法を開発し，これとMPS 法と合わせたMPS MAFL法によってメッシュ レスALE法を実現した． このMPS MAFL法を用いて大気圧下の水の プール核沸騰を計算した．バルクの水温は96˚C のサブクール状態で，底面の加熱壁は110˚Cと する．図7は気泡の成長と離脱の1サイクルに おけるx y 2次元計算の結果である．気泡は急 速に成長した後，浮力により加熱壁を離脱する． 上昇する気泡は加熱壁付近の高温水を引き上げて いて，これがサブクール核沸騰における主要な熱 伝達のメカニズムである．計算結果から得られた 熱伝達量は実験と定量的に一致している．沸騰は これまで実験的に得られた熱伝達率を設計等に用 いてきたが，MPS MAFL法では熱伝達率その ものを基礎方程式のシミュレーションから導出す ることができ，沸騰の直接シミュレーションに初 めて成功したと言えよう． この計算が難しい点は，気泡の成長と離脱に伴 い気液界面が大きく変形するのと同時に，きわ めて薄い温度境界層が加熱壁と気液界面に接し て現れることである．すなわち，大変形する界 面の近傍で高い空間解像度が必要となる．従っ て，核沸騰のシミュレーションは，格子を用いる 有限体積法だけでなく，完全ラグランジュ法の MPS法でも難しい．MPS MAFL法はALE法 であるため，流れとは独立に粒子を配置すること ができ，大変形に対応しつつ必要な場所の空間解 像度を高くすることができる．そのため，核沸騰 のシミュレーションをおこなうことができた． 3.5 蒸気爆発における溶融液滴の細粒化過程9, 16） 蒸気爆発は，高温の溶融物が低温の水に接触す ると生じる．原子炉では，溶融炉心が冷却水中 に落下すると発生する可能性が指摘されていた． 蒸気爆発の原因は，溶融液滴の細粒化により伝熱 面積が拡大し，急速な沸騰を生じるためであると 考えられている．しかしながら，極めてミクロか つ高速な現象であるため，溶融液滴の細粒化のメ カニズムについてはこれまで明らかにされてい なかった．MPS法を用いて，溶融液滴に複数の 水ジェットが衝突する計算をおこなったところ， 針状の多数の突起物が形成されつつ細粒化してい くという結果が得られた（図8）．これは，実験に おいて蒸気爆発の瞬間をX線で撮影した映像と 極めて良く一致している． 本計算から得られた知見を溶融炉心の場合に適 用したところ，水ジェットが溶融炉心液滴に衝突 した時点で液滴表面は固化してしまうため，蒸気 爆発が生じにくいこと，また，たとえ大規模な蒸 気爆発が生じたとしても熱エネルギーから機械エ ネルギーへの変換効率は1 %未満であるという結 論が得られた17）．これは原子炉の過酷事故時の 安全性にとって重要な成果である． 従来より原子力工学では，安全が非常に重視さ れていること，実規模の実験が不可能であること から，数値解析が大きな役割を果たしてきた．最 近では，原子炉が溶融するような過酷事故につい てのシミュレーションがおこなわれるようになっ

図6 円筒水槽の変形と液面振動の相互作用

図7 サブクール核沸騰における単一気泡の成長と離脱

粒子法による流れの数値解析 てきた．この場合，原子炉の溶融・凝固，冷却水 の沸騰・凝縮，鋼材やコンクリートの融解など， 非常に複雑な多相・多成分熱流動問題になる．本 稿では1例しか紹介しないが，MPS法はこうし た原子炉過酷事故に関する様々な問題に適用さ れ，成果を挙げている． 3.6 液滴の分裂18） 液滴がそれとは異なる流体中を進む場合，その ウエーバー数が臨界値を越えていると液滴の分 裂が生じる．液滴流における液滴径，および，気 泡流における気泡径は，この臨界ウエーバー数に よって決定される．液滴径や気泡径が決まれば気 液間の界面面積がわかるので，臨界ウエーバー数 は混相流解析にとって重要な値である．従来は実 験によって臨界ウエーバー数や平均液滴径などを 決めていた． ウエーバー数が13の場合の液滴の分裂過程の 計算結果を図9に示す．液滴は図中の左方向に運 動している．液滴は，周囲の流体との相対速度の ため，ベルヌーイの定理により側面の圧力が低下 し，進行方向と垂直な方向に引き伸ばされる．表 面張力はこれを元の形状に戻そうとする．ウエー バー数は，液滴を引きちぎろうとする慣性力と， これを防ごうとする表面張力の比である．ウエー バー数が13の場合には慣性力が勝っており，2 図9 ウエーバー数13における液滴の分裂 個の液滴に分裂した．一方，ウエーバー数が12 では分裂せずに1個の液滴が保たれた．従って， 臨界ウエーバー数として13が得られた．ただし， 本計算はx y 2次元であり，実際の液滴は3次 元である． 混相流の詳細解析は，複雑な界面挙動を扱う必 要があるため，数値流体力学における難問の1つ である．本稿で示した液滴の分裂の計算例は単純 な体系であり実験データも存在するが，他に，実 験データの無い場合として蒸気膜に覆われている 液滴の臨界ウエーバー数を求め，値として50を 得ている．MPS法を用いることで，混相流にお ける界面面積を基礎方程式から計算することがで きるようなってきたと言えよう． 4 おわりに 本稿では，粒子法のこれまでの研究の歴史を 概観し，次に，MPS法の原理を説明した．また， MPS法と他の手法との関係についても解説する とともに，現在の研究課題を紹介した．さらに， MPS法の計算例を示し，大変形を伴う自由液面 流れ，流体─構造連成挙動，沸騰，原子力安全に 係わるシミュレーション，混相流の詳細解析，に 対するMPS法の実績と可能性について述べた．

引 用 文 献

1）Harlow, F. H. : PIC and Its Progeny, Comput. Phys. Comm. 48（1988）1–10.

2）Monaghan, J. J. : An Introduction to SPH, Com-put. Phys. Comm. 48（1988）89–96.

3）Koshizuka, S., Tamako, H. & Oka, Y. : A Parti-cle Method for Incompressible Viscous Flow with Fluid Fragmentation, Comput. Fluid Dyn. J. 4 （1995）29–46.

4）Koshizuka, S. & Oka, Y. : Moving-Particle Semi-Implicit Method for Fragmentation of Incompress-ible Fluid, Nucl. Sci. Eng. 123（1996）421–434. 5）越塚誠一:数値流体力学（培風館, 1997）. 6）Koshizuka, S., Nobe, A. & Oka, Y. : Numerical

Analysis of Breaking Waves using the Moving Par-ticle Semi-implicit Method, Int. J. Numer. Meth. Fluids 26（1998）751–769.

7）後藤仁志,酒井哲郎,沖 和哉,芝原知樹:粒子法 による巻き波型砕波を伴う斜面遡上過程の数値 シミュレーション, 海岸工学論文集 45（1998） 181–185.

8）Yoon, H. Y., Koshizuka, S. & Oka, Y. : Direct Cal-culation of Bubble Growth, Departure, and Rise in Nucleate Pool Boiling, Int. J. Multiphase Flow 27 （2001）277–298.

9）Koshizuka, S., Ikeda, H. & Oka, Y. : Numerical Analysis of Fragmentation Mechanisms in Vapor Explosions, Nucl. Eng. Des. 189（1999）423–433. 10）Chikazawa, Y., Koshizuka, S. & Oka, Y. : A Parti-cle Method for Elastic and Visco-plastic Structures and Fluid-structure Interactions, Comput. Mech. 27（2001）97–106.

11）Chikazawa, Y., Koshizuka, S. & Oka, Y. : Nu-merical Analysis of Three-dimensional Sloshing in an Elastic Cylindrical Tank using Moving Particle Semi-implicit Method, Comput. Fluid Dyn. J. 9 （2001）376–383

12）Ghoniem, A. F. & Sherman, F. S. : Grid-free Sim-ulation of Diffusion using Random Walk Methods, J. Comput. Phys. 61（1985）1–37.

13）Koshizuka, S., Ohta, K. & Oka, Y. : Development of a 3-D Calculation Scheme using Moving Parti-cle Semi-implicit Method for Thermal Hydraulics,

Proc. 6th Int. Conf. Nucl. Eng.（ICONE-6）, San Diego, May 10-15, 1998, ICONE-6225.

14）梅垣菊男,高橋俊,三木一克:粒子法による非圧縮 性粘性流れの解析,数値解析43（1992）17–24. 15）Koshizuka, S., Chikazawa, Y. & Oka, Y. : Particle

Method for Fluid and Solid Dynamics, Proc. First

MIT Conf. on Computational Fluid and Solid Me-chanics, Boston, June 12-15, 2001, 1269–1271.

16）Koshizuka, S. & Oka, Y. : Application of Moving Particle Semi-implicit Method to Nuclear Reactor Safety, Comput. Fluid Dynamics J. 9（2001）366– 375.

17）Liu, J., Koshizuka, S. & Oka, Y. : Investigation on Energetics of Ex-vessel Vapor Explosion Based on Spontaneous Nucleation Fragmentation, J. Nucl. Sci. Technol. 39（2002）31–39.

18）Nomura, K., Koshizuka, S., Oka, Y. & Obata, H. : Numerical Analysis of Droplet Breakup Behavior using Particle Method, J. Nucl. Sci. Technol. 38 （2001）1057–1064.