u t x p u t q () p u C q (3) ls x x I I u 0 (4) t x p 0 q f (5) x xs q (6) * x xs 計 算 を,CIP 法 に 代 わりにTHINC/WLIC 法 を 用 いて 実 施 する. 固 相 にLagrange 粒 子 を 配

3次元固気液多相乱流数値モデルDOLPHIN-3D

への高精度界面捕獲法THINC/WLIC法の導入

川崎浩司

1

_{・松野哲弥}

2 1_{正会員 株式会社ハイドロソフト技術研究所（〒530-6126 大阪市北区中之島3-3-23中之島ダイビル26F）} E-mail:[email protected] 2_{正会員 パシフィックコンサルタンツ株式会社（〒451-0046 名古屋市西区牛島町2-5）} E-mail: [email protected] 本 研 究で は，異 相 界面 の捕獲 手 法の 高精度 化 と質 量保存 性 の向 上を目 指 し， 高精度 界 面捕 獲 法 THINC/WLIC法を3次元固気液多相乱流数値モデルDOLPHIN-3Dに導入し，その妥当性と有用性を検証す ることを目的としている．まず，THINC/WLIC法を方形波の移流計算，円に対するせん断流れ，球に対す る渦流れに適用し，精度検証を行った．その結果，同手法は十分な精度を有することを確認した．ついで， THINC/WLIC法を導入したDOLPHIN-3Dを水柱崩壊問題に適用し，非常に高い質量保存性を有することを 示した．さらに，固気液相の相互作用現象に対する本モデルの妥当性を検証するために，段波と矩形剛体 の衝突・漂流問題に適用した．その結果，同モデルは剛体に作用する圧力を良好に再現するとともに，段 波の矩形体衝突時における水塊の跳ね上がりも表現可能となり，異相界面の捕獲精度の向上を確認した．

Key Words : THINC/WLIC method, multiphase flow simulation, dam break, interface capturing

1. 序論 近年，東海・東南海・南海地震など，南海トラフ地 震に伴う津波災害に対して警戒を強めており，東北地 方太平洋沖地震以降，従来の防災対策を見直す必要に 迫られている．また，地球温暖化に伴う海面上昇，台 風勢力の増大化が懸念され，高潮・高波による水害を 防ぐために更なる対策を講じなくてはならない．高 潮・高波災害は，直接的な浸水被害のみならず，沿岸 部にある多数の木材やコンテナ等が漂流することがあ る．津波災害においても，コンテナ，木材，船舶等が 漂流し，漂流物の衝突による二次被害が発生する．こ れらの被害を最小限に抑えるためには，波浪変形や越 波・越流現象の詳細を把握し，さらには波の作用によ る漂流物の挙動・衝突力の検討が求められる． 以上の背景から，今後の沿岸防災を考える上で，波 浪－構造物の相互作用や気泡連行を伴う砕波現象など， 複雑な沿岸海象現象を把握・解明することが極めて重 要となる．そこで，著者ら1)_{は，これまで3次元場を対} 象とした圧縮性・非圧縮性流体を含む固気液多相乱流 場 を 高 精 度 に 解 析 可 能 な 数 値 モ デ ル DOLPHIN-3D （Dynamic numerical model Of muLti-Phase flow with Hydrody-namic INteractions -3 Dimension version）を開発してきた．

そして，同モデルを様々な水理現象に適用し，その有 用性を検証してきた．しかしながら，自由表面の捕獲 精度や質量保存性に改善の余地が残されていた．

本研究では，異相界面の捕獲手法の高精度化と質量 保存性の向上を目的とし，高精度界面捕獲法の一つで あるTHINC/WLIC（Tangent of Hyperbola for INterface Captur-ing / Weighed Line Interface Calculation）法2)を3次元固気液多 相乱流数値モデルDOLPHIN-3Dに導入し，その妥当性と 有用性を検証する． 2. 3次元固気液多相乱流数値モデルDOLPHIN-3D (1) 基礎方程式 3次元多相乱流数値モデルDOLPHIN-3Dの基礎方程式 は，Navier-Stokes方程式 (1)，質量保存式(2)，圧力方程式 (3)，異相間の割合を表す密度関数の移流方程式(4)，バ ロトロピー流体に関する状態方程式(5)から構成される． i i ij ij j si i i j i j i u D S x f g x p x u u t u       _{ }                    2 1 (1)

 

