How to Handle the Interface#

In a compressible two-fluid run, the fluxes are rarely what stalls the job. The stall is how a cell writes down the surface where the two fluids meet.

Three branches.

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

Flip the mode toggle among VOF, Level Set, and Diffuse and the same circle changes thickness and sign.

VOF locks color inside a straight cut of the cell. Level set draws a smooth zero contour. Diffuse spreads a band across a few cells.

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 line is the one that splits. Advection of $\alpha_1$ is non-conservative. That is the spot that keeps the interface pressure from ringing.

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 $\alpha_1 \nabla \cdot \mathbf{u}$ on the right is not zero in the discrete divergence, the pressure rings.

Expand the non-conservative advection $\partial\alpha/\partial t + \mathbf{u}\cdot\nabla\alpha = 0$ and it becomes $\partial\alpha/\partial t + \nabla\cdot(\alpha\mathbf{u}) - \alpha\nabla\cdot\mathbf{u} = 0$. Moving the last term across the equal sign is the equation above. On a contact, pressure and velocity are uniform, so $\nabla\cdot\mathbf{u}$ should be zero. If the discrete divergence misses that zero, the right-hand side survives. Call the leftover divergence $\delta$ in an interface cell with $\alpha=0.5$. The source is $0.5\delta$. Alpha moves by that amount, the mixture density changes, and the equation of state does not come back to the same pressure. That term is the seed of the oscillation.

The non-conservative form never builds that right-hand side. Alpha follows $\mathbf{u}$, and the face-flux difference does not increment alpha by itself. That is why the last line of the five-equation model is non-conservative. Mass is kept by the two lines above, $\partial(\alpha_k\rho_k)/\partial t + \nabla\cdot(\alpha_k\rho_k\mathbf{u})=0$. Alpha does not have to be written in conservation form for each partial density to stay honest.

The re-initialization formula already carries its own number. When $|\nabla\phi|=1$ breaks, $\mathrm{sign}(\phi_0)(|\nabla\phi|-1)$ pushes $\phi$ in pseudo-time $\tau$. If $|\nabla\phi|$ is above 1 the term trims $\phi$. If it is below 1 the term pushes the other way. The only stop is $|\nabla\phi|=1$. If the zero $\phi=0$ slips by even one cell, mass leaks. The integral of VOF alpha has no such slip. The normal $\mathbf{n}=\nabla\alpha/|\nabla\alpha|$ only has a slope in cells where alpha sits between 0 and 1, so curvature rides the noise in those cells. $\kappa=-\nabla\cdot(\nabla\phi/|\nabla\phi|)$ returns the same normal smoothly when $\phi$ is a distance. CLSVOF glues the two. Mass from $\alpha$, curvature from $\phi$.

7-equation Model (Baer-Nunziato type)#

Each phase can carry its own pressure and its own velocity. That route is already opened in the Baer–Nunziato seven-equation model.

$$\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#

A conservative alpha update rings the pressure at a contact#

In a compressible table, the rows that stay are the bottom two. The five-equation model keeps mass in $\partial(\alpha_k \rho_k)/\partial t + \nabla\cdot(\alpha_k \rho_k \mathbf{u})=0$ and leaves only the alpha advection non-conservative. If the discrete divergence does not kill $\alpha\nabla\cdot\mathbf{u}$, an interface cell at $\alpha=0.5$ pushes the mixture density and the pressure moves.

Whether to pull the smeared band back every step is a separate call. That fix lives in interface sharpening by anti-diffusion.