Example: a droplet next to a bead in a slit

A bead of radius \(R\) sits between two parallel walls \(z = \pm 0.5\), almost touching them (\(R = 0.48\): a gap of 0.02 at each wall). Beads repeat along \(x\) with spacing \(\ell\) (centre to centre), and a droplet wets the bead and both walls, reaching across to its neighbour’s droplet. By symmetry only a quarter cell is modelled: mirror planes \(x = 0\) and \(y = 0\) through the bead’s centre, and \(x = \ell/2\) halfway to the next bead, where the droplet meets its mirror image.

This builds on the liquid bridge: several walls at once, mirror planes that cut the problem down to a quarter, measured contact angles on each solid, and a liquid region closed face by face for a volume mesh.

Parameters:

• ell: the bead spacing, centre to centre (\(\ell \ge 1\); the beads would touch at \(\ell = 2R\));

• theta_bead, theta_wall: contact angles (degrees, measured inside the liquid) on the bead and on the walls;

• volume: the liquid volume in the quarter cell;

• radius: the bead radius (below 0.5; see Notes for why it must not touch).

Surface tension 1; lengths in units of the wall spacing.

import numpy as np
import matplotlib.pyplot as plt
import pyvista as pv
import pysurfaceevolver as pyse
from pysurfaceevolver import constraints as C

pv.set_jupyter_backend("static")    # "trame" in Jupyter for interactive views
pv.global_theme.window_size = [700, 450]

The model

The liquid–gas surface is the only meshed part. Its edges end on five constraints: the walls, the bead, and the mirror planes \(x = 0\) and \(x = \ell/2\) (the plane \(y = 0\) cuts through the liquid; the surface does not reach it).

The walls and the bead are not meshed either. Their wetted areas (energy \(-\cos\theta\) per area, Young’s law) and their part of the liquid’s volume enter as integrals along the contact lines, which pyse.constraints writes:

• Mirror planes (C.mirror) carry no energy (the surface meets them at 90° by itself) and, being vertical, add nothing to Evolver’s volume \(\int z\, n_z\, dA\).

• Walls (C.plane): the wetted area is \(\oint t\, ds\) in coordinates \(s = x\), \(t = \pm y\) on the wall. That integrand vanishes along \(y = 0\) and along the mirrors \(x = 0, \ell/2\), which close the wetted region.

• Bead (C.sphere): areas on a sphere are integrals over the azimuth about the \(z\) axis, and both mirror planes are meridians (no contribution). The bead’s poles lie in the wetted region, where the azimuth is undefined; the builder is told so (wet_poles) and the quarter (span), and knows what the integrals leave out: a volume (\(-\pi R^3/3\), the quarter bead), which goes into the body’s volconst since the body touches the bead, and the energy constant \(-\cos\theta_\text{bead}\,\pi R^2\), which ev.energy includes (see Energies).

The contact angles are given as the datafile parameters theta_wall and theta_bead, so they can be changed on a loaded surface.

The starting shape is a flat sheet of liquid at \(y = w\), from \(x = 0\) to \(\ell/2\) and wall to wall, with the bead cut out of it (\(w = 0.4\); relaxing meets the target volume). Its volume, areas and energy are known exactly, which checks all the signs.

def constraints(ell, radius):
    # normals point into the liquid; contact angles are datafile parameters
    return {1: C.plane((0, 0, 1), -0.5, contact_angle="theta_wall"),           # bottom wall
            2: C.plane((0, 0, -1), point=(0, 0, 0.5), contact_angle="theta_wall"),
            3: C.sphere((0, 0, 0), radius, contact_angle="theta_bead",
                        wet_poles=("north", "south"), span=np.pi/2),          # a quarter bead
            4: C.mirror("x"), 5: C.mirror("x", ell/2)}

def sheet_volume(ell, w, radius):
    # the liquid behind a flat sheet at y = w: a box minus the bead's part of it
    return ell/2*w - np.pi*(radius**2*w - w**3/3)/2

