本稿では OpenFOAM 6 に実装したSSFモデルの動作確認を行う。
確認内容は次の文献に準ずる。

Vachaparambil, K. J., Einarsrud, K. E. 2020. On sharp surface force model: Effect of sharpening coefficient

SSFモデルの実装については、次の記事を参照されたい。

OpenFOAMにおける表面張力モデルの実装

1.イントロダクション

Raeiniら(2012)が開発したSharp Surface Force(SSF)モデルは、一般的に使用されているCSFモデルと比較して、スプリアス・カレントを減らすことが知られている。
本稿では、OpenFOAM 6上で開発したSSFモデルで使用されるsharpening係数の効果を調査し、毛細管上昇のような動的ケースと、停滞気泡のような静的ケースをモデル化する。

2. ソルバー設定

2.1 最大タイムステップ

シミュレーションは緩和係数なしで実行される。
最大タイムステップ $\Delta t$ は次式を満たす値を設定する。

\Delta t \leq \frac{1}{2} \biggl(C_2 \tau_\mu + \sqrt{(C_2 \tau_\mu)^2 + 4 C_1 \tau_\rho^2} \biggr)\tag{1}

\Delta t \leq max(C_2 \tau_\mu, 10 C_1 \tau_\rho)\tag{2}

ここで、

C_1 = 0.01, C_2 = 10, \tau_\mu = \frac{\mu_{avg} \Delta x}{\sigma}, \tau_\rho = \sqrt{\frac{\rho_{avg}(\Delta x)^3}{\sigma}}

$\mu_{avg}$ および $\rho_{avg}$ は相間の平均粘性係数および密度として定義される。$\Delta x$ はメッシュ幅である。

2.2 代数方程式ソルバー

代数方程式ソルバーの設定を表2-1に示す。

表2-1 代数方程式ソルバーの設定

代数方程式	線形ソルバー	Smoother/Preconditioner	許容値
圧力補正方程式	PCG	GAMG	$10^{-20}$
運動方程式	smoothSolver	symGaussSeidel	$10^{-10}$
体積分率方程式	smoothSolver	symGaussSeidel	$10^{-10}$

代数方程式の求解で使用されるその他のパラメータを表2-2に示す。

表2-2 その他のパラメータ

パラメータ	設定値	メモ
nAlphaCorr	2	$\alpha_1$の補正回数
nAlphaSubCycles	1	$\alpha_1$方程式内のサブサイクル回数
cAlpha	1	方程式内の界面圧縮で使用される
MULESCorr	yes	準陰解法MULESを選択
nLimiterIter	3	制限を超えたMULESの反復回数
momentumPredictor	no	SIMPLE法で最初に運動方程式を解くかどうかの設定
minIter	1	運動方程式の最低反復計算回数
nOuterCorrectors	1	PIMPLE法の反復ループ回数
nCorrectors	3	PISO法の圧力補正ループ回数
nNonOrthogonalCorrectors	0	非直交補正ループ回数

2.3 離散化スキーム

離散化スキームの設定を表2-3に示す。

表2-3 離散化スキーム

モデル項	キーワード	スキーム
時間微分項	ddtSchemes	Euler
発散項	div(rhoPhi,U)	vanLeerV
	div(phi,alpha)	vanLeer
	div(phirb,alpha)	interfaceCompression
勾配項	gradSchemes	linear
ラプラシアン項	laplacianSchemes	linear corrected
その他	snGradSchemes	corrected
その他	interpolationSchemes	linear

3.シミュレーション

3.1 毛細管上昇

3.1.1 計算領域

計算領域は $1 mm \times 20 mm$ とし、メッシュの分割数は $20 \times 400$ とする。平行板の間隔は $1 mm$ とする。

3.1.2 境界条件

体積分率 (alpha.water) は出口 (atmosphere) でゼロ勾配、入口 (inlet) はディリクレ条件 (固定値1)、両壁面 (walls) は静的接触角45°のゼロ勾配とする。圧力 (p_rgh) は、出入口でディリクレ条件 (固定値0)、両壁面でノイマン条件とする。速度 (U) は両壁面でスリップなし、出入口は、pressure inlet outlet velocity 条件とする。この条件では、流出速度にはゼロ勾配条件が適用され、流入速度は内部セルのパッチフェイス法線成分から取得された速度を適用する。境界条件の一覧を表3-1に示す。

表3-1 境界条件

	inlet	atmosphere	walls
U	pressureInletOutletVelocity of (0 0 0)		noSlip
alpha.water	fixedValue of 1	zeroGradient	constantAlphaContactAngle of 45 deg
p_rgh	fixedValue of 0		fixedFluxPressure

3.1.3 物性値