q x u t i i   _      ₍₂₎       _          q x u C x p u t p i i ls i i 2  (3) 0       i I i I x u t   ₍₄₎

 

p f   (5)







s



s x x x x q q       0_* (6) ここで，は流体の密度，t は時間，uiは各方向の速度 ベクトル，xiは座標ベクトル，q は造波ソースのわき出 し強さであり，式(6)に示されるように造波側線 x = xsで q*の値をもつ．p は圧力，giは重力加速度ベクトル，fsiは

Brackbill et al. 3)が考案した CSF（Continuum Surface Force） モデルによる表面張力，ijは Salvetti and Banerjee 4)が開発

した dynamic 二変数混合モデル DTM（Dynamic Two-parameter Mixed model）に基づく LES（Large Eddy Simula-tion）モデルによる SGS（Sub Grid Scale）応力，は粘性 係数，Sijはひずみ速度テンソル，Diはエネルギー減衰領 域における係数を示す．また，Clsは局所音速，はあ る計算格子における I 相（I = 1 ~ 3；：固相，：液 相 ，  ： 気 相 ） の 密 度 関 数 で あ り ， 







の関係を満たす． (2) 計算アルゴリズム 本モデルは，Navier-Stokes方程式，圧力方程式を時間 分離解法により移流段階と非移流段階に分割し，計算 を行っている．計算フローチャートを図-1に示す．は じめに，初期条件・境界条件の設定を行う．ここで， 固相のみLagrange粒子を自動的に配置することが可能で ある．造波させたい場合，波高・周期の入力条件から 造波ソースを設定した位置におけるソースの湧き出し 強さを計算することができる．その後，流体の流速ui ， 圧力pの移流段階計算に，CIP（Constrained Interpolation Profile）法5), 6)を適用する．そして，算定された移流段階 後の各物理量を用いて，圧縮性・非圧縮性流体を同時 に解析可能な拡張SMAC（Simplified Marker And Cell）法 により非移流段階計算を実施する．なお，気液界面に 生じる表面張力の影響はCSFモデルにより算定し，乱流 量についてはdynamic二変数混合モデルDTMに基づく LESモデルにより評価する．非移流段階計算後，更新さ れた流速を使用し，密度をCIP-CSL2（Constrained Inter-polation Profile Conservative Semi-Lagragian of second-order）7) により計算する．前述の流動解析を固気液全相に対し て行った後，固相内部の圧力を使用して複数剛体の運 動解析を実施する．その後，剛体速度の流体場への影 響を考慮した流速を用いて，全相の密度関数の移流 計算を，CIP法に代わりにTHINC/WLIC法を用いて実施 する．固相にLagrange粒子を配置した場合，固相の流速 を固相に配置された計算粒子に内挿補間し，粒子によ るLagrange剛体解析を行い，剛体形状を修正する．最終 的に各物理量を更新し，次の時間ステップへ移行する． 上述した計算を時間ステップごとに繰り返し行うこ とにより，本モデルでは複数剛体の動的挙動を含む固 気液多相乱流場を精微に数値解析することができる． 3. THINC / WLIC法 (1) THINC / WLIC法の概要 THINC法は界面を高精度に解析することを目的とし， Xiaoら8)によって開発された高精度界面捕獲手法である． THINC法の特徴は，補間関数に双曲線正接関数を用い ることで，数値拡散や数値振動を抑えることができる． また，界面識別のためのステップ付加も必要なく，ア ルゴリズムが単純であり計算効率が良い．さらに， Yokoi 2)はTHINC法を改良し，法線方向を考慮した THINC/WLIC法を開発した．THINC/WLIC法は，界面の 法線方向に関するアルゴリズムが新たに付加されるが， 図-1 計算フローチャート

