Raising the time step 5x ran 1.9x faster; 10x ran no faster at all — the window implicit surface tension opens

Three days into one oscillating droplet#

A 2D oscillating droplet has been running for three days. The velocities are slow and the mesh is not large. The time step, however, is $10^{-6}$ s. The moment surface tension is treated explicitly, the time step stops being set by the flow and starts being set by capillary waves.

The usual suggestion at this point is to treat surface tension implicitly — with the interface at the new time instance inside the linear system. Break the constraint and the time step can go up by a factor of 5, or 10. Does that turn three days into one?

The answer is "up to about 5". The 2025 fully-coupled algorithm of Janodet, van Wachem and Denner measured both ends of that window at a density ratio of 1000. A stability limit closes it from above; the cost per step closes it from below. This article follows where those two walls stand, and why refining the mesh stops buying accuracy well before you expect it to.

What ties the time step is a capillary wave, not the flow speed#

With surface tension on the interface, there is a shortest capillary wave the mesh can resolve. Its wavelength is $\lambda_\sigma = 2\Delta x$. Use a time step larger than the time that wave needs to cross a cell, and an explicit surface-tension term blows up. Denner and van Wachem write the constraint as

where $\rho_A, \rho_B$ are the densities of the two fluids, $\sigma$ is the surface tension coefficient and $\Delta x$ is the cell size. The exponent $3/2$ is the problem. Halve the mesh spacing and the time step shrinks by 2.8. That closes in faster than the $\Delta x^1$ of an advective CFL condition. Diffusion terms can be solved implicitly and removed from the ledger; surface tension resisted that for a long time. Why the constraint exists and how the implicit treatment works is covered in an earlier post on the capillary time-step constraint.

Push dt/dt_sigma from 0.5 to 8. With the curvature error on, the blue amplitude barely leaves the grey dashed curve.

This is a capillary wave between two fluids at a density ratio of 1000, decaying under viscosity. The dashed grey line is Prosperetti's analytical solution, the blue line is the amplitude the discrete solver produces. Raise lambda/dx and the two curves close. With curvature dx^0.5 engaged, sweep dt/dt_sigma from 0.5 to 8 and watch how little the error moves — that is the subject of the next section but one.

A second ceiling stands where the first one fell#

Treating surface tension implicitly does let you pass $\Delta t_\sigma$. It does not give you an arbitrarily large time step. Following the analysis of Galusinski and Vigneaux, Denner et al. write the remaining limit as a competition between two time scales.

Here $\tau_{vc} = \hat\mu \lambda_\sigma / \sigma$ is the visco-capillary time scale and $\tau_\sigma = \sqrt{\hat\rho \lambda_\sigma^3 / \sigma}$ the capillary time scale, with $\hat\rho = \rho_A + \rho_B$ and $\hat\mu = \mu_A + \mu_B$. The constants $a_1, a_2$ are case dependent; $a_1 = 1/(16\pi)$ with $a_2 = 0$ recovers $\Delta t_\sigma$ exactly.

The ratio of the two scales is the mesh Ohnesorge number.

For $\mathrm{Oh}_{\Delta x} \ll 1$ inertia dominates and $\Delta t^{*} \propto \tau_\sigma$; for large $\mathrm{Oh}_{\Delta x}$ viscosity dominates and $\Delta t^{*} \propto \tau_{vc}$. The second regime is the one that pays: a large dynamic viscosity or a short capillary wave opens the limit considerably.

The painful number is the density ratio. For the stationary droplet (Laplace equilibrium) case, the limit in the $\mathrm{Oh}_{\Delta x} \ll 1$ regime was $1.5\,\Delta t_\sigma$ at a density ratio of 1000. The same family of algorithms reached $15\,\Delta t_\sigma$ at unit density ratio. A factor of ten disappeared. In the large-$\mathrm{Oh}_{\Delta x}$ regime the gap is an order of magnitude as well. Realistic gas-liquid density ratios narrow the window.

Halving the mesh once multiplies $\Delta t_\sigma$ by $2^{3/2}\approx 2.83$. Raising $\lambda/\Delta x$ from 25 to 200 cuts the spacing by eight. $8^{3/2}\approx 22.6$, so an explicit surface-tension run needs about 23 times as many steps for the same physical time. Three days becomes more than two months.

Going implicit does not reopen that. On a stationary droplet at density ratio 1000 the ceiling closes at $1.5\,\Delta t_\sigma$, one tenth of the $15\,\Delta t_\sigma$ the same family used at unit density ratio. $\Sigma=5$ sits outside that ceiling. Under inertia, five times diverges and ten times is farther out. Only when viscosity is large enough that $\Delta t^{*}\propto\tau_{vc}$ does a factor of five enter the stable range. With $a_2=0$ and $a_1=1/(16\pi)$, $\Delta t^{*}$ falls back to the original $\Delta t_\sigma$. The term that opens the window is $a_2\tau_{vc}$.