物性値を表3-2に示す。

表3-2 物性値

物理量	単位	設定値
水の密度 $\rho_1$	$kg/m^3$	$1000$
水の動粘性係数 $\nu_1$	$m^2/s$	$10^{-6}$
空気の密度 $\rho_2$	$kg/m^3$	$1$
空気の動粘性係数 $\nu_2$	$m^2/s$	$1.48 \times 10^{-5}$
表面張力 $\sigma$	$N/m$	$0.07$
重力加速度 $\vec{g}$	$m/s^2$	$(0\quad-10\quad0)$

3.1.4 その他の解析条件

初期水位高さは $8 mm$ とする。
最大タイムステップは、式(2)を満たす $3.5\mu s$ とする。
クーラン数は $0.1$ とする。
初期タイムステップは、最大タイムステップの10分の1の値とする。
計算終了時間は $1.5 s$ とする。
sharpening 係数 ($C_{sh}$) を $0.0$ から $0.5$ まで $0.1$ ずつ増加させたケースを実行する。

3.1.5 シミュレーション結果

毛細管現象の平衡点の高さ ($h_b$) の解析解は次の通りである。

h_b=\frac{2\sigma \cos \theta}{\Delta \rho |\vec{g}| a}\tag{3}

ここで、$\Delta \rho$ は液体と気体の密度差、$a$ は平行板の間隔である。
上述の物理量を代入すると、$h_b=9.91 mm$ となる。
シミュレーション結果の毛細管高さ $h_{b,calc}$ は次式で計算する。

h_{b,calc}=\frac{\int \alpha_1 dS}{a}\tag{4}

毛細管高さの誤差は次式で評価する。

E(h)=\frac{h_{b,calc}-h_b}{h_b}\tag{5}

スプリアス・カレント $U_{sc}$ は次式で計算する。

U_{sc}=max(|\vec{U}|)\tag{6}

スプリアス・カレントの大きさはシミュレーション終了時の値で評価する。

$C_{sh}$ を変化させたときの毛細管高さとその誤差、スプリアス・カレントの大きさを表3-3に示す。

表3-3 $C_{sh}$ を変化させたときの毛細管高さとその誤差、スプリアス・カレントの大きさ

$C_{sh}$	$U_{sc} (m/s)$	$h_{b,calc} (mm)$	$E(h)$	備考
$0.0$	$0.1735$	$9.49$	$0.042$
$0.1$	$0.0078$	$9.38$	$0.053$
$0.2$	$0.0066$	$9.47$	$0.045$
$0.3$	$0.0089$	$9.43$	$0.048$	振動している
$0.4$	$0.0038$	$9.43$	$0.049$
$0.5$	$0.0025$	$9.42$	$0.049$

$C_{sh}$ を変化させたときの毛細管高さの時系列をそれぞれ図3-1に示す。$C_{sh} = 0.3$ のとき毛細管高さがわずかに振動し、$1s$ から $1.5s$ の間で計算すると、近似的に $h_{b,calc} = 9.43±0.008 mm$ となった。

図3-1 $C_{sh}$ を変化させたときの毛細管高さの時系列

$C_{sh}$を変化させたときのスプリアス・カレントの時系列をそれぞれ図3-2に示す。$C_{sh} = 0.0$ のときのスプリアス・カレントが、他のケースよりも飛び抜けて高い。$C_{sh} = 0.3$ のときは毛細管高さと同様にスプリアス・カレントも振動する。

図3-2 $C_{sh}$を変化させたときのスプリアス・カレントの時系列

3.2 停滞気泡

3.2.1 計算領域

直径$2R$の気泡を、寸法$4R \times 4R$の正方形領域の中心に配置する。$R$ は $2.5mm$ とする。

3.2.2 境界条件

速度 (U) および体積分率 (alpha.water) はすべての境界でゼロ勾配とする。一方、圧力 (p_rgh) はすべての境界でディリクレ条件 (固定値$101,325Pa$)とする。境界条件の一覧を表3-4に示す。

表3-4 境界条件

	walls
U	zeroGradient
alpha.water	zeroGradient
p_rgh	fixedValue of 101325

3.2.3 物性値

物性値を表3-5に示す。

表3-5 物性値

物理量	単位	設定値
水の密度 $\rho_1$	$kg/m^3$	$1000$
水の動粘性係数 $\nu_1$	$m^2/s$	$10^{-6}$
空気の密度 $\rho_2$	$kg/m^3$	$1$
空気の動粘性係数 $\nu_2$	$m^2/s$	$1.48 \times 10^{-5}$
表面張力 $\sigma$	$N/m$	$0.07$
重力加速度 $\vec{g}$	$m/s^2$	$(0\quad0\quad0)$

