Un-smearing the Interface Every Step — Implementing Anti-Diffusion Sharpening (IST)

I ran a shock-bubble case with a 5-equation model. After 200 steps the helium bubble's edge had swollen from 6 cells to 20. The physical quantities were fine, but the interface was mush. That is the fate of the diffuse interface method (an approach that deliberately lets the interface spread over a finite width). Nguyen et al. (2021) undo that smearing with a single post-processing pass after every time step. Today we build that anti-diffusion technique from scratch and watch, with our own eyes, why the interface thickness locks to a constant.

Turn on numerical diffusion and the 6-cell edge swallows cells; turn on sharpening (IST) and that spread stops near ε/a.

The yellow band is the target thickness. Raising compression a shrinks it, raising regularization ε fattens it again, and the measured cell count tracks that ratio.

Paper information#

• Title: Numerical modeling of multiphase compressible flows with the presence of shock waves using an interface-sharpening five-equation model
• Authors: Van-Tu Nguyen, Thanh-Hoang Phan, Warn-Gyu Park (Pusan National University)
• Source: International Journal of Multiphase Flow, 2021, 103542
• DOI: 10.1016/j.ijmultiphaseflow.2020.103542
• One-line summary: Bolt interface sharpening (IST) onto a 5-equation compressible two-phase model as a post-processing step to keep interfaces at a constant thickness even next to shocks.

Why the interface fattens on its own#

The 5-equation model is a quasi-conservative system: conservation laws plus one color function. The volume fraction $\alpha_1$ (the share of fluid 1 inside a cell) obeys an advection equation like this.

Here $\alpha_1$ is the volume fraction of fluid 1, $\mathbf{u}$ is the mixture velocity, and the Kapila term on the right corrects for the difference in compressibility between the two phases.

The trouble starts when you solve this with a shock-capturing scheme. The instant you interpolate cell-face values smoothly, the step in $\alpha_1$ that should stay razor-sharp smears a little each step. That numerical diffusion accumulates. Over many steps the interface spreads from a few cells to dozens, and key features — bubble shape, shock reflection positions — vanish entirely.

High-order schemes like WENO reduce the diffusion, but they introduce oscillations and are expensive in multiple dimensions. So the authors turn sideways. Instead of touching the scheme itself, they fix up the $\alpha$ field separately after every step.

Anti-diffusion: the sharpening equation#

The core is the regularization equation of Shukla et al. (2010). We evolve $\alpha$ for a few iterations in a pseudo-time $\tau$, not physical time.

$\varepsilon$ is the diffusion coefficient that regulates thickness, $a$ is the compression strength that tightens the interface, and $\hat{\mathbf{n}} = \nabla\alpha / |\nabla\alpha|$ is the interface normal direction.

The two terms play tug-of-war. The first is ordinary diffusion, so it spreads $\alpha$ out. The second, with the opposite sign, is anti-diffusion: it compresses the interface along its normal. Thanks to the $\alpha(1-\alpha)$ factor, compression acts only at the interface ($0<\alpha<1$) and switches off inside the pure fluids ($\alpha=0$ or $1$).

Solving the steady state ($\partial_\tau \alpha = 0$) in one dimension, the two fluxes balance.

The interface converges to a hyperbolic-tangent shape, and its thickness scales like $\varepsilon/a$. In other words, you choose the thickness. This is the "constant thickness" property the authors emphasize. When the interface always keeps the same number of cells, shocks and contact discontinuities can be captured stably.

In 1D: the tug-of-war between smearing and sharpening#

The same tug-of-war, written out, is the script below. It mimics the smear the transport step leaves behind, then runs the sharpening equation in pseudo-time.

import numpy as np
 
N, dx = 201, 1.0 / 200
x = np.linspace(0, 1, N)
 
def band_width(a, lo=0.05, hi=0.95):
    # number of cells counted as interface (thickness measure)
    return int(np.sum((a > lo) & (a < hi)))
 
def smear_once(a, D=0.16):
    # the transport scheme fattening the interface each step
    lap = np.zeros_like(a)
    lap[1:-1] = a[2:] - 2 * a[1:-1] + a[:-2]
    return a + D * lap
 
def sharpen_sweep(a, comp_a, eps, iters=200):
    # flux-form anti-diffusion regularization (Shukla 2010, eqs. (38)-(39))
    dtau = 0.9 * min(dx * dx / (2 * eps), dx / comp_a)
    for _ in range(iters):
        af = 0.5 * (a[:-1] + a[1:])          # face-centered alpha
        dA = a[1:] - a[:-1]
        s = np.sign(dA)
        J = comp_a * af * (1 - af) * s - eps * dA / dx   # face flux
        a[1:-1] = np.clip(a[1:-1] - dtau / dx * (J[1:] - J[:-1]), 0, 1)
    return a
 
# start from a sharp step -> let it smear 20 times
a = 0.5 * (1 + np.tanh((x - 0.5) / 0.012))
for _ in range(20):
    a = smear_once(a)