Refining the mesh eightfold cut the error only in half#

The second validation case is a decaying capillary wave. Density and viscosity ratios both 1000, Laplace number $\mathrm{La} = \rho\lambda\sigma/\mu^2 = 300$, meshes $\lambda/\Delta x = \{25, 50, 100, 200\}$ and time steps $\Delta t/\Delta t_\sigma = \{0.5, 2, 8\}$. The distance from the analytical solution is measured as an $L_2$ norm of the amplitude.

What stands out in the resulting table is not the size of the error but the order of convergence. Most entries sit between 0.46 and 0.95. The same problem at unit density ratio converges at second order. Refining the mesh by a factor of eight cuts the error roughly in half.

The paper does not blame the temporal discretisation. It blames the interface transport, in two lines. The interface-capturing scheme in use is at best second-order accurate. Curvature is a second derivative of the colour function, so it loses two orders. Curvature is therefore at best zeroth order. On a sufficiently fine mesh the convergence order of the amplitude error goes to zero: the error settles on a constant and stops falling.

The owner of the convergence rate, in Python#

The argument reduces to a single damped oscillator. In the linear regime the wave amplitude obeys $A'' + 2\nu k^2 A' + \omega_0^2 A = 0$. What the solver sees is not $\omega_0$ but a frequency carrying the curvature error, $\omega_{num} = \omega_0\sqrt{1 + C(\Delta x/\lambda)^q}$. Vary $q$ alone, march with the trapezoidal rule, and read off the $L_2$ norm against the analytical solution together with its order.

import math
 
SIGMA, RHO_HAT, LAMBDA, K, LA = 1.0, 1.0, 2*math.pi, 1.0, 300.0
MU = math.sqrt(RHO_HAT * LAMBDA * SIGMA / LA)
NU = MU / RHO_HAT
A0, T_END = LAMBDA / 100.0, 25.0
 
 
def capillary_omega(dx, q, c_kappa=0.6):
    """the frequency the discrete solver actually sees, curvature error O(dx^q)"""
    w0 = math.sqrt(SIGMA * K**3 / RHO_HAT)
    return w0 * math.sqrt(1.0 + c_kappa * (dx / LAMBDA) ** q)
 
 
def analytic_amplitude(t):
    """exact solution of A'' + 2*nu*k^2*A' + w0^2*A = 0"""
    w0 = math.sqrt(SIGMA * K**3 / RHO_HAT)
    g = NU * K**2
    wd = math.sqrt(w0**2 - g**2)
    return A0 * math.exp(-g*t) * (math.cos(wd*t) + g/wd * math.sin(wd*t))
 
 
def march_amplitude(dt, w, n_steps):
    """march [A, A'] with the trapezoidal (Crank-Nicolson) rule"""
    g = NU * K**2
    a, v, hist = A0, 0.0, [A0]
    for _ in range(n_steps):
        h = 0.5 * dt
        rhs_a, rhs_v = a + h*v, v + h*(-w**2 * a - 2*g*v)
        det = (1 + 2*g*h) + h*h*w**2
        a = ((1 + 2*g*h) * rhs_a + h * rhs_v) / det
        v = (-h * w**2 * rhs_a + rhs_v) / det
        hist.append(a)
    return hist
 
 
def l2_amplitude(hist, dt):
    """L2 error norm of the amplitude (Eq. 61 of the paper)"""
    acc = 0.0
    for i, a in enumerate(hist):
        w = 0.5 if i in (0, len(hist)-1) else 1.0
        acc += w * (a - analytic_amplitude(i*dt))**2 * dt
    return math.sqrt(acc / (len(hist)-1) / dt) / A0
 
 
def order_of(e_coarse, e_fine):
    return math.log(e_coarse / e_fine) / math.log(2.0)
 
 
for label, q in [("curvature error ~ dx^2", 2.0), ("curvature error ~ dx^0.5", 0.5)]:
    print(f"\n{label}")
    print("lam/dx |  dt/dt_s=0.5        dt/dt_s=2          dt/dt_s=8")
    prev = {}
    for n in [25, 50, 100, 200]:
        dx = LAMBDA / n
        dt_sigma = math.sqrt(RHO_HAT * dx**3 / (2*math.pi*SIGMA))
        w = capillary_omega(dx, q)
        row = []
        for s in [0.5, 2.0, 8.0]:
            dt = s * dt_sigma
            e = l2_amplitude(march_amplitude(dt, w, int(T_END/dt)), dt)
            tag = "  (-- )" if s not in prev else f" ({order_of(prev[s], e):4.2f})"
            row.append(f"{e:.3e}{tag}")
            prev[s] = e
        print(f"{n:6d} | " + "  ".join(row))

