Transcript S-CLSVOF法

OpenFOAMに
実装したS-CLSVOF法
の検証
２０１３年１１月２３日
大阪大学大学院基礎工学研究科
物質創成専攻化学工学領域
修士２年 山本卓也
VOF法の界面再構築法
VOF
精度悪い
界面がなまりやすい
0
0
0
0
0
0
0
0
0.1
0.3
0
0
0.5
0.95
1.0
0
0.4
1.0
1.0
1.0
0
0.7
1.0
1.0
1.0
A. Albadawi et al. (2013)
S-CLSVOF (Simple Coupled Volume of Fluid with Level Set method)
OpenFOAM
VOF法(interFoam)
Air
S-CLSVOF法
誤差大
water
Air
実験
脱離時間
約35~50%
脱離体積
約25%
誤差小
実験
脱離時間
約2%
脱離体積
約3%
今回はS-CLSVOF法の実装とその検証
A. Albadawi et al., Int. J. Multiphase Flow, 53, 11-28 (2013).
interFoam (VOF法)
• 支配方程式
Navier-Stokes 式
v
 v  v  P   2 v  F  g
t
F  kn s
CSFモデル
流体率aの移流方程式
a
   av   0
t
液相領域 a    avl   0
t
a
   av   0
t
固相領域 a    1  a v g   0
:: liquid phase
a 1
0  a  1 :: interface
a 0
:: gas phase
r = arl + (1- a )rg
m = aml + (1- a )m g
小文字l, gはそれぞれ液相、気相を表す。
t
再定義
v  avl  1  a v g
v r  vl  v g
vr: 相関速度（圧縮速度）
interFoam (VOF法)
• 支配方程式
Navier-Stokes 式
v
 v  v  P   2 v  F  g
t
F  kn s
流体率aの移流方程式
¶a
+ Ñ × (av ) = 0
¶t
:: liquid phase
a 1
0  a  1 :: interface
a 0
:: gas phase
r = arl + (1- a )rg
m = aml + (1- a )m g
最終形
a
   av     1  a av r   0
t
alphaEqn.H中で設定
この項は界面上に働くもの
(1-a)aが入っているため
実際の物理現象では界面厚み
がないため数値計算のために
用いられる仮想的なもの
interFoam (VOF法)
a
   av     1  a avr   0
t
離散化（program中で解く形に変換）
a式の設定
d
adV  av  ndS 

S
dt V


 1  a av
S
r
イメージ
 ndS
Sf
avn
離散化
中点公式により近似
av  n f  S f
ここで、fはセル界面上を表す。
Sfは表面積
1  a av r  n f  S f
MULESでこのfluxを計算
interFoam (VOF法)
人工的に界面幅を圧縮する項
1  a a v  n  S
r
f
a場（赤:流体,
青:気体）
a
f
OpenFOAMのコード内

 uf Sf
uf Sf
n f min Ca
, max
 S

Sf
f


n f  n fv  S f
n fv =
( Ña ) f
( Ña ) f + dN
1.0e 8
N 
( Vi / N )1/ 3

N
解を安定化するもの
(nfvが無限大になるのを防ぐ)




