第 25 回数値流体力学シンポジウム
D11-4非構造格子における界面捕獲手法の構築
Development of an interface capturing method on an unstructured mesh
伊井 仁志, 東大, 文京区本郷 7-3-1, sii@fel.t.u-tokyo.ac.jp:
肖 鋒, 東工大, 横浜市緑区長津田町 4259, xiao@es.titech.ac.jp:

Satoshi Ii, The University of Tokyo, 7-3-1 Hongo Bunkyo-ku
Feng Xiao, Tokyo Institute of Technology, 4259 Nagatsuta, Midori-ku
An interface capturing method is newly developed in a philosophy of a smoothed volume-of-ﬂuid (VOF) approach,on the unstructured meshes, i.e., triangular mesh and tetrahedral mesh. The present method, based on the tangentof hyperbola interface capturing scheme (Xiao et al., Int. J. Numer. Meth. Fluids, 48 (2005) 1023-1040), employsa multi-dimensional hyperbolic tangent (or sigmoid) function with a discretized volume fraction data, thus the
VOF function keeps a smooth proﬁle over a few computational meshes, oﬀering that a derivative value of the
VOF function is directly obtained from itself. The numerical results show that the present method maintains an
initial VOF proﬁle well, and also infer a possibility to a multi-phase ﬂow simulation on the unstructured meshes.

はじめに
自由界面の大変形や大移動，またそれに伴うトポロジー
変化が頻繁に起こるような多相流問題では，それらを容易に扱える固定格子を用いた解析手法が有用である．ここで固定格子において界面を陰的に取り扱う界面捕獲手法として，例えば，流体の体積占有率を用いた VOF(Volumeof-ﬂuid) 法 (5) ，距離関数を用いた LS(Level-set) 法 (9,17) ，また VOF 法と LS 法と組み合わせた CLSVOF(Coupledlevel-set and volume-of-ﬂuid) 法 (18) ，密度関数を用いた密度関数法 (22,20) ，また相場関数を用いた Phase-Field 法(7)などが挙げられる．本研究では，質量保存性 (界面に囲まれた領域の体積保存性とは別に定義) に優れた VOF 法に関して新たな手法を提案する．VOF 法で幅広く用いられている PLIC(Piecewise linear interface calculation) 法(25,11,12,14,10)
では，離散変数である VOF 関数を用い，各計算要素において線形界面形状を持つヘビサイド関数が近似的に構築される．ヘビサイド関数の特異性により近似関数構築は幾何的な操作を必要とし複雑化するが，線形界面であれば解析かつ容易に関数が構築される (13) ．一方，
このような幾何的操作を避けるために，滑らかな VOF
関数 (あるいは各相を表す指標関数) を導入し，各計算要
素での指標関数を tanh 関数で近似した THINC(Tangentof hyperbola interface capturing) 法 (19) が提案されており，その発展手法 (24,21) も提案されている．また多次元において，界面を準二次曲面と仮定し，多次元 tanh 関数を用いた MTHINC(Multidimensional THINC) 法 (6) が近年著者らのグループで開発されている．
ここで，上記で述べた先行研究らは直交格子において定式化がなされているが，実用問題においては，複雑形状を精度良く扱うことができる非構造を用いた定式化も今後重要となると考えられる．PLIC 法で用いられる幾何的構築を基礎とした三角形あるいは四面体での定式化
は，例えば，文献 (15,23) などに見られるが，幾何的構築に伴うアルゴリズムの複雑化により線形界面を仮定している．そこで本研究では，三角形格子また四面体格子において，将来的に高次界面 (曲面) を用いた高精度界面捕獲手法の開発を目指し，直交格子において定式化された
MTHINC 法を非構造拡張に拡張する．その初段階として，線形界面を仮定し手法を構築し，いくつかの移流テストにより提案手法の妥当性を検証する．

1.

2.

非構造格子における界面捕獲手法

2.1 滑らかな指標関数により定義される VOF 関数
本研究では，気相と液相は混じり合わないとし，それらは連続体領域 (あるいは計算要素幅) で滑らかに遷移していると仮定する．その際，全領域 Ω を液相領域 Ωl ，気相領域 Ωg および遷移領域 Γ に分割し，各相とその遷移
1

領域を表す指標関数 (indicator function): χ(x) を以下のように導入する．




1,


χ(x) =x ∈ Ωl ,0,x ∈ Ωg ,
0 < χ(x) < 1, x ∈ Γ,

(1)

