{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": "# Example: draining liquid from a bead chain in a slit\n\nThis continues the [slit droplet example](slit_droplet.ipynb), and is the longest run\nin these pages: a quasi-static sequence of equilibria with events along the way, the\nmesh remeshed at every step (see [Mesh resolution](relaxing.md#mesh-resolution)). Beads of radius\n$R = 0.48$ sit midway between two walls $z = \\pm 0.5$, in a periodic chain along $x$ with\nspacing $\\ell$ (centre to centre; the beads touch at $\\ell = 2R = 0.96$). We start with the\nbeads fully immersed and remove liquid slowly, so the liquid passes through a sequence of\nequilibria. Contact angles: $\\theta = 40°$ on the beads, $60°$ on the walls. Lengths are in\nunits of the wall spacing, and the surface tension is 1, so the pressure jump across the\nfree surface is $\\Delta p = 2H$, with $H$ the mean curvature.\n\nWhat happens as the liquid drains:\n\n1. **Immersed.** The liquid fills the slit around the beads, bounded by a meniscus that is\n straight along $x$ and a circular arc across the gap. Its curvature is known exactly:\n $\\Delta p = -2\\cos\\theta_\\text{wall} = -1$, independent of the volume.\n2. **The bead emerges.** When the meniscus reaches the beads, a dry spot opens on each\n bead. That happens exactly when the meniscus reaches $y = R$ on the mid-plane. The\n pressure jumps.\n3. **The band.** The liquid forms a band along the chain on each side. It thins, and its\n pressure falls steeply.\n4. **The snap.** The two bands meet across the chain at the plane $y = 0$ and snap there.\n After that, the liquid around each bead no longer connects to the next bead's along the\n band, and we stop.\n\nBy symmetry we compute one eighth of a unit cell: $0 \\le x \\le \\ell/2$, $y \\ge 0$,\n$0 \\le z \\le 1/2$, with the bead at the origin and mirror planes $x = 0$, $x = \\ell/2$,\n$y = 0$ and $z = 0$. The code lives in\n[drainage_helpers.py](https://github.com/gautelinga/pySurfaceEvolver/blob/main/docs/drainage_helpers.py)\nnext to this notebook:\n\n* `Cell(ell)` defines the eighth cell. The bead and wall contact energies and the liquid\n volume under them are line integrals along the contact lines (Evolver constraint\n integrals). The bead's own volume is subtracted with a `volconst`, from the bead's\n constraint (given explicitly: `make_datafile` adds the `volconst` of the walls the\n surface touches, and the immersed meniscus doesn't touch the bead).\n* `Cell.meniscus(y0)` is the immersed start, and `Cell.sheet(w)` is a flat sheet with a dry\n spot on the bead, the restart after the bead emerges.\n* `drain(cell, v_start, v_end)` steps the volume down. After each step it relaxes the\n surface with gradient steps, equiangulation and vertex averaging, then a few Newton\n steps. It watches for the two events above. Near an event it takes 0.5% steps.",
"id": "cell-0"
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": "import time\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport pyvista as pv\nimport pysurfaceevolver as pyse\nimport drainage_helpers as dh\n\npv.set_jupyter_backend(\"static\")\npv.global_theme.window_size = [700, 450]",
"id": "cell-1"
},
{
"cell_type": "markdown",
"metadata": {},
"source": "## The immersed start\n\nThe meniscus sits at $y = y_0$ on the mid-plane. Its volume and energy are known in closed\nform (`Cell.meniscus_exact`), so the computed state can be checked directly. The\nenergy includes the bead's wetting energy, counted from a dry bead.",
"id": "cell-2"
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": "cell = dh.Cell(1.1)\ny0 = 0.55\nvolume, energy = cell.meniscus_exact(y0)\nev = cell.load(*cell.meniscus(y0), volume)\ncell.relax(ev)\nfor _ in range(2):\n ev.refine()\n cell.relax(ev)\nprint(f\"volume {ev.body(1).volume:.6f} exact {volume:.6f}\")\nprint(f\"energy {ev.energy:.6f} exact {energy:.6f}\")\nprint(f\"pressure {ev.body(1).pressure:+.6f} exact {-2*np.cos(np.radians(cell.theta_wall)):+.6f}\")",
"id": "cell-3"
},
{
"cell_type": "markdown",
"metadata": {},
"source": "## Draining at $\\ell = 1.1$\n\nThe beads are $0.14$ apart. We drain from the immersed state to a liquid volume of\n$0.02$ per eighth cell ($0.16$ per bead) on a coarse mesh, which takes a few minutes. Below\nthat the coarse mesh no longer resolves the band near the beads, and the stored run at\ntwice the resolution takes over.\n\nOn the band the mesh is remeshed by edge length at every step. As the dry spot grows, the\ncontact line on the bead stretches, and the surface beside it is squeezed into the corner\nwhere the contact line meets the mid-plane $z = 0$. Without remeshing, thin triangles pile\nup there and the pressure turns noisy.",
"id": "cell-4"
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": "cell = dh.Cell(1.1)\ncell.band_step = 0.03\nstart = time.perf_counter()\nstates, events = dh.drain(cell, cell.meniscus_exact(0.55)[0], 0.02, levels=2)\nprint(f\"{len(states)} states in {time.perf_counter() - start:.0f} s\")\nfor e in events:\n print(e)",
"id": "cell-5"
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": "def show_state(state, ell, radius=0.48):\n # the eighth cell mirrored to two neighbouring unit cells\n faces = np.hstack([np.full((len(state[\"facets\"]), 1), 3), state[\"facets\"]]).ravel()\n eighth = pv.PolyData(state[\"vertices\"], faces)\n surface = eighth.merge(eighth.reflect((0, 0, 1), point=(0, 0, 0)))\n surface = surface.merge(surface.reflect((0, 1, 0), point=(0, 0, 0)))\n surface = surface.merge(surface.reflect((1, 0, 0), point=(0, 0, 0)))\n surface = surface.merge(surface.reflect((1, 0, 0), point=(ell/2, 0, 0)))\n plotter = pv.Plotter()\n for cx in (0, ell):\n plotter.add_mesh(pv.Sphere(radius=radius, center=(cx, 0, 0), theta_resolution=60,\n phi_resolution=60), color=\"sandybrown\")\n for z in (-0.5, 0.5):\n plotter.add_mesh(pv.Plane(center=(ell/2, 0, z), i_size=ell + 1.6, j_size=1.9),\n color=\"gray\", opacity=0.12)\n plotter.add_mesh(surface, color=\"lightskyblue\", smooth_shading=True, opacity=0.9)\n plotter.camera_position = [(ell/2 + 1.0, -3.6, 2.4), (ell/2, 0, -0.05), (0, 0, 1)]\n plotter.show()\n\nshow_state(states[-1], cell.ell)",
"id": "cell-6"
},
{
"cell_type": "markdown",
"metadata": {},
"source": "The band at the end of the live run, mirrored to two unit cells: the dry spots on the\nbeads face the reader, and the band runs along the chain on either side of the beads.\n\n## Pressure and curvature\n\nThe pressure stays at $-1$ while the beads are immersed and jumps when they emerge. Then it\nfalls steeply as the band thins. The full run below (`drainage_runs.npz`, made with the same\nmethod on a mesh refined once more) continues to the snap. In the overlap the two agree to\n0.03% typically. Just after the emergence they differ by up to 1%, because the stored run\nfound the emergence by stepping (about 1% from the exact volume) instead of computing it.",
"id": "cell-7"
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": "runs = np.load(\"_static/drainage_runs.npz\")\nCELLS = 8 # eighth cells per bead\n\ndef pressure_axes(ax):\n ax.axhline(0, color=\"black\", lw=0.5)\n ax.set_xlabel(\"liquid volume per bead\")\n ax.set_ylabel(r\"pressure $\\Delta p$\")\n ax.secondary_yaxis(\"right\", functions=(lambda p: p/2, lambda h: 2*h)).set_ylabel(\n r\"mean curvature $H = \\Delta p/2$\")\n ax.invert_xaxis()\n\nfig, ax = plt.subplots(figsize=(8, 4.5))\nv, p, stage = runs[\"l1.1_volume\"], runs[\"l1.1_pressure\"], runs[\"l1.1_stage\"]\nfor s, label in ((0, \"immersed\"), (1, \"band\")):\n ax.plot(CELLS*v[stage == s], p[stage == s], \"-\", lw=2, color=f\"C{s}\",\n label=f\"{label} (stored, finer mesh)\")\nlive = np.array([(st[\"volume\"], st[\"pressure\"]) for st in states])\nax.plot(CELLS*live[:, 0], live[:, 1], \"k.\", ms=3, label=\"this run\")\nfor x, name in zip(runs[\"l1.1_events\"], (\"bead emerges\", \"snap\")):\n ax.axvline(CELLS*x, color=\"gray\", ls=\":\", lw=1)\n ax.annotate(name, (CELLS*x, 0.3), rotation=90, ha=\"right\", va=\"bottom\", fontsize=9)\npressure_axes(ax)\nax.set_title(r\"$\\ell = 1.1$\")\nax.legend(loc=\"lower left\")\nplt.show()",
"id": "cell-8"
},
{
"cell_type": "markdown",
"metadata": {},
"source": "## Three spacings\n\nThe same drainage, stored, for beads $0.04$ apart ($\\ell = 1$), $0.14$ apart ($\\ell = 1.1$)\nand well apart ($\\ell = 2$). Closer beads hold the band longer: it reaches lower volume and\nmuch lower pressure before it snaps.",
"id": "cell-9"
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": "fig, ax = plt.subplots(figsize=(8, 4.5))\nfor i, (key, ell) in enumerate(((\"l1\", 1.0), (\"l1.1\", 1.1), (\"l2\", 2.0))):\n v, p = runs[f\"{key}_volume\"], runs[f\"{key}_pressure\"]\n ax.plot(CELLS*v, p, \"-\", lw=2, color=f\"C{i}\", label=fr\"$\\ell = {ell:g}$\")\n snap = CELLS*runs[f\"{key}_events\"][1]\n ax.plot([snap], [p[-1]], \"o\", color=f\"C{i}\")\n emerge = CELLS*dh.Cell(ell).meniscus_exact(0.48)[0]\n print(f\"ell = {ell:g}: bead emerges at {emerge:.3f} (exact), \"\n f\"snap at {snap:.4f} per bead (Δp = {p[-1]:+.2f} just before)\")\npressure_axes(ax)\nax.legend(loc=\"lower left\")\nplt.show()",
"id": "cell-10"
},
{
"cell_type": "markdown",
"metadata": {},
"source": "## Movies\n\nEach frame is an equilibrium. Time in the movies is the volume removed, and they end at\nthe snap. The snap volume was found by stepping the volume in 0.5% steps near it, stopping\nwhen the band comes within $5\\cdot10^{-4}$ of the plane $y = 0$. At $\\ell = 2$ a mesh\nrefined once more moves it by 0.05%.\n\n$\\ell = 2$:\n\n\n\n$\\ell = 1$ (beads $0.04$ apart):\n\n\n\nLimits of the method:\n\n* **Quasi-static.** Every state is an equilibrium at the given volume, so the liquid moves\n between configurations only through the sequence of equilibria. Each event is\n instantaneous: the state jumps to the new branch at that volume.\n* **Symmetric.** The eighth cell imposes the mirror symmetries, so modes that break them\n (for example liquid gathering at one bead) are not followed.\n* **After the snap.** The liquid left sits around each bead, for example in rings at the\n bead–wall gaps. Following it needs other starting shapes, and it isn't computed here.\n\n## Related work\n\n* C. M. Cejas, L. A. Hough, C. Frétigny and R. Dreyfus, Effect of geometry on the\n dewetting of granular chains by evaporation, *Soft Matter* 14 (2018) 6994–7002,\n [doi:10.1039/C8SM01179F](https://doi.org/10.1039/C8SM01179F). Chains of beads drying,\n from loose to close packing, in experiments and with Surface Evolver: above a critical\n packing angle the liquid network breaks up into separate bridges.\n* J. L. Hilden and K. P. Trumble, Numerical analysis of capillarity in packed spheres:\n planar hexagonal-packed spheres, *J. Colloid Interface Sci.* 267 (2003) 463–474,\n [doi:10.1016/S0021-9797(03)00683-0](https://doi.org/10.1016/S0021-9797(03)00683-0).\n The capillary pressure at which liquid drains through a layer of packed spheres.\n\n## Next\n\nThat ends the worked examples. [Brakke's samples](brakke_films.ipynb) run the sample\ndatafiles that come with Evolver, each with his recipe and with `relax()`. The drainage\nband is one of two cases where `relax()`'s defaults alone fall short\n([where the defaults stop working](relaxing.md#where-the-defaults-stop-working)); the\n[stress cases](stress_cases.md) are the other hard problems it is checked against.",
"id": "cell-11"
}
],
"metadata": {
"kernelspec": {
"name": "python3",
"display_name": "Python 3",
"language": "python"
},
"language_info": {
"name": "python"
},
"mystnb": {
"execution_timeout": 900
}
},
"nbformat": 4,
"nbformat_minor": 5
}