A Motionless Droplet Melted to Half Its Mass — The Curvature Term Conservative Allen

The sharpening term in conservative Allen–Cahn is not a device for holding interface thickness — it is the term that cancels curvature-driven shrinkage exactly.

The velocity field was zero and the droplet still halved#

Should an interface with curvature but no flow move on its own? That question is easy to skip when a phase-field model (one that treats an interface as a transition layer of finite thickness) first goes into a code. Run a 2D droplet initialized at rest through the Allen–Cahn equation for 4000 steps and the answer shows up. The radius drops from 20.083 to 14.286, and 49.4% of the total volume fraction disappears. The velocity field stays exactly zero the whole time.

This post locates that loss term by term. Set classical Allen–Cahn beside the conservative form and the difference is a single curvature term. Then the same term is followed into a lattice Boltzmann (LB) code, down to the moment it has to sit in. By the end, the W in your code and the source term are two separate promises, not one.

Term-by-term comparison of the two Allen–Cahn forms#

Item	Classical Allen–Cahn	Conservative Allen–Cahn
Form	reaction–diffusion	flux divergence
Interface maintenance	double well $\phi(1-\phi)(1-2\phi)$	counter flux $\frac{4}{W}\phi(1-\phi)\hat{\mathbf n}$
$\int\phi\,dV$	not conserved	conserved at the discrete level
1D equilibrium	tanh	the same tanh
Curvature response	$v_n=-M\kappa$	zero

The two equations read as follows.

\frac{\partial \phi}{\partial t} = M\left(\nabla^2\phi - \frac{16}{W^2}\phi(1-\phi)(1-2\phi)\right)

\frac{\partial \phi}{\partial t} + \nabla\cdot(\phi\mathbf u) = \nabla\cdot\left[M\left(\nabla\phi - \frac{4\phi(1-\phi)}{W}\hat{\mathbf n}\right)\right]

Here $\phi$ is the volume fraction, $M$ the mobility (the coefficient that sets how fast the interface relaxes), $W$ the interface thickness, $\hat{\mathbf n}=\nabla\phi/|\nabla\phi|$ the interface normal, and $\kappa=\nabla\cdot\hat{\mathbf n}$ the curvature. The M and W sliders move the left half-disk (classical) and the right half-disk (conservative) on separate radii.

R0 = 20 luW = 4.0 luM = 0.050

t = 0  | classical R = 20.000, m/m0 = 1.0000 | conservative R = 20.000, m/m0 = 1.000000000

The left half-disk is the classical form, the right half the conservative one. Raise M and only the left half melts faster, tracking the dotted $\sqrt{R_0^2-2Mt}$ curve. Change W and the left-hand shrink rate barely moves — the culprit is curvature, not thickness.

An identity the tanh interface hands over for free#

Why $W$ enters the two equations in different shapes comes out of the 1D equilibrium. For a planar interface the solution is this tanh.

\phi(x) = \frac{1}{2}\left(1 + \tanh\frac{2x}{W}\right)

Differentiating it hands over two identities for free.

\frac{d\phi}{dx} = \frac{4}{W}\phi(1-\phi), \qquad \frac{d^2\phi}{dx^2} = \frac{16}{W^2}\phi(1-\phi)(1-2\phi)

The first says the slope can be written using $\phi$ alone. No interface coordinate is needed: a local value is enough to build the sharpening flux. The second is where the classical double-well coefficient $16/W^2$ comes from. Substitute both into either equation and the right-hand side is zero. In 1D the two models share the same equilibrium.

Does the curvature term survive or cancel?#

The difference appears the moment the interface bends. Take a normal coordinate $s$ and the Laplacian splits into two pieces.

\nabla^2\phi = \frac{\partial^2\phi}{\partial s^2} + \kappa\frac{\partial\phi}{\partial s}

Substitute the identities above. On the classical right-hand side, $\partial^2\phi/\partial s^2$ and the double-well term cancel exactly, and only the curvature piece survives.

\left.\frac{\partial\phi}{\partial t}\right|_{\text{classical}} = M\kappa\frac{\partial\phi}{\partial s} = M\kappa\frac{4\phi(1-\phi)}{W}

That one term is $v_n=-M\kappa$ , curvature-driven motion. For a droplet $\kappa=-1/R$ , so $dR/dt=-M/R$ , which integrates to $R^2=R_0^2-2Mt$ . With radius 20 and $M=0.05$ the lifetime is 4000 time units.