def load_slit(ell, theta_bead, theta_wall, volume, radius=0.48, w=0.4, ns=6, nt=8):
    rho = np.sqrt(radius**2 - w**2)        # the bead's cut in the plane y = w
    seg, arc = 0.5 - rho, np.pi*rho
    total = 2*seg + arc
    def left(t):
        # the sheet's edge at x ~ 0, by arc length: the x = 0 mirror below the
        # bead, around the bead, the mirror above it
        s = t*total
        if s <= seg:
            return 0.0, -0.5 + s, 4
        if s <= seg + arc:
            a = -np.pi/2 + (s - seg)/rho
            return rho*np.cos(a), rho*np.sin(a), 3
        return 0.0, rho + s - seg - arc, 4
    corners = [0, seg/total, (seg + arc)/total, 1]
    ts = np.unique(np.concatenate([np.linspace(corners[i], corners[i + 1], n + 1)
                                   for i, n in enumerate((2, nt, 2))]))
    ss = np.linspace(0, 1, ns + 1)
    vertices, on = [], []
    for t in ts:
        xl, zl, k = left(t)
        for s in ss:      # a ruled sheet from the left edge to the x = ell/2 mirror
            vertices.append(((1 - s)*xl + s*ell/2, w, (1 - s)*zl + s*(-0.5 + t)))
            c = set()
            if s == 0:
                c.add(k)
                if np.isclose(abs(zl), rho) and np.isclose(xl, 0):
                    c |= {3, 4}                   # where the bead meets the mirror
            if s == 1: c.add(5)
            if t == 0: c.add(1)
            if t == 1: c.add(2)
            on.append(c)
    n = len(ss)
    faces = [[i*n + j, (i + 1)*n + j, (i + 1)*n + j + 1, i*n + j + 1]   # normals +y
             for i in range(len(ts) - 1) for j in range(ns)]
    ev = pyse.Evolver()
    ev.load_arrays(np.array(vertices), faces, constraints=constraints(ell, radius),
                   vertex_constraints=on,            # the constraints of each vertex
                   bodies=[pyse.Body(volume=volume)],
                   parameters={"radius": radius, "ell": ell, "theta_bead": theta_bead,
                               "theta_wall": theta_wall})
    return ev

The check: a fine flat sheet (no relaxation), with no volume target, against the exact volume and energy. ev.energy includes the constant \(-\cos\theta_\text{bead}\,\pi R^2\) that Evolver’s ev.total_energy leaves out.

ell, w, R, tb, tw = 2.0, 0.3, 0.48, 40.0, 60.0
ev = load_slit(ell, tb, tw, None, R, w=w, ns=32, nt=32)
ev.recalc()
rho2 = R**2 - w**2
exact_energy = (ell/2 - np.pi*rho2/2) - 2*np.cos(np.radians(tw))*ell/2*w \
               - np.cos(np.radians(tb))*np.pi*R*w
print(f"volume {ev.bodies().volume[0]:.6f}, exact {sheet_volume(ell, w, R):.6f}")
print(f"energy {ev.energy:.6f}, exact {exact_energy:.6f}")

volume 0.205500, exact 0.205564
energy 0.132993, exact 0.132910

The small differences are the polygonal contact line on the bead; they shrink fourfold when its points double.

Relaxing

ev.relax(levels=3) relaxes each refinement level until the energy settles, equiangulating and averaging vertices every ten steps so facets at the contact lines don’t collapse; Newton steps finish it.

def slit_droplet(ell=2.0, theta_bead=40.0, theta_wall=60.0, volume=0.3, radius=0.48):
    ev = load_slit(ell, theta_bead, theta_wall, volume, radius)
    ev.relax(levels=3)
    return ev

ell, theta_bead, theta_wall, volume, radius = 2.0, 40.0, 60.0, 0.3, 0.48

ev = load_slit(ell, theta_bead, theta_wall, volume, radius)
view = ev.live_view(color="lightskyblue", walls=True)     # the walls and the bead
view.plotter.camera_position = [(2.2, 2.6, 1.6), (0.4, 0.2, 0), (0, 0, 1)]
ev.relax(levels=3, callback=view.update, every=10);

The quarter cell mirrored into the full picture: the droplet spans the slit and joins its neighbour’s halfway between the beads.

def show_slit(ev):
    # the quarter cell mirrored to two neighbouring cells, with both beads
    ell, radius = ev.parameters["ell"], ev.parameters["radius"]
    solids = constraints(ell, radius)
    ev.plot(mirror=["y", "x", ("x", ell/2)], color="lightskyblue", smooth_shading=True,
            opacity=0.85, show_edges=False,
            walls=[solids[1], solids[2], solids[3], C.sphere((ell, 0, 0), radius)])

show_slit(ev)

Checking the result

