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

非構造格子に対する HRIC 法を応用した自由界面捕獲スキームの開発と実現象への適用
Development of interface capturing scheme based on HRIC method for unstructured meshes and theapplication to actual phenomena⃝

八登 浩紀，(株) ソフトウェアクレイドル, 大阪市北区梅田 3-4-5, E-mail： yatou@cradle.co.jp
吉田 一，(株) ソフトウェアクレイドル, 大阪市北区梅田 3-4-5, E-mail： yoshida@cradle.co.jp

Hiroki YATOU, Software Cradle Co., Ltd., 3-4-5, Umeda Kita-ku, Osaka
Hajime YOSHIDA, Software Cradle Co., Ltd., 3-4-5, Umeda Kita-ku, OsakaWe develop a new interface capturing scheme based on HRIC method for free surface flow. Bounded-Downwind
scheme as a compressive scheme and SMART scheme as a high-resolution scheme are blended depending on theangle between mesh and free surface. Unlike HRIC method, the dependency of the spatial accuracy of convective
term by the Courant number is not considered. High accuracy is achieved in the validation of some free-surface
analyses by using this scheme.

えば 1 以上) においても自由界面を精度良く追跡することが可能となる.

序論
近年, コンピュータ環境の性能向上, とりわけ複数のコアや CPU を持つ計算機環境が一般化してきていることから, 数値流体力学 (CFD) の解析による設計の効率化がますます広く行われるようになっている. その応用範囲は, 自動車や建築、電子機器、ターボ機械など多岐にわたる.
そのなかでも自由界面を有する気液二相流れは, 産業分野において非常に重要なトピックスである. たとえば,航行中の船舶周りの流れ, 水力・波力発電, 化学プラントにおける液体の貯蔵や輸送, ヒートパイプ内の液体の蒸発といった応用例が挙げられる. また, 津波のような自然現象における気液二相流も注目される.
このような広い応用範囲を持つ自由界面問題に対して,さまざまな界面捕獲のための数値計算法が提案されてきた. その手法は主に, (1) 界面に粒子を配置するなどして界面の運動をラグランジュ的に追跡する手法, (2) 界面を識別するためのスカラー値に対する移流方程式をオイラー的に計算する手法, に大別される.
後者の代表的な手法に VOF 法 (1) が挙げられる. VOF法は VOF 値 (Volume Of Fluid: 流体体積率) という界面識別のためのスカラー値の時間発展式を解くことにより,自由界面の捕獲を行う. この手法は保存性が良いことで知られているが, 同時に界面が拡散して不鮮明になりやすいという特徴があるため, 実現象に適用するためには相応の工夫が必要となる. 具体的には, I. 要素内の界面形状を考慮することでより精度の高い移流計算を行う手法(たとえば PLIC 法 (2) ), II. VOF 法を Level Set 法 (3) と組み合わせることで適切な界面の再初期化を行い界面の拡散を抑える手法 (CLSVOF 法 (4) ）, III. 代数的な計算のみでも界面が拡散しないような移流スキームを用いる手法（たとえば HRIC 法 (5) ）, が挙げられる. I, II の手法は計算精度には優れているが, 非構造格子系に適用する場合に計算アルゴリズムは複雑なものとなり, 一般にその計算負荷は小さくない. このため, 実用的な解析においては III の手法が用いられることがほとんどである.
本稿では, III の代数的な手法である HRIC 法を応用した新しい計算スキームを提案する. このスキームでは,HRIC 法と同様に, 界面の拡散を抑えるスキームと, 高精度スキームを組み合わせる処理を行う. 界面の拡散を抑えるスキームとしては Bounded-Downwind 法 (6) , 高精度スキームとして SMART 法 (7) を用いる. 一方, HRIC法で行われる高クーラン数における移流スキームの一次風上化は行わない. これにより, 大きいクーラン数（たと1.

数値計算手法SCRYU/Tetra はハイブリッド要素（四面体, 五面体,六面体）からなる非構造格子を用いた汎用流体解析ソフトウェアである. ソルバーのほかに, 自動格子生成プログ
ラムと可視化プログラムを備える. ソルバーでは節点中心型の有限体積法を採用している. 基礎方程式には 3 次
元圧縮性/非圧縮性 Navier-Stokes 方程式にくわえ, エネ
ルギー保存式や乱流エネルギー式などを必要に応じて計算する. これらの式に対しては, 移流項の離散化には二次精度の風上差分（MUSCL）(8) , 圧力補正には SIMPLEC法 (9) を用いる. 時間進行は陰解法で行う.

2.

2.1 VOF 法
VOF 法は, 自由界面を有する二相流解析のために広く用いられる計算手法である. VOF 法では, VOF 値 (本稿では F で表す) とよばれるスカラー値を用いて二相を表現する. ここでは, VOF 値の値に応じて, 以下の状態を表すものとする.
相1
 F =1(1)0 < F < 1 界面
F =0相2
この VOF 値に対して, 以下の移流方程式を計算する.

∂F∂F+ ui=0∂t∂xi(2)

ここで, u は流速を表す. これにより F の時間発展が求まり, その結果として, 0 < F < 1 の領域に相当する界面の時間発展を追跡することができる.

2.2 HRIC 法
VOF 法では VOF 値の分布を用いて界面を表現するため, 界面は有限の厚みを持つ. しかし, 式 (2) の移流方程式を解くあいだに, 数値拡散によりこの界面厚みは次第に大きくなり, 界面は徐々に不鮮明となる, このような不鮮明な界面では精度良く界面を追跡できないため, この界面厚みの増大を抑制するような移流項の離散化スキームを用いる必要性が生じる.
そのような手法として考案されたものに HRIC 法や
CICSAM 法 (10) がある. これらの手法では, 界面厚みを1c 2013 by JSFM
Copyright ⃝第 27 回数値流体力学シンポジウム
C08-4

圧縮するようなスキーム (Compressive Scheme) と高精度スキーム (High-Resolution Scheme) を組み合わせることで界面幅が増大することを抑える（これらのスキームを本稿では混合型移流スキームと記述する）. HRIC 法では, Compressive Scheme として Bounded-Downwind 法,
High-Resolution Scheme として一次風上が用いられる.
具体的な計算手順は以下のとおりである. まず, 以下の記号を導入する.

F̃ =F − FU
FD − FU(3)

ここで, FD は下流側のコントロールボリューム (CV) のVOF 値, FU は上流側の VOF 値を表す.
VOF の移流計算を高精度に行うためには, CV 界面における VOF 値 (F̃f ) を精度よく評価する問題に帰着する.
HRIC 法においては, CV 界面の VOF 値は以下のように決定される.

F̃f = γ F̃f (comp) + (1 − γ)F̃f (HR)(4)

ここで, F̃comp は Compressive Scheme により評価した
CV 界面における VOF 値, F̃f (HR) は高精度スキームにより決定した CV 界面における VOF 値を表す.
重み γ は, 要素面の法線と自由界面の法線とのなす角を θ としたとき, 以下のように決定される.

γ=√cos(θ)

Fig. 1: Time evolution of free surface of dam breakanalysis.

(5)

HRIC 法ではさらに, 局所的なクーラン数がしきい値(Co = 0.7) を越えた場合にも一次風上に帰着するような処理が行われる. このように決定した CV 界面の VOF 値を用いて, CV 間の VOF 値の移流量が決定される.
2.3 HRIC 法を応用した混合型移流スキームの提案
HRIC 法は界面の拡散を抑えるための自由界面捕獲スキームであるが, クーラン数が大きい場合に一次風上に帰着することから, 界面付近でクーラン数が小さくなるように解析を行うことが必須となる.
一方, 非構造格子系における汎用解析ソフトウェアを用いた解析では, 解析モデルの複雑さから要素数が大きくなるケースも多く, 小さなクーラン数を維持した解析を現実的な解析時間で行うことが難しい場合もしばしばである.
そこで本稿では, HRIC 法をベースに, より汎用解析ソフトウェアに適した混合型移流スキームを提案する. 具
体的には, High-Resolution Scheme として SMART 法を採用する. 加えて, HRIC 法で行われる高クーラン数における一次風上への切り替えも行わない. 一方, Compressive Scheme については HRIC 法と同様に BoundedDownwind 法を用い, 重みについても同様に (5) 式により算出する. SMART 法は MUSCL と同様に二次精度の風上差分であり, 一次風上より高精度でかつ界面の拡散を抑えた解析を行うことができる. このスキームを BoundedDownwind 法と組み合わせることで, 精度を保ちながら高クーラン数においても VOF 値の拡散が起きにくい, 界面捕獲に適した移流計算を行うことが可能となる.

3.

2.4 その他の計算条件
表面張力は, CSF モデル (11) により運動方程式における体積力として考慮される. また, 界面での密度は VOF
値の流体体積率に従った二相の算術平均により求める. 粘性係数については二相の調和平均または算術平均により算出する.

4.

Fig. 2: Comparison between calculation and experiment results of dam break problem.

基礎検証
自由表面解析の検証問題としてよく知られているダム崩壊問題を用いて, 本稿で提案する混合型移流スキームの検証解析を行う. ダム崩壊問題では, 初期に静止した液柱が液柱自身の重力により崩壊する現象を扱う. 本解析
で液注は初期に高さ 29.2cm, 幅 14.6cm とする. クーラン数は Co = 1.0 に設定する.
Fig.1 は, 自由表面の時間発展を F = 0.5 の等値面を描くことによって可視化したものである. 計算初期において, 液柱の前線は時間とともに x 方向に移動する. この前線の x 座標と時間の関係について, 数値解析と実験の結果 (12) を比較したものが Fig.2 である. 図のように, 解析
結果は実験の結果とよく一致しており, 本稿で提案するスキームが高い精度で解析できていることを示している.

解析事例
本節では, 新しい混合型移流スキームを用いたいくつかの計算事例について述べる.

2c 2013 by JSFM
Copyright ⃝第 27 回数値流体力学シンポジウム
C08-4

Fig. 3: Contour map of wave height around vessel body.

Fig. 5: Computational mesh around moving object.
Blue line denotes water surface.
動方向を変えるような解析を行うことができる.
後者の解析の例として, 本節では物体の水面への落下の解析結果を示す. 物体の大きさは, 直径 10cm, 高さ 5cm であり, 物性としては木材を模擬し, 密度を ρ = 300kg/m3で与える. 物体は重力や流体力を受け, 任意の方向に加速されながら運動する.
Fig. 5 は, 初期状態での計算格子を表している. 四角形（五面体要素）で表示される格子は静止空間（独立領域）の格子であり, 三角形（四面体要素）で表される格子は落下物体周辺の領域（従属領域）の格子である. 独立領域の格子は移動せず, 従属領域の格子は移動物体の速度に追随して移動する. 本解析での計算格子の総要素数はおよそ 29 万要素である.
Fig. 6 は, 重力により自由落下した物体が水の中に入るとともに水からの浮力をうけ, 安定な体勢に移行したのち水面に浮かんだ状態になるまでを時系列で示している. 重合格子機能と混合型移流スキームを併用することで, 物体が自由界面を有する流体と相互作用する様子が定性的によく捉えられる. なお, 本解析における計算時間は 12 コアを用いて 40 分弱であり, このような複雑な現象についても, 自由表面の大きな散逸を伴わずに, 十分短い時間で解析することができる.
本解析では物体形状については簡略化されたものを用いているが, 形状をたとえば前節の船舶のモデルに置き換えて船舶が波に揺られる様子を解析するなど, より実用的な系への応用も期待できる.

Fig. 4: Comparison between experiment and calculation results of wave height near vessel body. The calculations are performed with HRIC method and the newmethod we propose in this article.

4.1 船舶の解析
航行中の船舶まわりの水面のゆらぎに関する検証解析を行う. 船舶のモデルとして DTMB5415 を用いる. このモデルは自由表面解析のベンチマークに広く用いられており, 解析結果との比較のための各種実測データも公開されている (13) .
本検証解析では, 船体自体は静止しており, 船体前方より一様な流れを与えた解析を行う. 船体の長さ（垂線間長）は Lpp = 5.72m, 流れの速度は x 方向に 2.1m/s とする.
Fig. 3 は自由界面形状がほぼ定常的になった段階での解析結果のスナップショットである. 水の流れが船体にぶつかることによる船体まわりでの水面の揺らぎが確認できる. Fig. 4 は, この波の高さについて, 上記の公開実測データと新手法による結果, ならびに HRIC 法による数値計算結果の比較を行ったものである. この結果は,y/Lpp = 1.72 の平面上での波の高さをプロットしている.
HRIC 法により得られる波高は, 中央のピーク付近において波の高さを小さく見積もっていることに対して, 本稿で提案するスキームでは, より実測に近い波高が得られており, 本手法の有効性が示されている.

結論
非構造格子系汎用流体解析ソフトウェア SCRYU/Tetraの自由表面解析機能に, HRIC 法を応用した混合型移流スキームを組み込んだ. この手法は従来手法と比較して高クーラン数においても精度を保ちかつ界面幅の増大を抑えることを可能とする. 自由表面流れの検証問題としてよく知られているダム崩壊について解析を行い, 本手法により実験に近い界面の時間発展が得られることが確認された. 加えて, 船体まわりの流れといった応用的な事例についても実測をよく近似する結果が得られた. また,本稿で例示したように, 重合格子との併用など, 移動物体を含む多様な流れの解析への応用が期待できる.

5.

4.2 浮体の解析
SCRYU/Tetra は重合格子機能を有し, 静止した空間と移動物体周辺の空間に異なる格子を配置し, それらを重ね合わせた解析を行うことができる. これにより, 形状再現性を損なわずに自由度の高い移動物体解析を行うことが可能となる. たとえば, 回転領域同士が重なっている複数枚の羽根車の解析, 物体が流体力を受けて不規則に運参考文献(1) Hirt, C. W. and Nichols, B. D., “Volume of Fluid(VOF) Method for the Dynamics of Free Boundaries,” J. Comp. Phys., 39 (1981), pp. 201-225.

(2) Youngs, D. L., “Time dependent multimaterial flow3c 2013 by JSFM
Copyright ⃝第 27 回数値流体力学シンポジウム
C08-4(9) Patankar, S. V., “Numerical Heat Transfer andFluid Flow,” Hemisphere Publishing, New York(1980).
(10) Ubbink, O. and Issa, R., “A Method for Capturing Sharp Fluid Interfaces on Arbitrary Meshes”,
J. Comp. Phys., 153 (1999), pp. 26-50.
(11) Brackbill, J. U., Kothe, D. B. and Zemach, C.,“A Continuum Method for Modelling Surface Tension”, J. Comp. Phys., 100 (1992), pp. 335-354.
(12) Martin, J. C. and Moyce, W. J., “An Experiment
Study of the Collapse of Liquid Columns on a Rigid
Horizontal Plane,” Philos. Trans. R. Soc. London
Ser. A, 244 (1952), pp. 312-324.
(13) Gothenburg 2000 A Workshop on CFD in Ship Hydrodynamics. http://www.nmri.go.jp/institutes/
fluid performance evaluation/cfd rd/cfdws05/gothenburg2000/index.html

Fig. 6: Time evolution of moving object falling to water
surface.
with large fluid distortion”, in Morton, K. M. andBaines, M. J. (eds), Numerical Methods for Fluid
Dynamics, Academic Press, New York (1982), pp.
27-39.
(3) Sussman, M. Smereka, P. and Osher, S., “A levelset approach for computing solutions to incompressible two-phase flows,” J. Comp. Phys., 114(1994), pp. 146-159.
(4) Sussman, M. and Puckett, E. G., “A coupled levelset and volume-of-fluid method for computing 3D
and axisymmetric incompressible two-phase flows,”
J. Comp. Phys., 162 (2000), pp. 301-337.
(5) Muzaferija S., and Peric M., “Computation of freesurface flows using interface-tracking and interfacecapturing methods” in Mahrenholtz O, Markiewicz
M. (eds.), Nonlinear Water Wave Interaction, Computational Mechanics Publications, Southampton(1998).
(6) Darwish, M., Moukalled, F., “Convective Schemes
for Capturing Interfaces of Free-Surface Flows on
Unstructured Grids,” Numerical Heat Transfer,Part B., 49 (2006), pp. 19-42.
(7) Darwish, M., Moukalled, F., “Normalized variable and space formulation methodology for highresolution schemes,” Numerical Heat Transfer, PartB., 26 (1994), pp. 79-96.
(8) van-Leer, B., “Towards the Ultimate ConservationDiﬀerence Scheme. 4, A New Approach to Numerical Convection,” J. Comp. Phys., 23 (1977), pp.
276-299.
4c 2013 by JSFM
Copyright ⃝