Two droplets coalesce: a sharp interface that changes topology
Two drops sit in a void, their rims 0.02 apart. Surface tension pulls them together; they touch, fuse, and the merged drop begins to round itself. The Cahn–Hilliard tutorial does the same physics the capturing way — the interface is a field, so topology costs nothing. This is the other route: the liquid is the mesh, its surface is a real boundary, and the topology changes because the mesh is re-triangulated on its own moving nodes.
with \(\mathbf w\) the mesh velocity and \(\mathbf D\) the symmetric gradient. There is no Dirichlet condition anywhere and no pressure pin: the traction sets the pressure level by itself.
The surface is a term, and so is its motion
The capillary traction is a boundary term; the kinematic condition is a geometry term. Both are ordinary
entries in the jno.fem([...]) list:
capillary = SIGMA * div_G(vs) # T n = -sigma H n, integrated by parts onto the test function
uf = u.bind(x=xs, y=ys).freeze(np.zeros((n_node, 2))) # the solved velocity, delivered each step
fem = jno.fem([momentum, continuity, capillary, u(x0, y0)[0] - 0.0, u(x0, y0)[1] - 0.0,
xs.d(ts) - uf[0], ys.d(ts) - uf[1]]) # the surface is Lagrangian: dX/dt = u
xs.d(ts) - u moves the mesh, and inside the weak form xi.d(ti) is \(\mathbf w\), so the convective
velocity is \(\mathbf u-\mathbf w\) — the ALE form. The nodal values ride with their vertices; nothing is
projected between meshes.
How the topology changes
traj = fem.solve(nonlinear=jno.solve.newton(direct=True),
adapt=jno.solve.remesh(alpha=1.2, every=1))
alpha= re-triangulates the nodes — a Delaunay triangulation filtered to the triangles whose
circumradius is below \(\alpha h\), the alpha shape of the Particle Finite Element Method. When the gap
between two bodies closes below about \(2\alpha h\), the bridging triangles survive the filter and two meshes
become one.
Re-triangulating is not enough on its own, though, and the reason is worth knowing: Delaunay already maximises the smallest angle for the points it is given, so it cannot repair a bad point distribution. A body that stretches thins its nodes; one that merges inherits whatever spacing the two surfaces happened to have. Left alone, this run degraded from a 41.5° smallest angle to 8.3°, with the largest cell 57× the smallest and the slivers lined up along the plane where the drops joined. So nodes are inserted and dropped as well — an edge longer than \(1.5h\) gains its midpoint (on a boundary edge that midpoint lies on the edge, so the liquid area is untouched), and an interior node crowding within \(0.55h\) of another is removed, never a boundary node. With that, the smallest angle stays above 18.6° for the whole run and no cell ends below 20°, for 8 % more nodes.
When the node set is unchanged the P1 state carries across by identity; when nodes are inserted or dropped it is interpolated from the old mesh, the same transfer a remesh uses.
What it produces
300 steps — 1.06 capillary times \(t_\sigma=\sqrt{\rho R^3/\sigma}\) — in about 10 minutes on this machine (CPU), with 130 reconnections and the mesh growing from 246 to 266 nodes.
| \(t/t_\sigma\) | bodies | area | width | height | neck | surface | nodes |
|---|---|---|---|---|---|---|---|
| 0.00 | 2 | 0.24972 | 0.820 | 0.400 | 0.078 | 2.509 | 246 |
| 0.27 | 1 | 0.25915 | 0.816 | 0.390 | 0.300 | 2.113 | 255 |
| 0.53 | 1 | 0.25898 | 0.804 | 0.387 | 0.387 | 2.044 | 258 |
| 1.06 | 1 | 0.25873 | 0.782 | 0.508 | 0.508 | 1.948 | 266 |
The neck opens from 0.078 to 0.508 — as tall as the drop itself — and the free surface shrinks from 2.509 to 1.948, which is what surface tension is for. Peak speed is 9.5, against the capillary velocity \(\sqrt{\sigma/\rho R}=7.07\). The surface never lengthens in a single step (worst change +0.00 %), so nothing is being bitten out of it.
Area is two statements, and conflating them would flatter the scheme. Merging adds liquid: the bridge fills the 0.02 gap the drops left between them, +3.6 % here and nearly all of it in the first few steps. That is the contact model, not an error in the solver. What the solver owes is conservation once the topology has settled, and there it holds to −0.16 % from just after the merge to the end, across 130 re-triangulations that each rebuild the problem.
Speed \(|\mathbf u|\) at \(t=0\), \(0.53\,t_\sigma\) and \(1.06\,t_\sigma\), each drawn on the mesh that step ran on. The drops start at rest; the flow is fastest near the neck as it opens.
Three things that are silently wrong if you get them wrong
The viscous term must be \(2\eta\,\mathbf D(\mathbf u)\!:\!\mathbf D(\mathbf v)\). The \(\eta\,\nabla\mathbf u\!:\!\nabla\mathbf v\) form — correct behind Dirichlet walls, and what the stabilised-flow and Cahn–Hilliard tutorials use — makes the natural condition a pseudo-traction. A rigidly rotating drop then loses its spin: \(\omega\) decays from 2.0 to 1.17 in 0.1 time units.
The SUPG/PSPG \(\tau\) must be scaled for this regime. A capillary drop is nearly inviscid and nearly stagnant, the opposite of advection-dominated flow. With the unscaled recipe a drop does not oscillate at all — the fitted frequency collapses from Lamb's 63.9 to 0.23 — while a spurious \(n=4\) mode grows sevenfold. Scaling \(\tau\) by \(10^{-4}\) recovers the physics.
Reconnection must run every step. At every=2 the mesh tangles (after 352 s at \(h=0.05\), 224 s at
\(h=0.04\)): the neck opens faster than the triangulation is refreshed.
Scope — what this is not
- One capillary time, not full relaxation. The merged drop is still a rounded dumbbell at \(1.06\,t_\sigma\) (width 0.783 against the final circle's 0.566). Relaxing all the way takes about \(5\,t_\sigma\) — roughly an hour here, because each reconnection re-assembles the problem (~2 s/step at 246 nodes). Cells are compile-time constants in the assembler; making this cheap needs connectivity as a runtime, padded argument.
- The time step is tighter than the usual estimate. \(\sqrt{\rho h^3/2\pi\sigma}\) suggests \(8\times10^{-4}\) at \(h=0.03\); \(5\times10^{-4}\) already tangles, and \(2\times10^{-4}\) tangles at \(h=0.05\) even reconnecting every step. This run uses \(10^{-4}\).
- Merging is a mesh-length contact model. Bodies join when their gap is comparable to the element size (\(\approx 2\alpha h\)), not by resolving film drainage. A larger \(\alpha\) merges earlier and adds area through wider bridges: at \(1.06\,t_\sigma\) the liquid has grown +3.3 % at \(\alpha=1.2\), +3.6 % at 1.6 and +6.4 % at 2.2 — most of it the bridge that fills the initial 0.02 gap.
- Two mechanisms keep the filter steady, and both are needed. Re-deciding every cell each step makes a triangle sitting near \(\alpha h\) drop out and return, so the surface loses and regains wedges: with one threshold and no node management, the worst one-step change in surface length was +25 %. Holding an existing cell to a wider threshold (hysteresis) brought that to +13 % over 2 steps of 300; adding node insertion and removal brought it to +0.00 %, because the edges that sat near the threshold are the over-long ones that now get split.
- 2-D and P1 only. Reconnection refuses P2 — a new edge bridging two bodies has its midpoint in the void — and 3-D, where Delaunay plus an alpha filter leaves sliver tetrahedra.
The same machinery is pinned by tests: tests/test_fem_free_surface_tension.py (the Laplace jump, per-drop
curvature, a rigid rotation, and the pseudo-traction trap), tests/test_fem_drop_oscillation.py (Lamb's
frequency to −0.4 % and the viscous decay to +4 %), tests/test_fem_droplets_merge.py and
tests/test_reconnect_alpha.py (the merge threshold and the alpha filter itself).
"""Two droplets, meshed as BODIES IN A VOID, touch and become one -- the sharp-interface route.
The liquid is the mesh. Its surface is a real boundary carrying the capillary traction, and it moves
with the fluid:
rho (du/dt|_X + ((u - w).grad)u) = -grad p + div(2 eta D(u)), div u = 0 (ALE, in the liquid)
T.n = -sigma H n on the whole surface (a traction, no BC)
dX/dt = u on the whole surface (it is Lagrangian)
with w the MESH velocity and D(u) the symmetric gradient. Two things make the topology change possible.
The surface rides with the fluid, written as geometry terms `xs.d(ts) - u`, so `jno.fem` marches the mesh
and the flow together. And `remesh(alpha=...)` re-decides which nodes form elements each step -- the alpha
shape of the Particle Finite Element Method -- so when the gap between two bodies closes below ~2 alpha h,
the bridging triangles survive the filter and two meshes become one.
Re-connecting alone is not enough: Delaunay already maximises the smallest angle FOR THE POINTS IT IS
GIVEN, so it cannot repair a point distribution that the motion has ruined. Without node management this
run fell from a 41.5 degree smallest angle to 8.3, the largest cell reaching 57x the smallest, the slivers
lined up along the plane where the drops joined. So `alpha=` also inserts a midpoint on any edge longer
than 1.5 h and drops interior nodes that crowd within 0.55 h -- never boundary nodes, which would cut
liquid away. The smallest angle then holds above 18.6 degrees for 8 % more nodes. Where the node set is
unchanged the P1 state carries by identity; where it changes it is interpolated from the old mesh.
Three details are not optional, and each is measured in the tests:
* the viscous term must be `2 eta D(u):D(v)`; `eta grad u : grad v` is a pseudo-traction that is correct
only behind Dirichlet walls, and at a free surface it silently destroys the motion;
* the SUPG/PSPG products must be written so each one has a single temporal order (a scalar times the
strong residual distributes; `inner(a, u.t + ...)` is refused). This file used to also say `tau` must be
scaled down, because the unscaled recipe over-damped the drop tenfold -- that was an assembly defect, and
the drop now rings at Lamb's frequency either way (tests/test_fem_drop_oscillation.py). TAU_SCALE is
kept at 1e-4 below only so this example's recorded output is unchanged;
* reconnection must run EVERY step (`every=1`). The neck opens fast, and at `every=2` the mesh tangles.
"""
import os
os.environ.setdefault("XLA_PYTHON_CLIENT_MEM_FRACTION", "0.6")
import jax
jax.config.update("jax_enable_x64", True)
import numpy as np # noqa: E402
import jno # noqa: E402
inner, symgrad = jno.np.inner, jno.np.symgrad
ddot = lambda a, b: inner(a, b, n_contract=2) # noqa: E731 A:B
RHO, ETA, SIGMA = 1.0, 0.01, 10.0 # density, viscosity, surface tension
NU = ETA / RHO
R, GAP, H = 0.20, 0.02, 0.04 # two drops of radius R, rims GAP apart, element size H
DT, N_STEPS = 1e-4, 300 # the capillary limit is tighter than sqrt(rho h^3 / 2 pi sigma) suggests
C_I, TAU_SCALE = 36.0, 1e-4 # not required for damping (see above); kept so the recorded output is unchanged
T_CAP = np.sqrt(RHO * R**3 / SIGMA) # capillary time: 0.0283
cx = R + GAP / 2.0
d = (jno.shape.disk(-cx, 0.0, R, size=H) | jno.shape.disk(cx, 0.0, R, size=H)).domain(time=(0.0, N_STEPS * DT, N_STEPS + 1))
n_node = int(np.asarray(d.mesh.points).shape[0])
u, v = d.fem_symbols(value_shape=(2,), names=("u", "v"), order=1) # P1/P1, stabilised
p, q = d.fem_symbols(names=("p", "q"), order=1)
xi, yi, ti = d.variable("interior", split=True)
xs, ys, ts, nx, ny = d.variable("boundary", normals=True, split=True) # the free surface, with normals
x0, y0, _t0 = d.variable("initial", split=True)
B = dict(x=xi, y=yi, t=ti)
ub, vv, pp, qq = u.bind(**B), v.bind(**B), p.bind(**B), q.bind(**B)
vs = v.bind(x=xs, y=ys)
D = lambda w: symgrad(w, [xi, yi]) # noqa: E731
ndv = lambda f, i: nx * f.x[i] + ny * f.y[i] # noqa: E731 (grad f_i).n
div_G = lambda f: f.x[0] + f.y[1] - (nx * ndv(f, 0) + ny * ndv(f, 1)) # noqa: E731 surface divergence
# The MESH velocity w is `coord.d(t)` read inside the weak form, so u - w is the ALE convective velocity.
# Componentwise: a vector cannot be assembled from two scalar expressions.
c0, c1 = ub[0] - xi.d(ti), ub[1] - yi.d(ti)
conv = lambda i: c0 * ub.x[i] + c1 * ub.y[i] # noqa: E731
G = d.cell_metric
gG = lambda a: inner(a, inner(G, a, n_contract=1), n_contract=1) # noqa: E731
tau = jno.lag(TAU_SCALE * ((2.0 / DT) ** 2 + gG(ub) + C_I * NU**2 * inner(G, G, n_contract=2)) ** -0.5)
r0, r1 = ub.t[0] + conv(0) + pp.x / RHO, ub.t[1] + conv(1) + pp.y / RHO # lap(u) vanishes on P1
momentum = (
RHO * (ub.t[0] * vv[0] + ub.t[1] * vv[1])
+ RHO * (conv(0) * vv[0] + conv(1) * vv[1])
+ 2.0 * ETA * ddot(D(ub), D(vv)) # the TRUE stress: grad:grad would be a pseudo-traction
- pp * (vv.x[0] + vv.y[1])
+ tau * ((c0 * vv.x[0] + c1 * vv.y[0]) * r0 + (c0 * vv.x[1] + c1 * vv.y[1]) * r1) # SUPG
)
continuity = -qq * (ub.x[0] + ub.y[1]) - tau * (qq.x * r0 + qq.y * r1) # PSPG: no pressure pin needed
capillary = SIGMA * div_G(vs) # T n = -sigma H n, integrated by parts onto the test function
uf = u.bind(x=xs, y=ys).freeze(np.zeros((n_node, 2))) # the solved velocity, delivered each step
fem = jno.fem(
[
momentum,
continuity,
capillary, # the traction also fixes the pressure LEVEL: there is no Dirichlet anywhere
u(x0, y0)[0] - 0.0, # both drops start at rest
u(x0, y0)[1] - 0.0,
xs.d(ts) - uf[0], # the surface is Lagrangian: dX/dt = u
ys.d(ts) - uf[1],
]
)
traj = fem.solve(
nonlinear=jno.solve.newton(direct=True), # a saddle step wants a sparse-direct Newton
adapt=jno.solve.remesh(alpha=1.2, every=1), # re-triangulate the moved nodes EVERY step
)
def bodies(k):
"""How many separate liquid bodies frame ``k`` has (connected components of the triangulation)."""
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import connected_components
pts, cells = np.asarray(traj.meshes[k][0]), np.asarray(traj.meshes[k][1])
r = np.concatenate([cells[:, 0], cells[:, 1], cells[:, 2]])
c = np.concatenate([cells[:, 1], cells[:, 2], cells[:, 0]])
used = np.zeros(len(pts), bool)
used[cells.reshape(-1)] = True
_n, lab = connected_components(coo_matrix((np.ones(r.size), (r, c)), shape=(len(pts),) * 2), directed=False)
return int(np.unique(lab[used]).size)
def boundary_loops(k):
"""How many closed boundary loops frame ``k`` has. A simply-connected drop has exactly one per body;
an interior hole adds another -- which is what the alpha filter leaves behind if it deletes stretched
interior triangles, and what a connected-component count alone will NOT notice."""
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import connected_components
from jno.utils.solver.fem_adapt import _boundary_edges_from_triangles
pts = np.asarray(traj.meshes[k][0])
e = np.asarray(_boundary_edges_from_triangles(np.asarray(traj.meshes[k][1])))
used = np.zeros(len(pts), bool)
used[e.reshape(-1)] = True
_n, lab = connected_components(coo_matrix((np.ones(len(e)), (e[:, 0], e[:, 1])), shape=(len(pts),) * 2), directed=False)
return int(np.unique(lab[used]).size)
def perimeter(k):
"""Total length of the free surface at frame ``k``."""
from jno.utils.solver.fem_adapt import _boundary_edges_from_triangles
pts = np.asarray(traj.meshes[k][0])
e = np.asarray(_boundary_edges_from_triangles(np.asarray(traj.meshes[k][1])))
return float(np.sum(np.linalg.norm(pts[e[:, 0]] - pts[e[:, 1]], axis=1)))
def shape_of(k):
"""``(area, x-extent, y-extent, neck width)`` of frame ``k``."""
pts, cells = np.asarray(traj.meshes[k][0]), np.asarray(traj.meshes[k][1])
t = pts[cells]
a, b = t[:, 1] - t[:, 0], t[:, 2] - t[:, 0]
area = float(np.abs(0.5 * np.sum(a[:, 0] * b[:, 1] - a[:, 1] * b[:, 0])))
mid = np.abs(pts[:, 0]) < 0.02
return area, np.ptp(pts[:, 0]), np.ptp(pts[:, 1]), (np.ptp(pts[mid, 1]) if mid.any() else 0.0)
last = len(traj.states) - 1
print(f"{n_node} nodes, {fem.dofs} dofs, {N_STEPS} steps = {N_STEPS * DT / T_CAP:.2f} capillary times")
for k in (0, last // 4, last // 2, last):
ar, ex, ey, neck = shape_of(k)
print(
f" t = {float(traj.times[k]) / T_CAP:4.2f} t_cap: {bodies(k)} body(ies), area {ar:.5f}, "
f"width {ex:.4f}, height {ey:.4f}, neck {neck:.4f}"
)
reconnects = sum(h["remeshed"] for h in fem.adapt_history)
area0, area1 = shape_of(0)[0], shape_of(last)[0]
print(f"{reconnects} reconnections; area {area0:.5f} -> {area1:.5f} ({100 * (area1 / area0 - 1):+.2f} %)")
assert bodies(0) == 2, "the drops did not start as two separate meshes"
assert bodies(last) == 1, "the drops never merged"
assert shape_of(last)[3] > 0.5 * shape_of(last)[2], "the neck never opened"
# Area is two statements, not one. MERGING ADDS LIQUID: the bridge fills the gap the drops left between
# them, which is what a mesh-length contact model does (+3.3 % here, most of it in the first few steps) --
# it is not a conservation error, and quoting it as one would flatter or damn the scheme by accident.
# What the scheme owes is conservation AFTER the topology settles, and that is the tight bound.
_a_merged = shape_of(last // 4)[0]
assert abs(area1 / area0 - 1.0) < 0.05, f"the merge added {100 * (area1 / area0 - 1):.1f} % of liquid"
assert abs(area1 / _a_merged - 1.0) < 0.01, f"area drifted {100 * (area1 / _a_merged - 1):.2f} % after the merge"
for _k in range(last + 1): # the liquid must stay SOLID: no voids punched by the alpha filter
assert boundary_loops(_k) == bodies(_k), f"frame {_k} has interior holes ({boundary_loops(_k)} loops)"
assert np.isfinite(np.asarray(traj.states[_k])).all(), f"frame {_k} carries non-finite values"
# Surface tension SHRINKS the surface, so a sudden LENGTHENING means the alpha filter bit a wedge out of
# it -- and a notch that reaches the surface adds no boundary loop, so the loop count above cannot see it.
# The filter is held steady by hysteresis (an existing cell survives to a wider threshold), which removes
# most of that: measured here, 2 steps of 300 still move the perimeter by more than 5 % (worst +13 %),
# because a Delaunay edge flip creates triangles that are NEW and so face `alpha` alone. Hence a loose
# catastrophe bound, plus the check that matters for the figure: the last frame must be representative of
# its neighbours rather than a bitten outlier.
_per = np.array([perimeter(k) for k in range(last + 1)])
_jump = float(np.max(_per[1:] / _per[:-1] - 1.0))
_median_tail = float(np.median(_per[-20:]))
assert _per[-1] < _per[0], "the free surface grew: surface tension does not do that"
assert abs(_per[-1] / _median_tail - 1.0) < 0.03, f"the last frame is an outlier ({_per[-1]:.4f} vs {_median_tail:.4f})"
assert _jump < 0.20, f"the alpha filter bit deep into the free surface (perimeter jumped {100 * _jump:.1f} % in a step)"
print(
f"perimeter {_per[0]:.4f} -> {_per[-1]:.4f}; worst one-step change {100 * _jump:+.2f} %; "
f"last frame {100 * (_per[-1] / _median_tail - 1):+.2f} % from its neighbours' median"
)
References. S. R. Idelsohn, E. Oñate & F. Del Pin, Int. J. Numer. Methods Eng. 61 (2004) 964–989 (PFEM). H. Edelsbrunner & E. P. Mücke, ACM Trans. Graph. 13 (1994) 43–72 (alpha shapes). H. Lamb, Hydrodynamics, 6th ed. (1932), §275 (drop oscillation). E. Bänsch, Numer. Math. 88 (2001) 203–235 (free-surface finite elements).