Why Droplet Simulations Get Trapped by Δt — Solving Surface Tension Implicitly

Double the mesh and the time step shrinks by about 2.8×. The physics did not get faster. On an interface with surface tension, the capillary time-step constraint is almost always the harshest limit.

Push the $\Delta t/\Delta t_\sigma$ slider past 0.9 and only the orange trace (explicit) blows its amplitude. The cyan trace (implicit) holds the same standing wave.

The dashed line is the analytic decay. Raise mesh Oh and both traces settle sooner. The 2025 coupling of Janodet, van Wachem and Denner puts that orange blow-up inside one matrix.

The paper under review is R. Janodet, B. van Wachem, F. Denner, "A fully-coupled algorithm with implicit surface tension treatment for interfacial flows with large density ratios", arXiv:2410.17757 (2024).

The capillary time-step shackle#

An incompressible interfacial-flow solver usually treats the diffusion term implicitly, erasing the viscous time-step limit. What remains as the strongest restriction is surface tension. In the form Brackbill first wrote in 1992,

$\rho_A, \rho_B$ are the densities of the two fluids, $\sigma$ is the surface-tension coefficient, and $\Delta x$ is the mesh spacing. The exponent $3/2$ is the problem. Halve $\Delta x$ to sharpen the interface and the time step drops by $2^{1.5}\approx 2.8$. The smaller the droplet, the more capillary-scale the flow, the more the step count explodes. With explicit surface tension there is no way out.

Why Δx^1.5 — the capillary-wave dispersion relation#

This restriction is really a CFL condition. A small sinusoidal perturbation of the interface oscillates as a capillary wave (a wave whose restoring force is surface tension). In the inviscid deep-water limit the dispersion relation is

where $\omega$ is the angular frequency and $k$ the wavenumber. The shortest wave the mesh can resolve is $k_{\max} = \pi/\Delta x$. If $\Delta t$ exceeds that wave's oscillation period, explicit integration cannot follow the oscillation and diverges. Computing $\omega_{\max}\Delta t_\sigma$ gives $\sqrt{(k_{\max}\Delta x)^3/(2\pi)} = \pi/\sqrt{2}\approx 2.22$. So $\Delta t_\sigma$ is exactly the CFL number of the fastest capillary wave.

Surface tension, made implicit — the core in a one-line oscillator#

You do not need the whole solver. A single capillary-wave mode reduces to a damped harmonic oscillator.

$\eta$ is the interface amplitude and $\gamma$ the viscous damping. Explicit surface tension is the same as integrating this oscillator with symplectic Euler, which diverges for $\omega\Delta t > 2$ — the very capillary CFL above. Evaluate the surface-tension term at the new time level (backward Euler) and it becomes unconditionally stable.

The orange trace up top breaks near 0.9 because that is the symplectic-Euler limit $\omega\Delta t = 2$ divided by the printed $\omega\Delta t_\sigma = 2.221$. And $2/2.221 \approx 0.900$. Above that line the explicit amplitude grows. Cyan still follows the dashed curve out to 3×. Raise mesh Oh (the viscous scale) and stronger damping settles both sooner.

A fake flow on a still interface shows up when this restoring force is discretised out of balance. That spot is parasitic current from a broken force balance.

The stability boundary, seen through the amplification factor#

The oscillator's amplification factor $|G|$ (the ratio by which the amplitude is multiplied each step) makes the boundary sharp. Thirty lines of numpy reproduce it. The toy problem is not a shock tube but the decay of a standing capillary wave.

import numpy as np
 
def capillary_omega(sigma, k, rho_hat):
    # dispersion relation: inviscid deep-water capillary wave omega^2 = sigma k^3 / (rhoA+rhoB)
    return np.sqrt(sigma * k**3 / rho_hat)
 
def dt_capillary(sigma, dx, rho_hat):
    # Eq.(4): Brackbill capillary time step ~ dx^{3/2}
    return np.sqrt(rho_hat * dx**3 / (2*np.pi*sigma))
 
def run_standing_wave(omega, gamma, dt, steps, implicit):
    eta, v = 1.0, 0.0
    peak = 0.0
    for _ in range(steps):
        if implicit:        # backward Euler: surface tension at the new level (unconditionally stable)
            v = (v - dt*omega**2*eta) / (1 + dt**2*omega**2 + 2*gamma*dt)
            eta = eta + dt*v
        else:               # symplectic Euler: explicit surface tension (omega*dt<2 limit)
            v = v - dt*omega**2*eta - dt*2*gamma*v
            eta = eta + dt*v
        peak = max(peak, abs(eta))
    return peak
 
sigma, rho_hat, L = 0.07, 1001.2, 1.0      # water/air
for N in (64, 128, 256):
    dx = L/N; k = np.pi/dx                  # shortest resolvable capillary wave
    om = capillary_omega(sigma, k, rho_hat)
    dts = dt_capillary(sigma, dx, rho_hat)
    print(f"N={N:4d}  dt_sigma={dts:.2e}  omega*dt_sigma={om*dts:.3f}")
