Example: draining liquid from a bead chain in a slit

This continues the slit droplet example, and is the longest run in these pages: a quasi-static sequence of equilibria with events along the way, the mesh remeshed at every step (see Mesh resolution). Beads of radius \(R = 0.48\) sit midway between two walls \(z = \pm 0.5\), in a periodic chain along \(x\) with spacing \(\ell\) (centre to centre; the beads touch at \(\ell = 2R = 0.96\)). We start with the beads fully immersed and remove liquid slowly, so the liquid passes through a sequence of equilibria. Contact angles: \(\theta = 40°\) on the beads, \(60°\) on the walls. Lengths are in units of the wall spacing, and the surface tension is 1, so the pressure jump across the free surface is \(\Delta p = 2H\), with \(H\) the mean curvature.

What happens as the liquid drains:

1. Immersed. The liquid fills the slit around the beads, bounded by a meniscus that is straight along \(x\) and a circular arc across the gap. Its curvature is known exactly: \(\Delta p = -2\cos\theta_\text{wall} = -1\), independent of the volume.

2. The bead emerges. When the meniscus reaches the beads, a dry spot opens on each bead. That happens exactly when the meniscus reaches \(y = R\) on the mid-plane. The pressure jumps.

3. The band. The liquid forms a band along the chain on each side. It thins, and its pressure falls steeply.

4. The snap. The two bands meet across the chain at the plane \(y = 0\) and snap there. After that, the liquid around each bead no longer connects to the next bead’s along the band, and we stop.

By symmetry we compute one eighth of a unit cell: \(0 \le x \le \ell/2\), \(y \ge 0\), \(0 \le z \le 1/2\), with the bead at the origin and mirror planes \(x = 0\), \(x = \ell/2\), \(y = 0\) and \(z = 0\). The code lives in drainage_helpers.py next to this notebook:

• Cell(ell) defines the eighth cell. The bead and wall contact energies and the liquid volume under them are line integrals along the contact lines (Evolver constraint integrals). The bead’s own volume is subtracted with a volconst, from the bead’s constraint (given explicitly: make_datafile adds the volconst of the walls the surface touches, and the immersed meniscus doesn’t touch the bead).

• Cell.meniscus(y0) is the immersed start, and Cell.sheet(w) is a flat sheet with a dry spot on the bead, the restart after the bead emerges.

• drain(cell, v_start, v_end) steps the volume down. After each step it relaxes the surface with gradient steps, equiangulation and vertex averaging, then a few Newton steps. It watches for the two events above. Near an event it takes 0.5% steps.

import time
import numpy as np
import matplotlib.pyplot as plt
import pyvista as pv
import pysurfaceevolver as pyse
import drainage_helpers as dh

pv.set_jupyter_backend("static")
pv.global_theme.window_size = [700, 450]

The immersed start

The meniscus sits at \(y = y_0\) on the mid-plane. Its volume and energy are known in closed form (Cell.meniscus_exact), so the computed state can be checked directly. The energy includes the bead’s wetting energy, counted from a dry bead.

cell = dh.Cell(1.1)
y0 = 0.55
volume, energy = cell.meniscus_exact(y0)
ev = cell.load(*cell.meniscus(y0), volume)
cell.relax(ev)
for _ in range(2):
    ev.refine()
    cell.relax(ev)
print(f"volume     {ev.body(1).volume:.6f}   exact {volume:.6f}")
print(f"energy     {ev.energy:.6f}   exact {energy:.6f}")
print(f"pressure  {ev.body(1).pressure:+.6f}   exact {-2*np.cos(np.radians(cell.theta_wall)):+.6f}")

volume     0.105276   exact 0.105276
energy     -0.177349   exact -0.177354
pressure  -1.000041   exact -1.000000

Draining at \(\ell = 1.1\)