The measured contact angles should match the parameters on all three solids. The pressure is the volume constraint’s Lagrange multiplier (negative: the meniscus is concave and pulls the walls together). ev.eigen_counts().negative counts unstable directions (0: stable). The last line is a pinch-off check: the free surface must stay clear of the \(y = 0\) mirror, or the droplet would really meet its mirror image and split.

# the angle between the surface's normal out of the liquid and the wall's out of
# the solid, at each contact-line vertex (corners, on two walls, left out)
for k, name in ((1, "bottom wall"), (2, "top wall"), (3, "bead")):
    angles = ev.contact_angles(k)
    print(f"{name:12s} {angles.mean():6.2f}° ± {angles.std():.2f}°")
print(f"pressure {ev.body(1).pressure:.4f}, energy {ev.energy:.5f}, "
      f"unstable directions {ev.eigen_counts().negative}")
print(f"closest approach to the y = 0 mirror: {ev.mesh().vertices[:, 1].min():.3f}")

bottom wall   60.36° ± 0.05°
top wall      60.36° ± 0.05°
bead          40.16° ± 0.23°
pressure -1.0424, energy -0.02675, unstable directions 0
closest approach to the y = 0 mirror: 0.359

Varying the wall contact angle

The less the liquid wets the walls, the less concave its meniscus between them and the less negative its pressure. Cross-sections of the free surface in the mid-plane \(z = 0\) (from above) and on the plane \(x = \ell/2\) (where the droplet meets its neighbour):

fig, (top, side) = plt.subplots(1, 2, figsize=(10, 4), width_ratios=[2.2, 1])
s = np.linspace(0, np.pi/2, 50)
top.fill_between(radius*np.cos(s), 0, radius*np.sin(s), color="sandybrown", alpha=0.5)
pressures = []
for tw, color in zip((40, 50, 60, 70, 80), plt.cm.viridis(np.linspace(0, 0.9, 5))):
    ev = slit_droplet(ell, theta_bead, tw, volume, radius)
    pressures.append(ev.body(1).pressure)
    m = ev.mesh()
    mid = m.slice((0, 0, 1))[0]                       # the mid-plane z = 0
    top.plot(mid[:, 0], mid[:, 1], color=color, label=f"θ_wall = {tw}°")
    edge = m.slice((1, 0, 0), (ell/2, 0, 0))[0]       # the mirror x = ell/2
    side.plot(edge[:, 1], edge[:, 2], color=color)
top.set_xlabel("x"); top.set_ylabel("y"); top.set_aspect("equal")
top.set_title("mid-plane z = 0"); top.legend(fontsize="small")
side.set_xlabel("y"); side.set_ylabel("z"); side.set_aspect("equal"); side.set_xlim(0, 0.6)
side.set_title("plane x = ℓ/2")
for z in (-0.5, 0.5):
    side.axhline(z, color="gray")
plt.tight_layout()
print("pressures:", np.round(pressures, 3))

pressures: [-1.65  -1.367 -1.042 -0.685 -0.305]

Exporting for finite elements

The liquid region is bounded by the free surface and by pieces of the walls, the bead and the mirror planes. Its free surface ends on all of them, so its single boundary loop runs along six arcs: the bottom wall, the \(x = 0\) mirror in the gap below the bead, the bead, the gap above it, the top wall, and the \(x = \ell/2\) mirror. Each of the seven other faces of the liquid region is bounded by some of those arcs and by known edges: straight lines where two planes meet, circles on the bead.

So the region is closed face by face. Each face is meshed with Gmsh as a 2D polygon, its boundary points fixed: the flat faces in their plane, the bead’s after a stereographic projection (from the point of the bead opposite the liquid), its new points mapped back onto the sphere. Faces share their edges point for point, so the result is watertight, and volume_mesh fills it with tetrahedra.

import gmsh

def boundary_loop(facets):
    # the vertex rows along a surface's (single) boundary loop, in order
    edges = np.sort(np.vstack([facets[:, [0, 1]], facets[:, [1, 2]], facets[:, [2, 0]]]), axis=1)
    unique, count = np.unique(edges, axis=0, return_counts=True)
    neighbours = {}
    for p, q in unique[count == 1].tolist():
        neighbours.setdefault(p, []).append(q)
        neighbours.setdefault(q, []).append(p)
    loop = [next(iter(neighbours))]
    while True:
        nxt = [v for v in neighbours[loop[-1]] if len(loop) < 2 or v != loop[-2]][0]
        if nxt == loop[0]:
            return np.array(loop)
        loop.append(nxt)

