23【建設工学研究所論文報告集第55号〔論文〕2014年11月】

粒子法をもとにした気液2相流計算法の乱流への適用性の検討
ApplicabilityofTwo−PhaseFlowAnalysisBasedonParticle
MemodtoTurbulentFlows

中山昭彦

久末信

幸

Akihiko Nakayama Nobuyuki Hisasue

1．はじめに乱流の数値計算は、空間に固定された計算格子点での変数の値を微分方程式の差分近似により解く手法が主流であるが、界面が大きく変動する自由水面流や混相流などの場合には、流速など流れの変数を流体粒子とともに移動する点で定義し、これらの値より空間微分項を内挿により求める、いわゆる粒子法（SmoothedParticleHydrodymmics（SPH）
1）やMovingParticle Semi・imphcit（MPS）2））が有用であることが示され応用も進んでいる3）。これらの手法は、Navier−StokeS式原型のみならず、レイノルズ平均運動方程式や空間フィルタ平均運動方程式にも適用できるので乱流の計算も可能である。ColagmSSi＆LandriniO．Ⅵ01eau＆IsSa6）、後藤ら6）などは通常のオイラー乱流統計量の輸送方程式に適用し、乱流計算を試みている。粒子法は流体運動に沿った座標を用いるので乱流のラグランジュ平均量の予測法にもなる。Holm7）はラグランジュ座標系での運動方程式を厳密に平均した式を導出し、ラグランジュ平均には粒子の分散による効果があることを指摘し、ラグランジュ平均とオイラー平均の違いや、移動速度の平均と運動量の平均は異なることを指摘している。著者ら部ま粒子内挿操作自体に空間平均効果があるので、小スケール乱れのある乱流の瞬時流れ場に粒子内指を施し得られる運動方程式をモデル化する方法を提案した。この場合変動する粒子に厳密に追随するラグランジュ記述は不可能で、粒子平均した軌道に沿う座標での運動方程式を扱うことになる。この粒子軌跡変動の影響はMonaghanl）が計算を安定させる目的で導入した項に関連している。しかし、分散効果やサブグリッド応力に相当する付加応力など乱流計算において発生する粒子軌跡変動の効果については十分研究されていない3）。
ここでは、SPH粒子法を移流拡散方程式の離散手法として、固定座標上で空間フィルタ平均されたLarge−EddySimulation（LES）の基礎式に適用する。この場合移流速度は空間平均流速であるので、サブフィルタースケール効果は固定座標モデルと同様のものを適用できるが、その適用法や境界条件の設定法に注意を要する。ここでは乱流のシミュレーション法としての一手法を提案し、基礎的流れで妥当性や効果を検証する。また水面が大きく変化し、空気混入や液滴の飛散などのある流れへの適用の可能性も調べる。

2．SPH法2．1 SPH基礎方程式
粒子法は、例えば流れの変数Aの任意の点′での値は周辺の値を内挿することにより評価できることに基づいている。
すなわちAの内指値AIはAI（′）＝臨申・h函）dr・（1）24とカーネル積分で表わされる。ここでmr−r男は距離中一′●iがおよそhの以内の領域で値をもつ内挿カーネル関数で脚Wい，hV＝1・悪W（r−γ・・h）＝6（′−T・）（2）の正規条件と距離ゼロの極限でその点での変数の値になる条件を満たす。ここで糾−′，）はD血Cのデルタ関数である。内指値を数値計算で評価するには流れ場空間に固定された計算格子は必要でなく、任意の離散点でAが与えられていればよい。AOが有限個の粒子で代表される離散点れ（添字bは粒子を表わす指標）で与えられている場合、（1）式の積分は

4（・）＝∑偽筆W（′一両）（3）と粒子bについての和で近似できる。ここでmb、竹、A小は点れ近傍の流体を代表する粒子bの質量、密度およびA（ので、比mb／励ま体積になる。これよりAsの勾配ベクトルは甑′）＝∑mb舎∇W（′一朝￠）とWの勾配を重みとする平均で表わせることになる。離散内指値As、∇AsなどをA、∇Aなどの近似値と見なせば運動方程式、質量保存式中の項も粒子内指値で表わすことが出来、これらについての基礎式になる。ただし、Mona孤独1）によれば（4）式を圧力勾配項に適用する場合計算の安定性などの面から、Aに直接圧力pを代入するのではなく、