The conservative form builds the same term a second time inside the divergence. Expanding the counter-flux divergence gives

\nabla\cdot\left[\frac{4\phi(1-\phi)}{W}\hat{\mathbf n}\right] = \frac{16}{W^2}\phi(1-\phi)(1-2\phi) + \frac{4\phi(1-\phi)}{W}\kappa

so it delivers the double-well term and a curvature term together. The classical form only had the first piece. Adding the second erases $M\kappa\,\partial\phi/\partial s$ . The place where an absent phase leaked into the interface in three-phase Allen–Cahn had the same structure. Which term you drop decides in advance which conserved quantity breaks.

Same tanh, two update rules, 4000 steps#

Both update rules run on an axisymmetric 1D grid. Solving only the radial direction brings the curvature in as $1/r$ automatically. The Laplacian and the sharpening flux are both written as face fluxes, so the divergence part is conservative at the discrete level in both models. Every bit of the classical form's mass loss therefore comes from the double-well term alone.

import math
 
NR, DR, R0, W, M, DT, NSTEP = 256, 0.5, 20.0, 4.0, 0.05, 0.5, 4000
rc = [(i + 0.5) * DR for i in range(NR)]          # cell centers
rf = [(i + 1.0) * DR for i in range(NR - 1)]      # interior faces
 
def tanh_disk(w):
    return [0.5 * (1.0 - math.tanh(2.0 * (r - R0) / w)) for r in rc]
 
def disk_mass(f):
    return sum(r * v for r, v in zip(rc, f)) * DR
 
def disk_radius(f):
    return math.sqrt(2.0 * disk_mass(f))
 
def peak_slope(f):
    return max(abs(f[i + 1] - f[i]) / DR for i in range(NR - 1))
 
def well_drive(v, w):
    return (16.0 / w ** 2) * v * (1.0 - v) * (1.0 - 2.0 * v)
 
def flux_div(F):
    out = [0.0] * NR
    for i in range(NR):
        lo = rf[i - 1] * F[i - 1] if i > 0 else 0.0
        hi = rf[i] * F[i] if i < NR - 1 else 0.0
        out[i] = (hi - lo) / (rc[i] * DR)
    return out
 
def grad_faces(f):
    return [(f[i + 1] - f[i]) / DR for i in range(NR - 1)]
 
def sharp_faces(f, w, eps=1e-12):
    F = []
    for i in range(NR - 1):
        g = (f[i + 1] - f[i]) / DR
        pf = 0.5 * (f[i] + f[i + 1])
        F.append(g - (4.0 / w) * pf * (1.0 - pf) * g / (abs(g) + eps))
    return F
 
def step_classical(f, w, dt):
    lap = flux_div(grad_faces(f))
    return [v + dt * M * (lap[i] - well_drive(v, w)) for i, v in enumerate(f)]
 
def step_conservative(f, w, dt):
    div = flux_div(sharp_faces(f, w))
    return [v + dt * M * div[i] for i, v in enumerate(f)]
 
f0 = tanh_disk(W)
m0 = disk_mass(f0)
fk, fc = list(f0), list(f0)
print("    t   R_classic  R_conserv   mass_classic  mass_conserv")
for s in range(NSTEP + 1):
    if s % 1000 == 0:
        print("%5.0f  %9.3f  %9.3f  %13.9f  %13.9f"
              % (s * DT, disk_radius(fk), disk_radius(fc),
                 disk_mass(fk) / m0, disk_mass(fc) / m0))
    fk = step_classical(fk, W, DT)
    fc = step_conservative(fc, W, DT)
print("curvature law  sqrt(R0^2-2*M*t) = %.3f" % math.sqrt(R0 ** 2 - 2 * M * NSTEP * DT))
print("interface width  init %.3f  classic %.3f  conserv %.3f"
      % (1 / peak_slope(f0), 1 / peak_slope(fk), 1 / peak_slope(fc)))
print("mass drift  classic %+.3e  conserv %+.3e"
      % (disk_mass(fk) / m0 - 1, disk_mass(fc) / m0 - 1))