R3 = {(x, y, z) | x ∈ Ω = Ωl ∪ Ωg ∪ Γ} .
この時，検査領域 δV (x) における液相の体積占有率(volume-fraction function: VOF 関数): φ(x) は以下のように定義される.

φ(x) ≡1|δV (x)|∫δV (x)

χ(x0 )dx0 = hχ(x0 )iδV (x)(2)

2.2 VOF 関数の半離散時間発展式
本研究では，気液界面を VOF 関数により陰的に表現
する．各相が混じり合わない場合，χ(x) の時間発展式は
以下のように与えられる．
∂χDχ
= 0, →+ ∇ · (vχ) = χ∇ · v,
Dt
∂t

(3)

ここで各計算要素 i (i = 1, ..., Ne ) での検査領域を δViとすると，定義式 (2) および式 (3) より，VOF 関数 φi の時間発展式は，∂φi1+∂t|δVi |

Iφi(vn χ)dS =
|δV
i|
δSi

I
vn dS,(4)

δSiと与えられる．ここで，vn = v · n は法線方向 n の射影H速度，また δS dSi は検査領域 δVi での周回積分を表す．
式 (4) より，非圧縮条件 ∇ · v = 0 では φ の保存性は保証されるが，通常，非圧縮条件は数値的な誤差を伴い満たされるため右辺項は残しておく．計算要素を構成する面要素 (二次元三角形では 3 つの直線，三次元四面体では 4 つの三角形) での面積分から周回積分を評価し，数値積分を適用すると空間方向に離散化された時間常微分方程式が得られる．本研究では，面要素の数値積分に四次精度ガウス求積法 (2) また時間積分には三次精度 TVD
ルンゲ・クッタ法 (16) を適用する．

Copyright c 2011 by JSFM

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

2.3 近似関数の構築
VOF 関数を用いた界面捕獲法においては，各計算要素領域において式 (4) の左辺第 2 項に含まれる χ(x) を離散
変数 φi (i = 1, .., Ne ) から如何に精度良く近似するかが重要な点である．例えば PLIC 法 (25,11,12,14,10) では，各計算要素領域 δVi において線形界面を仮定し，特異点を含む多次元ヘビサイド関数 Hi (x) により χ(x) を近似する．すなわち，{x ∈ δVi : χ(x)} = Hi (x)．この際，ヘビサイド関数の特異性により近似関数構築に幾何的な操作を必要とするが，直交格子での解析的表現 (13) ，また三角形および四面体格子での解析的表現 (23) が提案されており，線形界面での近似関数が容易に構築される．ただし曲面を考慮した際，幾何的操作によりアルゴリズムはとたんに複雑になり (4) ，一般的に解析的構築はほぼ不可能となる．一方，MTHINC 法 (6) では連続関数である多次元シグモイド関数を用いるため幾何的な操作を必要とせず，曲面に関しても同様の手順で近似関数の構築を行うことができる．よって，将来的に非構造格子において曲面を持った界面を取り扱う際に有用と考えられる．
本研究では初段階として，各計算要素において界面は線形であると仮定する．その際，各計算要素における近似関数 Ĥi (x) を以下のように与える．これ以降，簡単化のため要素インデックス i を省略して記述する．
1
Ĥ(x) = (1 + tanh(β(nx x + ny y + nz z + c))) . (5)2
ここで，n = (nx , ny , nz ) は界面の単位法線ベクトルで後述の評価法より求める．また β = β 0 /have と定義され，
have は以下のように定義した平均格子幅である．
e1 ∑(6|δVi |)1/3 .
Ne i=1N

have =(6)

Ne は全要素数，また β 0 は VOF 関数の遷移領域幅を決定するパラメータであり，本研究では文献 (21) およびい
くつかの数値テストにより β 0 = 2.3 とする．
単位法線ベクトル n は要素中心での値 nc により近似する．これらは Youngs らの方法 (25) を非構造格子に適用
し，要素頂点での VOF 関数の勾配を隣接する離散 VOF関数から最小二乗法により求め，要素中心での平均勾配√∇φc から nc = ∇φc / |∇φc |2 + 10−16 と評価する．本手法では，離散 VOF 関数が格子幅に対してある程度の遷移領域を持ち滑らかに分布していることを利用し，VOF関数の勾配から直接法線ベクトルを求めているが，空間収束性に問題があることが知られているため (3) ，高精度な法線ベクトル評価法が今後の課題と考えられる．
近似関数の構築にあたり，四面体格子において局座標
系 (ξ, η, ζ) を導入する (Fig.1)．基準点を x0 とした時，