幸＝∇（れ薯と変形して得られる

軸訊講）∇弥（5）を用いるのが有効であるとされている。ここで、左辺は粒子aの位置での圧力勾配、pかpb、は粒子a、bの圧力、旋＝欣亙れh）、∇。耽，は耽あのaの座標についての勾配である。速度の2階微分を含む粘性項にも（4）式を再度微分した式を用いるのでなく、カーネル関数の1階微分のみで表わされるモデル式が用いられる。本研究では画調＆M00a如m9）
による式∇・備V＋（∇↓刑＝−∑硯b∇脇。

（6）b

現b二一五一Vb）高一の（7）pJも（仏十のhを用いる。〟は流体粘性係数で、仏は粒子aの位置での粘性係数、Gは次元に依存する定数で、3次元流れの場合K如紬
＆Monaghan9）は112／15としている。また分母の…01h2はT。＝rbの特異性を避けるものである。以上の関係を用いると粒子aの運動方程式および質量保存式は

誓一言mb（講叫b）∇拙g＋傭および（8）25誓＝∑扉。高∇脇（9）bと離散化される。ここで関り、p。（t‖ま、粒子aの位置判りでの流速と圧力、gは重力加速度、凧は重力以外の体積力で、下記に述べる境界条件壁面モデルにより生ずる。
粒子aの位置互t‖ま速度関りとの関係

告＝Va（10）を積分することにより求められる。
SPH法では非圧縮流も弱圧縮性で近似し、圧力は密度のみの関数とするモデル状態方程式から陽的に求める。その場
合音速が流体の流れより十分大きいようにすれば圧縮性は流れに殆ど影響しないので通常次のようなモデル式

pa＝B（（到（11）を仮定する。ここでpoは密度の基準値、γは7としBは、音速が予測される最大流速の10倍程度になるよう設定される定数であるlr3）。Bの値のとり方で、計算安定性や圧力分布に影響を及ぼす可能性があるので、計算結果を見て判断する。
MPS法では密度は一定とし、式（9）の右辺副より、あるいは粒子密度が変化しないとする式より圧力を求める。本研究では、カーネル関数耽bには次に示す3次スプライン関数を用いる。

k2−。）3−4（1−。）31，0≦。≦1，1′へ

Wb（q）＝＿、3（2−q）
0．

1≦q≦2，（12）q＞2

ここでq＝厄一心／hで、勾配∇。Wbは∇頗q）＝誓∇q（13）から計算される。

2．2乱流モデル

上の式は層流あるいは粒子内挿距離h以下のスケールの乱れのない流れについてであり、高レイノルズ数流れで粘性スケールがh以下の場合適用できない。「はじめに」にも述べたように、乱流を対象とする場合、レイノルズ平均あるいは空間フィルタ平均された運動方程式を粒子法により離散化する方法と、粒子荷開式を平均前の瞬時流れ場に適用する方法がある。後者の場合、粒子運動に適した乱流モデルを導入することが出来、解析精度も向上できる可能性があるが、ここでは、空間フィルタ平均された運動方程式に乱流モデルとして渦粘性モデルを仮定した基礎式に上記の粒子法による離散化を行うとする。この場合分子粘性項の分子粘性係数に渦粘性係数を付加することで乱流の扱いが可能である。
すなわち粘性項脇bを

26

1㌔b二一q（仏＋祐）（偽＋侮） 。一高危一物）（14）

pJみ（仏＋掬＋仏＋侮）h 山一扉＋0．01h2と通常の分子粘性にサブフィルタ渦粘性係数仏、塩を付加する。これらにSma印hsbモデルを適用すると、たとえば仏は也＝qh2is。I（15）であるo鎚S皿gO叫定数、録はひずみ速度テンソル缶担叩周で、軌の値は∇V。＝∑m虐土山。耽
pb

（16）

