[Paper Review] The Capillary Shackle Came Off and CFL 0.05 Was Still There — The Real Ceiling on VOF Interface Advection

The value a compressive scheme uses to stand the interface up is divided by the Courant number. Raise the time step and that quotient is the first thing to vanish.

There is a sentence in the conclusions of Janodet et al. (2025): "the proposed algorithm allows simulating realistic gas-liquid flows with time steps larger than the capillary time-step constraint, as long as other time-step constraints are satisfied." The emphasis is mine. One shackle came off by making surface tension implicit, and the largest CFL number the paper could actually run at was 0.05. The scheme carrying the colour function had taken hold of the time step instead. This post shows where that ceiling comes from, on the NVD diagram, and then puts a number on what it costs using a vortex advection experiment.

The capillary constraint itself, and the implicit treatment that breaks it, were covered in the post on implicit surface tension. What follows is only the sequel.

The value that stands the interface up comes from downwind#

In algebraic VOF — advecting the colour function directly instead of reconstructing the interface — there is exactly one reason interfaces thicken. Upwind face values always smear them. So compressive schemes pull the face value toward the downwind cell. Take the downwind value outright and the step is confined to a single cell.

The trouble is that the downwind value carries no boundedness guarantee. Once the colour function pokes below 0 or above 1, the density goes negative and the run is over. Some rule has to say how far the pull may go. That rule takes the Courant number as an argument, and that is the whole story here.

Try it in the simulation below.

Courant C0.05

blend gamma_f0.60transition cells0slab width0

Start at C = 0.05: the amber window fills almost the whole box, the blue HYPER-C curve pins itself to phi~_f = 1, and the slab keeps a two-cell edge forever. Drag C toward 0.9 and the ceiling min(1, phi~_D/C) folds down onto the dashed upwind diagonal — the window is a sliver, the pink face states have nowhere to sit, and the slab bleeds out over a dozen cells.

Leave Courant C at 0.05 and the amber admissible region on the left nearly fills the box, while the slab on the right keeps a two-cell edge. Push C to 0.9 and the ceiling folds down onto the dashed diagonal, which is upwind. The pink dots run out of places to sit, and the slab bleeds outward as time passes. Same scheme, same mesh. Only the time step grew.

On the NVD box, the Courant number lowers the ceiling#

The normalised variable diagram (NVD) takes one face and rescales the upwind cell $U$ , donor cell $D$ , and acceptor cell $A$ like this.

\tilde{\phi}_D = \frac{\phi_D - \phi_U}{\phi_A - \phi_U}, \qquad \tilde{\phi}_f = \frac{\phi_f - \phi_U}{\phi_A - \phi_U}

$\tilde{\phi}_D$ says where the donor cell sits between upwind and acceptor; $\tilde{\phi}_f$ says where the face value lands. $\tilde{\phi}_f = \tilde{\phi}_D$ is upwind and $\tilde{\phi}_f = 1$ is downwind.

Leonard's Convection Boundedness Criterion (CBC) nails down where the face value has to be under explicit time marching.

\tilde{\phi}_D \le \tilde{\phi}_f \le \min\!\left(1,\ \frac{\tilde{\phi}_D}{C}\right), \qquad 0 \le \tilde{\phi}_D \le 1

Here $C = u\,\Delta t/\Delta x$ is that face's Courant number. The meaning of the ceiling $\tilde{\phi}_D/C$ is direct. One step drains a volume $C$ out of the donor cell, and if the colour function carried in that volume exceeds what the cell held to begin with, the cell goes negative. That condition is $C \cdot \tilde{\phi}_f \le \tilde{\phi}_D$ , which rearranges into the inequality above.

At $C = 0.05$ the ceiling is $20\,\tilde{\phi}_D$ . A donor barely past 0.05 already lets the face value climb to 1. At $C = 0.8$ the ceiling is $1.25\,\tilde{\phi}_D$ , a thin ribbon just above the diagonal. The available compression scales as $1/C$ .

Why CICSAM sits at 0.01 and THINC/QQ at 0.05#

