界面をどのように扱うか

圧縮性の二流体で止まる場所は、流束そのものではありません。二つの流体が接する面を、セルの中へどう書き残すかです。

分かれ道は三つです。

1. Volume of Fluid (VOF)
2. Level Set
3. Diffuse Interface (拡散界面)

モードのトグルを VOF、Level Set、Diffuse に切り替えると、同じ円の厚みと符号が変わります。

VOF はセル内の直線片で色を閉じ込めます。Level Set はゼロ等高線が滑らかで、Diffuse は帯が数セルに広がります。

1. Volume of Fluid (VOF)#

原理

各セルにおいて特定の流体が占める 体積分率 (volume fraction) $\alpha$ を追跡します。

$$\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = 0$$

$\alpha = 1$ であればそのセルは流体Aで満たされており、 $\alpha = 0$ であれば流体B、$0 < \alpha < 1$ であれば界面がセル内部を通過していることを意味します。

界面再構成 (Interface Reconstruction)#

VOFの核心は、$\alpha$ の値からセル内部の界面位置を再構成することにあります。

• SLIC (Simple Line Interface Calculation): 界面を座標軸に平行に近似します。
• PLIC (Piecewise Linear Interface Calculation): 界面を任意の傾きの直線として近似します。

PLICにおいて界面の法線ベクトル $\mathbf{n}$ は、$\alpha$ フィールドの勾配から推定されます。

$$\mathbf{n} = \frac{\nabla \alpha}{|\nabla \alpha|}$$

長所と短所

• 長所: 質量が正確に保存されます（保存型）。
• 短所: 界面再構成が複雑で計算コストが高いです。3Dへの拡張が難しく、曲率の計算が不正確になりがちです。

2. Level Set Method#

原理

界面を 符号付き距離関数 (signed distance function) $\phi$ のゼロ等値面として定義します。

$$\frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi = 0$$

• $\phi > 0$: 流体A
• $\phi < 0$: 流体B
• $\phi = 0$: 界面

再初期化 (Re-initialization)#

移流 (advection) の過程で、$\phi$ は符号付き距離関数の性質 ($|\nabla\phi| = 1$) を失っていきます。 これを復元するために、再初期化方程式を解きます。

$$\frac{\partial \phi}{\partial \tau} + \text{sign}(\phi_0)(|\nabla\phi| - 1) = 0$$

ここで $\tau$ は仮想時間 (pseudo-time) です。

曲率と表面張力

Level Setの大きな利点は、幾何学的な量の計算が自然に行えることです。

$$\kappa = -\nabla \cdot \left(\frac{\nabla\phi}{|\nabla\phi|}\right)$$

表面張力はCSF (Continuum Surface Force) モデルによって体積力へと変換されます。

$$\mathbf{f}_s = \sigma \kappa \delta(\phi) \nabla\phi$$

長所と短所

• 長所: 曲率や法線の計算が容易です。トポロジーの変化 (topology change) を自然に扱えます。
• 短所: 質量保存に失敗します。再初期化の過程で界面の位置が移動する可能性があります。

CLSVOF: 二つの手法の結合#

VOFの質量保存性とLevel Setの幾何学的な利点を組み合わせた Coupled Level Set and VOF (CLSVOF) 法が広く用いられています。

3. Diffuse Interface Method (拡散界面)#

哲学の転換

これまでの二つの手法が界面を 鋭い (sharp) ものとして捉えるのに対し、拡散界面法は界面が本質的に 有限の厚み を持つものと考えます。

このアプローチでは、別途界面追跡方程式を必要とせず、体積分率が保存則システムの一部として組み込まれます。

5-equation Model (Allaire et al., 2002; Kapila et al., 2001)#

最も広く使われている低次モデル (reduced model) です。

$$\frac{\partial (\alpha_1 \rho_1)}{\partial t} + \nabla \cdot (\alpha_1 \rho_1 \mathbf{u}) = 0$$ $$\frac{\partial (\alpha_2 \rho_2)}{\partial t} + \nabla \cdot (\alpha_2 \rho_2 \mathbf{u}) = 0$$ $$\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u} + p\mathbf{I}) = 0$$ $$\frac{\partial E}{\partial t} + \nabla \cdot ((E + p)\mathbf{u}) = 0$$ $$\frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 = 0$$