′ I●から求められる。
粒子法の基礎方程式（（8）、（9）、（10））は個々の粒子の時間についての常微分方程式で、その解には空間座標についての境界条件は必要ない。とくに自由表面での拘束条件は必要としないので、複雑な界面追跡などなしに自由表面流の計算が可能である。気体と接する液体の流れも、それぞれの物性値をもつ粒子を挿入することで、混相流の計算が可能である。上述の式は全て密度、粘性係数は変数として扱えるため、特別な操作なしにそのまま適用できる。ただし、大きく異なる密度の粒子の音速の設定、それに伴う時間ステップ幅の設定など注意を要する点もある。
他の境界での拘束条件は支配方程式の粒子補間項を境界近傍で評価する際に考慮される。不透面近傍では、反発力を付加したり、仮想粒子を壁面あるいは壁面外に配置し、これらの仮想粒子の運動あるいは位置に境界条件を設定する方法と、粒子補間項を評価する際に境界面での値を導入する方法10）がある。乱流の場合、とくに壁面近傍乱れを解像できない場合の境界条件設定の手法を提案する。

2．3

壁面モデル

境界での条件であるが、完全滑り面、完全粘着面あるいは流速に限らず、温度などの物理量についてのデリクレーあるいはノイマン条件の設定法は種々の堤案がなされ、良い結果も得られることも示されている。しかし乱流LESの場合、単純に速度やその勾配を設定するのでなく、複雑な関係を含む壁面モデルの導入が必要になる。本報では次のような壁面処理を提案する。
乱流壁面法則は固定格子で解く場合でも壁面上でなく流れ内の点で適用されるが、粒子法では壁面近傍に壁面層（厚さhb）を設定し、この壁面層に入った粒子に壁面条件を与える。下記の計算ではhb＝h＝1．3ro（小は初期粒子間隔）としている。不透過面近傍では法線方向流速成分には式（8）を適用せず非圧縮連続式を満たすとする。その際壁面では法線方向成分はゼロとする。
壁面接線方向流速は式（8）を用いて求めるが、右辺凡は壁法則で定められる壁面せん断応力とする。これは壁粒子などを設定する手法の壁面相互作用に相当する。壁面までの最短距離を右そこでの接線方向流速を巧とすると、塙くhbとなる場合、接線方向勾配壁面せん断応力による接線方向体積力

旦＝五二一CD航hppとする。V。二回でGは3層壁面則

（17）

27

2．51n吐五十5．5，吐五＞30
巧＿5．01n吐五一3．05，5≦吐五≦30

〟＊

V（18）

V王立≦5V

より求められる摩擦速度がおよびq＝￠拘2から定める。この方法は固定座標を用いるLES法に使われる壁モデル11）を粒子法に適用させたものである。

2．4

数値計算手法の概要

本手法の支配方程式（8）、（9）、（10）式は個々の粒子についての時間発展常微分方程式で、既往の種々の解法を適用できる。工夫を要するのはこれらの式の右辺が全粒子についての和で粒子数に比例して計算量が増加することであるい近傍粒子の探索法を用いるべきで本研究ではHockney＆Eastwood12）によるLinkCcll法を適用する。時間進行は1次精度
Euler法を用いる。LinkCellの幅は粒子内挿距離hの2倍とし、各LinkCell内の粒子の番号を各時間ステップの始めにあらかじめ検索しておくことにより、全ての粒子の和をとる計算を省略できることになる。

3．検証計算3．1 層流開水路流計算はまず十分発達した平坦開水路流に適用する。
目的は乱流のLESであるが、まず粘性項と境界条件の扱い法を確認するため、平行壁面間流れである
平面PoiBeuille流の片側半分に相当する開水路流について計算する（図−1）。この流れでは重力と底面せん断がバランスし、フルード数の低い乱流の場合、圧力差で流れる平行平板間流れとほぼ同等の流れである。層流の場合、流速分布は放物曲線となり、水図−1斜面上層流開水路流面での最大流速肋間肌は斜面勾配により決まるので、粘性モデルと粘着条件の検証になる。水深をH、斜面角度を0、重力加速度をgとすると底面せん断応力品はpg椛hOで動粘性係数をVとすると、最大流速構m歓は互群hO〟Vである。

60

向速度成分Vはゼロで、請の分布は放物曲線である。