CICSAM blends two curves inside this box. One rides the ceiling itself — HYPER-C — and the other is the gentler ULTIMATE-QUICKEST. The blending weight $\gamma_f$ comes from the angle between the interface normal and the face vector. Normal to the face sends it toward HYPER-C; oblique sends it toward UQ. Compressing an oblique interface produces artificial staircase wrinkles, which is what the blend avoids.

The CICSAM blend button and the blend gamma_f slider above are that mixture. Lower $\gamma_f$ and the curve drops off the ceiling, thickening the slab immediately. You can also see that a one-dimensional aligned interface has $\gamma_f = 1$ , so CICSAM effectively collapses onto HYPER-C there.

That is where CICSAM's practical CFL ceiling near 0.01 comes from. On a real three-dimensional interface, $\gamma_f$ swings between 0 and 1, and only the HYPER-C component satisfies CBC on its own — the UQ part has to be re-limited separately. Keeping the blended result under the ceiling requires a small $C$ . To dodge that constraint, the paper uses THINC/QQ instead of CICSAM.

One tanh redraws the interface inside the cell#

THINC (Tangent of Hyperbola for INterface Capturing) skips the choice of a face value and draws the in-cell distribution outright. On a cell rescaled to $[0,1]$ with coordinate $\tilde{x}$ , it sets

\Phi(\tilde{x}) = \frac{1}{2}\left[1 + \gamma \tanh\!\big(\beta(\tilde{x} - \tilde{x}_c)\big)\right]

Here $\beta$ is the interface sharpness (usually near 2), $\gamma = \pm 1$ is the interface orientation read from the neighbours, and $\tilde{x}_c$ locates the tanh jump. $\tilde{x}_c$ is fixed by requiring the cell average $\bar{\phi}$ to be reproduced exactly, and it solves in closed form.

\tilde{x}_c = \frac{1}{\beta}\,\mathrm{artanh}\!\left(\frac{\cosh\beta - e^{\beta\gamma(2\bar{\phi}-1)}}{\sinh\beta}\right)

The amount crossing the face follows from integrating this curve over the departure region.

F_{i+1/2} = \int_{1-C}^{1} \Phi(\tilde{x})\,\mathrm{d}\tilde{x}

THINC/QQ adds a quadratic surface reconstruction on top, which handles curved interfaces better. The same tanh used on the shock-capturing side appears in TENO-THINC reconstruction.

The point is that this integral carries $C$ explicitly. As $C$ grows the integration window approaches the full cell width, and the scheme ends up transporting nothing more than the cell average. Because tanh never uses the downwind value, CBC holds automatically — but the decay of compressive power with $C$ is identical.

The time-step budget has three line items#

Now look at the whole budget. Resolving capillary waves explicitly costs

\Delta t_\sigma = \sqrt{\frac{(\rho_A + \rho_B)\,\Delta x^3}{2\pi\sigma}}

and the flow's CFL constraint $\Delta t_{\text{CFL}} = C_{\max}\Delta x / U$ sits alongside it. What the paper did was strike the first item from the budget. What remains is the $C_{\max}$ the interface advection scheme will allow.

Race two solvers to the same physical time below.

D / dx120U (flow speed)2.0e-2

interface CFL cap0.05breach factor S8.0

Leave U small — capillary-driven flow — and lane A is bound by dt_sigma while lane B runs away: that is the paper’s selling point. Now drag the interface CFL cap down to 0.01, the CICSAM value: lane B collapses back onto lane A even though surface tension is still implicit. Push the cap up to 0.5 instead and the volume-error readout is what pays for it.

Keep U small — capillary-driven flow — and lane A is pinned to $\Delta t_\sigma$ while lane B pulls ahead. That is the paper's selling point. Now drag interface CFL cap down to 0.01: lane B falls back beside lane A even though surface tension is still implicit. Push it to 0.5 instead and the volume-error readout underneath tells you the price.

One vortex, and the volume error it measures#

