毛細管上昇の流体解析
どうやら財団版OpenFOAM7ではnoSlip境界条件に問題があるらしく, 計算精度を上げるほどおかしくなるらしい
(1)：Numerical Simulation of Capillary Rise Using interFoam
(2)：A comparative study of transient capillary rise using direct numerical simulations
財団版OpenFOAMはESI版よりも積極的に新機能を投入する思想らしい
debianとubuntuの関係性みたいだ
このnoSlip境界条件の問題がESI版にも存在するのか検証するために論文の解析ケースをOpenFOAMv2012で真似てみる
物性値
毛細管上昇って数値によっては釣り合う位置から慣性でオーバーシュートするみたい
オーバーシュートするか否かは以下の式でわかる
$\sigma$：表面張力
$\eta$：粘度
$\theta$：接触角
$\rho$：密度
$g$：重力加速度
$R$：管の半径
$\Omega > 2$：過減衰
$\Omega = 2$：臨界減衰
$\Omega < 2$：振動（オーバーシュート）
それで、流体解析に使う物性値は$\Omega = 1$になるように設定するらしい
| 名称 | 文字 | 値 | 単位 |
|---|---|---|---|
| 管直径（半径） | $d(R＝d/2)$ | $0.01(0.005)$ | $\mathrm{m}$ |
| 液体密度 | $\rho$ | $83.1$ | $\mathrm{kg/m^3}$ |
| 液体粘度 | $\eta$ | $0.01$ | $\mathrm{Pa\cdot s}$ |
| 重力加速度 | $g$ | $4.17$ | $\mathrm{m/s^2}$ |
| 表面張力 | $\sigma$ | $0.04$ | $\mathrm{N/m}$ |
| 接触角 | $\theta$ | $30$ | $\mathrm{deg}$ |
| 管高さ | - | $0.04$ | $\mathrm{m}$ |
管の直径の単位は、(1)では$\mathrm{mm}$になっているが$\mathrm{m}$が正解なので注意
計算すると…ほぼ１
解析解
$$ h(t)=\sqrt{h^2_c+\exp(-\zeta\omega_0 t)\cos(\sqrt{1-\zeta^2}\omega_0 t+\varphi)} $$導出は難しくてわかんね
初期水面高さを任意の値にして、オーバーシュートのモデリングを含めるとこうなるよう
excelでグラフ化した
(1)の解析解のグラフと同じグラフが出たので、これでよさそう
橙の破線と水色の実線は、力が釣り合って静止したときの水面高さ
本来の式は水色の実線に収束するような式だが、(1)では簡略化のために橙の破線に収束するような式になっている
解析上の理論値は実線に収束する
流体解析ケースの作成
(1)と同じように、二次元解析で計算する
OpenFOAMのtutorialにCapillaryRiseがあるので少し変えるだけ
↓
$FOAM_TUTORIALS/multiphase/interFoam/laminar/capillaryRise
blockMeshDict
管は直径0.01, 高さ0.04で、解析上では図の感じで管の半分だけ計算する
0番が原点(0 0 0)で、順に指定
二次元解析ってどうやるんだろうと思ってたが、z軸方向にメッシュを切らないようにすればよいらしい
1回目は$4[\mathrm{cell/R}]$の分割数で解析したいのと、メッシュは正方形に分割したとのことなので、
x方向に$4$分割
高さは直径の8倍なので、y方向は$4\times 8=32$分割
2次元解析なので、z方向は$1$分割
グレーディングはしない
scale 0.001;
vertices
(
(0 0 0)
(5 0 0)
(5 40 0)
(0 40 0)
(0 0 1)
(5 0 1)
(5 40 1)
(0 40 1)
);
blocks
(
hex (0 1 2 3 4 5 6 7) (4 32 1) simpleGrading (1 1 1)
);
boundary内のwallsを消して、2つに分ける
wallLeft
{
type wall;
faces
(
(0 3 7 4)
);
}
wallMiddle
{
type patch;
faces
(
(1 5 6 2)
);
}
こうなりました
controlDict
最大時間ステップは表面張力波の位相が1cellを超えないように設定する
$\Delta x$は、1セルの辺の長さで、先程$0.005[\mathrm{m}]$の幅に4分割と決めたため、
$\Delta x=\frac{0.005}{4}=0.00125[\mathrm{m}]$
計算すると$\Delta t_{max}=2.84\times 10^{-4}[\mathrm{s}]$
分割数を変えるごとにこの値も変わる
| cell/R | $\Delta t_{max}$ |
|---|---|
| 4 | $2.84\times 10^{-4}$ |
| 8 | $1.0\times 10^{-4}$ |
| 16 | $3.55\times 10^{-5}$ |
| 32 | $1.26\times 10^{-5}$ |
| 64 | $4.44\times 10^{-6}$ |
| 128 | $1.57\times 10^{-6}$ |
クーラン数は$0.4$
明記されていないが、おそらく$\Delta t_{max}=\Delta t$
グラフ化用に、functionsで水面位置をデータとして出力する
directionの向きをもとに、locationsで指定した点からの距離を出力する
locationsの点は複数個指定可能
ノルムではなく、各軸で返されるので、斜め方向にdirectionを設定すると、2つの軸での距離が出力される
今回は座標軸に沿っているので、同じ数値が2列出力される
application interFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 1;
deltaT 2.84e-04;
writeControl adjustable;
writeInterval 0.01;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable yes;
adjustTimeStep yes;
maxCo 0.4;
maxAlphaCo 0.4;
maxDeltaT 2.84e-04;
functions
{
interfaceHeight1
{
type interfaceHeight;
libs (fieldFunctionObjects);
locations ((0.005 0 0));
alpha alpha.water;
liquid true;
direction (0 -1 0);
interpolationScheme cellPoint;
log false;
writeControl runTime;
writeInterval 0.01;
}
}
constant/g
value (0 -4.17 0);
constant/transportProperties
気相は変更なし
water
{
transportModel Newtonian;
nu 1.2e-04; /* ν = η/ρ */
rho 83.1;
}
air
{
transportModel Newtonian;
nu 1.48e-05;
rho 1;
}
sigma 0.04;
0.orig/U
wallsを置き換え
wallLeft
{
type noSlip;
}
wallMiddle
{
type zeroGradient;
}
0.orig/alpha.water
wallsを置き換え
wallLeft
{
type constantAlphaContactAngle;
theta0 30;
limit gradient;
value uniform 0;
}
wallMiddle
{
type zeroGradient;
}
0.orig/p_rgh
wallsを置き換え
wallLeft
{
type fixedFluxPressure;
}
wallMiddle
{
type zeroGradient;
}
仮説
財団版OpenFOAM7での問題がESI版には存在しなかった場合、メッシュ解像度を上げていけば普通に解析解に近づいていく
存在した場合、ここでの壁面計算は[1]の2倍なので、誤差はより大きくなる.はず
メッシュ解像度を上げていくと、解析解から離れていく
結果
$4[\mathrm{cell/R}]$のアニメーション
両面あるものもやってみた
$120\mathrm{cell/R}$の解析は長いのでやめた
(1)より精度はあるものの、やはり解像度を上げていくごとにオーバーシュートしなくなる
定常値は全ケースにおいて$h_{apex}$に届かず
$4\sim16\mathrm{cell/R}$までは定常値がほぼ同じところにあるが、それ以上は下がっていく
$4, 8 ,16 \mathrm{cell/R}$を$5[s]$解析した
メッシュが細かくなるごとに発振の波長が伸びている
$4\mathrm{cell/R}$の最大時間ステップを$2.84\times 10^{-4}$と$10$で10秒まで解析してみた
7秒あたりでクーラン数が仕事してる
ESI版でもこの問題は起こるようだ
OpenFoamがっていうより流体解析そのものの問題な気がする
noSlip境界条件は表面張力が支配的になるような小さい空間でかなり不安定
接触線上で散逸が積分不可能だとかなんとか
ナビエ滑り境界条件を適応することメッシュ精度を上げても定常値の変化が緩和するようだが、OpenFoam自体にナビエ滑り境界条件はなく、論文独自の実装だった
論文と同じような実装ができれば表面張力を正確にモデリングできるやもしれぬ