きの鉛直分布の計算結果を示す。理論解の底面垂直方80

図一2に層流の場合の速度成分録、∴V及び圧力p

粒子数3200で奥行き方向には4列の初期粒子配置であるが、流速分布は理論値によく合う結果となっている。圧力は計算初期水面での密度をp。とし、p。gHで無次元化した値をプロットしている。理論では密度一定の静水庄分布であるが、計算では弱圧縮を仮定し、圧力は（11）式を満たすので、正確な静水

0020．40」60．8

11．21．4〟かma嶋〆鵬H

圧分布にはならない。本計算では底面近傍でやや大きくなる結果になっているが、）船0．1では静水圧

図−2斜面上層流開水路流の流速および圧力分布

280．00．51．01．53．03．54．0

胸図一3平坦開水路乱流中鉛直断面内瞬時粒子位置と流速ベクトル0）乱れ強度分布（3 平均流速分布

持雪ぎ計ミ至言もき

252015105○均のSロ

Vの調

△男mS

口のロ含量難曲曲田

01000．2

0．4

0．60．811．2yH

図一4平坦開水路乱流中鉛直断面内平均流速分布と乱れ強度分布分布に近い値になっている。

3．2 完全発達開水路乱流
次に同じ斜面上流れであるが、レイノルズ数が大きく完全発達した乱流の場合をとりあげ、摩擦速度がと平均水深Hで定義される摩擦レイノルズ数の値で約800、平均流速叫Vを基にしたRe三脚Vは16000の乱流についての計算を行った。計算領域は流れ方向に〟1奥行き方向にHで、合計62，720の粒子を流れ方向、鉛直方向、奥行き方向にそれぞれ112、28、20個等間隔に配置し、初期流速分布は固定格子の計算結果13）を内挿して求めた。計算領域外の流れは計算領域内と同一のものが周期的に繰り返すとし、粒子が計算領域外に出た場合、反対側境界内に同じ位置と流速をもつ粒子を挿入している。この粒子数でも、10倍程度の固定格子を用いる場合の計算時間を要する。しかしGPUなど並列演算そ用いることにより計算時間は大幅に改善できる。
図−3は鉛直断面内粒子の瞬時位置とその流速分布の一例である。計算開始後0．5耽経過したときもので、乱れの小さい上の方の粒子は初期位置から平行移動したような位置であるが、底面近傍の乱れの大きい領域では不規則になっている。分布密度は全体であまり変わらないほぼ一様性を保っている。底面近傍の乱れは主流内に伝わって行く様子も見られる。定量的特性は次に示す平均量で評価する。時間平均は、瞬時流れ粒子計算結果を60露0露0の固定格子点青上で内挿し、同一点で時間平均し求めた。速度成分の内指値は（3）式の粒子補間A（扉∑耽重昭（19）

29

実験（Koshizukaetal．2））！≒0．1sec！≒0．2sec

I〇回！≒0．4sec￣

l＝0．8sec′

手車
r l上p l i 直子
一lih！

1・一一 一一 一一一 二二° ．
言 i

1，−−ー
−
一一−ー

計算結果
404030301010．0 20．0 30．0 40．0 50．0．両．In）

〇ーき．
︶
−
︵520亨．

1000．010．0 20．0 30．0 40．0 50．0時（周）00．010．020．030．0 40．0 50．0重（，nl0．010．020．030．0 40．0 50．0．†（（−I）il

図−5障害物なし水槽内水柱崩壊流れを用いて求められた。粒子の間隔は壁面粘性単位で30程あり粘性底層やバッファ層は解像できていない。時間平均のための内挿点の間隔もこれとほぼ同じである。水面外の点でも内拝は可能であるので表示しているが、内挿距離hの外の点ではゼロになる。
図一4（a）は平均流速分布を壁面則座標6＋＝撮り旬〟＋『諦＊）で示す。滑面標準壁法則にほぼ従う結果となっている。
図−4（b）はRMS乱れ強度成分（流れ方向的関、底面垂直方向γ関、奥行き方向W，棚）の鉛直分布を示す．通常〟棚は最も大きく，1㍍鳩の倍程度であるが、本計算では底面近傍以外の領域で壁面垂直方向乱れV恥が流れ方向成分〟朋に近い値で、より等方的な結果になっている。これは壁面境界の設定法、モデルにも影響されていると考えられるが、平均時間も十分でない可能性もあり、用いられた少ない粒子数にも関わらず基本的分布特性は、これまでの実験、オイラー座標上固定格子でのシミュレーション結果（久末ら13））と整合していると言える。