What actually degrades as $C$ grows? In directionally split advection, each sweep cannot see the solenoidal velocity field, so a dilatation correction is required.

\phi^{*}_{i} = \phi^{n}_{i} - \left(F_{i+1/2} - F_{i-1/2}\right) + \phi^{n}_{i}\left(C_{i+1/2} - C_{i-1/2}\right)

This term keeps a uniform $\phi = 1$ region from breaking apart under a single sweep. But in interface cells, $\phi^n_i$ differs from the value actually in play mid-sweep, and that difference survives as a volume error. I measured it by putting THINC advection on the Rider–Kothe single vortex with time reversal, $T=2$ .

from math import atanh, cos, cosh, exp, log, log1p, pi, sin, sinh
 
N, BETA, EPS = 40, 2.0, 1e-6
H = 1.0 / N
 
 
def lncosh(z):
    a = abs(z)
    return a + log1p(exp(-2.0 * a)) - log(2.0)
 
 
def thinc_slab(pbar, g, a, b):
    """Integrate the donor cell's tanh reconstruction over [a, b]."""
    s = g * (2.0 * pbar - 1.0)
    r = max(-0.999999, min(0.999999, (cosh(BETA) - exp(BETA * s)) / sinh(BETA)))
    xc = atanh(r) / BETA
    return 0.5 * ((b - a) + (g / BETA) * (lncosh(BETA * (b - xc)) - lncosh(BETA * (a - xc))))
 
 
def face_flux(pm, p0, pp, c):
    """Colour fraction crossing the face. p0 is the donor cell, c its Courant number."""
    if abs(c) < 1e-14:
        return 0.0
    g = 1.0 if pp > pm else (-1.0 if pp < pm else 0.0)
    if g == 0.0 or p0 < EPS or p0 > 1.0 - EPS:
        return c * p0
    return thinc_slab(p0, g, 1.0 - c, 1.0) if c > 0 else -thinc_slab(p0, g, 0.0, -c)
 
 
def line(col, vel, k):
    """One periodic 1-D sweep, dilatation correction included."""
    n, out = len(col), [0.0] * len(col)
    for i in range(n):
        cw, ce = vel[i] * k, vel[i + 1] * k
        fw = face_flux(col[(i - 2) % n], col[(i - 1) % n], col[i], cw) if cw > 0 else \
            face_flux(col[(i - 1) % n], col[i], col[(i + 1) % n], cw)
        fe = face_flux(col[(i - 1) % n], col[i], col[(i + 1) % n], ce) if ce > 0 else \
            face_flux(col[i], col[(i + 1) % n], col[(i + 2) % n], ce)
        out[i] = col[i] - (fe - fw) + col[i] * (ce - cw)
    return out
 
 
def run(courant, tend=2.0):
    uf = [[-sin(pi * i * H) ** 2 * sin(2 * pi * (j + .5) * H) for i in range(N + 1)] for j in range(N)]
    vf = [[sin(pi * j * H) ** 2 * sin(2 * pi * (i + .5) * H) for i in range(N)] for j in range(N + 1)]
    nstep = max(1, int(tend * max(abs(x) for r in uf for x in r) / (courant * H)))
    dt = tend / nstep
    f = [[1.0 if ((i + .5) * H - .5) ** 2 + ((j + .5) * H - .75) ** 2 < .15 ** 2 else 0.0
          for i in range(N)] for j in range(N)]
    f0, m0 = [r[:] for r in f], sum(sum(r) for r in f)
    for n in range(nstep):
        w = cos(pi * (n + .5) * dt / tend)          # Rider-Kothe time reversal
        for ax in ((0, 1) if n % 2 == 0 else (1, 0)):
            if ax == 0:
                f = [line(f[j], [x * w for x in uf[j]], dt / H) for j in range(N)]
            else:
                cols = [[f[j][i] for j in range(N)] for i in range(N)]
                vv = [[vf[j][i] * w for j in range(N + 1)] for i in range(N)]
                cols = [line(cols[i], vv[i], dt / H) for i in range(N)]
                f = [[cols[i][j] for i in range(N)] for j in range(N)]
    lo = min(min(r) for r in f)
    hi = max(max(r) for r in f)
    dm = (sum(sum(r) for r in f) - m0) / m0
    err = sum(abs(f[j][i] - f0[j][i]) for j in range(N) for i in range(N)) / m0
    return nstep, lo, hi, dm, err
 
 