d = β(n · x0 + c),

(11)

(ξ)

Ĥ

(1/4, 1/4) = φ,Case2 : aη = max(aξ , aη , aζ ),
(12)

(η)

Ĥ

(1/4, 1/4) = φ,Case3 : aζ = max(aξ , aη , aζ ),
(13)

(ζ)

Ĥ

(1/4, 1/4) = φ.

ここで，Case1 に注目すると以下の手順より d を逆問題
(11) から求めることができる．
()
1
1 − AQd = ln,2
AB(Q − A)(7)

A = exp(aξ /2),

(14)

B = exp((aη + aζ )/2),となり，それに伴い近似関数 (5) は局座標系における近似関数 Ĥ(ξ, η, ζ) に置き換えられる．

aξ = βn · xξ , aη = βn · xη , aζ = βn · xζ ,VOF 関数の定義式 (2) より，近似関数 (8) の体積積分値から逆問題として d が決定されるが，tanh(ξ, η, ζ) の
多次元積分が初等関数では表現できないため，解析積分と１点数値積分を併用する．ここで，∫ 1−η−ζ(ξ)
1Ĥ (η, ζ) ≡Ĥ(ξ, η, ζ)dξ,
1−η−ζ 0∫ 1−ξ−ζ(η)1(10)
Ĥ (ξ, ζ) ≡Ĥ(ξ, η, ζ)dη,
1−ξ−ζ 0∫ 1−ξ−η(ζ)
1Ĥ(ξ, η, ζ)dζ,Ĥ (ξ, η) ≡1−ξ−η 0
Case1 : aξ = max(aξ , aη , aζ ),(0 ≤ ξ, η, ζ ≤ 1, ξ + η + ζ ≤ 1),

0 ≤ ξ, η, ζ ≤ 1, ξ + η + ζ ≤ 1 :1Ĥ(ξ, η, ζ) = (1 + tanh(aξ ξ + aη η + aζ ζ + d)) .
2
ここでそれぞれ，

Fig. 1: Schematic of a local coordinate system on atetrahedral mesh.

を定義し，以下のような場合分けにより d を決定する．

xξ = x1 − x0 , xη = x2 − x0 , xζ = x3 − x0 ,x = ξxξ + ηxη + ζxζ ,である．

Q = exp(aξ (φ − 1/2)).
他の場合 (12)，(13) に関しても，係数を順次入れ替えてd を求めることができる．
式 (14) より，φ = 0 あるいは φ = 1 においては定義が成り立たないため，本研究では微小数  = 10−8 を導入
し，以下のように計算要素 i での指標関数を近似する．
Ĥ (x), for 0 ≤ φ ≤ 1,ii{x ∈ δVi : χ(x)} =
(15)
φ ,
otherwise.

(8)

(9)

i

2

Copyright c 2011 by JSFM

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

移流テスト
本研究では，計算要素での CFL(Courant-FriedrichsLewy) 数を以下のように定義する．
 max(|ui |,|vi |)∆t ,in 2D,(2|δSi |)1/2
CFLi =
 max(|ui |,|vi |,|wi |)∆t , in 3D,(16)(6|δV |)1/3

3.

i(i = 1, ..., Ne ).
い く つ か の 数 値 テ ス ト よ り，CFLmax=
max(CFLi )i∈[1,Ne ] が 0.25 以 下 に お い て 十 分 な 精 度を保つこと確認しており，この基準を満足する時間増分∆t を用いる．

(a) Ne = 1282 × 2

3.1 三角形格子
Delaunay 分割により作成された三角形格子において数値テストを行う．ここで，全要素数を直交格子と対比させるため Ne = N 2 × 2 と表記する．
スロットを持つ円盤の剛体回転 (26)
計算要素 Ne = 712 × 2 ≈ 10000 において CFLmax =0.25 のもと計算を行う．初期結果，1 回転後また 5 回転後の計算結果を Fig.2 に示す．ここでは剛体部を拡大して表示している．計算結果より，1 回転後において界面形状を良く保持していることが確認できる．5 回転後の結果では，スロット部において回転方向に依存する誤差が見られるが，VOF 関数の数値拡散が抑えられているため体積保存性は良く保たれている．

(b) Ne = 2562 × 2

