Why a Droplet at Rest Starts to Flow — Parasitic Currents and Well-Balanced Surface Tension

The 300 Pa quoted for a 1 mm drop is the sphere's 2/R. The code's 2D cylinder is half of that, 144 Pa, and a 3% error leaves 4.3 Pa of it unmatched.

The droplet was at rest. The forces matched. The run started, and small vortices grew on the surface. Nobody pushed. The fluid moved anyway. That spurious velocity field is a parasitic current. Tallois et al. (2025) shut it off with a well-balanced surface tension.

Drop curvature error ε to zero and the velocity field vanishes. Raise the same slider and the vortices come back along the interface.

curvature error ε: 0.35 surface tension σ: 0.072 N/mdroplet radius R: 0.50 mm

Laplace jump ΔP = σκ = σ/R = 144.0 Pa

kinetic energy Σ½|u|² = 0.00e+0

The frame at ε = 0 is the balanced state. Switch to 1D · face and the same error locks the flow onto the grid axes. multiD · nodal spreads that same error into a much weaker field. The kinetic energy on the right is the number for that difference.

The Laplace Law — Pressure Jump from Curvature#

A curved interface, the boundary where two fluids meet, makes a pressure jump. That is the Laplace law.

\Delta P = \sigma \kappa

$\Delta P$ is the pressure difference across the interface. $\sigma$ is the surface tension coefficient, force per unit length. $\kappa$ is the curvature, the reciprocal of the radius.

For a 2D droplet, $\kappa=1/R$ . A smaller $R$ means a larger jump. A diameter of 1 mm, so $R=0.5\,\mathrm{mm}$ , read as a sphere, has $\kappa=2/R$ and $\Delta P=\sigma\cdot 2/R$ . With water's $\sigma=0.072\,\mathrm{N/m}$ that is $0.072\times 4000=288\,\mathrm{Pa}$ , the "about 300 Pa" in the usual line. The code is a 2D cylinder, so $\kappa=1/R=2000\,\mathrm{m^{-1}}$ and the jump is $0.072/0.0005=144\,\mathrm{Pa}$ . Half the sphere.

Interior pressure pushes out. Surface tension pulls in. When the two match, the droplet stays put. That state is the equilibrium.

Why a Droplet at Rest Flows#

The computer does not hit that balance on the cells.

Numerical methods turn surface tension into a body force. That is CSF, continuum surface force. The interface is smeared across a few cells and the force is spread over the band. The force scales with $\kappa$ . Curvature on a grid goes wrong easily. Pushing the same force out to a wall contact angle is a separate closure. What breaks here is earlier. Curvature and pressure stop canceling.

A small curvature error and the surface-tension force misses the pressure gradient. The leftover pushes the fluid.

\mathbf{r} = \nabla P - \sigma \kappa_{\text{num}} \nabla z

$\mathbf{r}$ is the residual force. It is zero when the balance holds. $\kappa_{\text{num}}$ is the numerical curvature. $\nabla z$ is the gradient of the volume fraction, the liquid share of a cell.

Exact $\kappa_{\text{num}}$ gives $\mathbf{r}=0$ . A few percent of error and $\mathbf{r}$ is not zero. That residual builds vortices along the interface. The ε on the figure above is those few percent.

Well-Balanced — Matching the Discrete Balance Exactly#

If the discretization keeps the static solution as a solution, it is well-balanced.

Tallois et al. treat surface tension as a non-conservative product and put the term inside the Riemann solver. Pressure gradient and surface-tension force, computed with the same discrete rule, cancel cell by cell. When a problem that should stay quiet instead jumps in pressure, the test for that break is in where a non-conservative term ruins a static solution. On the droplet, that term is the curvature.

The curvature stencil splits the result too. A face-based 1D stencil sees only neighboring faces. The residual lines up with the grid axes, the velocity spikes that way, and the interface shakes. A node-based multidimensional stencil sees around the node. The residual spreads and the parasitic current weakens. Surface tension is multidimensional, so the stencil has to be too. The kinetic-energy gap between 1D · face and multiD · nodal on the figure above is that choice.

The Droplet That Moves Again — Rayleigh Oscillation#

A parasitic current is a spurious flow to remove. Surface tension also makes a real motion.

An ellipse in a gas oscillates. It starts where the surface energy is highest. Surface tension returns it to a circle. Kinetic energy peaks there. Inertia overshoots and the ellipse comes back. With no damping, the exchange continues.

The period follows a modified Rayleigh formula.

\omega^2 = (n^3 - n)\, \frac{\sigma}{(\rho_l + \rho_g) R^3}

$\omega=2\pi/T$ is the angular frequency. $n$ is the mode. An ellipse is $n=2$ . $\rho_l,\rho_g$ are the liquid and gas densities. $R$ is the mean radius at rest.

For $n=2$ , $n^3-n=8-2=6$ . With the code's $\sigma=0.072$ and $R=5\times 10^{-4}$ , $R^3=1.25\times 10^{-10}$ and $6\sigma/R^3=3.456\times 10^{9}$ . $\omega^2$ is that number divided by $\rho_l+\rho_g$ . This post does not give the two densities as numbers. A larger denominator lengthens the period. Raising $\sigma$ grows the numerator and shortens the period.