Radius holds at 20.083, mass error at 1e-14#

    t   R_classic  R_conserv   mass_classic  mass_conserv
    0     20.083     20.083    1.000000000    1.000000000
  500     18.803     20.083    0.876607772    1.000000000
 1000     17.429     20.083    0.753192944    1.000000000
 1500     15.936     20.083    0.629689171    1.000000000
 2000     14.286     20.083    0.506056483    1.000000000
curvature law  sqrt(R0^2-2*M*t) = 14.142
interface width  init 4.021  classic 4.032  conserv 4.000
mass drift  classic -4.941e-01  conserv -2.631e-14

The classical 14.286 sits 1.0% off the analytic 14.142. That gap comes from the grid spacing of 0.5 and the area-equivalent radius definition. The trend follows $R^2=R_0^2-2Mt$ exactly. The conservative form does not budge from 20.083 to three decimals.

The mass numbers are blunter. Classical gives $-4.941\times10^{-1}$ , conservative $-2.631\times10^{-14}$ . The latter is accumulated double-precision round-off and will not shrink when the grid is refined. It is not an error that refinement reduces — it is an error that was never there. Interface thickness stays near 4 in both models, so thickness is not what separates them.

At t=0 both radii are 20.083 and both mass ratios are 1.000000000. At t=500 only the classical side has moved, to 18.803 and 0.876607772. At t=1000 it is 17.429 and 0.753192944. At t=1500, 15.936 and 0.629689171. At t=2000, 14.286 and 0.506056483. The conservative column stays on 20.083 and 1.000000000. The analytic radius at that same time, $\sqrt{R_0^2-2Mt}$ , is 14.142. The gap 14.286 minus 14.142 is 0.144, the 1.0% above. The mass drift $-4.941\times10^{-1}$ is 1 minus 0.506056483. Set beside $-2.631\times10^{-14}$ , the printed scales are 1e-1 and 1e-14. Interface widths are 4.021 at the start, 4.032 classical, 4.000 conservative. Radius 20 and $M=0.05$ give a lifetime of 4000 time units because $R^2=R_0^2-2Mt$ reads $400=0.1\,t$ . The printed t=2000 is half that life, and 14.142 is the analytic radius there.

In LB the term sits in the source, not the equilibrium#

The standard way to solve the conservative form in LB is to carry a second distribution function $h_i$ for $\phi$ . On a D2Q9 lattice the update reads

h_i(\mathbf x+\mathbf c_i\Delta t,\ t+\Delta t) = h_i(\mathbf x,t) - \frac{h_i-h_i^{eq}}{\tau_\phi} + F_i

h_i^{eq} = \phi\,w_i\left(1+\frac{\mathbf c_i\cdot\mathbf u}{c_s^2}\right), \qquad M = c_s^2\left(\tau_\phi-\frac{1}{2}\right)\Delta t

where $h_i$ is the distribution carrying $\phi$ , $\tau_\phi$ the dimensionless relaxation time, $w_i$ the lattice weights, and $c_s^2=1/3$ . The equilibrium holds advection only. The sharpening term is supplied by the source $F_i$ , not by the equilibrium.

F_i = \Delta t\,w_i\left(1-\frac{1}{2\tau_\phi}\right)\frac{4\phi(1-\phi)}{W}\frac{\mathbf c_i\cdot\hat{\mathbf n}}{c_s^2}

The Chapman–Enskog expansion (sorting small departures from equilibrium order by order in the lattice-spacing ratio) asks for two moment conditions. The zeroth moment must satisfy $\sum_i F_i = 0$ . Since $\phi=\sum_i h_i$ , a nonzero zeroth moment means the source creates or destroys volume fraction. The first moment must reproduce the counter flux exactly. Because $\sum_i w_i(\mathbf c_i\cdot\hat{\mathbf n})\mathbf c_i = c_s^2\hat{\mathbf n}$ , the $c_s^2$ in $F_i$ cancels and the target value comes out. This is the same species of constraint as why the D2Q9 equilibrium is a polynomial, except here the list of lines to respect ends at two.

Break the zeroth moment and mass leaks#

A common accident is writing $F_i$ as if it were an equilibrium — putting a 1 inside the bracket, as $w_i\,\mathrm{amp}\,(1+\mathbf c_i\cdot\hat{\mathbf n}/c_s^2)$ . The first moment still comes out right. Interface thickness therefore stays healthy and the code runs quietly.

