VOF法における表面張力の記述

• math
• fluidMechanics
• wip
• CFD

CSF Model

CSFモデルでは, 表面張力$\mathbf{f}_c$を次式で計算する

$\sigma$：表面張力係数
$k(=\nabla\cdot \mathbf{n}_s)$：曲率
$\mathbf{n}_s(=\frac{\nabla\alpha}{|\nabla\alpha|})$：単位法線ベクトル
$\delta_s$：デルタ関数

ヤング・ラプラスの式より,

$R$：曲率半径

デルタ関数

レベルセット法においては, 界面からの符号付き距離$\phi$（符号付きとは, 正負が存在するということ）を用い, 単位法線ベクトルにも$\phi$が使用される$\left(\frac{\nabla\phi}{|\nabla\phi|}\right)$.

意味合いとしては, 表面張力は面積力であるため, 界面にのみ働くように式を一般化する必要がある. デルタ関数は界面以外をゼロ, 界面を非ゼロとすることで場合分けする役割を担っている.

しかしVOF法では, 界面を体積として計算するため, 界面はメッシュ解像度に依存した, “滲んた界面"となる（セル境界面を界面に合わせて変形させ方法も存在する（界面追跡法など）が, 静的で変化の少ない場合に限られる）

そこで, $\delta_s\approx|\nabla\alpha|$と近似することで, $0<\alpha<1$の各セルに表面張力を分散させている.
よって式(1)は以下の通りに変形される.

$\nabla\alpha$は体積分率の勾配であり, 変化が激しいほど無限大に発散していく
体積分率$\alpha$の移流方程式により, メッシュ解像度を細かくするほど界面の拡散が抑えられる. だがその反面, 隣接するセルの$\alpha$変化率は大きなものとなる.

セル長を$\Delta x$とすると, $\nabla\alpha$は$\mathrm{O}(\Delta x^{-1})$ で増大していく. 先にも示した通り, 表面張力は$\mathbf{f}_c=\sigma k\nabla\alpha$で計算されるため, セルが小さいほど表面張力は極端に大きくなる.
単位を考慮すると, $\mathbf{f}_c$は即ち体積力密度（$[\rm{N/m^3}]$）であり, $1[\rm{m^3}]$における圧力を表すが, 物理的には表面張力による曲面の形成はメートルスケールで発生しない. ただし, セル内で発生する表面張力は最終的に$[\mathrm{N}]$で記述される. 単位変換では体積力密度に対してセルの体積$\Delta x^{3}$が乗算されるため, 力自体は$\mathrm{O}(\Delta x^{2})$である（理想的な立方格子の場合） .

しかし加速度を求める上では, メートルスケールが使用されるため, メッシュ解像度によっては非常に大きな値を計算することになる.

$a[\rm{m/s^2}]$：加速度
$F_c[\rm{N}]$：表面張力
$\rho[\rm{kg/m^3}]$：密度
$V[\rm{m^3}]$：セルの体積

ここで問題となるのは数値誤差で, 例えば外力の無い無重力空間に存在する, 直径が毛管長より小さい液滴をシミュレーションしたとき, 力の均衡は以下のように記述される.

即ち, $\nabla p$と$\sigma k\nabla\alpha$は等しくなければならない.
解析上の誤差は避けられないが, ここで扱われるスケールは$[\rm{m}]$であるため$\mathrm{O}(\Delta x^{-1})$ であり, 致命的な誤差に発展する.
OpenFOAMではセルのパラメータを共有するなど整合をとる設計を組み込んではいるものの, 曲率$k$も係数として存在する故にこちらの誤差も物理的不正に大きく影響する.
多くの論文では, 時間刻みを細かくて誤差のスケールを減少させるなどの対策例を挙げている.
また, SSFモデルなどの曲率平滑化, 鋭敏化モデルを用いることで, 偽流の発生をいくらか抑えることができる.
近年のOpenFOAMではsmoothedCSFモデルのLS平滑化が標準で実装されている.