# output: omega*dt_sigma = 2.221 (= pi/sqrt(2)) for every N -- the CFL IS the capillary limit
 
N = 128; dx = L/N; k = np.pi/dx
om = capillary_omega(sigma, k, rho_hat); dts = dt_capillary(sigma, dx, rho_hat)
dt = 1.5*dts; gamma = 0.1*om                # forced to 1.5x the capillary limit
pe = run_standing_wave(om, gamma, dt, 200, implicit=False)
pi_ = run_standing_wave(om, gamma, dt, 200, implicit=True)
print(f"dt=1.5*dt_sigma  explicit peak={pe:.1e}  implicit peak={pi_:.3f}")
# example output: explicit peak=3.6e+07  implicit peak=1.000

The listing sets $\sigma=0.07$, a density sum of $1001.2$, and a length of $1$. That $1001.2$ is water plus air. $0.07$ is the surface tension. At $N=64$, $\Delta x=1/64$ and the same formula gives $\Delta t_\sigma \approx 9.32\times 10^{-2}$. $N=128$ divides by $2^{3/2}\approx 2.828$, about $3.29\times 10^{-2}$. $N=256$ divides again, about $1.16\times 10^{-2}$. From 64 to 256 the spacing drops by 4 and the step count for the same physical time rises by 8, because $4^{3/2}=8$. omega*dt_sigma stays 2.221 on every mesh. The comment's example, at $\Delta t=1.5\Delta t_\sigma$, 200 steps, and $\gamma=0.1\omega$, lists an explicit peak of $3.6\times 10^{7}$ and an implicit peak of $1.000$. There $\omega\Delta t=1.5\times 2.221=3.3315$, past the symplectic-Euler limit of 2.

The horizontal axis of the $|G|$ curve is $\Delta t/\Delta t_\sigma$. The mesh and $\sigma$ sliders move $\Delta t_\sigma$ itself.

Explicit (orange) holds $|G|=1$ up to $\Delta t/\Delta t_\sigma \approx 0.9$, then crosses 1. Implicit (cyan) keeps $|G|<1$ at every step. Raise $\sigma$ by 4 and $\Delta t_\sigma$ falls by $\sqrt{1/4}=1/2$. The shackle is not only a mesh slope.

Full coupling is the point — colour function, pressure, velocity in one matrix#

The hard part of going implicit is that the surface-tension force $\mathbf{f}_\sigma = \sigma\kappa\nabla\psi$ is nonlinearly entangled with the colour function (the VOF variable marking fluid type from 0 to 1) $\psi$. The curvature $\kappa$ is itself set by derivatives of $\psi$. The paper applies a Newton linearisation of this term with respect to the new-level $\psi$, and solves the interface-advection equation together with the continuity and momentum equations as a single linear system $\mathbf{A}\boldsymbol{\phi}=\mathbf{b}$, in which $\boldsymbol{\phi}=(p, u, v, w, \psi)$ is coupled implicitly all at once. Denner et al.'s earlier algorithm was locked to a density ratio of 1. Janodet's contribution is to write continuity and momentum in conservative form and treat the density implicitly with respect to $\psi$ in the transient terms, holding the coupling together even at a 1000:1 water/air density ratio. The interface advection is discretised with THINC/QQ instead of CICSAM, keeping the interface sharp at large CFL. Unlike OpenFOAM's segregated PIMPLE, this only pays off in a block-coupled solver.

What I ran into reproducing it#

"Unconditionally stable" comes with a caveat. Even when implicit surface tension breaks the capillary limit, the largest stable $\Delta t$ is still finite. The Galusinski–Vigneaux visco-capillary upper limit sets it. At a density ratio of 1000 this limit is an order of magnitude tighter than at unit density ratio, and the paper itself reports stability only up to about $1.5\,\Delta t_\sigma$ in the inviscid limit. The oscillator model above simplifies that visco-capillary coupling, so it makes the implicit side look more generous than reality. Backward Euler is also stable at the cost of artificially damping fast capillary waves. In energy-conserving problems this numerical damping can be mistaken for physical decay, which is why the paper verifies force balance and energy conservation separately.

Stop at 1.5 Δt_σ, or drop the density ratio of 1000#

Implicit time stepping can walk past the capillary CFL and still hit a ceiling. In the inviscid limit the paper's stable window is about $1.5\,\Delta t_\sigma$. The comment's case at 1.5× and 200 steps, with an explicit peak of $3.6\times 10^{7}$ and an implicit peak stuck at 1.000, is the oscillator, not the solver. At a density ratio of 1000 the Galusinski–Vigneaux bound is an order of magnitude tighter than at unit ratio. The oscillator's generosity does not survive that far.

Backward Euler pays for the stability by numerically killing fast capillary waves. If the energy has to stay put, that damping gets mistaken for physics. The window after pushing $\Delta t$ further is open at five times and shut at ten.