VOF 法を用いた自由表面流れ解析における高精度解析手法の構築

VOF 法を用いた自由表面流れ解析における高精度解析手法の構築

中央大学大学院 学生員 ○ 弘崎 聡 日本工営株式会社 正会員 桜庭 雅明 中央大学 正会員 樫山 和男

1. はじめに

自由表面を有する複雑な流れを把握することは工学上重 要であり，数多くの数値解析手法がこれまで提案されてい る．その中で，砕波などを含む複雑な自由表面流れを解析 することに優位である代表的な手法として界面捕捉法の1 つであるVOF(Volume of Fluid)法¹⁾がある．VOF法は 流れ場の流速から体積占有率を表す界面関数（VOF関数）

を移流させることにより，自由表面位置を決定させる手法 である．しかし，VOF法は安定かつ高精度な解析を行うこ とに大きな問題点を抱えている．

そこで本研究は，従来のVOF法に対して流れ場の計算 には安定化有限要素法，移流方程式の解法にはCIVA法お よび新しい体積補正手法を導入して，自由表面流れ解析を 安定かつ高精度に行う手法を構築した．数値解析例として 3次元矩形貯水槽内スロッシング解析を取り上げ，本手法の 有効性を検討した．

2. 数値解析手法

非圧縮性粘性流体の基礎方程式であるNavier-Stokesの 運動方程式と連続式は以下のようになる．

ρ µ∂ui

∂t +ujui,j−fi

¶

−σij,j= 0 in Ω (1)

ui,i= 0 in Ω (2)

ここに，uiは流速成分，ρは密度，fiは物体力を表してい る．また，σijは応力テンソルであり，式(3)のように表す．

σ_ij =−pδ_ij+µ µ∂ui

∂xj +∂uj

∂xi

¶

(3)

p は 圧 力 ，µ は 粘 性 係 数 を 表 す ．基 礎 方 程 式 (1)，(2) に 対 し て SUPG/PSPG(Streamline Upwind Petrov Galerkin/Pressure Stabilizing Petrov Galerkin) 法 に 基 づく安定化有限要素法²⁾を適用し，P1/P1（流速・圧力1 次）要素を用いて空間方向の離散化を行った．また，時間 方向の離散化には運動方程式に対してCrank-Nicolson法，

圧力および連続式の項に対して陰的オイラー法を適用し ており，連立１次方程式の解法にはElement-by-Element Bi-CGSTAB法を用いた．

(2) 自由表面位置の計算

自由表面位置を追跡する手法として，界面捕捉法の1つ であるVOF法¹⁾を用いた．VOF法は，自由表面位置を表 現するVOF関数φ（液体であれば1，気体であれば0，自 由表面上であれば0.5）に関する移流方程式を解くことによ

り自由表面位置の履歴を追跡する手法である．

∂φ

∂t +uiφ,i= 0 in Ω (4) ここで，uiは流速成分である．なお，各節点における気体，

液体の密度と粘性係数は計算されたVOF関数値を用いて，

式(5)のように決定できる．

ρ=ρlφ+ρg(1−φ), µ=µlφ+µg(1−φ) (5)

ここで，ρl, ρg, µl, µg はそれぞれ液体，気体の密度および 粘性係数である．また，式 (4) の解法にはCIVA(Cubic Interporation with Volume/Area Coordinate)法³⁾を用い た．CIVA法は CIP法⁴⁾を3 角形要素または4 面体要素 へ拡張した手法であり，移流方程式の厳密解が式(6)とな ることから，上流側の三角形要素内のVOF関数を高次の 補間を行う手法である．3次元の場合の補間式は体積座標 (L1, L2, L3, L4)を用い式(7)のような3次関数で表すこと ができる．

φⁿ⁺¹(x, y, z, t) =φⁿ(x−u∆t, y−v∆t, z−w∆t, t−∆t) (6) φ(L1, L2, L3, L4) =

X4 i=1

αiLi+1 2

X4

j,k=1 (j6=k)

βjkLjLk(1+Lj−Lk)

(7) ここで，係数α, βは上流側の三角形要素での節点のVOF 関数φとその空間微係数を用いることにより決定され，以 下のようになる．

αi=φi, βjk=φk−φj+ X3

l=1

(xlj−xlk)∂φk

∂xl (8) (3) 体積補正手法

本体積補正手法は，界面近傍(0< φ <1)のVOF関数を 単純に増減させることで補正を行う．まず，自由表面近傍 (0< φ <1)であることを判定する関数として式(9)を用い て，解析領域内における自由表面近傍の体積を式(10)を用 いて計算する．

D(φ) = 1 + cos{2π(φ−0.5)} (9) A(t) =

Z Z Z

D(φ)dxdydz (10)

ここで，D(φ)は自由表面位置では2，自由表面から遠ざか るにつれて0に近づくような関数である．式(10)で算出し たA(t)を用いて，界面近傍に補正するVOF関数φerrは 次式のように計算することができる．

φerr= Verr(t)