Fig. 3: Numerical results in a 2D time-dependent single
vortex ﬂow at t = 0 (left), t = T /2 (middle) and t = T(right) in case where T = 8. The contour lines of the
VOF function are shown at 0.05 (blue), 0.5 (black) and
0.95 (red), respectively.

Tab. 1: Numerical errors and convergence rates for the
spatial resolution in a 2D time-dependent single vortex
ﬂow.

Fig. 2: Numerical results in a slotted-disk rotation on atriangular mesh. Black line: initial, red line: 1 rotation,
blue line: 5 rotations.

Ne

642 × 21282 × 2

2562 × 2

Error

1.59E-2

3.13E-3

6.03E-4

Order

–

2.34

2.37

二次元の時間反転渦度場における形状再現性 (1,12)
計算要素 Ne = 642 × 2，1282 × 2 また 2562 × 2 において CFLmax ≈ 0.17 のもと計算を行う．ここで，T = 8
を 1 周期とする時間反転速度場を与え，初期形状の再現性を検証する．Fig.3 に計算結果を示す．t = T /2 において界面間の厚みは非常に薄くなり空間解像度が厳しくなるため，t = T において Ne = 1282 × 2 要素では初期円
形状が若干崩れている．Ne = 2562 × 2 要素を用いた計算結果 (Fig.3(b)) より，これらは空間解像度を上げることにより改善できていることが確認できる．また，Fig.3
において φ = 0.05，0.5 また 0.95 での等高線を示しており，本手法は VOF 関数の時間的な散逸を抑制できていることが確認できる．ここで，空間解像度に対する L1 ノルムによる初期解との相対誤差を Tab.1 に示す．本問題においては，空間二次精度の収束性を示している．

Fig. 4: The regular tetrahedral meshes divided from the23 cubic meshes.

3.2 四面体格子
次に四面体格子において数値テストを行う．ここで，四面体格子は三次元直交格子を 5 分割した規則性のある格子を用いる (Fig.4)．その際，全要素数は Ne = N 3 × 5と表記する．
3

Copyright c 2011 by JSFM

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

三次元の時間反転渦渦度場における形状再現性 (8)計算要素 Ne = 303 × 5，603 × 5 また 1203 × 5 において CFLmax ≈ 0.12 のもと計算を行う．本テストでは，
T = 3 を 1 周期とする時間反転速度場を与え，初期形状の再現性を検証する．Fig.5 に t = 0，T /2 また T での計算結果を示す．二次元問題と同様に，t = T /2 において界面間の厚みは最も薄くなり空間解像度が厳しくなるが，
空間解像度を上げることで改善できており，それに伴い初期球形状の再現性も増している．更に，空間解像度に対する L1 ノルムによる初期解との相対誤差 (Tab.2) より，本問題においても空間二次精度が実現されている．

(4) S.V. Diwakar, S.K. Das, T. Sundararajan, A
Quadratic Spline based Interface (QUASI) reconstruction algorithm for accurate tracking of twophase ﬂows, J. Comput. Phys., 228 (2009) 91079130.
(5) C.W. Hirt, B.D. Nichols, Volume of ﬂuid (VOF)
method for the dynamics of free boundaries, J.
Comput. Phys., 39 (1981) 201-225.
(6) S. Ii, K. Sugiyama, S. Takeuchi, S. Takagi, Y. Matsumoto, F.Xiao, An interface capturing methodwith a continuous function: the THINC methodwith multi-dimensional reconstruction, J. Comput.
Phys., revised (2011).
(7) D. Jacqmin, Calculation of two-phase NavierStokes ﬂows using phase-ﬁeld modeling, J. Comput.
Phys., 1455 (1999) 96-127.
(8) R. LeVeque, High-Resolution Conservative Algorithms For Advection In Incompressible Flow,
SIAM J. Numerical Analysis, 33 (1996) 627-665.

(a) Ne = 603 × 5(9) S.J. Osher, J.A. Sethian, Fronts propagating withcurvature dependent speed. Algorithms based on
Hamilton-Jacobi formulations, J. Comput. Phys.,
79 (1988) 12-49.
(10) J.E. Pilliod Jr., E.G. Puckett, Second-order accurate volume-of-ﬂuid algorithms for tracking material interfaces, J. Comput. Phys., 199 (2004) 465502.

(b) Ne = 1203 × 5(11) E.G. Puckett, A.S. Almgren, J.B. Bell, D.L. Marcus and W.J. Rider, A High-Order Projection
Method for Tracking Fluid Interfaces in Variable
Density Incompressible Flows, J. Comput. Phys.,
130 (1997) 269-282.