The beads are \(0.14\) apart. We drain from the immersed state to a liquid volume of \(0.02\) per eighth cell (\(0.16\) per bead) on a coarse mesh, which takes a few minutes. Below that the coarse mesh no longer resolves the band near the beads, and the stored run at twice the resolution takes over.

On the band the mesh is remeshed by edge length at every step. As the dry spot grows, the contact line on the bead stretches, and the surface beside it is squeezed into the corner where the contact line meets the mid-plane \(z = 0\). Without remeshing, thin triangles pile up there and the pressure turns noisy.

cell = dh.Cell(1.1)
cell.band_step = 0.03
start = time.perf_counter()
states, events = dh.drain(cell, cell.meniscus_exact(0.55)[0], 0.02, levels=2)
print(f"{len(states)} states in {time.perf_counter() - start:.0f} s")
for e in events:
    print(e)

52 states in 70 s
{'kind': 'bead emerges', 'volume': 0.08602600789914044, 'before': -1.0000528011717804, 'after': -0.7292784201370199}

def show_state(state, ell, radius=0.48):
    # the eighth cell mirrored to two neighbouring unit cells
    faces = np.hstack([np.full((len(state["facets"]), 1), 3), state["facets"]]).ravel()
    eighth = pv.PolyData(state["vertices"], faces)
    surface = eighth.merge(eighth.reflect((0, 0, 1), point=(0, 0, 0)))
    surface = surface.merge(surface.reflect((0, 1, 0), point=(0, 0, 0)))
    surface = surface.merge(surface.reflect((1, 0, 0), point=(0, 0, 0)))
    surface = surface.merge(surface.reflect((1, 0, 0), point=(ell/2, 0, 0)))
    plotter = pv.Plotter()
    for cx in (0, ell):
        plotter.add_mesh(pv.Sphere(radius=radius, center=(cx, 0, 0), theta_resolution=60,
                                   phi_resolution=60), color="sandybrown")
    for z in (-0.5, 0.5):
        plotter.add_mesh(pv.Plane(center=(ell/2, 0, z), i_size=ell + 1.6, j_size=1.9),
                         color="gray", opacity=0.12)
    plotter.add_mesh(surface, color="lightskyblue", smooth_shading=True, opacity=0.9)
    plotter.camera_position = [(ell/2 + 1.0, -3.6, 2.4), (ell/2, 0, -0.05), (0, 0, 1)]
    plotter.show()

show_state(states[-1], cell.ell)

The band at the end of the live run, mirrored to two unit cells: the dry spots on the beads face the reader, and the band runs along the chain on either side of the beads.

Pressure and curvature

The pressure stays at \(-1\) while the beads are immersed and jumps when they emerge. Then it falls steeply as the band thins. The full run below (drainage_runs.npz, made with the same method on a mesh refined once more) continues to the snap. In the overlap the two agree to 0.03% typically. Just after the emergence they differ by up to 1%, because the stored run found the emergence by stepping (about 1% from the exact volume) instead of computing it.

runs = np.load("_static/drainage_runs.npz")
CELLS = 8      # eighth cells per bead

def pressure_axes(ax):
    ax.axhline(0, color="black", lw=0.5)
    ax.set_xlabel("liquid volume per bead")
    ax.set_ylabel(r"pressure $\Delta p$")
    ax.secondary_yaxis("right", functions=(lambda p: p/2, lambda h: 2*h)).set_ylabel(
        r"mean curvature $H = \Delta p/2$")
    ax.invert_xaxis()

fig, ax = plt.subplots(figsize=(8, 4.5))
v, p, stage = runs["l1.1_volume"], runs["l1.1_pressure"], runs["l1.1_stage"]
for s, label in ((0, "immersed"), (1, "band")):
    ax.plot(CELLS*v[stage == s], p[stage == s], "-", lw=2, color=f"C{s}",
            label=f"{label} (stored, finer mesh)")
