[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.

Janodet et al. (2025) close with this line: "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 added here. Implicit surface tension took one shackle off. The largest CFL the paper actually ran was 0.05. The scheme carrying the color function had grabbed the time step. The NVD box shows where that ceiling folds, and a vortex advection puts a number on the volume it bills.

The capillary formula itself shows up later. What follows is only the $C_{\max}$ left after that line is struck.

The value that stands the interface up comes from downwind#

In algebraic VOF — advecting the color 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 color 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.

Raise Courant C from 0.05 to 0.9 and the admissible ceiling folds down onto the upwind diagonal.

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 color 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 integral carries $C$ in the limits, in plain sight. 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.

The gap between the two lanes shows up again in the race right after the printed table. U and interface CFL cap are the handles there.

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

Drop U on the capillary-driven side and lane A stays tied to $\Delta t_\sigma$ while lane B pulls ahead. Pull interface CFL cap down to 0.01 and B falls back beside A. Push it to 0.5 and the volume error underneath grows.

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.

At a cap of 0.05 the table shows 1595 steps and a volume error of $2.94\times 10^{-3}$ . Slide the cap toward 0.5 and the step count falls, while the ledger, 0.058 times $C$ , grows in step with it.

Boundedness holds. Undershoots sit at $10^{-6}$ , so THINC stayed inside 0 and 1. Volume error scales with $C$ . The fifth column, the fourth divided by $C$ , is nailed near 0.058. From $C = 0.05$ to $0.80$ , a factor of 16, the coefficients read 0.0588, 0.0587, 0.0580, 0.0583, 0.0577. The spread is inside 2%.

The integral window makes that proportion ordinary. The face amount is an integral whose lower limit is $1 - C$ , so the window width is $C$ itself. At $C = 0.05$ only the last 5% of the cell is carried. At $C = 0.8$ , 80% of it goes. The step count falls with it. At $T = 2$ and $N = 40$ , $C = 0.05$ takes 1595 steps, 0.1 takes 797, 0.2 takes 398, 0.4 takes 199, and 0.8 takes 99. Nearly a halving each time. Double $C$ and the window doubles, the step count halves, and the volume error per unit physical time doubles anyway.

$2.94\times 10^{-3}$ at $C = 0.05$ is 0.294%. $4.62\times 10^{-2}$ at $C = 0.8$ is 4.62%. Their ratio matches the ratio of the Courant numbers, 16. On a two-dimensional area that is about 2.3% in droplet diameter. Curvature is the inverse of the radius, so the Laplace pressure jump leans by that same 2.3%. With $\beta = 2$ the tanh jump is closed by the cell average alone, and a wider window still keeps the undershoot under $10^{-6}$ . Sharpness stays. Only the volume ledger follows $C$ .

The CBC ceiling reads the same $C$ . At $C = 0.05$ the ceiling is $20\,\tilde{\phi}_D$ , so a donor barely past 0.05 can already push the face value to 1. At $C = 0.8$ the ceiling is $1.25\,\tilde{\phi}_D$ . 20 divided by 1.25 is 16, the same factor the volume error grew by. Compression room shrinks as $1/C$ while the ledger error grows as $C$ .

Shape error in the last column wanders through 0.199, 0.212, 0.185, 0.177, 0.221, independent of $C$ . That floor belongs to the $40^2$ mesh.

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$ .

Stack that bill on a curvature error and the ledger is wrong twice. A volume off by 0.5% shifts the radius by about that much, and curvature, being the inverse of the radius, leans by the same fraction. An unbalanced surface tension dirties the velocity, and the dirty velocity shakes $C$ on the next sweep. The opening left by implicit surface tension, where a larger $\Delta t$ still closes stays shut if interface advection keeps $C$ at 0.05.

What has to change to turn 0.05 into 0.5#

The paper grabs the robustness of the height function first. On an under-resolved interface a collapsed height function kills the curvature outright. Advection is next. In the paper's words, a scheme that allows a larger CFL is where the performance still has room. Make the advection itself implicit and step out of the CBC's $1/C$ ceiling, switch to unsplit geometric VOF (PLIC) and drop the dilatation correction, or pick the face value apart from the transport the way MUSCL-THINC-BVD reconstruction does. All three hand back some of algebraic VOF's cheapness.

What is left is not a new ceiling but a new bottleneck. The capillary constraint vacated the seat and the interface-advection CFL sat down. $\Delta t_\sigma$ shrinks as $\Delta x^{3/2}$ . This one shrinks only as $\Delta x$ . A finer mesh gives it relatively more air. Whether to accept ten times the 0.294% volume error at $C = 0.05$ , or to keep the 1595 steps, is the next call.