3.2.4 メッシュおよび最大タイムステップ

メッシュおよび最大タイムステップを表3-6に示す。

表3-6 メッシュおよび最大タイムステップ

メッシュ	メッシュ解像度$(mm^2)$	総セル数	$2R/\delta x$	最大タイムステップ
$M1$	$0.25 \times 0.25$	$1600$	$20$	$3 \times 10^{-5}$
$M2$	$0.125 \times 0.125$	$6400$	$40$	$1 \times 10^{-5}$
$M3$	$0.083 \times 0.083$	$14400$	$60$	$6 \times 10^{-6}$

3.2.5 その他の解析条件

クーラン数は $0.1$ とする。
緩和係数は $1$ とする。
初期タイムステップは最大タイムステップの10分の1の値とする。
計算終了時間は $0.05s$ とする。
sharpening 係数 ($C_{sh}$) を $0.0$ から $0.5$ まで $0.1$ ずつ増加させたケースを実行する。

3.2.6 シミュレーション結果

2次元気泡の場合、ラプラス圧力はヤング・ラプラス方程式を用いて次のように計算できる。

\acute{\Delta p_c} = \frac{\sigma}{R}\tag{7}

気泡内のラプラス圧力は次のように計算される。

\Delta p_c = \frac{\int \alpha_2 p dV}{\int \alpha_2 dV} - p_0\tag{8}

ここで、$p_0$ は大気圧 $(101,325 Pa)$、$\alpha_2$ は空気の体積分率で、$\alpha_1$ を水の体積分率とすると、$\alpha_1 + \alpha_2 = 1$ の関係がある。ラプラス圧力に関連する平均誤差は次のように定義される。

\bar{E} (\Delta p_c) = \frac{\bar{\Delta p_c} - \acute{\Delta p_c}}{\acute{\Delta p_c}}\tag{9}

ここで、バーは時間平均値を表す。

シミュレーション終了時の速度(U)の大きさ(M3メッシュ)を図3-3に示す。

図3-3 シミュレーション終了時の速度(U)の大きさ(M3メッシュ)

平均スプリアス・カレント、平均ラプラス圧力およびその誤差を表3-7に示す。

表3-7(1/3) 平均スプリアス・カレント、平均ラプラス圧力およびその誤差(メッシュ$M1$)

$C_{sh}$	$\bar{U_{sc}}$	$\bar{\Delta p_c}$	$E(\bar{\Delta p_c})$
$0.0$	$0.112$	$24.129$	$-0.1383$
$0.1$	$0.078$	$24.266$	$-0.1334$
$0.2$	$0.071$	$24.416$	$-0.1280$
$0.3$	$0.066$	$24.515$	$-0.1245$
$0.4$	$0.061$	$24.573$	$-0.1224$
$0.5$	$0.057$	$24.608$	$-0.1211$

表3-7(2/3) 平均スプリアス・カレント、平均ラプラス圧力およびその誤差(メッシュ$M2$)

$C_{sh}$	$\bar{U_{sc}}$	$\bar{\Delta p_c}$	$E(\bar{\Delta p_c})$
$0.0$	$0.075$	$25.257$	$-0.0980$
$0.1$	$0.050$	$25.210$	$-0.0997$
$0.2$	$0.044$	$25.260$	$-0.0979$
$0.3$	$0.038$	$25.278$	$-0.0972$
$0.4$	$0.035$	$25.271$	$-0.0975$
$0.5$	$0.033$	$25.294$	$-0.0966$

表3-7(3/3) 平均スプリアス・カレント、平均ラプラス圧力およびその誤差(メッシュ$M3$)

$C_{sh}$	$\bar{U_{sc}}$	$\bar{\Delta p_c}$	$E(\bar{\Delta p_c})$
$0.0$	$0.060$	$25.488$	$-0.0897$
$0.1$	$0.053$	$25.481$	$-0.0900$
$0.2$	$0.048$	$25.498$	$-0.0894$
$0.3$	$0.042$	$25.496$	$-0.0894$
$0.4$	$0.037$	$25.497$	$-0.0894$
$0.5$	$0.034$	$25.501$	$-0.0893$

より大きな sharpening 係数を使用すると、シミュレーションにおけるラプラス圧力の計算の誤差とスプリアス・カレントが減少する傾向がある。また、メッシュサイズを小さくしても、必ずしもスプリアス・カレントが増加するわけではない。これは、CSF モデルを使用した場合とは異なる。

4. まとめ

「OpenFOAMにおける表面張力モデルの実装」で実装した表面張力モデル（SSFモデル）の計算結果は、文献とほぼ同じ結果であり、実装に問題がないことを確認できた。