live = np.array([(st["volume"], st["pressure"]) for st in states])
ax.plot(CELLS*live[:, 0], live[:, 1], "k.", ms=3, label="this run")
for x, name in zip(runs["l1.1_events"], ("bead emerges", "snap")):
    ax.axvline(CELLS*x, color="gray", ls=":", lw=1)
    ax.annotate(name, (CELLS*x, 0.3), rotation=90, ha="right", va="bottom", fontsize=9)
pressure_axes(ax)
ax.set_title(r"$\ell = 1.1$")
ax.legend(loc="lower left")
plt.show()

Three spacings

The same drainage, stored, for beads \(0.04\) apart (\(\ell = 1\)), \(0.14\) apart (\(\ell = 1.1\)) and well apart (\(\ell = 2\)). Closer beads hold the band longer: it reaches lower volume and much lower pressure before it snaps.

fig, ax = plt.subplots(figsize=(8, 4.5))
for i, (key, ell) in enumerate((("l1", 1.0), ("l1.1", 1.1), ("l2", 2.0))):
    v, p = runs[f"{key}_volume"], runs[f"{key}_pressure"]
    ax.plot(CELLS*v, p, "-", lw=2, color=f"C{i}", label=fr"$\ell = {ell:g}$")
    snap = CELLS*runs[f"{key}_events"][1]
    ax.plot([snap], [p[-1]], "o", color=f"C{i}")
    emerge = CELLS*dh.Cell(ell).meniscus_exact(0.48)[0]
    print(f"ell = {ell:g}: bead emerges at {emerge:.3f} (exact), "
          f"snap at {snap:.4f} per bead (Δp = {p[-1]:+.2f} just before)")
pressure_axes(ax)
ax.legend(loc="lower left")
plt.show()

ell = 1: bead emerges at 0.584 (exact), snap at 0.0435 per bead (Δp = -7.51 just before)
ell = 1.1: bead emerges at 0.688 (exact), snap at 0.0566 per bead (Δp = -5.95 just before)
ell = 2: bead emerges at 1.630 (exact), snap at 0.3298 per bead (Δp = -1.82 just before)

Movies

Each frame is an equilibrium. Time in the movies is the volume removed, and they end at the snap. The snap volume was found by stepping the volume in 0.5% steps near it, stopping when the band comes within \(5\cdot10^{-4}\) of the plane \(y = 0\). At \(\ell = 2\) a mesh refined once more moves it by 0.05%.

\(\ell = 2\):

\(\ell = 1\) (beads \(0.04\) apart):

Limits of the method:

• Quasi-static. Every state is an equilibrium at the given volume, so the liquid moves between configurations only through the sequence of equilibria. Each event is instantaneous: the state jumps to the new branch at that volume.

• Symmetric. The eighth cell imposes the mirror symmetries, so modes that break them (for example liquid gathering at one bead) are not followed.

• After the snap. The liquid left sits around each bead, for example in rings at the bead–wall gaps. Following it needs other starting shapes, and it isn’t computed here.

Related work

• C. M. Cejas, L. A. Hough, C. Frétigny and R. Dreyfus, Effect of geometry on the dewetting of granular chains by evaporation, Soft Matter 14 (2018) 6994–7002, doi:10.1039/C8SM01179F. Chains of beads drying, from loose to close packing, in experiments and with Surface Evolver: above a critical packing angle the liquid network breaks up into separate bridges.

• J. L. Hilden and K. P. Trumble, Numerical analysis of capillarity in packed spheres: planar hexagonal-packed spheres, J. Colloid Interface Sci. 267 (2003) 463–474, doi:10.1016/S0021-9797(03)00683-0. The capillary pressure at which liquid drains through a layer of packed spheres.

Next

That ends the worked examples. Brakke’s samples run the sample datafiles that come with Evolver, each with his recipe and with relax(). The drainage band is one of two cases where relax()’s defaults alone fall short (where the defaults stop working); the stress cases are the other hard problems it is checked against.