3．3 水柱崩壊流れ
次に水面変動の大きい流れの例として、水柱崩壊流れの数値解析結果と実験結果との比較を示す。実験はKoshizl．Lkaetal．2）により行われたもので、長さ59cm深さ40cm以上の水槽の左端に幅14．5cm高さ29cmの水柱を板で支え、時間FOに支持板を取り除いた後の水柱の崩壊過程流れの様子を可視化したものである。図−5の上段に実験結果を下段に本計算結果を示す。この流れは静止から加速していき、初期は乱れの影響は小さいと考えられ、数値解法の検証のみにも使われているが14）、崩壊水塊が壁面に衝突するまでには流速は毎秒数10cmとなり、乱流遷移が起こり、石壁衝突後は上方い跳ね上がり、そのご落下しながら左方向に戻って行く複雑な流れである。
解析結果は全体の流れの傾向を良く捉えている。従来用いられた過大な疑似粘性とすべり境界条件では水柱先端の角度が小さく進行速度が大きく見積もられていたのが、本計算法の声0．2S∞の結果では実験に近い先端形状と移動速度になっているのが分かる。また固体壁面境界近傍で反発力を与えるMona如ml）の方法では石壁衝突後水が壁面から離れる傾向になり、跳ね上げられた水滴が水槽中ほどに落下するが、壁面近傍法線方向速度を与える本計算ではそれがなく、戻り流れの巻き波の形状なども実験に近い結果となっているのが分かる。

3．4 空気混入を伴う障害物のある水柱崩壊流れ
上の水柱崩壊計算は少なくとも時間声0．8secまででは空気と水は混合せず液相のみの計算で流れの再現ができている。
しかし、空気が水の中に閉じ込められる場合、空気の運動はともかく水と壁の間に閉じ込められる空気の影響を考慮す30

実験結果
三0．2sec（「）（＝￣値．蛤「（h）！ニ0．18倍

FO．4seci＝0．38cc（｛販＝0：撮t＝0．5sec

Iり（ニ0了肌（Ol＝0疑cc

液相のみの計算

01020∞

4050010青Illml2000lIl405omI0102000
1IllnII

405001020

鞠I ・00∫Irml

500102030∴∴40
1IlWIl50

2相流計算、こ ●・、、JL∴二ー責（品 −i

図−6底面に障害物の設置された水柱崩壊流れる必要がある。ここでは、やはりKo血柵k訓こより行われた実験で上記の実験条件と同一であるが、底面に障害物が置かれている場合の水柱崩壊過程を計算する。この場合、水柱は障害物により跳ねあげられ、その後落下するがその際水塊と底面、側壁の間に空気が閉じ込められる。
図−6に実験結果と本計算手法による液相のみの計算と水柱周辺の空気の運動も計算する2相流計算を行った結果を示す。空気を考慮する2相流計算では、計算負荷を考慮し、水内に閉じ込められる位置にある底面から水柱高さの1／3の深さの位置にのみ空気粒子を配置しそれ以外では水粒子空気粒子とも配置していない。また音速差による計算不安定を考慮し、空気相の初期密度は水の1／20としている。このため空気の速度は実際より小さくなっているが、本計算では単相計算では無祝されていた空気の存在の効果を考慮するのが主目的で、その運動の再現性は今後の課題である。
液相のみの計算結果はFO．3S∞までは実験結果に合った結果になっているが、創．3secでは飛び跳ねた水塊が底面に落下し、実験と大きく異なっている。《これに対し空気を考慮した計算では戸0・4sec以降飛び跳ねた水塊がそれと底面の間の空気により落下せず石壁に接するようなかたちになっている。仁0．5S∞では閉じ込められた空気塊が二つに分かれ、障害物先端からはねた水が底面に落ちていく様子など実験に近い結果になっているのが分かる。