def mesh_polygon(uv, size):
    # triangles filling a 2D polygon with its points kept (no new boundary
    # points); returns the new interior points and triangles indexing
    # [polygon points, new points]
    gmsh.model.add("polygon")
    geo = gmsh.model.geo
    corners = [geo.addPoint(u, v, 0, size) for u, v in uv]
    sides = [geo.addLine(corners[i], corners[(i + 1) % len(corners)]) for i in range(len(corners))]
    for side in sides:
        geo.mesh.setTransfiniteCurve(side, 2)
    geo.addPlaneSurface([geo.addCurveLoop(sides)])
    geo.synchronize()
    gmsh.model.mesh.generate(2)
    index = {int(gmsh.model.mesh.getNodes(0, c)[0][0]): i for i, c in enumerate(corners)}
    tags, xyz, _ = gmsh.model.mesh.getNodes(2, -1)
    for t in tags.tolist():
        index[int(t)] = len(index)
    triangles = np.vectorize(index.get)(gmsh.model.mesh.getElementsByType(2)[1].reshape(-1, 3))
    gmsh.model.remove()
    return xyz.reshape(-1, 3)[:, :2], triangles

def close_slit(ev):
    ell, R = ev.parameters["ell"], ev.parameters["radius"]
    m = ev.mesh()
    points = [m.vertices]
    def add(xyz):        # new points; their rows
        points.append(np.atleast_2d(np.asarray(xyz, float)))
        start = sum(len(p) for p in points[:-1])
        return list(range(start, start + len(points[-1])))

    # the rim, cut into its six arcs at the corners (vertices on two constraints)
    rim = boundary_loop(m.facets)
    on = np.array([ev.on_constraint(k) for k in range(1, 6)]).T
    rim = np.roll(rim, -np.flatnonzero(on[rim].sum(1) == 2)[0])
    cuts = list(np.flatnonzero(on[rim].sum(1) == 2)) + [len(rim)]
    arcs = {}
    for a, b in zip(cuts[:-1], cuts[1:]):
        arc = list(rim[a:b]) + [rim[b % len(rim)]]
        k = int(np.flatnonzero(on[arc[1]])[0]) + 1      # its constraint
        name = {1: "bottom", 2: "top", 3: "bead", 5: "far"}.get(k)
        if k == 4:
            name = "below" if m.vertices[arc, 2].mean() < 0 else "above"
        arcs[name] = arc
    def corner(a, b):
        return (set(arcs[a]) & set(arcs[b])).pop()
    def run(name, start):    # an arc from one of its ends
        return arcs[name] if arcs[name][0] == start else arcs[name][::-1]

    # edges that are not on the rim, with points about as far apart as the rim's
    h = np.linalg.norm(np.diff(m.vertices[rim], axis=0), axis=1).mean()
    def line(p, q):
        P = np.vstack(points)
        n = max(1, int(np.ceil(np.linalg.norm(P[q] - P[p])/h)))
        t = np.linspace(0, 1, n + 1)[1:-1, None]
        return [p] + (add(P[p] + t*(P[q] - P[p])) if len(t) else []) + [q]
    def bead_arc(p, q):      # a great-circle arc on the bead
        P = np.vstack(points)
        a, b = P[p], P[q]
        angle = np.arccos(np.clip(a @ b/R**2, -1, 1))
        n = max(1, int(np.ceil(R*angle/h)))
        t = np.linspace(0, 1, n + 1)[1:-1]
        inner = [(np.sin((1 - s)*angle)*a + np.sin(s*angle)*b)/np.sin(angle) for s in t]
        return [p] + (add(inner) if inner else []) + [q]
    b0, b1, t0, t1, pb, pt, eq = add([(0, 0, -0.5), (ell/2, 0, -0.5), (0, 0, 0.5), (ell/2, 0, 0.5),
                                      (0, 0, -R), (0, 0, R), (R, 0, 0)])
    c_bb, c_bf = corner("bottom", "below"), corner("bottom", "far")
    c_tf, c_ta = corner("top", "far"), corner("top", "above")
    c_ba, c_bb2 = corner("bead", "above"), corner("bead", "below")
    floor, roof, back = line(b0, b1), line(t1, t0), line(b1, t1)
    axis_b, axis_t = line(pb, b0), line(t0, pt)
    meridian = bead_arc(pt, eq)[:-1] + bead_arc(eq, pb)        # bead in y = 0
    gap_b, gap_t = bead_arc(c_bb2, pb), bead_arc(pt, c_ba)       # bead in x = 0
    wall_b0, wall_b1 = line(c_bb, b0), line(b1, c_bf)
    wall_t0, wall_t1 = line(c_ta, t0), line(t1, c_tf)

    def loop(*chains):       # closed polygon from chains that join end to start
        return [i for c in chains for i in c[:-1]]
    flat = lambda d: (lambda p: p[:, d])
    sphere_to_plane = lambda p: R*p[:, [0, 2]]/(p[:, [1]] + R)
    def plane_to_sphere(uv):
        s = 2*R**2/((uv**2).sum(1) + R**2)
        return np.column_stack([s*uv[:, 0], -R + s*R, s*uv[:, 1]])
    faces = {   # name: (polygon, to 2D, back to 3D, direction out of the liquid)
        "bottom wall": (loop(run("bottom", c_bf), wall_b0, floor, wall_b1), flat([0, 1]), None, (0, 0, -1)),
        "top wall": (loop(run("top", c_ta), wall_t1[::-1], roof, wall_t0[::-1]), flat([0, 1]), None, (0, 0, 1)),
        "mirror x = ell/2": (loop(run("far", c_tf), wall_b1[::-1], back, wall_t1), flat([1, 2]), None, (1, 0, 0)),
        "mirror x = 0, below": (loop(run("below", c_bb), gap_b, axis_b, wall_b0[::-1]), flat([1, 2]), None, (-1, 0, 0)),
        "mirror x = 0, above": (loop(run("above", c_ba), wall_t0, axis_t, gap_t), flat([1, 2]), None, (-1, 0, 0)),
        "mirror y = 0": (loop(floor, back, roof, axis_t, meridian, axis_b), flat([0, 2]), None, (0, -1, 0)),
        "bead": (loop(run("bead", c_bb2), gap_t[::-1], meridian, gap_b[::-1]),
                 sphere_to_plane, plane_to_sphere, None),
    }
    out = m.facet_bodies[:, 0] == 1
    cells = [np.where(out[:, None], m.facets, m.facets[:, ::-1])]     # free surface
    tags, names = [np.zeros(len(m.facets), int)], {0: "free surface"}
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    gmsh.option.setNumber("Mesh.MeshSizeMax", h)
    for k, (name, (polygon, to_2d, to_3d, outward)) in enumerate(faces.items(), start=1):
        P = np.vstack(points)
        uv, triangles = mesh_polygon(to_2d(P[polygon]), h)
        if to_3d is None:        # a flat face: put back its constant coordinate
            d = int(np.flatnonzero(np.abs(outward))[0])
            xyz = np.insert(uv, d, P[polygon[0], d], axis=1)
        else:
            xyz = to_3d(uv)
        rows = np.array(polygon + (add(xyz) if len(xyz) else []))
        T = rows[triangles]
        P = np.vstack(points)
        normal = np.cross(P[T[:, 1]] - P[T[:, 0]], P[T[:, 2]] - P[T[:, 0]])
        want = -P[T].mean(1) if outward is None else np.asarray(outward, float)   # bead: inward
        flip = (normal*want).sum(1) < 0
        T[flip] = T[flip][:, ::-1]
        cells.append(T)
        tags.append(np.full(len(T), k))
        names[k] = name
    gmsh.finalize()
    T = np.vstack(cells)
    surface = pyse.BodySurface(1, np.vstack(points), T, "triangle", np.zeros(len(T), int),
                              pyse.is_watertight(T), 0, cap_ids=np.concatenate(tags))
    return surface, names

