1. はじめに
「OpenFOAMにおける表面張力モデルの実装」で実装した表面張力モデル（Smoothed CSFおよびSSFモデル）の妥当性確認および解析を行う。ここでも次の文献を参考にする。
Vachaparambil, K. J., Einarsrud, K. E. 2019a. Comparison of surface tension models for the volume of fluid method. Processes, 7: 542.
2. ソルバー設定
2.1 最大タイムステップ
次の３つの最大タイムステップが提案されている。
- Brackbillらが提案した最大タイムステップ
\Delta t < \sqrt{\frac{\rho_{avg}(\Delta x)^3}{2\pi\sigma}}\tag{1}
ここで、$\Delta x$、$\sigma$および$\rho_{avg}$はそれぞれグリッド幅、表面張力および平均密度である。
- GalusinskiとVigneauxが提案した最大タイムステップ
\Delta t \leq \frac{1}{2}(C_2 \tau_\mu+\sqrt{(C_2 \tau_\mu)^2+4C_1 \tau_\rho^2})\tag{2}
\tau_\mu = \frac{\mu_{ave} \Delta x}{\sigma}, \tau_\rho = \sqrt{\frac{\rho_{avg} (\Delta x)^3}{\sigma}}
ここで、$\mu_{avg}$は平均粘性係数である。
- Deshpandeらが提案した最大タイムステップ
\Delta t \leq max(C_2 \tau_\mu, 10 C_1 \tau_\rho)\tag{3}
Deshpandeらは、interFoamの$C_1$と$C_2$の値がそれぞれ$0.01$と$10$であると計算した。
2.2 離散化スキームの設定
離散化スキームの設定を表2-1に示す。
表2-1 離散化スキームの設定
| モデル項 | キーワード | スキーム |
|---|---|---|
| 時間微分 | ddtSchemes | Euler |
| 発散項 | div(rhoPhi,U) | vanLeerV |
| div(phi,alpha) | vanLeer | |
| div(phirb,alpha) | interfaceCompression | |
| 勾配項 | gradSchemes | linear |
| ラプラシアン項 | laplacianSchemes | linear corrected |
| その他 | snGradSchemes | corrected |
| interpolationSchemes | linear |
2.3 代数方程式ソルバーの設定
代数方程式ソルバーの設定を表2-2に示す。
表2-2 代数方程式ソルバーの設定
| 代数方程式 | 線形ソルバー | Smoother/Preconditioner | 許容値 |
|---|---|---|---|
| 圧力補正方程式 | PCG | DIC | $10^{-20}$ |
| 運動方程式 | smoothSolver | symGaussSeidel | $10^{-12}$ |
| alpha.water方程式 | smoothSolver | symGaussSeidel | $10^{-12}$ |
代数方程式の求解で使用されるその他のパラメータを表2-3に示す
表2-3 その他のパラメータ
| パラメータ | 設定値 | コメント |
|---|---|---|
| nAlphaCorr | 2 | alpha.waterの補正回数 |
| nAlphaSubCycles | 1 | alpha.water方程式内のサブサイクル回数 |
| cAlpha | 1 | 界面圧縮係数 |
| MULESCorr | yes | 準陰解法MULESを選択 |
| nLimiterIter | 3 | 制限を超えたMULESの反復回数 |
| momentumPredictor | no | SIMPLE法で最初に運動方程式を解くかどうかの設定（解かない） |
| minIter | 1 | 運動方程式の最低反復計算回数 |
| nOuterCorrectors | 1 | PIMPLE法の反復ループ回数（1の場合PISO法と同じ） |
| nCorrectors | 3 | PISO法の圧力補正ループ回数 |
| nNonOrthogonalCorrectors | 0 | 非直交補正ループ回数 |
3. 妥当性確認
3.1 気泡上昇
3.1.1 テストケース
テストケースは次式で定義されるレイノルズ数$Re$、エトベス数$Eo$およびキャピラリー数$Ca$に基づいて区別する。
Re = \frac{U_{air} L}{\nu_{water}}\tag{4}
Eo = \frac{\rho_{water} U_{air}^2 L}{\sigma}\tag{5}
Ca = \frac{\mu_{water} U_{air}}{\sigma}=\frac{Eo}{Re}\tag{6}
$Ca$が小さいケースをテストケースTC1、$Ca$が大きいケースをテストケースTC2とする。検証にはHysingらによって報告された$160 \times 320$の均一グリッドのデータを使用する。
3.1.2 計算領域
計算領域は $1m \times 2m$ とし、分割数は $160 \times 320$ とする。気泡の直径は $0.5m$ で、気泡の中心は底面および側面から $0.5m$ の位置にある。
3.1.3 境界条件
境界条件を表3.1-1に示す。
表3.1-1 境界条件
| bottom | top | walls | |
|---|---|---|---|
| U | noSlip | slip | |
| alpha.water | zeroGradient | ||
| p_rgh | zeroGradient | fixedValue of 0 | zeroGradient |
3.1.4 物理量
各テストケースで使用する物理量を表3.1-2に示す。
表3.1-2 物理量
| 物理量 | TC1 | TC2 |
|---|---|---|
| 水の密度$\rho_{water}$ | $1000$ | $1000$ |
| 水の動粘性係数$\nu_{water}$ | $10^{-2}$ | $10^{-2}$ |
| 空気の密度$\rho_{air}$ | $100$ | $1$ |
| 空気の動粘性係数$\nu_{air}$ | $10^{-2}$ | $10^{-1}$ |
| 表面張力$\sigma$ | $24.5$ | $1.96$ |
| 重力加速度$\vec{g}$ | $(0\quad-0.98\quad0)$ | $(0\quad-0.98\quad0)$ |
| $Re$ | $35$ | $35$ |
| $Eo$ | $10$ | $125$ |
| $Ca$ | $0.286$ | $3.571$ |
3.1.5 その他
- 式(3)を満たす $1.4 \times 10^{-2} s$ を最大タイムステップとする。
- クーラン数は $0.01$ とし、初期タイムステップは最大タイムステップの10分の1とする。
- 出力間隔は $0.5s$ とし、計算終了時間は $3s$ とする。
3.1.6 解析結果
テストケースTC1
このケースはキャピラリー数が小さいため、表面張力の効果が支配的なケースに相当する。
各モデルのTC1の気泡の時間的変化を図3.1-1および図3.1-2に示す。
図3.1-1 TC1の気泡の時間的変化: (a) t = 0.5s (b) t = 1.5s (c) t =2.5s
図3.1-2 TC1の計算終了時(t = 3s)の気泡形状: (a)気泡形態および(b)詳細
界面は楕円形の気泡に変形する。
各モデルのTC1の気泡上昇速度の時間的変化を図3.1-3に示す。
図3.1-3 TC1の気泡上昇速度: (a)気泡上昇速度の時間的変化および(b)詳細
各モデルのTC1の最大気泡上昇速度とそれが発生した時間を表3.1-3に示す。
表3.1-3 TC1の最大気泡上昇速度とそれが発生した時間
| パラメータ | CSF | Smoothed CSF | SSF | ベンチマーク |
|---|---|---|---|---|
| $V_{max}$ | $0.2381$ | $0.2392$ | $0.2386$ | $0.2421$ |
| $t(V_{max})$ | $0.9556$ | $0.9424$ | $0.9245$ | $0.9313$ |
ここで、最大気泡上昇速度は気泡の平均速度である。
最大気泡上昇速度については、Smoothed CSFモデルが最もベンチマーク近く、最大気泡上昇速度が発生する時間については、SSFモデルが最もベンチマークに近い結果となった。
テストケースTC2
このケースはキャピラリー数が大きいため、表面張力の影響が低いケースに相当する。
各モデルのTC2の気泡の時間的変化を図3.1-4および図3.1-5に示す。
図3.1-4 TC2の気泡の時間的変化: (a) t = 0.5s (b) t = 1.5s (c) t =2.5s
図3.1-5 TC2の計算終了時(t = 3s)の気泡形状: (a)気泡形態および(b)詳細
時間が進むにつれ界面の変形が大きくなり、最終的には細いフィラメントを持つスカート状の気泡が形成され、それが小さな液滴に分解される。
各モデルのTC2の気泡上昇速度の時間的変化を図3.1-6に示す。
図3.1-6 TC2の気泡上昇速度: (a)気泡上昇速度の時間的変化および(b)詳細
各モデルのTC2の最大気泡上昇速度とそれが発生した時間を表3.1-4に示す。
表3.1-4 TC2の最大気泡上昇速度とそれが発生した時間
| パラメータ | CSF | Smoothed CSF | SSF | ベンチマーク |
|---|---|---|---|---|
| $V_{max1}$ | $0.2499$ | $0.2493$ | $0.2490$ | $0.2514$ |
| $t(V_{max1})$ | $0.7483$ | $0.7466$ | $0.7354$ | $0.7281$ |
| $V_{max2}$ | $0.2320$ | $0.2309$ | $0.2284$ | $0.2440$ |
| $t(V_{max2})$ | $1.9809$ | $1.9781$ | $1.9772$ | $1.9844$ |
TC2では最大気泡上昇速度の極値が2つある。
最初のピークでの最大気泡上昇速度は、どのモデルでも誤差が小さくほぼ一緒であるが、ピークのタイミングについては、SSFモデルの誤差が最も小さい。2番目のピークでは、最大気泡上昇速度の誤差はどのモデルでも大きいが、ピークのタイミングについては、どのモデルも小さい。
3.2 毛細管上昇
液体が細い管内または2枚の平行板の間を上昇する現象は、壁面が液体で濡れることにより発生し、毛細管上昇として知られている。液体が上昇するにつれて、表面張力によって及ぼされる力の垂直成分が上昇した液柱に作用する重力とバランスが取れ、平衡点に達する。その平衡点の高さ $h_b$ の解析解は次の通りである。
h_b = \frac{2 \sigma \cos \theta}{\Delta \rho |\vec g| a}\tag{7}
ここで、$\Delta \rho$ は液体と気体の密度差、$a$ は平行板の間隔である。
後述する物理量を代入すると、$h_b$ の値は $9.91 mm$ となる。
3.2.1 計算領域
計算領域は $1mm \times 20mm$ とし、分割数は $20 \times 400$ とする。平行板の間隔 $a$ は $1mm$ とする。
3.2.2 境界条件
境界条件は下表の通りとする。
表3.2-1 境界条件
| inlet | atmosphere | walls | |
|---|---|---|---|
| U | pressureInletOutlet of (0 0 0) | noSlip | |
| alpha.water | fixedValue of 1 | zeroGradient | constantAlphaContactAngle of 45 degree |
| p_rgh | fixedValue of 0 | fixedFluxPressure | |
3.2.3 物理量
物理量は下表の通りとする。
表3.2-2 物理量
| 物理量 | 設定値 |
|---|---|
| 水の密度$\rho_{water}$ | $1000$ |
| 水の動粘性係数$\nu_{water}$ | $10^{-6}$ |
| 空気の密度$\rho_{air}$ | $1$ |
| 空気の動粘性係数$\nu_{air}$ | $1.48 \times 10^{-5}$ |
| 表面張力$\sigma$ | $0.07$ |
| 重力加速度$\vec g$ | $(0\quad-10\quad0)$ |
3.2.4 その他
- 初期水位高さは $8mm$ とする。
- 最大タイムステップは式(3)を満たす $3.5 \mu s$ とする。
- クーラン数は $0.1$ とする。
- 初期タイムステップは、最大タイムステップの10分の1の値とする。
- 計算終了時間は $1.5s$ とする。
3.2.5 解析結果
毛細管高さは$h_{b,calc}$は次式で計算する。
h_{b,calc} = \frac{\int \alpha_1 dS}{a}\tag{8}
各表面張力モデルでの毛細管高さの計算結果を図3.2-1に示す。
図3.2-1 各表面張力モデルでの毛細管高さの比較
各表面張力モデルの計算終了時の毛細管高さの誤差を表3.2-1に示す。
表3.2-1 各表面張力モデルの計算終了時の毛細管高さの誤差
| 表面張力モデル | 毛細管高さ | 誤差 |
|---|---|---|
| CSF | 9.26 | -0.0066 |
| Smoothed CSF | 計算時間内に毛細管高さは静定せず | |
| SSF | 9.38 | -0.0053 |
ここで、誤差は次式で計算した。
E(h_{b,calc}) = \frac{(h_{b,calc} - h_b)}{h_b}\tag{9}
図3.2-1および表3.2-1より、SSFモデルはCSFモデルよりも毛細管高さを正確に予測できることを示している。
Smoothed CSFモデルでは、水柱の振動により毛細管高さを算出できなかった。
スプリアス・カレントを用いてこの状態を説明する。
スプリアス・カレント $U_{sc}$ を次式のように定義する。
U_{sc} = max(|\vec U|)\tag{10}
各表面張力モデルでのスプリアス・カレントの計算結果を図3.2-2に示す。
図3.2-2 各表面張力モデルでのスプリアス・カレントの比較
Smoothed CSFモデルにおけるスプリアス・カレントの周期的な増加と減少により、界面の非現実的な動きが生じるが、スプリアス・カレントの値がはるかに大きい CSFモデルでは周期性がはるかに高くなり、液体と気体の界面の正味の動きが減少する。CSFモデルおよび Smoothed CSFモデルと比較すると、SSFモデルでのスプリアス・カレントの変化は2桁近く低減される。
4. 解析
4.1 停滞気泡
4.1.1 計算領域
直径 $2R$ の気泡を、寸法 $4R \times 4R$ の正方形領域の中心に配置する。
4.1.2 境界条件
境界条件を表4.1-1に示す。
表4.1-1 境界条件
| walls | |
|---|---|
| U | zeroGradient |
| alpha.water | |
| p_rgh | fixedValue of 101325 |
4.1.3 物理量
物理量を表4.1-2 に示す。
| 物理量 | 設定値 |
|---|---|
| 水の密度$\rho_{water}$ | $1000$ |
| 水の動粘性係数$\nu_{water}$ | $10^{-6}$ |
| 空気の密度$\rho_{air}$ | $1$ |
| 空気の動粘性係数$\nu_{air}$ | $1.48 \times 10^{-5}$ |
| 表面張力$\sigma$ | $0.07$ |
| 重力加速度$\vec g$ | $(0\quad0\quad0)$ |
4.1.4 その他
- 緩和係数を $0.9$ とする。
- 気泡の直径 $2R$ を $5mm$ とする。
- 最大タイムステップは式(3)を満たす値とする。
- クーラン数は $0.1$ とする。
- 初期タイムステップは、最大タイムステップの10分の1の値とする。
- 計算終了時間は $0.05s$ とする。
4.1.5 精度
表面張力モデルの精度は、次のパラメータに基づいて計算される。
- ラプラス圧力
- スプリアス・カレントの大きさ
2次元気泡の場合、ラプラス圧力はヤング・ラプラス方程式を用いて次のように計算できる。
\acute{\Delta p_c} = \frac{\sigma}{R}\tag{11}
気泡内のラプラス圧力は次のように計算される。
\Delta p_c = \frac{\int \alpha_2 p dV}{\int \alpha_2 dV} - p_0\tag{12}
ここで、$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{13}
ここで、バーは時間平均値を表す。
4.1.6 テストケース
メッシュ解像度の異なる次の4ケースの計算を行う。
表4.1-3 テストケース一覧
| メッシュ | メッシュ解像度$(mm^2)$ | 総セル数 | $R/\delta x$ | 最大タイムステップ |
|---|---|---|---|---|
| $M0$ | $0.5 \times 0.5$ | $400$ | $5$ | $9 \times 10^{-5}$ |
| $M1$ | $0.25 \times 0.25$ | $1600$ | $10$ | $3 \times 10^{-5}$ |
| $M2$ | $0.125 \times 0.125$ | $6400$ | $20$ | $1 \times 10^{-5}$ |
| $M3$ | $0.083 \times 0.083$ | $14400$ | $30$ | $6 \times 10^{-6}$ |
M3メッシュでの計算終了時のスプリアス・カレントを図4.1-1に示す。
図4.1-1 M3メッシュでの計算終了時のスプリアス・カレント
この計算で使用している表面張力モデルでは、スプリアス・カレントは界面の周囲に発生し、その大きさは気泡の外側よりも内側の方が大きくなる。スプリアス・カレントを定量化するために、スプリアス・カレントの大きさとラプラス圧力の時間平均値を表4.1-4にまとめる。
表4.1-4 スプリアス・カレントの大きさとラプラス圧力の時間平均値 (緩和係数 $0.9$)
| 表面張力モデル | メッシュ | $\bar{U_{sc}}$ | $Ca = \frac{\rho_1 \nu_1 \bar{U_{sc}}}{\sigma}$ | $\bar{\Delta p_c}$ | $\bar{E}(\Delta p_c)$ |
|---|---|---|---|---|---|
| $CSF$ | $M0$ | $0.1430$ | $2.043 \times 10^{-3}$ | $21.55$ | $-0.2302$ |
| $M1$ | $0.1743$ | $2.490 \times 10^{-3}$ | $22.79$ | $-0.1862$ | |
| $M2$ | $0.1735$ | $2.479 \times 10^{-3}$ | $23.96$ | $-0.1442$ | |
| $M3$ | $0.1890$ | $2.700 \times 10^{-3}$ | $24.51$ | $-0.1245$ | |
| $Smoothed CSF$ | $M0$ | $0.1045$ | $1.493 \times 10^{-3}$ | $23.68$ | $-0.1543$ |
| $M1$ | $0.0891$ | $1.273 \times 10^{-3}$ | $24.70$ | $-0.1179$ | |
| $M2$ | $0.0613$ | $8.754 \times 10^{-4}$ | $25.36$ | $-0.0944$ | |
| $M3$ | $0.0482$ | $6.892 \times 10^{-4}$ | $25.66$ | $-0.0834$ | |
| $SSF$ | $M0$ | $0.0386$ | $5.519 \times 10^{-4}$ | $24.00$ | $-0.1438$ |
| $M1$ | $0.0713$ | $1.019 \times 10^{-3}$ | $24.65$ | $-0.1195$ | |
| $M2$ | $0.0293$ | $4.190 \times 10^{-4}$ | $25.60$ | $-0.0856$ | |
| $M3$ | $0.0373$ | $5.322 \times 10^{-4}$ | $25.49$ | $-0.0897$ |
表面張力モデルによって生成されるスプリアス・カレントは、SSFおよびSmoothed CSFモデルの両方において、メッシュが細かくなるにつれて減少する傾向がある。また、ラプラス圧力についてもSSFおよびSmoothed CSFモデルはCSFモデルよりも良い予測を提供している。
4.1.7 タイムステップの効果
前掲のM3メッシュケースでは、最大タイムステップをDeshpandeらが提案した値(式(3))である $6 \mu s$ とした。ここでは最大タイムステップをBrackbillらが提案した値(式(1))である $25 \mu s$ としてタイムステップ効果を確認する。各表面張力モデルでのスプリアス・カレントの時間変化を図4.1-2に示す。
図4.1-2 タイムステップの効果
CSFモデルでは、Deshpandeらが提案したタイムステップ制約を使用した場合、Brackbillらが提案したタイムステップ制約を使用した場合よりもスプリアス・カレントは半分以下に減少する。一方、Smoothed CSFおよびSSFモデルではスプリアス・カレントの絶対差はそれぞれ約4%、約3%程度である。
M3メッシュおよび表面張力モデルに基づくタイムステップ制約での平均スプリアス・カレントの比較を表4.1-5に示す(緩和係数 $0.9$ を使用)。
表4.1-5 タイムステップ制約での平均スプリアス・カレントの比較 (緩和係数 $0.9$)
| 表面張力モデル | Brackbillらの制約 | Deshpandeらの制約 |
|---|---|---|
| $CSF$ | $0.402775$ | $0.188972$ |
| $Smooehtd CSF$ | $0.050295$ | $0.048244$ |
| $SSF$ | $0.036170$ | $0.037256$ |
4.1.8 緩和係数の効果
緩和係数の影響を理解するために、緩和係数を $1$ としたケースを実施した。
スプリアス・カレントの大きさとラプラス圧力の時間平均値を以下にまとめる。
表4.1-6 スプリアス・カレントの大きさとラプラス圧力の時間平均値 (緩和係数 $1$)
| 表面張力モデル | メッシュ | $\bar{U_{sc}}$ | $Ca = \frac{\rho_1 \nu_1 \bar{U_{sc}}}{\sigma}$ | $\bar{\Delta p_c}$ | $\bar{E}(\Delta p_c)$ |
|---|---|---|---|---|---|
| $CSF$ | $M0$ | $0.1701$ | $2.429 \times 10^{-3}$ | $21.57$ | $-0.2295$ |
| $M1$ | $0.2854$ | $4.077 \times 10^{-3}$ | $22.87$ | $-0.1831$ | |
| $M2$ | $0.5128$ | $7.326 \times 10^{-3}$ | $24.23$ | $-0.1347$ | |
| $M3$ | $0.7250$ | $1.036 \times 10^{-2}$ | $24.23$ | $-0.1312$ | |
| $Smoothed CSF$ | $M0$ | $0.1664$ | $2.378 \times 10^{-3}$ | $23.46$ | $-0.1620$ |
| $M1$ | $0.1244$ | $1.778 \times 10^{-3}$ | $24.65$ | $-0.1198$ | |
| $M2$ | $0.1036$ | $1.480 \times 10^{-3}$ | $25.33$ | $-0.0953$ | |
| $M3$ | $0.0751$ | $1.073 \times 10^{-3}$ | $25.66$ | $-0.0835$ | |
| $SSF$ | $M0$ | $0.0567$ | $8.100 \times 10^{-4}$ | $23.68$ | $-0.1544$ |
| $M1$ | $0.0600$ | $8.578 \times 10^{-4}$ | $24.63$ | $-0.1202$ | |
| $M2$ | $0.0350$ | $4.997 \times 10^{-4}$ | $25.57$ | $-0.0868$ | |
| $M3$ | $0.0350$ | $5.004 \times 10^{-4}$ | $25.62$ | $-0.0851$ |
緩和係数を$1$とした場合と$0.9$とした場合を比較すると、$1$とした方がCSFおよびSmoothed CSFモデルでスプリアス・カレントが大きくなることがわかる。SSFモデルは両方の緩和係数で最小のスプリアス・カレントとなる。
緩和係数を$1$とした場合のM3メッシュおよび表面張力モデルに基づくタイムステップ制約での平均スプリアス・カレントの比較を以下に示す。
表4.1-7 タイムステップ制約での平均スプリアス・カレントの比較 (緩和係数 $1$)
| 表面張力モデル | Brackbillらの制約 | Deshpandeらの制約 |
|---|---|---|
| $CSF$ | $0.727061$ | $0.725033$ |
| $Smooehtd CSF$ | $0.079818$ | $0.075083$ |
| $SSF$ | $0.035506$ | $0.035026$ |
タイムステップ制約に基づいたスプリアス・カレントの進展では、Smoothed CSF および CSF モデルと比較したとき、SSF モデルは最も小さいスプリアス・カレントを生成する。
5. まとめ
「OpenFOAMにおける表面張力モデルの実装」で実装した表面張力モデル（Smoothed CSFおよびSSFモデル）の計算結果は、文献とほぼ同じ結果であり、実装に問題がないことを確認できた。