従来のスキームの単純性を保持したまま高精度な解析 が可能である． 式(7)に示す移流方程式を対象に，一次元のTHINC法 について詳述する．

 

0          x u x u t    ₍₇₎ ここで，tは時間，uは移流速度，は流体率，xは空間座 標である． 図-2は，THINC法の概念図である．ここで，i-1，i， i+1は格子点位置を表している．また，i は補間関数の ジャンプの中心座標である．図-2のように流体率を双 曲線正接関数によって補間する．補間関数i は式(8)で 表される．

 

_                   _    i i i x x x x      1/2 tanh 1 2 1 ₍₈₎







1 1



1 1 1 1            i i i i      (9) ここで，は関数の勾配を制御する変数，は界面勾配 の方向であり，式(9)によって決められる．なお，は通 常3.5が用いられる．勾配の決定方法に関しては後述す る．また，i は未知数であり，式(10)を満たすように決 められる．

 

n i x x i i i i dx x x   



   2 / 1 2 / 1 1 ₍₁₀₎ そして，式(7)における流束fの算定にTHINC法が用い られ，式(11)から求められる．

 



      t u x x iup i i i i dx x f 1/2 1/2 2 / 1 2 / 1 (11)

ここで，iupはui+1/2<0のときi+1，ui+1/2



0のときはiである．

また，THINC/WLIC法では，重み付け関数を使用して 次式で流束が計算される． z z y y x x f f f f    (12) ここで，は重み付け関数，fx，fy，fzは各軸と平行な界 面である． (2) THINC/WLIC法の精度検証 まず，1次元矩形波の移流計算を行った．計算領域は 500，計算格子間隔xは0.01m，計算時間間隔tは0.001s， 流速uは0.80m/sとした．5000ステップ後の計算結果を図-3に示す．図-3から，計算結果は数値拡散，数値振動も なく初期形状と精度良く一致していることがわかる． 次に，円のせん断流れに本手法を適用する．計算領 域は400×400×3，計算格子間隔はx =y =z = 0.01m， 計算時間間隔tは0.005sとした．(u, v, w)=(sin(x)cos(y), -cos(x)sin(y), 0.0) の流速により円は時間経過とともに引き 延ばされ，2000ステップ以降は流速を逆にすることで 変形された円を元の位置，形状へと戻る．なお，円の 半径は0.3m，中心座標は（0.76，0.76）である．図-4に それぞれx-y断面における初期形状から4000ステップま での計算結果を示す．同図から，円がせん断流れによ って細長く変形され，4000ステップ後では元の円形に 戻る様子がわかる．同図(d)から，界面は数値拡散がな く初期の円形状が保持されていることが確認できる． 式(8)における勾配の値による計算精度について検討 する．図-5はそれぞれ=1.0，3.5，20.0における4000ステ ップ後の計算結果である．同図(b)から，=1.0の場合， 円の形状は保持されているものの，異相界面で数値拡 散が生じている．また，= 20.0の場合では，形状が保 持されていない．一方で，図-5(c)から，= 3.5の場合， 数値拡散もほとんどなく形状も保持されている．以上 から，本論では= 3.5を用いることとする． 3 次元流動場における渦流れによって球体を変形させ る計算を実施した．計算領域は 200×200×200，計算格 子間隔はx =y =z =0.005m，計算時間間隔tは 0.0001s とした． u = 2sin2 (x) sin (2y) sin (2z) cos (t / T ) , v = -sin

(2x) sin2 (y) sin (2z) cos (t / T ) , w = -sin (2x) sin (y) sin2

図-2 THINC法の概念図

