{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b00",
   "metadata": {},
   "source": "# Example: a liquid bridge between two beads\n\nA droplet wets two nearly touching spherical beads and forms a **capillary\nbridge**: the shape of minimal surface energy at fixed liquid volume, meeting\neach bead at a prescribed **contact angle** $\\theta$. Bridges like this hold wet\ngranular media together; the force they exert depends on the gap between the\nbeads and on $\\theta$.\n\nThis is the drop on a table of [Getting started](getting_started.md) with two spheres\nfor the table: contact angles on curved solids, the force between them, and a volume\nmesh of the liquid. It uses the [conventions](getting_started.md#conventions) set out\nthere, and `relax()` as in the [tutorial](tutorial.ipynb).\n\nThe beads have radius 1 and the liquid–gas surface tension is 1, so lengths are\nin bead radii, energies in surface tension × radius², and forces in surface\ntension × radius. The adjustable parameters:\n\n* `theta`: the contact angle in degrees, measured inside the liquid\n  (below 90°: wetting);\n* `gap`: the distance between the bead surfaces (0: touching);\n* `volume`: the liquid volume."
  },
  {
   "cell_type": "code",
   "id": "b01",
   "metadata": {},
   "source": "import numpy as np\nimport matplotlib.pyplot as plt\nimport pyvista as pv\nimport pysurfaceevolver as pyse\nfrom pysurfaceevolver import constraints as C\n\npv.set_jupyter_backend(\"static\")    # \"trame\" in Jupyter for interactive views\npv.global_theme.window_size = [600, 450]",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "b02",
   "metadata": {},
   "source": "## The model\n\nOnly the liquid–gas surface is made of facets. The beads are **constraints**:\nthe two contact lines live on the spheres\n$x^2 + y^2 + (z \\pm c)^2 = 1$, with $c = 1 + \\mathrm{gap}/2$, and slide freely on\nthem. The wetted part of each bead is not meshed; it enters as integrals along\nthe contact line, which Evolver adds to the edges on a constraint (as in the\n`mound.fe` sample, there on a plane):\n\n* **Energy.** Wetting a unit area of bead replaces solid–gas by solid–liquid\n  surface, a change of $\\gamma_{SL} - \\gamma_{SG} = -\\cos\\theta$ (Young's law).\n  So the energy is the liquid area plus $-\\cos\\theta$ times the wetted area of\n  each bead.\n* **Volume.** Evolver computes a body's volume from its facets as\n  $\\int z\\, n_z\\, dA$; the wetted caps close the liquid's boundary, so their part of\n  that integral is added too.\n\n[`pyse.constraints.sphere`](api.md) writes both integrals. Over the azimuth $\\phi$\nabout the $z$ axis, with $h$ the height above a bead's centre, the wetted area is\n$\\oint R\\,h\\,d\\phi$ up to the pole facing the gap (the integrand is singular at the\npoles), and the bead's share of $\\int z\\, n_z\\, dA$ is $\\oint (c_z h^2/2 + h^3/3)\\, d\\phi$.\nTold which pole is wet (`wet_poles`), the builder knows what those integrals\nleave out: the volume over the pole, which it takes into the volume integral, and a\nconstant energy $-\\cos\\theta \\cdot 2\\pi R^2$ per bead, which `ev.energy` includes\nand Evolver's own `ev.total_energy` leaves out (see\n[Energies](getting_started.md#energies)). It doesn't change with the gap, so forces\ndon't see it. The contact angle is the datafile parameter `theta` and the bead\ncentres follow the parameter `gap` ($c$ = `radius + gap/2` in the formulas), so both\ncan be changed on a loaded surface."
  },
  {
   "cell_type": "code",
   "id": "b03",
   "metadata": {},
   "source": "def beads(radius=1.0):\n    # the two beads, wetted around the poles that face each other; their\n    # centres follow the datafile parameter gap, the contact angle theta\n    return {1: C.sphere((0, 0, \"-(radius + gap/2)\"), radius, contact_angle=\"theta\",\n                        wet_poles=(\"north\",)),\n            2: C.sphere((0, 0, \"radius + gap/2\"), radius, contact_angle=\"theta\",\n                        wet_poles=(\"south\",))}\n\ndef load_bridge(gap, theta, volume, radius=1.0, a=0.5):\n    # a cylinder of liquid of radius a between the beads, its ends on the bead\n    # surfaces (relaxing meets the volume); normals point out of the liquid\n    zb = radius + gap/2 - np.sqrt(radius**2 - a**2)\n    tube = pyse.shapes.tube(a, -zb, zb, nz=4, nt=12)\n    ev = pyse.Evolver()\n    ev.load_arrays(tube.vertices, tube.faces, constraints=beads(radius),\n                   vertex_constraints={1: tube.rims[\"bottom\"], 2: tube.rims[\"top\"]},\n                   bodies=[pyse.Body(volume=volume)],      # volconst: from the beads\n                   parameters={\"radius\": radius, \"gap\": gap, \"theta\": theta})\n    return ev",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "b04",
   "metadata": {},
   "source": "## Relaxing\n\n`ev.relax(levels=3)` refines the mesh three times (to about 6000 facets), with\ngradient steps in between, and ends with Newton steps. Every ten steps it\nequiangulates and averages vertices: without that, facets at the contact line\nget thin and collapse when the contact line moves. On the coarse starting mesh\nthe contact line can still collapse within a few steps; `relax()` notices the\njump in the residual, undoes those steps and moves on\n([Relaxing a surface](relaxing.md) has the details).\n\nVery low contact angles at larger gaps make the contact line travel far over\nthe beads; there the bridge is relaxed at 40° first and `theta` is then lowered\nin steps of 10°, relaxing after each."
  },
  {
   "cell_type": "code",
   "id": "b05",
   "metadata": {},
   "source": "def liquid_bridge(gap=0.2, theta=40.0, volume=0.3, callback=None):\n    ev = load_bridge(gap, max(theta, 40.0), volume)\n    ev.relax(levels=3, callback=callback, every=10)\n    while ev.parameters[\"theta\"] > theta + 1e-9:\n        ev.parameters[\"theta\"] = max(theta, ev.parameters[\"theta\"] - 10)\n        ev.relax()\n    return ev",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "b06",
   "metadata": {},
   "source": "Here is a bridge relaxing from its starting cylinder. In Jupyter (trame backend)\nthe live view redraws during the run; these docs show its first frame, and the\nrelaxed bridge below."
  },
  {
   "cell_type": "code",
   "id": "b07",
   "metadata": {},
   "source": "gap, theta, volume = 0.2, 40.0, 0.3\n\nev = load_bridge(gap, theta, volume)\nview = ev.live_view(color=\"lightskyblue\", opacity=0.9, walls=True)   # the beads too\nview.plotter.camera_position = [(3.2, -3.2, 1.2), (0, 0, 0), (0, 0, 1)]\nev.relax(levels=3, callback=view.update, every=10);",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "id": "b08",
   "metadata": {},
   "source": "ev.plot(color=\"lightskyblue\", smooth_shading=True, show_edges=False, walls=True)\nev",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "b09",
   "metadata": {},
   "source": "## Checking the result\n\nThe contact angle is not imposed on the mesh: it comes out of the energy\nbalance at the contact line. Measuring it is a check of the whole model: at each\ncontact-line vertex, the angle between the liquid surface's normal (out of the\nliquid) and the bead's outward normal.\n\nThe pressure is the Lagrange multiplier of the volume constraint. The force\nbetween the beads follows from the neck, the narrowest circle in the mid-plane:\nthere the surface is vertical by symmetry, so the force transmitted through\nthat plane is $F = 2\\pi r_n \\gamma - \\pi r_n^2\\, \\Delta p$ (positive: the beads\nattract). It must also equal the change of energy with the gap at fixed volume,\n$F = dE/d\\,\\mathrm{gap}$, here from bridges built at the gap plus and minus 0.01\n(`ev.derivative(name)` does this for a parameter the constraint formulas use; the\nbead centres are numbers in them)."
  },
  {
   "cell_type": "code",
   "id": "b10",
   "metadata": {},
   "source": "def neck_radius(ev):\n    # mean radius where the liquid surface crosses the mid-plane z = 0\n    ring = ev.mesh().slice((0, 0, 1))[0][:-1]            # a closed curve\n    return np.hypot(ring[:, 0], ring[:, 1]).mean()\n\ndef force(ev):\n    rn, p = neck_radius(ev), ev.body(1).pressure\n    return 2*np.pi*rn - np.pi*rn**2*p\n\n# (the mesh is rotationally symmetric, so every contact-line vertex gives the same angle;\n#  it is first order in the edge length, hence a little above theta)\nangles = np.concatenate([ev.contact_angles(1), ev.contact_angles(2)])\nprint(f\"contact angle {angles.mean():.2f}° ± {angles.std():.2f}° (theta = {theta}°)\")\nprint(f\"pressure {ev.body(1).pressure:.4f}, neck radius {neck_radius(ev):.4f}, \"\n      f\"force {force(ev):.4f}, unstable directions {ev.eigen_counts().negative}\")\n\nh = 0.01\ndE = (liquid_bridge(gap + h, theta, volume).energy\n      - liquid_bridge(gap - h, theta, volume).energy) / (2*h)\nprint(f\"dE/dgap {dE:.4f}\")",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "b11",
   "metadata": {},
   "source": "The measured angle is within about a degree of `theta` (it is estimated from\nflat facets, and is less accurate at very low and very high angles), and the two\nforces agree to a fraction of a percent. `ev.eigen_counts().negative` counts the\nnegative eigenvalues of the Hessian: none, so this bridge is stable.\n\n## Profiles for several contact angles\n\nThe more the liquid wets the beads, the further it spreads over them and the\nnarrower its neck. Cross-sections through the plane $y = 0$:"
  },
  {
   "cell_type": "code",
   "id": "b12",
   "metadata": {},
   "source": "def profile(ev):\n    # (x, z) where the liquid surface crosses the half-plane y = 0, x > 0\n    curve = max(ev.mesh().slice((0, 1, 0)), key=lambda c: c[:, 0].mean())\n    return curve[:, 0], curve[:, 2]\n\nfig, ax = plt.subplots(figsize=(6, 4))\nc = 1 + gap/2\ns = np.linspace(0, 2*np.pi, 200)\nfor z in (-c, c):\n    ax.fill(np.cos(s), z + np.sin(s), color=\"sandybrown\", alpha=0.5, lw=0)\nfor th, color in zip((20, 40, 60, 80), plt.cm.viridis(np.linspace(0, 0.9, 4))):\n    x, z = profile(liquid_bridge(gap, th, volume))\n    ax.plot(x, z, color=color, label=f\"θ = {th}°\")\n    ax.plot(-x, z, color=color)\nax.set_aspect(\"equal\")\nax.set_xlim(-0.9, 0.9)\nax.set_ylim(-0.55, 0.55)\nax.set_xlabel(\"x\")\nax.set_ylabel(\"z\")\nax.legend(loc=\"center right\", fontsize=\"small\");",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "b13",
   "metadata": {},
   "source": "## Force against gap\n\nWetting liquids pull the beads together, the more strongly the smaller the\ncontact angle. Non-wetting liquids ($\\theta \\gtrsim 90°$) push them apart at small\ngaps, but there the symmetric bridge is **unstable**: the Hessian has a pair of\nnegative eigenvalues, sideways shifts in $x$ and $y$, so the drop would rather\nslide out of the gap. Such points are drawn hollow. (Each point is a fresh\nbridge, about half a second each.)"
  },
  {
   "cell_type": "code",
   "id": "b14",
   "metadata": {},
   "source": "import warnings\n\ngaps = np.linspace(0.0, 0.6, 7)\nfig, ax = plt.subplots(figsize=(6, 3.5))\nfor th, color in zip((20, 60, 100), (\"tab:blue\", \"tab:green\", \"tab:red\")):\n    forces, stable = [], []\n    for g in gaps:\n        with warnings.catch_warnings():   # the unstable bridges: see above\n            warnings.simplefilter(\"ignore\", pyse.UnstableEquilibriumWarning)\n            b = liquid_bridge(g, th, volume)\n        # every Evolver is a handle to the one engine: use b before the next build\n        forces.append(force(b))\n        stable.append(b.eigen_counts().negative == 0)\n    forces, stable = np.array(forces), np.array(stable)\n    ax.plot(gaps, forces, \"-\", color=color, label=f\"θ = {th}°\")\n    ax.plot(gaps[stable], forces[stable], \"o\", color=color)\n    ax.plot(gaps[~stable], forces[~stable], \"o\", mfc=\"white\", color=color)\nax.axhline(0, color=\"gray\", lw=0.5)\nax.set_xlabel(\"gap\")\nax.set_ylabel(\"force (attractive > 0)\")\nax.set_title(f\"liquid volume {volume}\")\nax.legend();",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "b14a",
   "metadata": {},
   "source": "## How far a bridge stretches\n\nPull the beads apart and a bridge of fixed volume breaks at some gap. The beads follow\nthe parameter `gap`, so `recipes.arclength` can pull them, as on the\n[bridge page](bridge_rupture.md) it shrinks the liquid, until the bridge's first\nunstable mode. What ends the symmetric bridge is one of two things. At its **fold** no\nsymmetric bridge of that volume spans a wider gap: the neck thins and pinches in the\nmiddle. Or, before that, at a **branch point**, bridges whose neck sits off the\nmid-plane branch off and the symmetric one turns unstable: the liquid moves toward one\nbead. The exact axisymmetric bridges (`stress_helpers.py`, as on the bridge page) say\nwhich comes first: here, the volume that breaks at gaps 0.4 and 0.7. Pulling a bridge\nof that volume apart should then break it at that gap.\n\nA rule widely used in simulations of wet granular matter is the rupture distance of Lian,\nThornton and Adams, $s_c \\approx (1 + \\theta/2)\\,V^{1/3}$ ($\\theta$ in radians, lengths in\nbead radii)."
  },
  {
   "cell_type": "code",
   "id": "b14b",
   "metadata": {},
   "source": "from pysurfaceevolver import recipes\nfrom stress_helpers import bridge_spheres_branch, bridge_spheres_fold\n\ndef breaking_volume(gap, theta):\n    # the least volume a stable symmetric bridge holds at this gap (exact), and why\n    necks = (0.005, 0.6)\n    fold = bridge_spheres_fold(gap, theta, necks=necks)\n    try:\n        branch = bridge_spheres_branch(gap, theta, necks=necks)\n    except ValueError:\n        branch = None\n    if branch is not None and branch[\"neck\"] > fold[\"neck\"]:\n        return branch[\"volume\"], \"branch point\"\n    return fold[\"volume\"], \"fold\"\n\nresults = {}\nfor th in (20.0, 60.0):\n    for g in (0.4, 0.7):\n        v, why = breaking_volume(g, th)\n        # a bridge of that volume at a gap 0.07 narrower (built the usual way and\n        # shrunk to it), pulled apart until it turns unstable\n        b = liquid_bridge(g - 0.07, th, 0.15)\n        stable = None\n        with warnings.catch_warnings():\n            warnings.simplefilter(\"ignore\", pyse.UnstableEquilibriumWarning)\n            for step in recipes.continuation(b, b.body(1), np.geomspace(0.15, v, 8)[1:]):\n                pass\n            for step in recipes.arclength(b, \"gap\", 0.01, 60):\n                if step.result.negative_modes:\n                    break\n                stable = step.value                 # the widest stable gap so far\n        lian = (1 + np.radians(th)/2)*v**(1/3)\n        results[th, g] = stable, v\n        print(f\"θ = {th:.0f}°, volume {v:.4f}: stable up to gap {stable:.3f} \"\n              f\"(exact {g}, at the {why}); Lian et al. {lian:.3f}\")\n\nfig, ax = plt.subplots(figsize=(6, 3.5))\nvs = np.geomspace(0.01, 0.3, 100)\nfor th, color in ((20.0, \"tab:blue\"), (60.0, \"tab:red\")):\n    ax.plot(vs, (1 + np.radians(th)/2)*vs**(1/3), color=color, lw=1, label=f\"Lian et al., θ = {th:.0f}°\")\n    for g in (0.4, 0.7):\n        stable, v = results[th, g]\n        ax.plot(v, g, \"o\", mfc=\"white\", color=color)\n        ax.plot(v, stable, \".\", color=color)\nax.set_xscale(\"log\")\nax.set_xlabel(\"liquid volume\")\nax.set_ylabel(\"gap at which the bridge breaks\")\nax.set_title(\"open: exact; dots: Evolver's widest stable bridge\", fontsize=\"small\")\nax.legend(fontsize=\"small\");",
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "id": "b14c",
   "metadata": {},
   "source": "At 20° the bridge reaches its fold: the neck pinches in the middle, and `arclength`\nturns back there. At 60° the branch point comes first, well before the fold (at gap 0.4\nthe fold is at volume 0.0136), and the liquid would go to one bead. Evolver's widest\nstable bridge is within 0.01 of the exact gap (the step), on a mesh of about 6000\nfacets. Lian et al.'s rule gives the gap to within 5% in all four, without telling the\ntwo ways of breaking apart."
  },
  {
   "cell_type": "markdown",
   "id": "b15",
   "metadata": {},
   "source": "## Exporting for finite elements\n\nFor a simulation of the liquid (flow inside the bridge, say), the liquid region\nneeds a **closed surface** and a **volume mesh**. Evolver's mesh is only the\nliquid–gas interface: the wetted parts of the beads are integrals, not facets,\nso the body's surface has two openings. `body_surfaces(cap=True)` closes each\nopening with a cap on the constraint its rim lies on (Evolver projects the cap's\npoints, so this works for any constraint formula), and tells the caps apart.\n\nEvery `Evolver` object is a handle to the one engine, and the sweep above left\nits last bridge there; build the one to export again:"
  },
  {
   "cell_type": "code",
   "id": "b16",
   "metadata": {},
   "source": "ev = liquid_bridge(gap, theta, volume)\nsurface = ev.body_surfaces(cap=True)[1]          # body 1, the liquid\n\narea = surface.areas()                           # the liquid-gas surface, cap 1, cap 2\nprint(\"watertight:\", surface.watertight, \"| caps on constraints:\", surface.cap_constraints)\nprint(f\"enclosed volume {surface.volume:.5f} (target {volume})\")\nprint(f\"liquid-gas area {area[0]:.5f} (Evolver {ev.total_area:.5f}), \"\n      f\"wetted areas {area[1]:.4f}, {area[2]:.4f}\")",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "b17",
   "metadata": {},
   "source": "`volume_mesh` fills the closed surface with tetrahedra using\n[Gmsh](https://gmsh.info) (`pip install gmsh`), keeping the surface triangles as\nthey are. The physical groups: 1 for the liquid–gas facets, 2 and 3 for the caps\n(cap 1 and 2), and the body id for the tetrahedra. A `.msh` path keeps the group\nnames; that is what `dolfinx.io.gmshio.read_from_msh` reads."
  },
  {
   "cell_type": "code",
   "id": "b18",
   "metadata": {},
   "source": "import os, tempfile\nimport meshio\n\nfolder = tempfile.mkdtemp()\nmesh = surface.volume_mesh(size=0.04, path=os.path.join(folder, \"bridge.msh\"))\n\ntets, P = mesh.cells_dict[\"tetra\"], mesh.points\nv = np.einsum(\"ij,ij->i\", P[tets[:, 1]] - P[tets[:, 0]],\n              np.cross(P[tets[:, 2]] - P[tets[:, 0]], P[tets[:, 3]] - P[tets[:, 0]]))\nprint(f\"{len(tets)} tetrahedra, volume {np.abs(v).sum()/6:.5f}\")\ngroups = mesh.cell_data_dict[\"gmsh:physical\"]\nprint(\"boundary triangles per group:\", dict(zip(*np.unique(groups[\"triangle\"], return_counts=True))))",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "b19",
   "metadata": {},
   "source": "For other FEM codes, meshio writes the tetrahedra and the tagged boundary\ntriangles separately, for example as XDMF:"
  },
  {
   "cell_type": "code",
   "id": "b20",
   "metadata": {},
   "source": "meshio.write(os.path.join(folder, \"bridge.xdmf\"),\n             meshio.Mesh(P, [(\"tetra\", tets)], cell_data={\"region\": [groups[\"tetra\"]]}))\nmeshio.write(os.path.join(folder, \"bridge_boundary.xdmf\"),\n             meshio.Mesh(P, [(\"triangle\", mesh.cells_dict[\"triangle\"])],\n                         cell_data={\"boundary\": [groups[\"triangle\"]]}))\nprint(sorted(f for f in os.listdir(folder)))\n\n# a cut through the tetrahedra\ngrid = pv.from_meshio(meshio.Mesh(P, [(\"tetra\", tets)]))\nplotter = pv.Plotter()\nplotter.add_mesh(grid.clip(normal=\"y\", invert=False, crinkle=True), color=\"lightskyblue\",\n                 show_edges=True)\nfor zc in (-1 - gap/2, 1 + gap/2):\n    plotter.add_mesh(pv.Sphere(radius=1.0, center=(0, 0, zc)), color=\"sandybrown\", opacity=0.3)\nplotter.camera_position = [(0, -2.6, 0.4), (0, 0, 0), (0, 0, 1)]\nplotter.show()",
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "id": "b21",
   "metadata": {},
   "source": "## Notes and limits\n\n* **Changing parameters on a relaxed bridge.** `ev.parameters[\"theta\"] = ...` (or\n  `\"gap\"`) followed by `ev.relax()` works for small changes and for lowering the contact\n  angle; large jumps, or raising the angle so that the contact line has to\n  retreat, can leave folded facets at the contact line (the energy is still right,\n  but the measured angle is not). `liquid_bridge(...)` builds afresh each time,\n  which takes about half a second.\n* **Stability.** A relaxed bridge is an equilibrium, not necessarily a stable\n  one: check `ev.eigen_counts().negative`. Non-wetting bridges\n  at small gaps are saddles (see above).\n* **Rupture.** Past the volumes and gaps of\n  [how far a bridge stretches](#how-far-a-bridge-stretches) a bridge is unstable and\n  would break; the mesh cannot follow a pinch-off, so relaxing fails there (for example\n  volume 0.05 at gap 0.5).\n* **Unequal beads, gravity, …** follow the same pattern: another sphere radius or\n  centre in the constraint formulas, `Body(volume=..., density=...)` and `ev.gravity`.\n\n## Further reading\n\n* G. Lian, C. Thornton and M. Adams, A theoretical study of the liquid bridge forces\n  between two rigid spherical bodies, *J. Colloid Interface Sci.* 161 (1993) 138–147,\n  [doi:10.1006/jcis.1993.1452](https://doi.org/10.1006/jcis.1993.1452). The rupture\n  distance used above.\n\nA bridge between two beads is the smallest liquid structure in a wet bead pack. With more\nliquid, neighbouring bridges merge into clusters that wet three or more beads. Surface\nEvolver has computed these shapes in:\n\n* M. Scheel, R. Seemann, M. Brinkmann, M. Di Michiel, A. Sheppard, B. Breidenbach and\n  S. Herminghaus, Morphological clues to wet granular pile stability, *Nature Materials*\n  7 (2008) 189–193, [doi:10.1038/nmat2117](https://doi.org/10.1038/nmat2117). X-ray\n  tomography of wet bead packs, with bridges, trimers, pentamers and filled tetrahedra\n  computed in Surface Evolver: the pack's cohesion barely depends on the amount of liquid.\n* C. Semprebon, M. Scheel, S. Herminghaus, R. Seemann and M. Brinkmann, Liquid\n  morphologies and capillary forces between three spherical beads, *Phys. Rev. E* 94\n  (2016) 012907, [doi:10.1103/PhysRevE.94.012907](https://doi.org/10.1103/PhysRevE.94.012907).\n  Two or three merged bridges between three beads: their shapes, forces and stability\n  limits.\n* J. Wang, E. Gallo, B. François, F. Gabrieli and P. Lambert, Capillary force and rupture\n  of funicular liquid bridges between three spherical bodies, *Powder Technology* 305\n  (2017) 89–98, [doi:10.1016/j.powtec.2016.09.060](https://doi.org/10.1016/j.powtec.2016.09.060).\n  The force and rupture of liquid between three spheres pulled apart.\n* K. Melnikov, R. Mani, F. K. Wittel, M. Thielmann and H. J. Herrmann, Grain-scale\n  modeling of arbitrary fluid saturation in random packings, *Phys. Rev. E* 92 (2015)\n  022206, [doi:10.1103/PhysRevE.92.022206](https://doi.org/10.1103/PhysRevE.92.022206).\n  A model of whole random packings, from dry to saturated, that takes its bridges'\n  pressure and volume from Semprebon et al.'s Surface Evolver tables.\n\n## Next\n\n[Solder joints](solder_joints.ipynb) computes forces with `ev.derivative` (for a\nparameter that moves a pad) and lets a loaded chip find its own height. The\n[slit droplet](slit_droplet.ipynb) puts a liquid against several walls and mirror\nplanes at once. More drops against walls, from Brakke's datafiles:\n[drops and capillarity](brakke_capillarity.ipynb)."
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3",
   "language": "python"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}