A droplet refined twice (fewer elements than above), closed and checked: the enclosed volume against the target, and the energy recomputed from the face areas (\(A_\text{free} - \cos\theta_\text{wall} A_\text{walls} - \cos\theta_\text{bead} A_\text{bead}\)) against Evolver’s.

ev = load_slit(ell, theta_bead, theta_wall, volume, radius)
ev.relax(levels=2)
surface, names = close_slit(ev)

area = surface.areas()                     # per face: index k is names[k]
by_face = {name: area[k] for k, name in names.items()}
energy = (by_face["free surface"]
          - np.cos(np.radians(theta_wall))*(by_face["bottom wall"] + by_face["top wall"])
          - np.cos(np.radians(theta_bead))*by_face["bead"])
print("watertight:", surface.watertight, f"| {len(surface.cells)} triangles")
print(f"enclosed volume {surface.volume:.5f} (target {volume})")
# the energy is a small difference of larger areas: the faces (flat caps on the walls)
# and Evolver's line integrals differ at second order, 4e-4 here and 1e-4 one
# refinement further
print(f"energy from the faces {energy:.5f}, Evolver's {ev.energy:.5f}")
for name, A in by_face.items():
    print(f"  {name:20s} area {A:.4f}")

watertight: True | 11904 triangles
enclosed volume 0.30016 (target 0.3)
energy from the faces -0.02631, Evolver's -0.02672
  free surface         area 0.9267
  bottom wall          area 0.4994
  top wall             area 0.4994
  mirror x = ell/2     area 0.4020
  mirror x = 0, below  area 0.0469
  mirror x = 0, above  area 0.0469
  mirror y = 0         area 0.6383
  bead                 area 0.5922