(z) cos (t / T ) の流速を与え，周期 T=3.0s で変形 された球体は元の形状へと戻る．なお，球は中心座 標（0.35，0.35，0.35），半径は 0.15m とした．計 算結果を図-6 に示す．図-6(b)と図-6(d)を比較する と，表面に若干の乱れがみられるが，形状に大きな 差は確認されない．また，図-6 から，初期形状と 8T/8s 後を比較すると，球から分離した微小な塊が 存在するものの，複雑な流動場においても形状が保 持されていることがわかる． 以上から，THINC/WLIC法に関してそれぞれ1次 元・2次元・3次元の物理量の移流計算を実施した結 果，その妥当性が確認された． 4. DOLPHIN-3DへのTHINC/WLIC法の導入 (1) 水柱崩壊問題への適用 前 節 に て 精 度 が 確 認 さ れ た THINC/WLIC 法 を DOLPHIN-3Dに導入し，水柱崩壊問題に適用する． 計算領域は150×10×80，計算格子間隔はx=z =0.005m，y=0.01m，計算時間間隔tは0.005sとし た．なお，水柱サイズは0.10×0.10×0.20mである． 段波の先端挙動に関して，Martin and Moyce 9)_と

Hu and Sueyoshi 10)による実験結果と比較検証した． 図-7に，実験結果，CIP法とTHINC/WLIC法による 計算結果を示す．計算結果はいずれも実験結果を良 好に再現している．しかし，質量保存の時系列を示 す図-8をみると，CIP法の場合，質量は初期質量か ら94%まで減少するが，THINC/WLIC法では質量は ほぼ100%のまま高精度に保持されている． (2) 段波と矩形剛体の衝突・漂流問題への適用 本モデルの妥当性を検証するために，著者ら11)_が 図-7 水塊先端の移動距離 (a) 初期状態 (b) 1000 ステップ後 (c) 2000 ステップ後 (d) 4000 ステップ後 図-4 円に対するせん断流れの計算結果 (a) 初期状態 (b) =1.0 (c) =3.5 (d) =20.0 図-5 による計算結果の比較（4000ステップ後） (a) 初期状態 (b) 2T/8s (c) 4T/8s (d) 6T/8s (e) 8T/8s 図-6 球に対する 3次元渦流れの計算結果

実施した水理模型実験に本モデルを適用する．実験 模型の概略は図-9のとおりであり，段波が矩形体に 衝突し発生する自由表面挙動と矩形体の作用波圧 に関して検討を行った．計算条件を表-1に示す． 剛体を固定した場合の段波の伝播過程に関する実 験結果とCIP法，THINC/WLIC法をそれぞれ用いた 計算結果の比較を図-10に示す．同図(a)から，時刻 t=0.4sで段波が矩形体に衝突する．図-10(b)をみる と，時刻t=0.6sでは，水塊が上方へ高く打ち上がっ ている．THINC/WLIC法を用いた場合，CIP法によ る計算結果に比べ，水塊の打ち上げ高さなど水理模 型実験を良好に再現している． ついで，矩形体前面に作用する波圧の比較を図-11に示す．同図から，すべての計測地点における 瞬間的な波圧の発生とその後の圧力低下に関して， 両計算結果は実験結果と概ね一致しているといえる． また，両計算結果は，実験結果と比べ僅かな時間差 がみられるが，これは，計算においてゲートの急開 動作を考慮していないことによるものと考えられる． 矩形体を非固定とした場合を対象に，固気液相の 相互作用現象に対する本モデルの妥当性を検証した． 段波の伝播過程と矩形体の漂流挙動に関する実験結 果，CIP法とTHINC/WLIC法をそれぞれ用いた計算 結 果 の 比 較 を 図 -12 に 示 す ． 同 図 (a) を み る と ， t=0.5sでは，実験結果において水塊が弓なりに高く 打ち上がっている．THINC/WLIC法による計算結果 は，CIP法による計算結果では再現できなかった水 塊の弓なりとなって打ち上がる様子を良好に表現し ている．その後，図-12(b)に示す時刻t=0.8sには， 実験結果，計算結果ともに水塊が矩形体を越流し激 しく打ち上がっている． (a) P1 地点 (b) P2 地点 (c) P3 地点 図-11 矩形体に作用する波圧の時系列変化 図-8 質量保存性 図-9 実験模型の概略図 CIP 法 CIP 法 THINC/WLIC 法 THINC/WLIC 法 (a) t=0.5s (b) t=0.6s 図-10 水柱崩壊に伴う段波と矩形体の衝突解析 表-1 計算条件 計算時間間隔 0.00004s 計算格子間隔 x=_z=0.005m y =0.01m 計算領域 230×20×80 水柱サイズ 0.15×0.20×0.20m 矩形体サイズ 0.04×0.20×0.04m