An ideal droplet oscillates without damping. A diffusive method shaves the amplitude down to a circle within a few periods. The paper holds that numerical damping down with a low-Mach correction.

Seeing the Balance Error in Python#

On a radial cut, smear the volume fraction with a tanh and compute the residual force.

import numpy as np
 
sigma = 0.072            # N/m, water-air
R = 0.5e-3               # droplet radius (m)
kappa_exact = 1.0 / R    # 2D cylinder curvature
 
def volume_fraction(r, radius, width):
    # smeared liquid fraction across the interface (0=gas, 1=liquid)
    return 0.5 * (1.0 - np.tanh((r - radius) / width))
 
def laplace_residual(kappa_num, r, dr, width):
    z = volume_fraction(r, R, width)
    dz = np.gradient(z, dr)                         # d z / d r
    f_st = sigma * kappa_num * dz                   # CSF surface tension force
    p = np.cumsum(sigma * kappa_exact * dz) * dr    # Laplace equilibrium pressure
    dp = np.gradient(p, dr)
    return dp - f_st                                # residual force (0 means well-balanced)
 
def parasitic_energy(residual, dt=1e-6, steps=200):
    # the leftover force accelerates the fluid: u <- u + dt * residual
    u = np.zeros_like(residual)
    for _ in range(steps):
        u += dt * residual
    return 0.5 * np.sum(u * u)
 
r = np.linspace(0.2e-3, 0.8e-3, 400)
dr = r[1] - r[0]
width = 3 * dr
 
res_ok  = laplace_residual(kappa_exact,        r, dr, width)   # exact curvature
res_bad = laplace_residual(kappa_exact * 1.03, r, dr, width)   # 3% error
 
print("exact curvature : max|r| = %.2e,  KE = %.2e"
      % (np.abs(res_ok).max(),  parasitic_energy(res_ok)))
print("3%% error       : max|r| = %.2e,  KE = %.2e"
      % (np.abs(res_bad).max(), parasitic_energy(res_bad)))

The run prints this.

exact curvature : max|r| = 4.1e-10,  KE = 2.3e-19
3% error : max|r| = 4.8e+05,  KE = 1.7e+05

surface tension σ: 340 N/mliquid density ρ_l: 100 kg/m³numerical damping: 0.40

Rayleigh period T = 2π/ω = 81.2 ms

ω = √[(n³−n)σ / ((ρ_l+ρ_g)R³)] , n = 2

Exact curvature has a peak residual of $4.1\times 10^{-10}$ . The kinetic energy $2.3\times 10^{-19}$ is effectively zero. The 3% case, $\kappa$ multiplied by 1.03, has a peak residual of $4.8\times 10^{5}$ and a kinetic energy of $1.7\times 10^{5}$ . Their ratio is $1.7\times 10^{5}/2.3\times 10^{-19}\approx 7.4\times 10^{23}$ , about 23.9 orders. That is the "24 orders" line.

The $4.8\times 10^{5}$ is the 3% of the pressure jump, turned into a gradient. $\sigma\cdot 0.03\cdot\kappa=0.072\times 60=4.32\,\mathrm{Pa}$ is the pressure error left on the interface. Adding $dt=10^{-6}$ two hundred times gives $2\times 10^{-4}$ . Times the peak residual, the velocity scale is $2\times 10^{-4}\times 4.8\times 10^{5}=96$ . The force a 3% error leaves piles that scale of speed onto the interface.

Set numerical damping to zero and the ellipse keeps its amplitude. Raise it and the collapse to a circle speeds up. Raise σ or lower ρ_l and the period $T$ shortens. The numerator above is $6\sigma/R^3$ , so doubling $\sigma$ doubles that numerator and multiplies $\omega$ by $\sqrt{2}$ .

The radius 0.5 mm sits in the middle of the samples. The code places 400 points from 0.2 mm to 0.8 mm. The span is 0.6 mm over 399 gaps, so one gap is about 1.50 μm. The interface is three gaps wide, about 4.51 μm. The radius is about 111 times that width, and the smeared band is thin beside the drop.

Sphere curvature 2/R with the same surface tension gives a jump of 288 Pa, and 3% of it is 8.64 Pa. The cylinder in the code has half the curvature, so the leftover pressure is half, 4.32 Pa. Feed the sphere formula into this two-dimensional residual and the error force reads twice as large.

Doubling the radius multiplies the cube by 8, so the angular frequency falls by the square root of 8, about 2.83, and the period grows by that factor. Doubling surface tension only multiplies the frequency by the square root of 2, about 1.41. Radius pushes the period harder. The density sum is still not a number. The period in seconds does not close here.

Leave the 3% in place and a larger time step still fails#

If the droplet at rest is not quiet, match the curvature before measuring a period. While the residual $4.8\times 10^{5}$ is alive, a larger $\Delta t$ only stacks the velocity scale of 96 faster. Well-balanced is what drops that residual to $4.1\times 10^{-10}$ and earns the right to talk about the time step. Leave the stencil on faces and the axial flow takes that right back.

The capillary time-step limit is still there after the residual dies. How far $\Delta t$ can be pushed is not answered by this residual. How wide the time step actually opens once surface tension is treated implicitly carries the next failure.