A(t) = (V(t)−Vinit)

A(t) (11)

KeyWords： VOF法，自由表面流れ解析，安定化有限要素法，CIVA法，体積補正手法 連絡先： 〒112-8551東京都文京区春日1-13-27 TEL 03-3817-1815 FAX 03-3817-1803

土木学会第58回年次学術講演会（平成15年9月）

‑483‑

CS9‑024

ここで，V，Verrは各ステップにおける液体の体積および初 期値と任意ステップにおける液体の体積の差であり，Vinit

は初期の液体の体積である．なお，φerrはその値が負であ れば液体の体積が減少，正であれば液体の体積が増加して いることを表現している．このφerr を用いて，式(12)の ように界面近傍のVOF関数φを単純に増減させることに より体積補正を行う．ただし，体積を補正する範囲はVOF 関数φの値によって異なり，φerrが正のときは界面より気 相側に，負のときは液相側に補正を施した．

φ(t) =φ(t)−2φerr

(0< φ <0.5, φerr>0

0.5< φ <1, φ_err<0 (12)

ここで，係数2は界面近傍(0< φ <1)の1/2の領域に補 正を施すことから，導き出される係数である．また，式(12) を用いるとVOF関数φが1以上の値もしくは0以下の値 になることがある．この場合，1以上の値であれば1を，0 以下の値であれば0とすることにより，最終的に体積補正 を行ったVOF関数を決定する．本体積補正手法を総括す ると，VOF関数値が変化をする領域では図-1のような概念 図で表せる．図より，本手法はVOF関数が液体であれば 1，気体であれば0に近づける働きを有しているため，界面 鋭敏化の効果を同時に発揮することができる方法であるこ とがわかる．

f=0 1

0.5

0

(fe r r )

2f

2f

f=1

(fe r r )

図– 1 本体積補正手法の概念図

3. 数値解析例

本手法の有効性を検討するために，数値解析例として，

図-2に示すような3次元矩形貯水槽内スロッシング問題を 取り上げた．なお，起振力として式(13)で表わされる水平 加速度を加えた．

1m

1m

0.5m 0.1m

X Y Z

50 5 50 15606

75000

D =0.001t s

A=9.3 10^-3m w=5.311rad/s

図– 2 矩形スロッシング解析モデルと解析条件

f =Aω²sinωt (13)

ここで，液体は水と仮定し，ρl= 998.0kg/m³, µl= 1.01× 10⁻³Ns/m²，気体は空気と仮定し，ρg= 1.205kg/m³, µg= 1.81×10⁻⁵Ns/m²を与えた．

3.0[s] 6.0[s] 9.0[s]   

図– 3 矩形貯水槽内スロッシング解析結果

0 2 4 6 8 10

0.4 0.6 0.8

waterlevelonleftwall[m]

time [s]

図– 4 左壁における時刻水位暦

図-3に計算結果の時系列を示す．この結果より本手法に よる計算は安定に計算ができていることが確認される．ま た，左壁における時刻水位暦を図-4に示す．なお，本手法 の精度を検証するために，岡本による実験⁵⁾およびALE法

6)（節点数4305）と比較を行った．図より，本手法による

解析結果は実験値，ALE法と良い一致を示しており，計算 精度の観点で有効性が確認できた．

4. おわりに

本報告では，VOF法を用いた自由表面流れ解析を安定か つ高精度に行う手法を構築した．数値解析例として, 3次元 矩形貯水槽内スロッシング解析を取り上げ本手法の有効性 を検討し，以下の結論を得た．流れ場の計算に安定化有限 要素法，移流方程式の解法にCIVA法および体積補正手法 を導入したことにより，3次元解析の結果は実験値および ALE法の解析値とほぼ一致を示し，安定かつ高精度な計算 が行え，本手法の有効性が示された．

今後は本手法の数値波動問題への適用を予定している．

参考文献

1) C.W.Hirt and B.D.Nichols: Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries, J. Comput.

Phys., Vol.39, pp.201-225，1981   

2) T.Tezduyar: Stabilized finite element formulations for incompressible flow computations, Advan. Appl. Mech., Vol.28, pp.1-44, 1991   

3) N.Tanaka: The CIVA Method for Mesh-Free Approaches:

Improvement of the CIP Method for n-Simplex, Comput.

Fluid Dynamics J., Vol.8, no.1, pp121-127, 1999    4) T.Yabe, T,Aoki: A universal solver for hyperbolic equa-

tion by cubic-polynomial intrpolation . One-dimentional solver, Comput. Phys. Commun., 66, pp233-242, 1991 5) 岡本隆: 任意ラグランジェ・オイラー有限要素法による大規模

スロッシング解析に関する研究，中央大学学位論文，1992 6) 桜庭雅明，田中聖三，樫山和男: PCクラスタを用いたALE

並列有限要素法による非線形自由表面流れ解析，応用力学論文 集，Vol.4, pp.113-120, 2001

土木学会第58回年次学術講演会（平成15年9月）

‑484‑

CS9‑024