Fig. 5: Numerical results in a 3D time-dependent vortex
ﬂow at t = 0 (left), t = T /2 (middle) and t = T (right)
in case where T = 3.

(12) W.J. Rider, D.B. Kothe, Reconstructing volume
tracking, J. Comput. Phys., 141 (1998) 112-152.

Tab. 2: Numerical errors and convergence rates for the
spatial resolution in a 3D time-dependent vortex ﬂow.
Ne

303 × 5603 × 51203 × 5

Error
Order

7.14E-3–

2.18E-3
1.71

5.31E-4
2.04(13) R. Scardovelli, S. Zaleski, Analytical relations connecting linear interfaces and volume fractions inrectangular grids, J. Comput. Phys., 164 (2000)
228-237.
(14) R. Scardovelli, S. Zaleski, Interface reconstruction
with least-square ﬁt and split Eulerian Lagrangian
advection, Int. J. Numer. Meth. Fluids, 41 (2003)
251-274.

まとめ
本研究では，格子サイズで分布する滑らかな離散 VOF関数を導入し，連続関数である多次元 tanh 関数により計算要素での指標関数を近似する多次元 THINC 法(MTHINC 法) を非構造格子に拡張した．本提案手法は
界面捕獲手法で重要な点である; 1)VOF 関数の時間的な散逸の抑制; 2) 体積保存性，に対し優れた計算結果を示すことができた．今後は，流体解析と連成させその有用性を検討していく．

4.

(15) K. Shahbazi, M. Paraschivoiu, J. Mostaghimi,
Second order accurate volume tracking based on
remapping for triangular meshes J. Comput. Phys.,
188 (2003) 100-122.
(16) C.W. Shu, Total-variation-diminishing time discretizations, SIAM J. Sci. Stat. Comput., 9 (1988)
1073-1084.
(17) M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible
two-phase ﬂows, J. Comput. Phys., 114 (1994)
146-159.

参考文献(1) J.B. Bell, P. Colella, H.M. Glaz, A second-orderprojection method for the incompressible NavierStokes equations, J. Comput. Phys., 85 (1989) 257283.

(18) M. Sussman, E.G. Puckett, A Coupled Level
Set and Volume-of-Fluid Method for Computing
3D and Axisymmetric Incompressible Two-Phase
Flows, J. Comput. Phys., 162 (2000) 301-337.

(2) G.R. Cowper, Gaussian Quadrature Formulas for
Triangle, Int. J. Numer. Meth. Eng., 7 (1973)
405-408.

(19) F. Xiao, Y. Honma, K. Kono, A simple algebraic interface capturing scheme using hyperbolic tangent
function, Int. J. Numer. Meth. Fluids, 48 (2005)
1023-1040.

(3) S.J. Cummins, M.M. Francois, D.B. Kothe, Estimating curvature from volume fractions, Comput.
Structures, 83 (2005) 425-434.
4

Copyright c 2011 by JSFM

第 25 回数値流体力学シンポジウム
D11-4(20) F. Xiao, A. Ikebata, An eﬃcient method for capturing free boundary in multi-ﬂuid simulations, Int.
J. Numer. Meth. Fluids, 42 (2003) 187-210.
(21) F. Xiao, S. Ii and C.G. Chen, Revisit to the THINC
scheme: a simple algebraic VOF algorithm, J.
Comput. Phys., 230 (2011) 7086-7092.
(22) T. Yabe, F. Xiao, Description of Complex and
Sharp Interface during Shock Wave Interaction
with Liquid Drop, J. Phys. Soc. Jpn., 62 (1993)
2537-2540.
(23) X. Yang, A.J. James, Analytic relations for reconstructing piecewise linear interfaces in triangular and tetrahedral grids, J. Comput. Phys., 214(2006) 41-54.
(24) K. Yokoi, Eﬃcient implementation of THINC
scheme: A simple and practical smoothed VOF algorithm, J. Comput. Phys., 226 (2007) 1985-2002.
(25) D.L. Youngs, Time-dependent multi-material ﬂow
with large ﬂuid distortion, in K.W. Morton and
M.J. Baines (eds), Numerical Methods for Fluid
Dynamics Academic, New York, 1982, 273-285.
(26) S.T. Zalesak, Fully multidimensional ﬂux-correctedtransport algorithms for ﬂuids, J. Comput. Phys.,
31 (1979) 335-362.

5

Copyright c 2011 by JSFM