print("band after smearing:", band_width(a), "cells")     # -> wider
 
a = sharpen_sweep(a, comp_a=1.0, eps=0.004)
print("band after sharpening:", band_width(a), "cells")    # -> shrinks to target

The band is wide after the smear and drops to the few cells $\varepsilon/a$ sets. The scheme stays untouched. sharpen_sweep only sees the array. It never learns which Riemann solver filled it.

Defaults are comp_a=1.0 and eps=0.004, so the thickness scale $\varepsilon/a$ is 0.004. The steady profile's slope factor is $a/(2\varepsilon)=125$. The tanh that starts the script has denominator 0.012, so the target is steeper than the initial step. $dx=1/200=0.005$. The pseudo-time cap is the smaller of $dx^2/(2\varepsilon)=0.000025/0.008=0.003125$ and $dx/a=0.005$. Diffusion wins. Multiply by 0.9 and $d\tau$ is 0.0028125. Past that cap the anti-diffusion flips sign and the interface rings.

The compressive piece of the face flux scales with $\alpha(1-\alpha)$. That product peaks at 0.25 when $\alpha=0.5$ and is zero in a pure fluid, so compression switches off there. The band counter only tallies cells between 0.05 and 0.95. Double $\varepsilon$ alone and the target thickness doubles; double $a$ alone and it halves. The 200 steps in the opening are physical time. Six cells becoming 20 is about a factor of 3.3. D=0.16 inside smear_once and iters=200 inside the sharpener are different layers. The paper's one to three passes per step assumes the physical step barely diffuses. The Kapila source $\alpha_1 K\nabla\cdot\mathbf{u}$ pushes volume fraction directly wherever a shock spikes the divergence. That source is not numerical diffusion, and IST does not erase it. Unlike the separate pressure relaxation in the Saurel (2009) procedure, this pass only puts the cell count back.

What it means to hold thickness constant#

Making the interface thin is not the same as pinning its thickness to a constant. This is exactly where the authors draw a line against earlier approaches (Tiwari et al. 2013). Adding a sharpening function to the source term reduces the diffusion error, but the thickness wanders in space and time. Fixing $\varepsilon/a$ as a post-processing step, by contrast, yields the same thickness everywhere.

In compressible flow there is one more layer. If you fix up $\alpha$, the mixture density and energy must change with it to keep the thermodynamics consistent. So the authors do not merely correct $\alpha$; they redistribute the entire conservative-variable vector ($\alpha_1\rho_1$, $\alpha_2\rho_2$, momentum, energy) according to the mixture rules. Thanks to this "mixture-consistent regularization," velocity, pressure, and temperature equilibrium at the interface are not broken. The paper's 1D pure-interface advection test shows exactly this: when water and air move at the same velocity, turning on sharpening keeps pressure and temperature intact.

In 2D: keeping the bubble interface round#

It worked in 1D, so on to a 2D bubble. The paper's shock-bubble figures keep a circle only a few cells thick. Turn on numerical diffusion and the rim goes blurry; turn on sharpening and raise compression a, and the rim tightens back into a ring.

The teal ring tightens and the bubble center stays put. It neither grows nor shrinks. Only the cell count in the ring drops. Compression follows $\hat{\mathbf{n}}$, so the curvature does not crumple.

What I hit while reproducing it#

In pure fluid $|\nabla\alpha|\to 0$, so the normal divides by zero. The paper skips that because $a\,\alpha(1-\alpha)$ is zero there, but the code still needs a small number in the denominator or it writes NaN. A large $d\tau$ makes the interface ring. A CFL-type cap is what keeps it quiet. The paper's one to three passes per step assume the physical step barely diffuses. A coarse grid or a strong shock can need more.

One more thing. The elegance of "post-processing" has a price. Sharpening changes $\alpha$ independently of the governing equations, so at that instant the advection is not a solution of the original PDE. In exchange for constant thickness, local conservation near the interface can drift slightly. The authors minimize this with the mixture-consistent redistribution, but it is not exactly zero. OpenFOAM's interfoam uses a compression term of the same family ($\nabla\cdot(\alpha(1-\alpha)\mathbf{u}_r)$), which makes a useful yardstick when porting this technique to a compressible solver.

Leave ε/a free and the cell count drifts again#

A sharpening source cuts the diffusion error and still lets the thickness wander in space and time. Fixing $\varepsilon/a$ after the step is what keeps the same cell count beside a shock. Correct $\alpha$ alone, without splitting the conserved variables back by the mixture rules, and pressure and temperature equilibrium go first.

If $d\tau$ exceeds $dx^2/(2\varepsilon)$, anti-diffusion rings the interface instead of tightening it. The price of a constant thickness is that, for that instant, the advection is not a solution of the original PDE. The small conservation drift sits in the same family as the nonconservative terms caught by the Abgrall criterion.