The volume mesh, with a physical group per face (named after it) and the body’s for the tetrahedra; .msh for dolfinx.io.gmshio, XDMF for others:

import os, tempfile
import meshio

folder = tempfile.mkdtemp()
mesh = surface.volume_mesh(size=0.05, names=names, path=os.path.join(folder, "slit.msh"))
tets, P = mesh.cells_dict["tetra"], mesh.points
v = np.einsum("ij,ij->i", P[tets[:, 1]] - P[tets[:, 0]],
              np.cross(P[tets[:, 2]] - P[tets[:, 0]], P[tets[:, 3]] - P[tets[:, 0]]))
print(f"{len(tets)} tetrahedra, volume {np.abs(v).sum()/6:.5f}")
groups = mesh.cell_data_dict["gmsh:physical"]
for g, n in zip(*np.unique(groups["triangle"], return_counts=True)):
    print(f"  group {g}: {names[g - 1]:20s} {n} boundary triangles")
meshio.write(os.path.join(folder, "slit.xdmf"),
             meshio.Mesh(P, [("tetra", tets)], cell_data={"region": [groups["tetra"]]}))
meshio.write(os.path.join(folder, "slit_boundary.xdmf"),
             meshio.Mesh(P, [("triangle", mesh.cells_dict["triangle"])],
                         cell_data={"boundary": [groups["triangle"]]}))

grid = pv.from_meshio(meshio.Mesh(P, [("tetra", tets)]))
plotter = pv.Plotter()
plotter.add_mesh(grid.clip(normal="z", invert=True, crinkle=True), color="lightskyblue",
                 show_edges=True)
plotter.camera_position = [(1.4, 1.6, 1.3), (0.5, 0.2, -0.25), (0, 0, 1)]
plotter.show()

66016 tetrahedra, volume 0.30016
  group 1: free surface         4608 boundary triangles
  group 2: bottom wall          1220 boundary triangles
  group 3: top wall             1218 boundary triangles
  group 4: mirror x = ell/2     1391 boundary triangles
  group 5: mirror x = 0, below  131 boundary triangles
  group 6: mirror x = 0, above  131 boundary triangles
  group 7: mirror y = 0         1677 boundary triangles
  group 8: bead                 1528 boundary triangles

Notes and limits

• A touching bead. With \(R = 0.5\) the bead touches the walls at single points, and a wetting liquid is drawn into the wedge around them, which closes to zero thickness; the mesh can’t follow it there. A small gap (\(R = 0.48\)) ends the wedge in a thin but finite meniscus.

• Pinch-off. Between walls a distance 1 apart, the meniscus sags inward by about \(0.5\,(1 - \sin\theta_\text{wall})/\cos\theta_\text{wall}\) (0.13 at 60°, 0.23 at 40°). With too little liquid for that, the free surface reaches the \(y = 0\) mirror: physically the droplet splits into a top and a bottom part, which the mesh can’t represent. Check the closest approach to the mirror; give more volume (or a larger wall angle).

• Non-wetting beads (\(\theta_\text{bead} > 90°\)) push the liquid away until it pinches in the same way.

• Stability. A relaxed droplet is an equilibrium; the relax result’s stable (or ev.eigen_counts().negative) tells whether it is a stable one.

• Export. The liquid’s free surface ends on several solids and mirrors at once, so body_surfaces(cap=True) (one cap per opening, on one constraint) cannot close it; close_slit above closes it face by face for this geometry.

Next

Draining the slit takes the same geometry through a long sequence of equilibria, removing liquid step by step until the liquid snaps.