以上から，THINC/WLIC法を導入した本モデルは， 異相界面の捕獲に優れ，非常に高い質量保存性を有 することが検証された． 5. 結論 本研究では，3次元固気液多相乱流数値モデル DOLPHIN-3DにTHINC/WLIC法を導入し，その妥当 性を検証した．実験結果，従来のモデルとの比較か ら質量保存性と界面捕獲の精度向上が確認された． したがって，本研究で構築した新たなDOLPHIN-3D は，固相－気相－液相の相互干渉を解析可能な高精 度かつ汎用性に優れた数値モデルであるといえる． 参考文献 1) 川崎浩司，袴田充哉：3 次元固気液多相乱流数値モ デル DOLPHIN-3D の開発と波作用下での漂流物の動 的解析，海岸工学論文集，第 54 巻，pp.31-35，2007. 2) Yokoi, K. : Efficient implementation of THINC scheme:

A simple and practical smoothed VOF algorithm, J. Com-putational Physics, Vol.226, pp.1985-2002, 2007. 3) Brackbill, J.U., Kothe, D.B. and C. Zemach : A contiuum

method for modeling surface tension, J. Computational Physics, Vol.100, pp.335-354, 1992.

4) Salvetti M.V. and Banerjee S. : A priori tests of a new dynamic subgrid-scale model for finite-diffrence

large-eddy simulations, Phys. Fluid, Vol.7, No.11, pp.2831-2847, 1995.

5) Yabe T. and Aoki T. : Universal solver for hyperbolic equations by cubic-polynominal interpolation One-dimensional solver, Computer Physics Communications, Vol.66, pp.219-232, 1991.

6) Yabe T. : Unified solver CIP for solid, liquid and gas, Computational Fluid Dynamics Review 1997, 1997. 7) Nakamura, T., Tanaka, R., Yabe, T. and Takizawa, K. :

Exactly conservative semi-Lagrangian scheme for multi-dimensional hyperbolic equations with directional splitting technique, J. Computational Physics, Vol.174, pp.171-207, 2001.

8) Xiao, F., Honma, Y. and Kono, T. : A simple algebraic interface capturing scheme using hyperbolic tangent func-tion, Int. J. Numer. Meth. fluids, Vol.48, pp.1023-1040, 2005.

9) Martin, J.C. and Moyce, W.J. : An experimental study of the collapse of liquid columns on a rigid horizontal plane, Philos. Trans. Roy. Soc. London Ser. A., Vol.224, pp.312-324, 1952.

10) Hu, C. and Sueyoshi, M. : Numerical simulation and ex-periment on dam break problem, J. Marine. Sci. Appl., Vol.9, pp.12-49, 2010.

11) 川崎浩司，山口 聡，袴田充哉，水谷法美，宮島正 悟：段波と矩形物体の衝突・漂流過程における作用 波圧特性，海岸工学論文集，第 53 巻，pp.786-790， 2006.

INTRODUCTION OF THINC/WLIC METHOD TO THREE-DIMENSIONAL

MULTIPHASE FLOW NUMERICAL MODEL "DOLPHIN-3D"

Koji KAWASAKI and Tetsuya MATSUNO

The purpose of this study is to improve mass conservation and accuracy for interface capturing in a three-dimensional multiphase flow numerical model "DOLPHIN-3D" by introducing a THINC/WLIC method. The validity of the THINC/WLIC method was firstly verified through 1D, 2D and 3D advection tests.The utility of the multiphase flow model developed here was comfirmed by applying it to dam break problem and collision problems between water collapse-induced bore and a rigid body. The results re-vealed that the model has sufficient mass conseravation and high accuracy for interface capturing.

実験結果 CIP 法 THINC/WLIC 法

(a) t=0.5s

実験結果 CIP 法 THINC/WLIC 法

(b) t=0.8s