はじめに
OpenFOAMのソルバの一つであるinterFoamのチュートリアルケースに毛管上昇(capillary rise)がある。
Foundation版もESI版もリリースのある時点から界面の振動が止まらなくなった。
果たして、振動が止まらないのが正しいのかどうか、調査を行った。
まだ、ESI版ではv2212からsmoothingの回数を指定できるようになったので、その回数が毛管上昇に与える影響を調査した。
毛管上昇動特性
Jurin's height
自由曲面の曲率を無視した重力と平衡となる毛細管での高さは次式で計算できる。
h_0=
\frac{2\sigma cos \theta_0}{\rho g R}
ここで、
| 名前 | 説明 | 単位 |
|---|---|---|
| $h_0$ | 液相の平衡高さ (Jurin's height) | $[m]$ |
| $\sigma$ | 表面張力 | $[N/m]$ |
| $\theta_0$ | 平衡接触角 | $[deg]$ |
| $\rho$ | 液相の密度 | $[kg/m^3]$ |
| $g$ | 重力加速度 | $[m/s^2]$ |
| $R$ | 毛管の半径 | $[m]$ |
力のバランス式
液相の高さを$h=h(t)$とした場合の力のバランス式は次の通り。
2 \pi R \sigma cos \theta_0 - \pi R^2 \rho g h - 8 \pi \mu h \frac{dh}{dt} = \pi R^2 \rho \frac{d}{dt}(h \frac{dh}{dt} )
この式の左辺は外力で、第一項が表面張力、第二項が重力、第三項が粘性力を表している。
この式を次の量で無次元化する。
H=\frac{h}{h_0}
s=\frac{t}{\tau}
\tau = \sqrt{\frac{h_0}{g}}
すると、無次元化した力のバランス式は次のようになる。
\frac{d}{ds}(H \frac{H}{ds}) + \Omega H \frac{dH}{ds} + H -1 = 0
ここで、$\Omega$ は、
\Omega = \sqrt{\frac{128 \mu^2 \sigma cos \theta_0}{R^5 \rho^3 g^2}}
状態が定常に近いと仮定して、
H(s) = 1 + \epsilon (s)
これを力のバランス式に代入して、微分を含む $\epsilon$ の自乗項を無視すると、
\frac{d^2 \epsilon}{ds^2} + \Omega \frac{d \epsilon}{ds} + \epsilon = 0
この式を振動の微分方程式とみなし、$\epsilon = e^{\lambda s}$ とすると、特性方程式は、
\lambda^2 + \Omega \lambda + 1 = 0
この方程式の解は、
\lambda = \frac{-\Omega ± \sqrt{\Omega^2 - 4}}{2}
この解は次の3つに分類できる。
- $\Omega^2 - 4 < 0$ の場合 ($\Omega < 2$)
この解は不足減衰で、振幅が時間とともに指数関数的に小さくなりながら振動する減衰振動となる。 - $\Omega^2 - 4 = 0$ の場合 ($\Omega = 2$)
この解は臨界減衰で、振動せずに釣り合いの位置に収束する非周期的減衰を示す。 - $\Omega^2 - 4 > 0$ の場合 ($\Omega > 2$)
この解は過減衰で、振動せずに釣り合いの位置に収束する非周期的減衰を示す。
まとめると、$\Omega < 2$ の場合は水位が振動する可能性が高いことを意味している。
OpenFOAMのチュートリアルケース(capillaryRise)
ESI版OpenFOAM v2406の設定値および解析値は次の通りである。
| 物理量 | 単位 | 値 | 備考 |
|---|---|---|---|
| 密度 $\rho$ | $[kg/m^3]$ | 1000 | |
| 動粘性係数 $\nu$ | $[m^2/s]$ | 1e-6 | |
| 粘性係数 $\mu$ | $[Pa s]$ | 1e-3 | $\mu = \nu \rho$ OpenFOAMの設定値ではない |
| 静的接触角 $\theta_0$ | $[deg]$ | 45 | |
| 毛管半径 $R$ | $[m]$ | 1e-3 | チュートリアルケース(capillaryRise)の幅に相当 |
| 重力加速度 $g$ | $[m/s^2]$ | 10 | 絶対値を記載 |
| 表面張力 $\sigma$ | $[N/m]$ | 0.0707106 | |
| 平衡水位 $h_0$ | $[m]$ | 0.01 | 初期水位は8e-3$[m]$ |
| $\Omega$ | $[-]$ | 0.253 | $\Omega < 2$ |
計算結果を次に示す。
図1 OpenFOAMのチュートリアルケース(capillaryRise)の計算結果
$\Omega$が2よりもかなり小さいためか、水位が大きく振動している。なお、チュートリアルケースの元々の計算終了時間は$0.5[s]$であったが、$2.0[s]$に延長した。
チュートリアルケースの改善
毛管半径に該当する幅を小さくすると、$\Omega > 2$ とすることができ、水位が振動しなくなることが予想される。
設定値および解析値は次の通りである。
| 物理量 | 単位 | 値 | 備考 |
|---|---|---|---|
| 密度 $\rho$ | $[kg/m^3]$ | 1000 | |
| 動粘性係数 $\nu$ | $[m^2/s]$ | 1e-6 | |
| 粘性係数 $\mu$ | $[Pa s]$ | 1e-3 | $\mu = \nu \rho$ OpenFOAMの設定値ではない |
| 静的接触角 $\theta_0$ | $[deg]$ | 45 | |
| 毛管半径 $R$ | $[m]$ | 4e-4 | 1e-3$[m]$から変更 |
| 重力加速度 $g$ | $[m/s^2]$ | 10 | 絶対値を記載 |
| 表面張力 $\sigma$ | $[N/m]$ | 0.0707106 | |
| 平衡水位 $h_0$ | $[m]$ | 0.025 | 初期水位を1e-2$[m]$に変更 |
| $\Omega$ | $[-]$ | 2.50 | $\Omega > 2$ |
計算結果を次に示す。
図2 改善したOpenFOAMのチュートリアルケース(capillaryRise)の計算結果
$\Omega$が2よりも大きいため、水位は振動していない。計算終了時の水位は$0.022[m]$で$h_0$に近い値となっている。
smoothing
改善したOpenFOAMのチュートリアルケース(capillaryRise)を対象にして、smoothingを1から3回行うケースを実施した。
計算結果を次に示す。
図3 smoothing計算結果
nAlphaSmoothCurvatureがsmoothingの回数の設定値で、system/fvSolutionで設定する。
下表にsmoothing計算結果の誤差を示す。
smoothingの回数が1および2の場合は、誤差$E(h)$が大きく減って行くが、3になると大きな変化はない。
| nAlphaSmoothCurvature | $h_{last}$ | $E(h)=\frac{(h_{last}-h_0)}{h_0}$ |
|---|---|---|
| 0 | 0.0222 | -0.1119 |
| 1 | 0.0230 | -0.0794 |
| 2 | 0.0241 | -0.0359 |
| 3 | 0.0242 | -0.0306 |
おわりに
参考文献によると、毛管上昇の振動は$\Omega$だけではなく、"contact line friction"の寄与もあるらしい。そうすると、$\Omega < 2$ ならば必ず振動するとは言い難くなる。
"contact line friction"は実際に観測しないとわからない量のようで、水の毛管上昇でのデータは見つからなかった。
参考文献
[1] M. Fricke, E. Ouro-Koura, S. Raju, R. von Klitzing, J. De Coninck, and
D. Bothe. An analytical study of capillary rise dynamics: Critical conditions and hidden oscillations 2023, DOI: 10.1016/j.physd.2023.133895.