4．階段状水路流れ
次に、気液2相流計算の適用例として、空気気泡混入、水滴飛散の顕著な流れで、実験計測などの行われている階段状水路流れを取り上げる。この流れはダム洪水吐や減勢工周りの流れに関連しており、このような流れの解析への適用性を考察することができる。
安田ら15）は実験により、階段状水路の流れの特性を調べている。こういった流れの定量的計測は非常に難しいので、本計算もその適用性の評価が目的で、定量的検証は今後の課題とする。図一7は安田ら15）の実験水路の模式図と実験風景を示す。計算はこの水路形状に合わせるが、計算領域は4段の階段のみとし、計算条件は水面形状が水路床形状の詳細に依存しないsHmmhgHowに分類される勾配15度の階段水路について行った。
水柱崩壊流れの解析と同様に、初期状態として水面上4列の空気粒子を配置する2相流計算である。水粒子は上流端

31−
二十￣
、i・
妻 豪 雪輔
「アヤ：
禦
＝＝隷
後
釜
勘，、遅〔ミ三三三三i雲暮雪l雲＝＝＝ご．
＿二∴
…
…
三二

上

薫 雲 霞雲：‡∴∴−書稜 二三三二＿

一

題闇 ￣−

ie

二一
埼 二心・一・三や−∴ 、ニ三一＿＿己∴ 「㌣￣￣ ￣ ∵i
愛妻 薬 曇雲 霞 襲薄 幸￣

図一7安田ら15）による階段状水路流れ（b）内拝による流速ベクトル値（b 粒子位置へSo︶へ
一−

5 0 5 0 5 0 5 0 5 02 2 1 1ーl 1 2

一一一一￣−● −● 一、−−−−−−−−● 一一−一一￣一一●一一￣−●一一ヽ一一一一一一一・、、●一一−一一●一一●一一一一一一一一●一一一＼●一・一
一

ー ＿一二
＿−

・一

・＿−
． −
＿． −
． 一
、 −
． −
． −
． −
．

一
事． ・−
． 一

一一
÷、

ご 一二二
ミ÷… ニ ー−、一二

竺一・三 三−

100102030405060可am）
大きい点：水粒子、小さい点：空気粒子

実線：等粒子密度比による水面位置
図−8階段状水路流れの計算結果

から摩擦速度0．2m／Sの対数則分布をもつ流れが流入する条件である。本計算では計算領域の上流端及び下流端のすぐ外に境界ゾーンを設定し、この境界ゾーンにある粒子に自由流出境界条件を与える。また対数則流速分布をもつ粒子列を一定の粒子間隔で導入している。下流端では境界ゾーンから流出する粒子は粒子リストから削除し粒子番号を振り直す操作を行う。
図一8に階段状水路流れへの2相流の適用例を示す。図一8（a）は粒子の瞬時位置を示す。大きい点が水粒子、小さい点が空気粒子の位置を示す。水面での空気の混入、階段隅角部に残る空気塊など2相流計算で可能になる現象が再現されている。図−8（b）はこの自由に移動する粒子の位置での流速結果を粒子補間により、等間隔固定格子上に内挿した結果である。水面位置は、内挿された密度の値が水粒子の値と空気粒子の値の中間値をとる面を検索したものである。飛散する液滴、水面下に閉じ込められた空気泡が再現されているのが分かる。実験でも流速分布など定量的結果が得られなく、詳細の検討は十分でないが、粒子の空気混入率分布など定量的評価も含め今後の課題とする。

5，結論
水面の大きく変化する流れや、空気、水の混合する2相流の解析に有用な粒子法のひとつSPH法を、工学流れの実用シミュレーション法に発展させるため、基礎的手法と必要なモデルの提案、さらに基礎的流れでの検証を行った。とくに、高いレイノルズ数流れの乱流効果、それの反映に不可欠な壁面モデルと境界条件の設定法を検討した。サブフィ
ルタモデルに渦粘性モデルを採用すると、分子粘性項と同様の扱いが可能で、分子粘性モデルの拡張で対処可能である。
壁面モデルとして、壁面近傍に粒子間隔の約1．3倍の厚さの境界域を設定し、その層内の粒子に壁法則を適用する手法を提案した。まず理論解のある層流で粘性と壁面境界条件の設定が確認され、乱流のシミュレーションではサブフィル