curvature error ~ dx^2
lam/dx |  dt/dt_s=0.5        dt/dt_s=2          dt/dt_s=8
    25 | 5.763e-04  (-- )  5.692e-04  (-- )  1.680e-02  (-- )
    50 | 1.516e-04 (1.93)  6.731e-05 (3.08)  2.023e-03 (3.05)
   100 | 3.885e-05 (1.96)  2.548e-05 (1.40)  2.345e-04 (3.11)
   200 | 9.833e-06 (1.98)  8.085e-06 (1.66)  2.506e-05 (3.23)
 
curvature error ~ dx^0.5
lam/dx |  dt/dt_s=0.5        dt/dt_s=2          dt/dt_s=8
    25 | 7.565e-02  (-- )  7.483e-02  (-- )  6.051e-02  (-- )
    50 | 5.449e-02 (0.47)  5.439e-02 (0.46)  5.266e-02 (0.20)
   100 | 3.897e-02 (0.48)  3.896e-02 (0.48)  3.874e-02 (0.44)
   200 | 2.775e-02 (0.49)  2.775e-02 (0.49)  2.773e-02 (0.48)

Take S from 1 to 5 and the blue lane finishes early. Push it to 10 and the lane barely moves.

Drop Oh_dx and the red wall $\Delta t^{*}$ slides left until the fast lane stops. That is a different wall from the lower table, where stretching the time step sixteenfold leaves the error unchanged.

The upper table shows second order at small time steps. The third order in the $\Delta t/\Delta t_\sigma = 8$ column is not a bonus: since $\Delta t \propto \Delta x^{3/2}$, a second-order temporal error falls as $\Delta x^3$.

The lower table is the paper's situation. The order locks near 0.5. More telling, the three columns carry essentially the same values. Cutting the time step by sixteen leaves the error untouched. What sets the floor on accuracy is the curvature, not the temporal discretisation. The measured orders of 0.46 to 0.95 in the paper land exactly on this picture.

The lam/dx = 200 row reads 2.775e-02, 2.775e-02, 2.773e-02. Stretching dt/dt_s from 0.5 to 8, sixteen times, moves only the last digit, and the orders are 0.49, 0.49, 0.48, inside that 0.46 to 0.95 band. The upper table's dt/dt_s = 0.5 column falls from 5.763e-04 to 9.833e-06 with orders 1.93, 1.96, 1.98. Dropping the curvature-error exponent q from 2 to 0.5 is what changes the floor. Swapping the time-step column is not.

How that curvature enters as a body force continues in contact angle and CSF. Refining the mesh still leaves a curvature that is at best zeroth order.

5x gave 1.9x, 10x gave nothing#

The third case is the damped oscillation of a 2D elliptical droplet. It starts with a major axis of 0.15 m and a minor axis of 0.1 m, oscillates in the $n=2$ mode and is damped by viscous stresses. The applied time step is the smaller of two constraints.

$\Sigma$ is the factor by which the capillary constraint is breached. The paper ran $\Sigma \in \{2, 5, 10\}$ with the maximum CFL number held at 0.05.

On accuracy first: the oscillation-frequency error was about 3% for $\Sigma = 2$ and $\Sigma = 5$, smaller than the roughly 4.5% obtained with an explicit surface-tension treatment at the same resolution. $\Sigma = 10$, by contrast, failed to follow the kinetic-energy decay. With the surface-tension-driven interface motion no longer resolved in time, the paper notes, the formal second-order accuracy of the temporal scheme cannot be expected.

The cost numbers are this article's title. Raising $\Sigma$ from 2 to 5 — a factor of 2.5 — cut the total wall clock time by a factor of 1.9. Pushing on to 10 did not keep that gain. The reduced computational time per step rose substantially. The reason is single: the larger the time step, the more slowly the nonlinear procedure converges within each step. Fewer steps, more work in each. Raising the step by 2.5 cut the wall clock only by 1.9. The missing 0.6 was eaten by extra Newton iterations on every step.

The walls stand in different places#

One case carries five separate limits and optima, and no two are alike.

An algorithm that breaks $\Delta t_\sigma$ erases exactly one row of that table. The others stay. How the force balance survives on a stationary droplet is covered in the post on parasitic currents.

$\Sigma=10$ stays stable and still misses the decay#

If $\mathrm{Oh}_{\Delta x}$ is far below 1, the stability window closes near $1.5\,\Delta t_\sigma$. The stationary droplet at density ratio 1000 is that case. Starting at $\Sigma=10$ is already outside the window. Count how many steps cut one period of the mode that matters. Frequency error at $\Sigma=2$ and $5$ was about 3%, smaller than the roughly 4.5% of the explicit treatment. $\Sigma=10$ did not diverge, and still failed to follow the kinetic-energy decay.

When the nonlinear iteration count grows in proportion to $\Sigma$, the wall clock has already passed its minimum. From 2 to 5, a factor of 2.5, time dropped by 1.9. At 10 that gain came back. Fewer steps, and the Newton iterations filled the gap. CFL is held at 0.05 in this run. On some meshes advection binds the time step before $\Sigma$ can.