print("   C   steps    min(f)     max(f)-1     dM/M    (dM/M)/C   shape err")
for c in (0.05, 0.1, 0.2, 0.4, 0.8):
    ns, lo, hi, dm, err = run(c)
    print(f"{c:5.2f} {ns:6d}  {lo:9.2e}  {hi - 1.0:9.2e}  {dm:8.2e}  {dm / c:8.4f}   {err:8.3e}")

   C   steps    min(f)     max(f)-1     dM/M    (dM/M)/C   shape err
 0.05   1595   2.23e-29  -6.03e-07  2.94e-03    0.0588   1.988e-01
 0.10    797  -6.51e-07  -6.33e-07  5.87e-03    0.0587   2.119e-01
 0.20    398  -1.62e-06  -5.27e-07  1.16e-02    0.0580   1.850e-01
 0.40    199  -3.87e-06   1.38e-07  2.33e-02    0.0583   1.771e-01
 0.80     99  -3.20e-06   1.43e-06  4.62e-02    0.0577   2.214e-01

Two things to read off. First, boundedness is intact. Undershoots sit at the $10^{-6}$ level, so THINC kept its promise. Second, the volume error is exactly proportional to $C$ . The fifth column, the fourth divided by $C$ , is nailed near 0.058. It holds to within 2% while $C$ grows by a factor of 16.

A volume error of 0.3% at $C = 0.05$ becomes 4.6% at $C = 0.8$ . Since this is a two-dimensional area, that converts to 2.3% in droplet diameter. For a surface-tension calculation this is fatal: curvature is the inverse of the radius, so the Laplace pressure jump is off by the same 2.3%.

The last column, shape error, wanders between 0.18 and 0.22 regardless of $C$ . That one is set by mesh resolution.

Does refining the mesh fix it?#

To find out where the 0.058 comes from, I reran $C = 0.2$ on $60^2$ instead of $40^2$ . The coefficient drops from 0.0580 to 0.0384. The ratio 0.66 is essentially the mesh spacing ratio $40/60$ . In other words,

\frac{\Delta M}{M} \approx 2.3\,C\,\Delta x \propto U\,\Delta t

The volume error is first order in time. Refine the mesh while holding $C$ and the error falls in proportion to $\Delta x$ . Refine the mesh while holding $\Delta t$ and $C$ grows by the same factor, leaving the error where it was. Enlarging the time step in interface advection is not free, and the bill scales precisely with $\Delta t$ .

That bill compounds with the parasitic currents problem. A volume off by 0.5% gives a curvature that is off, and a wrong curvature becomes an unbalanced surface-tension force that pollutes the velocity field all over again.

What has to change to turn 0.05 into 0.5#

The paper names two things in its own conclusions. The first is the robustness of the implicit height function: when the height function fails on an under-resolved interface, the curvature collapses wholesale. The second is the interface advection scheme — in the paper's words, improvements "should enable simulations with larger CFL numbers, bearing the potential to greatly improve the performance of the proposed algorithm."

Three routes look plausible. Make the advection itself implicit and escape the CBC's $1/C$ ceiling; switch to unsplit geometric VOF (PLIC) and delete the dilatation correction entirely; or separate interface reconstruction from advection the way anti-diffusion sharpening does. All three give up some of algebraic VOF's cheapness.

So what this paper delivers is not a new ceiling but a new bottleneck. The capillary constraint vacated its seat and the interface advection CFL sat down in it — and that one shrinks as $\Delta x$ , not as $\Delta x^{3/2}$ the way $\Delta t_\sigma$ does. Which means it gets relatively kinder as the mesh refines. That is useful information when choosing the next bottleneck to attack.