32

タ応力の導入効果を確認した。使用した粒子数が十分でないこともあり、固定格子のLES結果に比べやや劣る結果であるが、平均流速及び乱流強度の分布特性は正しく得られた。水面が大きく変化し、液相の飛散のある水柱崩壊流れの計算では、これまでの乱流モデル、壁面モデルの結果に比べ実験値に近い結果が得られた。今後粒子配置を多くすることで、さらに一般的自由水面乱流のシミュレーションに適用できることが期待される。

参考文献

1）Mona虫an，J・J・：ParticlemedlOdsfbrhydrodyni的Dics，anZputel・PhviαRqorw，Vol．3，PP．71−124，1985．
2）Koshizuka，S・，Tm調，H．andOka，YApaJticlememodfbrincompressiblevis∞uSnOWWidhnuid龍gmentation，（九mpu幼aonalF弛め砂場amics．kmmal，l硯．4，m．2946，1995．

3）Monaghm，J・J・：Smoo血edparticlehydodynamics，Rq．mpg．恥．Ⅶ1．68，Pp．1703−1759，2005．
4）Colagrossi，A・andLm血ini，M．：Numehcalsimulationofinterfh。ialnowsbysmoo血edparticlehydrodynalnics，JCbnzl融Phys．Vol．191，Pp．448475，2003．
5）Ⅵ01eau，D．andIssa，R言Nmericalmodellingof∞mplextuhulmtをee−Surfわenowswi心血eSPHme血od：anOVe重心ewJ証．JNme雛Me脱F被i Vbl．83，m．277−304，2007．

6）後藤仁志，酒井哲郎，芝原知樹：SPS乱流モデルの導入による新しい粒子法の展開，水工学論文集，第44巻，pp．575−580，2000．

7）Holm，D・D−Fluctuatione能ctson3DLagrangianmeanandEJehanmeallHuidmotion，DDISicaD，Vbl．133，
pp．215−269，1999．
8）hokuma，H・，Nakayama，A．andmenaga，K．：OnCalcuationofTuhulmtFree−Su魚ceFlowswi血ParticleMemod，Proc．
Seven血SymposiumonnlhulenceandShe紬FlowPhenomena，Ⅵ）1．1，3B2P2011．
9）KajtaJJ・andMonaghan，J・J言SPHsimAebionofswimninglinkedbodies，JαnlPut．Phv．，Ⅶ1．227，PP．8568−8587，
2008．
10）Ferrand．M・，Ⅵ01eau，D・，Rogers，B．D．andL細山rence，D．，Consistentwallbomdはytreatmentfbrl血ninarandtuJbulent
nowsinSPH，PnOe・qf協e6碗hte棚naOnal皮，hericwoI心h呼，Hbmbw電Gema妙PP．275−282，2011．
11）Kasai，TandNakayama，A−StudyofWdlModelinLESCalct山bionofTubulentFlowOverConplexBomdary，Jbumald

助成vscienceandtD肋uticE71gineering．Ⅶ1．28，Nol，pp．117−125，2010．

12）Hockney；R・WandEasw00d，J・W：a規則肋Simuh，蘭nZkingLhruehs，lns瓜uteofPhysicsmbli血ing，Bhstol，U．K．，
Philadelphia，U．S．A，1988．

13）久末信幸，中山昭彦，横嶋哲：実スケール開水路乱流解析のためのLESサブグリッドモデルの検討，建設工学研
究所論文報告集第55号，pp．47−55，2013．
14）番藤敦史，村上拓、中山司：有限要素法をベースとしたフリーラグランジュ法による粘性自由表面流れの数値解析，ながれ，Ⅵ1．28，m．399408，2009．
15）安田陽一，高橋正行，大津岩夫：階段状水路におけるSkimingflowの内部特性，水工学論文集，第48巻，pp．499−504，2004．

著

者

中山 昭彦

所員，m．D．，流体工学，水理学

久末 信幸

関西電力（株）博士（工学）水理学