最後の行が分かれ目です。$\alpha_1$ の移流が 非保存型 (non-conservative) です。界面の圧力を揺らさない場所がここです。

なぜ非保存型が必要なのか？

$\alpha_1$ を保存型で解くと以下のようになります。

$$\frac{\partial \alpha_1}{\partial t} + \nabla \cdot (\alpha_1 \mathbf{u}) = \alpha_1 \nabla \cdot \mathbf{u}$$

右辺の $\alpha_1 \nabla \cdot \mathbf{u}$ が離散発散で 0 にならないと、圧力が揺れます。

非保存の移流 $\partial\alpha/\partial t + \mathbf{u}\cdot\nabla\alpha = 0$ を展開すると $\partial\alpha/\partial t + \nabla\cdot(\alpha\mathbf{u}) - \alpha\nabla\cdot\mathbf{u} = 0$ です。最後の項を右へ移したのが上の式です。接触面で圧力と速度が一定なら $\nabla\cdot\mathbf{u}$ は 0 であるはずです。離散発散がその 0 を出せないと右辺が残ります。$\alpha=0.5$ の界面セルで残差発散を $\delta$ とすると、右辺は $0.5\delta$ です。$\alpha$ がその分動くと混合密度が変わり、状態方程式は同じ圧力へ戻れません。振動の種はこの項です。

非保存形はその右辺を作りません。$\alpha$ は $\mathbf{u}$ に沿って動き、セル面の流束差分が $\alpha$ を直接増減させません。5方程式の最後の行が非保存である理由はここです。質量は上の二行、$\partial(\alpha_k\rho_k)/\partial t + \nabla\cdot(\alpha_k\rho_k\mathbf{u})=0$ が別に守ります。$\alpha$ を保存形で書かなくても、各相の質量はその行が受け持ちます。

再初期化の数字も式の中にあります。$|\nabla\phi|=1$ が崩れると $\mathrm{sign}(\phi_0)(|\nabla\phi|-1)$ が仮想時間 $\tau$ で $\phi$ を押します。$|\nabla\phi|$ が 1 より大きいとその項が $\phi$ を削り、1 より小さいと逆へ押します。止まる場所は $|\nabla\phi|=1$ だけです。その過程でゼロ点 $\phi=0$ が 1 セルでも滑ると質量が漏れます。VOF の $\alpha$ 積分にはその滑りがありません。代わりに $\mathbf{n}=\nabla\alpha/|\nabla\alpha|$ で取る法線は、$\alpha$ が 0 と 1 の間に広がったセルでしか傾きを持ちません。曲率はそのセルのノイズに乗ります。$\kappa=-\nabla\cdot(\nabla\phi/|\nabla\phi|)$ は $\phi$ が距離関数のとき、同じ法線を滑らかに返します。だから CLSVOF が二つを結びます。質量は $\alpha$、曲率は $\phi$ です。

7-equation Model (Baer-Nunziato type)#

各相に圧力と速度を別々に持たせられる完全非平衡モデルです。その道は Baer–Nunziato の 7 方程式がすでに開いています。

$$\frac{\partial \alpha_k \rho_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k \mathbf{u}_k) = 0, \quad k = 1, 2$$

緩和 (relaxation) 過程を通じて圧力と速度の平衡へと収束させます。 数学的に最も完全ですが、計算コストは高いです。

比較のまとめ

αを保存形で解くと接触面の圧力が揺れる

圧縮性で表から残すのは下の二行です。5方程式は質量を $\partial(\alpha_k \rho_k)/\partial t + \nabla\cdot(\alpha_k \rho_k \mathbf{u})=0$ で守り、$\alpha$ の移流だけを非保存に残します。離散発散が $\alpha\nabla\cdot\mathbf{u}$ を 0 にできないと、$\alpha=0.5$ の界面セルが混合密度を押し、圧力が揺れます。

広がった帯を毎ステップ戻すかどうかは別の決定です。その操作は 反拡散で界面を立て直す実装へ続きます。