CX = [0, 1, -1, 0, 0, 1, -1, -1, 1]
CY = [0, 0, 0, 1, -1, 1, 1, -1, -1]
WT = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]
CS2 = 1.0 / 3.0
W, TAU, PHI, NX, NY = 4.0, 0.8, 0.5, 0.6, 0.8
 
def sharp_source(phi, nx, ny, w, tau, spill=False):
    amp = (1.0 - 0.5 / tau) * 4.0 * phi * (1.0 - phi) / w
    out = []
    for k in range(9):
        proj = (CX[k] * nx + CY[k] * ny) / CS2        # projection for the first moment
        out.append(WT[k] * amp * (proj + 1.0 if spill else proj))
    return out
 
amp = (1.0 - 0.5 / TAU) * 4.0 * PHI * (1.0 - PHI) / W
print("target            sum c_i F_i = [%+.6f %+.6f]" % (amp * NX, amp * NY))
for tag, spill in (("moment-safe", False), ("0th-moment spill", True)):
    F = sharp_source(PHI, NX, NY, W, TAU, spill)
    sx = sum(CX[k] * F[k] for k in range(9))
    sy = sum(CY[k] * F[k] for k in range(9))
    print("%-17s sum F_i = %+.3e   sum c_i F_i = [%+.6f %+.6f]" % (tag, sum(F), sx, sy))

target            sum c_i F_i = [+0.056250 +0.075000]
moment-safe       sum F_i = +8.674e-19   sum c_i F_i = [+0.056250 +0.075000]
0th-moment spill  sum F_i = +9.375e-02   sum c_i F_i = [+0.056250 +0.075000]

W_init = 3.0 luW_model = 6.0 lu

t = 0  |  W_eff = 3.000 lu  |  max |net flux| = 0.00e+0  | mass / mass0 = 1.000000

Turn source on and the zeroth moment breaks. Give W_init and W_model different values and the red net-flux curve spikes, then falls onto the dashed tanh. W_eff stays locked to W_model while the mass curve walks away from 1. A thickness reading alone never shows that hole.

The first moments of the two variants agree to six decimals. Only the zeroth moment moves, from zero to $9.375\times10^{-2}$ . That value is exactly $(1-1/2\tau_\phi)\cdot4\phi(1-\phi)/W$ , which means volume fraction is being injected at every interface node on every step.

The printed first moment is +0.056250 and +0.075000 in both modes. With $\phi=0.5$ , $\tau_\phi=0.8$ and $W=4$ , $(1-1/(2\tau_\phi))\cdot 4\phi(1-\phi)/W = 0.09375$ . Times $n_x=0.6$ that is 0.056250. Times $n_y=0.8$ that is 0.075000. The zeroth moment splits into +8.674e-19 and +9.375e-02. 9.375e-02 is 0.09375 written with an exponent. A width of 4.000 can still be right while that sum is not zero, and then mass leaks.

A table for choosing between the two models#

Situation	Classical	Conservative
Planar interface relaxation	suitable	suitable
Holding a static droplet for long runs	unsuitable	suitable
Modeling grain growth or curvature-driven motion	suitable	unsuitable
Reporting total volume fraction	needs correction	as is
LB implementation cost	no source term	needs $\hat{\mathbf n}$

The third row is the one that matters. The conservative form is not always the right answer. If curvature-driven motion is what you are trying to model, the classical form is the correct physics. That term is only an error when the goal is advecting a fluid interface. Either way $\hat{\mathbf n}$ has to be computed, so the conservative form leaves you guarding a division where $|\nabla\phi|$ approaches zero in the bulk.

What you sign for when you type W into an interface capturer#

Typing W = 4 into a code fixes more than a thickness. The double-well coefficient $16/W^2$ and the counter-flux coefficient $4/W$ both descend from that one number. Change only one of them and the equilibrium stops being a tanh. Whether both coefficients are derived from the same $W$ is the first thing to check.

The next check is the conserved quantity. Print the total volume fraction every step and 1e-14 separates from 1e-1 at a glance. Thickness looks healthy in both models, so watching thickness filters nothing. A conservative form whose interface thickness stayed put when mobility went up by 10 already split thickness from mass. In an LB implementation the check is one line: print $\sum_i F_i$ and confirm it is zero.

A Motionless Droplet Melted to Half Its Mass — The Curvature Term Conservative Allen–Cahn Erases