界面をどのように扱うか

圧縮性多相流シミュレーションにおいて、最も根本的な問いの一つは 「二つの流体間の界面 (interface) を数値的にどのように表現するか」です。

大きく分けて三つの系統の方法が存在します。

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

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}$ 項が離散化の過程で適切に処理されないと、EOSの混合による非物理的な圧力振動が発生します。

非保存型の移流は、$\alpha_1$ を界面でシャープに維持しつつ、圧力と速度の平衡を満足させます。

7-equation Model (Baer-Nunziato type)#

各相が 独立した圧力と速度 を持つことができる完全非平衡モデルです。

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

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

比較のまとめ

結論

圧縮性多相流においては、Diffuse Interface (特に 5-equation model) が最も自然で堅牢 (robust) な選択肢です。 衝撃波と界面の相互作用、キャビテーション、水中爆発などの問題において卓越した性能を発揮します。

次の記事では、これらの手法を実際にコードで実装する方法について解説します。

VOF / Level Set / Diffuse を切り替えて同じ円形界面の表現の違いを比較。