如何处理界面

可压缩两相算例里，卡住的往往不是通量。是两种流体相接的那一面，在单元里怎么记下来。

分叉有三条。

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$。五方程最后一行写成非守恒，原因在这里。质量由上面两行 $\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$ 哪怕滑过一个单元，质量就漏了。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 七方程已经写开了。

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

通过松弛 (relaxation) 过程收敛到压力/速度平衡。 在数学上最完整，但计算成本高。

比较摘要

把α按守恒型来更新，接触面上的压力会晃

可压缩问题里，表上留下的是底下两行。五方程用 $\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$ 的界面单元会推动混合密度，压力跟着晃。

要不要每步把抹开的带子拉回来，是另一个决定。那一步写在用反扩散把界面重新立起来里。