S-CLSVOF法
• 支配方程式
Navier-Stokes 式
¶v
+ v × Ñv = -ÑP + nÑ2v + Fs + r g
¶t
流体率aの移流方程式
¶a
+ Ñ × (av ) = 0
¶t
:: liquid phase
a 1
0  a  1 :: interface
a 0
:: gas phase
r = arl + (1- a )rg
m = aml + (1- a )m g
Level-Set関数 f
f0 = (2a -1)× G
G ; 無次元数
G = 0.75Dx
x ; 無次元数
Re-initialization equation
イメージ
¶f
= S(f0 ) (1- Ñf )
¶t
a
f ( x, 0) = f0 ( x )
反復回数fcorr
e
fcorr =
Dt
界面幅e
e =1.5Dx
Ñf
S-CLSVOF法
• 支配方程式
Navier-Stokes 式
¶v
+ v × Ñv = -ÑP + nÑ2v + Fs + r g
¶t
流体率aの移流方程式
¶a
+ Ñ × (av ) = 0
¶t
:: liquid phase
a 1
0  a  1 :: interface
a 0
:: gas phase
r = arl + (1- a )rg
m = aml + (1- a )m g
CSFモデル
Fs = s kdÑf
曲率
æ ( Ñf ) ö
f
÷
k = -Ñ × n f = -Ñ × ç
ç ( Ñf ) + d ÷
è
ø
f
ダイラック関数
d (f ) = 0
d (f ) =
æ pf öö
1æ
1+
cos
ç ÷÷
ç
è e øø
2e è
f >e
f £e
ヘビサイド関数H
H (f ) = 0
1 æ f 1 æ pf öö
H (f ) = ç1+ + sin ç ÷÷
2 è e p è e øø
H (f ) = 1
S-CLSVOF法
• 支配方程式
Navier-Stokes 式
H (f ) = 0
¶v
+ v × Ñv = -ÑP + nÑ2v + Fs + r g
¶t
1 æ f 1 æ pf öö
H (f ) = ç1+ + sin ç ÷÷
2 è e p è e øø
流体率aの移流方程式
¶a
+ Ñ × (av ) = 0
H (f ) = 1
:: liquid phase
a 1
0  a  1 :: interface
a 0
:: gas phase
r = arl + (1- a )rg
m = aml + (1- a )m g
m = H ml + (1- H )m g
¶t
ヘビサイド関数H
f < -e
f £e
f >e
r = H rl + (1- H )rg
一般には物性値をヘビサイド関数で更新
しかし、A. Albadawi et al. (2013)では
物性値はヘビサイド関数で更新せず
検証１（Laplace圧の測定）
• A. Albadawi et al.(2013)で行われている一つの検
元々は M. M. Francois et al., J. Comput. Phys., 213, 141-173 (2006).
証問題
Laplace圧
２種の流体を分ける表面を横切るときに生ずる静水圧の
ジャンプp (表面張力の物理学(p.8)より 著 ドゥジェ
ンヌ、ブロシャール−ヴィアール、ケレ 訳 奥村剛)
æ1 1 ö
Dp = g ç + ÷
è R R' ø
Dp = p - p
in
0
out
¥
p0in
気泡中心部の圧力
p¥out
壁境界での圧力
数値計算による圧力差と理論値の比較
検証１（Laplace圧の測定）
• 計算条件
0.05 m
liquid
物性値
g 0.01 N/m
g 1 kg/m3
mg 10-5 kg/(ms)
l 1000 kg/m3
ml 10-3 kg/(ms)
0.01 m
gas
無重力条件
（静置条件）
計算時間
0.1 sec. (t = 5x10-5 sec.)
理論値のラプラス圧
æ1
èR
D pexact = g ç +
0.05 m
等間隔格子
DX = 0.001 m (Fine)
= 0.0005 m (Coarse)
1ö
÷ = 20
R' ø
計算によるラプラス圧
D p = p0in - p¥out
p0in 気泡中心部の圧力
p¥out 壁境界での圧力
相対圧力誤差E0
D p - D pexact
E0 =
D pexact
検証１（Laplace圧の測定）
• 計算結果
0.05 m
liquid
0.01 m
gas
手法
S-CLSVOF (F)
S-CLSVOF (C)
VOF (F)
VOF (C)
理論値のラプラス圧
æ1
èR
D pexact = g ç +
0.05 m
E0
1.21
1.67
19.29
25.23
相対圧力誤差E0
1ö
÷ = 20 E = D p - D pexact ´100
0
R' ø
(％) D pexact
計算によるラプラス圧
D p = p0in - p¥out
p0in 気泡中心部の圧力
p¥out 壁境界での圧力
まとめ
S-CLSVOF法の実装に成功した
S-CLSVOF法の方が確実に精度が高い
計算時間少し上昇
他のケースにも実装予定
物性値の更新も実装予定
コードが見せれるようになればOpenにする予定(綺麗に
書き直し・コード検証)
• CLSVOF法も時間があれば挑戦予定
• LS関数の境界条件の実装は怪しい
• ダイナミックな変形問題への検証
•
•
•
•
•
•
計算例１（キャビティ内液滴）
0.5 m/s
目的
liquid 1
liquid 2
0.02 m
y
x
0.1 m
剪断によって液滴がどのように変形するか
（浮力差による影響は無視するため物性は
liquid 1とliquid 2で同じで表面張力のみ発生
するという系）
物性値
動粘度 1.0 x 10-3 m2/s
表面張力 10 mN/m
0.1 m
計算１
interFoam (VOF法)
計算２
sclsVOFFoam(S-CLSVOF法)
計算格子
200 x 200 (x, y direction)
計算例１（キャビティ内液滴）
VOF法
S-CLSVOF法
計算例１（キャビティ内液滴）
VOF法
S-CLSVOF法
界面の数値的な振動
界面がスムーズ
計算例2（ダムブレイク）
0.584 m
0.1461 m
phase 1
動粘度 1 x 10-6 m2/s
密度 1000 kg/m3
phase 2
0.584 m
phase 1
0.292 m
0.048 m
0.292 m
phase 2
動粘度 1.48 x 10-5 m2/s
密度 1 kg/m3
表面張力 70 mN/m
計算例2（ダムブレイク）
VOF法
S-CLSVOF法
計算時間 約1.36倍
計算例2（ダムブレイク）
VOF法
S-CLSVOF法
0.2 s
0.3 s
0.2 s
0.4 s
0.5 s
0.4 s
0.3 s
0.5 s