How to Handle the Interface#

One of the most fundamental questions in compressible multiphase flow simulation is "how to numerically represent the interface between two fluids."

There are three main categories of methods:

1. Volume of Fluid (VOF)
2. Level Set
3. Diffuse Interface

1. Volume of Fluid (VOF)#

Principles#

VOF tracks the volume fraction $\alpha$ occupied by a specific fluid in each cell.

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

If $\alpha = 1$, the cell is full of fluid A; if $\alpha = 0$, it's full of fluid B; and if $0 < \alpha < 1$, the interface passes through the cell.

Interface Reconstruction#

The core of VOF is reconstructing the interface position within the cell from the $\alpha$ values:

• SLIC (Simple Line Interface Calculation): Approximates the interface parallel to the coordinate axes.
• PLIC (Piecewise Linear Interface Calculation): Approximates the interface as a line with an arbitrary slope.

In PLIC, the interface normal vector $\mathbf{n}$ is estimated from the gradient of the $\alpha$ field:

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

Pros and Cons#

• Pros: Exact mass conservation (conservative form).
• Cons: Interface reconstruction is complex and computationally expensive. Difficult to extend to 3D. Inaccurate curvature calculation.

2. Level Set Method#

Principles#

The interface is defined as the zero-level contour of a signed distance function $\phi$.

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

• $\phi > 0$: Fluid A
• $\phi < 0$: Fluid B
• $\phi = 0$: Interface

Re-initialization#

During the advection process, $\phi$ loses its property as a signed distance function ($|\nabla\phi| = 1$). To restore this, a re-initialization equation is solved:

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

where $\tau$ is pseudo-time.

Curvature and Surface Tension#

A major advantage of the Level Set method is that geometric quantities are naturally calculated:

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

Surface tension is converted into a body force using the CSF (Continuum Surface Force) model:

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

Pros and Cons#

• Pros: Easy calculation of curvature/normals. Handles topology changes naturally.
• Cons: Fails mass conservation. The interface position can shift during the re-initialization process.

CLSVOF: Combining the Two#

The Coupled Level Set and VOF (CLSVOF) method, which combines the mass conservation of VOF with the geometric advantages of Level Set, is widely used.

3. Diffuse Interface Method#

Shift in Philosophy#

While the previous two methods view the interface as sharp, the diffuse interface method treats it as having an inherently finite thickness.

In this approach, no separate interface tracking equation is required, and the volume fraction is included as part of the system of conservation laws.

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

The most widely used 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$$

The last equation is key. Note that the advection equation for $\alpha_1$ is non-conservative. This is the secret to preventing pressure oscillations at the interface.

Why is Non-conservative Form Needed?#

If $\alpha_1$ is solved in conservative form:

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

If the $\alpha_1 \nabla \cdot \mathbf{u}$ term on the right-hand side is not handled properly during discretization, unphysical pressure oscillations caused by EOS mixing occur.

Non-conservative advection maintains $\alpha_1$ as sharp at the interface while satisfying pressure/velocity equilibrium.

7-equation Model (Baer-Nunziato type)#

A fully non-equilibrium model where each phase can have independent pressure and velocity:

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

Convergence to pressure/velocity equilibrium is achieved through a relaxation process. It is mathematically the most complete but has high computational costs.

Comparison Summary#

Conclusion#

In compressible multiphase flows, the Diffuse Interface method (especially the 5-equation model) is the most natural and robust choice. It shows outstanding performance in problems such as shock-interface interaction, cavitation, and underwater explosions.

In the next article, we will cover how to actually implement these techniques in code.

Toggle VOF / Level Set / Diffuse to compare how each represents the same circular interface.