{ "cells": [ { "cell_type": "markdown", "id": "s00", "metadata": {}, "source": "# Example: a droplet next to a bead in a slit\n\nA bead of radius $R$ sits between two parallel walls $z = \\pm 0.5$, almost touching\nthem ($R = 0.48$: a gap of 0.02 at each wall). Beads repeat along $x$ with spacing\n$\\ell$ (centre to centre), and a droplet wets the bead and both walls, reaching across to its\nneighbour's droplet. By symmetry only a quarter cell is modelled: mirror planes\n$x = 0$ and $y = 0$ through the bead's centre, and $x = \\ell/2$ halfway to the next\nbead, where the droplet meets its mirror image.\n\nThis builds on the [liquid bridge](liquid_bridge.ipynb): several walls at once, mirror\nplanes that cut the problem down to a quarter, measured contact angles on each solid,\nand a liquid region closed face by face for a volume mesh.\n\nParameters:\n\n* `ell`: the bead spacing, centre to centre ($\\ell \\ge 1$; the beads would touch at\n $\\ell = 2R$);\n* `theta_bead`, `theta_wall`: contact angles (degrees, measured inside the liquid)\n on the bead and on the walls;\n* `volume`: the liquid volume in the quarter cell;\n* `radius`: the bead radius (below 0.5; see *Notes* for why it must not touch).\n\nSurface tension 1; lengths in units of the wall spacing." }, { "cell_type": "code", "id": "s01", "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 = [700, 450]", "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "s02", "metadata": {}, "source": "## The model\n\nThe liquid–gas surface is the only meshed part. Its edges end on five constraints:\nthe walls, the bead, and the mirror planes $x = 0$ and $x = \\ell/2$ (the plane\n$y = 0$ cuts through the liquid; the surface does not reach it).\n\nThe walls and the bead are not meshed either. Their wetted areas (energy\n$-\\cos\\theta$ per area, Young's law) and their part of the liquid's volume enter\nas integrals along the contact lines, which\n[`pyse.constraints`](api.md) writes:\n\n* **Mirror planes** (`C.mirror`) carry no energy (the surface meets them at 90° by\n itself) and, being vertical, add nothing to Evolver's volume $\\int z\\, n_z\\, dA$.\n* **Walls** (`C.plane`): the wetted area is $\\oint t\\, ds$ in coordinates $s = x$,\n $t = \\pm y$ on the wall. That integrand vanishes along $y = 0$ and along the\n mirrors $x = 0, \\ell/2$, which close the wetted region.\n* **Bead** (`C.sphere`): areas on a sphere are integrals over the azimuth about the\n $z$ axis, and both mirror planes are meridians (no contribution). The bead's\n poles lie in the wetted region, where the azimuth is undefined; the builder is\n told so (`wet_poles`) and the quarter (`span`), and knows what the integrals\n leave out: a volume ($-\\pi R^3/3$, the quarter bead), which goes into the body's\n `volconst` since the body touches the bead, and the energy constant\n $-\\cos\\theta_\\text{bead}\\,\\pi R^2$, which `ev.energy` includes (see\n [Energies](getting_started.md#energies)).\n\nThe contact angles are given as the datafile parameters `theta_wall` and\n`theta_bead`, so they can be changed on a loaded surface.\n\nThe starting shape is a flat sheet of liquid at $y = w$, from $x = 0$ to $\\ell/2$ and\nwall to wall, with the bead cut out of it ($w = 0.4$; relaxing meets the target\nvolume). Its volume, areas and energy are known exactly, which checks all the signs." }, { "cell_type": "code", "id": "s03", "metadata": {}, "source": "def constraints(ell, radius):\n # normals point into the liquid; contact angles are datafile parameters\n return {1: C.plane((0, 0, 1), -0.5, contact_angle=\"theta_wall\"), # bottom wall\n 2: C.plane((0, 0, -1), point=(0, 0, 0.5), contact_angle=\"theta_wall\"),\n 3: C.sphere((0, 0, 0), radius, contact_angle=\"theta_bead\",\n wet_poles=(\"north\", \"south\"), span=np.pi/2), # a quarter bead\n 4: C.mirror(\"x\"), 5: C.mirror(\"x\", ell/2)}\n\ndef sheet_volume(ell, w, radius):\n # the liquid behind a flat sheet at y = w: a box minus the bead's part of it\n return ell/2*w - np.pi*(radius**2*w - w**3/3)/2\n\ndef load_slit(ell, theta_bead, theta_wall, volume, radius=0.48, w=0.4, ns=6, nt=8):\n rho = np.sqrt(radius**2 - w**2) # the bead's cut in the plane y = w\n seg, arc = 0.5 - rho, np.pi*rho\n total = 2*seg + arc\n def left(t):\n # the sheet's edge at x ~ 0, by arc length: the x = 0 mirror below the\n # bead, around the bead, the mirror above it\n s = t*total\n if s <= seg:\n return 0.0, -0.5 + s, 4\n if s <= seg + arc:\n a = -np.pi/2 + (s - seg)/rho\n return rho*np.cos(a), rho*np.sin(a), 3\n return 0.0, rho + s - seg - arc, 4\n corners = [0, seg/total, (seg + arc)/total, 1]\n ts = np.unique(np.concatenate([np.linspace(corners[i], corners[i + 1], n + 1)\n for i, n in enumerate((2, nt, 2))]))\n ss = np.linspace(0, 1, ns + 1)\n vertices, on = [], []\n for t in ts:\n xl, zl, k = left(t)\n for s in ss: # a ruled sheet from the left edge to the x = ell/2 mirror\n vertices.append(((1 - s)*xl + s*ell/2, w, (1 - s)*zl + s*(-0.5 + t)))\n c = set()\n if s == 0:\n c.add(k)\n if np.isclose(abs(zl), rho) and np.isclose(xl, 0):\n c |= {3, 4} # where the bead meets the mirror\n if s == 1: c.add(5)\n if t == 0: c.add(1)\n if t == 1: c.add(2)\n on.append(c)\n n = len(ss)\n faces = [[i*n + j, (i + 1)*n + j, (i + 1)*n + j + 1, i*n + j + 1] # normals +y\n for i in range(len(ts) - 1) for j in range(ns)]\n ev = pyse.Evolver()\n ev.load_arrays(np.array(vertices), faces, constraints=constraints(ell, radius),\n vertex_constraints=on, # the constraints of each vertex\n bodies=[pyse.Body(volume=volume)],\n parameters={\"radius\": radius, \"ell\": ell, \"theta_bead\": theta_bead,\n \"theta_wall\": theta_wall})\n return ev", "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "s04", "metadata": {}, "source": "The check: a fine flat sheet (no relaxation), with no volume target, against the\nexact volume and energy. `ev.energy` includes the constant\n$-\\cos\\theta_\\text{bead}\\,\\pi R^2$ that Evolver's `ev.total_energy` leaves out." }, { "cell_type": "code", "id": "s05", "metadata": {}, "source": "ell, w, R, tb, tw = 2.0, 0.3, 0.48, 40.0, 60.0\nev = load_slit(ell, tb, tw, None, R, w=w, ns=32, nt=32)\nev.recalc()\nrho2 = R**2 - w**2\nexact_energy = (ell/2 - np.pi*rho2/2) - 2*np.cos(np.radians(tw))*ell/2*w \\\n - np.cos(np.radians(tb))*np.pi*R*w\nprint(f\"volume {ev.bodies().volume[0]:.6f}, exact {sheet_volume(ell, w, R):.6f}\")\nprint(f\"energy {ev.energy:.6f}, exact {exact_energy:.6f}\")", "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "s06", "metadata": {}, "source": "The small differences are the polygonal contact line on the bead; they shrink\nfourfold when its points double.\n\n## Relaxing\n\n`ev.relax(levels=3)` relaxes each refinement level until the energy settles,\nequiangulating and averaging vertices every ten steps so facets at the contact lines\ndon't collapse; Newton steps finish it." }, { "cell_type": "code", "id": "s07", "metadata": {}, "source": "def slit_droplet(ell=2.0, theta_bead=40.0, theta_wall=60.0, volume=0.3, radius=0.48):\n ev = load_slit(ell, theta_bead, theta_wall, volume, radius)\n ev.relax(levels=3)\n return ev\n\nell, theta_bead, theta_wall, volume, radius = 2.0, 40.0, 60.0, 0.3, 0.48\n\nev = load_slit(ell, theta_bead, theta_wall, volume, radius)\nview = ev.live_view(color=\"lightskyblue\", walls=True) # the walls and the bead\nview.plotter.camera_position = [(2.2, 2.6, 1.6), (0.4, 0.2, 0), (0, 0, 1)]\nev.relax(levels=3, callback=view.update, every=10);", "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "s08", "metadata": {}, "source": "The quarter cell mirrored into the full picture: the droplet spans the slit and\njoins its neighbour's halfway between the beads." }, { "cell_type": "code", "id": "s09", "metadata": {}, "source": "def show_slit(ev):\n # the quarter cell mirrored to two neighbouring cells, with both beads\n ell, radius = ev.parameters[\"ell\"], ev.parameters[\"radius\"]\n solids = constraints(ell, radius)\n ev.plot(mirror=[\"y\", \"x\", (\"x\", ell/2)], color=\"lightskyblue\", smooth_shading=True,\n opacity=0.85, show_edges=False,\n walls=[solids[1], solids[2], solids[3], C.sphere((ell, 0, 0), radius)])\n\nshow_slit(ev)", "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "s10", "metadata": {}, "source": "## Checking the result\n\nThe measured contact angles should match the parameters on all three solids. The\npressure is the volume constraint's Lagrange multiplier (negative: the meniscus is\nconcave and pulls the walls together). `ev.eigen_counts().negative` counts unstable\ndirections (0: stable). The last line is a pinch-off check: the free surface must stay clear\nof the $y = 0$ mirror, or the droplet would really meet its mirror image and split." }, { "cell_type": "code", "id": "s11", "metadata": {}, "source": "# the angle between the surface's normal out of the liquid and the wall's out of\n# the solid, at each contact-line vertex (corners, on two walls, left out)\nfor k, name in ((1, \"bottom wall\"), (2, \"top wall\"), (3, \"bead\")):\n angles = ev.contact_angles(k)\n print(f\"{name:12s} {angles.mean():6.2f}° ± {angles.std():.2f}°\")\nprint(f\"pressure {ev.body(1).pressure:.4f}, energy {ev.energy:.5f}, \"\n f\"unstable directions {ev.eigen_counts().negative}\")\nprint(f\"closest approach to the y = 0 mirror: {ev.mesh().vertices[:, 1].min():.3f}\")", "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "s12", "metadata": {}, "source": "## Varying the wall contact angle\n\nThe less the liquid wets the walls, the less concave its meniscus between them and\nthe less negative its pressure. Cross-sections of the free surface in the mid-plane\n$z = 0$ (from above) and on the plane $x = \\ell/2$ (where the droplet meets its\nneighbour):" }, { "cell_type": "code", "id": "s13", "metadata": {}, "source": "fig, (top, side) = plt.subplots(1, 2, figsize=(10, 4), width_ratios=[2.2, 1])\ns = np.linspace(0, np.pi/2, 50)\ntop.fill_between(radius*np.cos(s), 0, radius*np.sin(s), color=\"sandybrown\", alpha=0.5)\npressures = []\nfor tw, color in zip((40, 50, 60, 70, 80), plt.cm.viridis(np.linspace(0, 0.9, 5))):\n ev = slit_droplet(ell, theta_bead, tw, volume, radius)\n pressures.append(ev.body(1).pressure)\n m = ev.mesh()\n mid = m.slice((0, 0, 1))[0] # the mid-plane z = 0\n top.plot(mid[:, 0], mid[:, 1], color=color, label=f\"θ_wall = {tw}°\")\n edge = m.slice((1, 0, 0), (ell/2, 0, 0))[0] # the mirror x = ell/2\n side.plot(edge[:, 1], edge[:, 2], color=color)\ntop.set_xlabel(\"x\"); top.set_ylabel(\"y\"); top.set_aspect(\"equal\")\ntop.set_title(\"mid-plane z = 0\"); top.legend(fontsize=\"small\")\nside.set_xlabel(\"y\"); side.set_ylabel(\"z\"); side.set_aspect(\"equal\"); side.set_xlim(0, 0.6)\nside.set_title(\"plane x = ℓ/2\")\nfor z in (-0.5, 0.5):\n side.axhline(z, color=\"gray\")\nplt.tight_layout()\nprint(\"pressures:\", np.round(pressures, 3))", "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "s14", "metadata": {}, "source": "## Exporting for finite elements\n\nThe liquid region is bounded by the free surface and by pieces of the walls, the\nbead and the mirror planes. Its free surface ends on all of them, so its single\nboundary loop runs along six arcs: the bottom wall, the $x = 0$ mirror in the gap\nbelow the bead, the bead, the gap above it, the top wall, and the $x = \\ell/2$ mirror.\nEach of the seven other faces of the liquid region is bounded by some of those\narcs and by known edges: straight lines where two planes meet, circles on the bead.\n\nSo the region is closed face by face. Each face is meshed with Gmsh as a 2D\npolygon, its boundary points fixed: the flat faces in their plane, the bead's\nafter a stereographic projection (from the point of the bead opposite the liquid),\nits new points mapped back onto the sphere. Faces share their edges point for point,\nso the result is watertight, and `volume_mesh` fills it with tetrahedra." }, { "cell_type": "code", "id": "s15", "metadata": {}, "source": "import gmsh\n\ndef boundary_loop(facets):\n # the vertex rows along a surface's (single) boundary loop, in order\n edges = np.sort(np.vstack([facets[:, [0, 1]], facets[:, [1, 2]], facets[:, [2, 0]]]), axis=1)\n unique, count = np.unique(edges, axis=0, return_counts=True)\n neighbours = {}\n for p, q in unique[count == 1].tolist():\n neighbours.setdefault(p, []).append(q)\n neighbours.setdefault(q, []).append(p)\n loop = [next(iter(neighbours))]\n while True:\n nxt = [v for v in neighbours[loop[-1]] if len(loop) < 2 or v != loop[-2]][0]\n if nxt == loop[0]:\n return np.array(loop)\n loop.append(nxt)\n\ndef mesh_polygon(uv, size):\n # triangles filling a 2D polygon with its points kept (no new boundary\n # points); returns the new interior points and triangles indexing\n # [polygon points, new points]\n gmsh.model.add(\"polygon\")\n geo = gmsh.model.geo\n corners = [geo.addPoint(u, v, 0, size) for u, v in uv]\n sides = [geo.addLine(corners[i], corners[(i + 1) % len(corners)]) for i in range(len(corners))]\n for side in sides:\n geo.mesh.setTransfiniteCurve(side, 2)\n geo.addPlaneSurface([geo.addCurveLoop(sides)])\n geo.synchronize()\n gmsh.model.mesh.generate(2)\n index = {int(gmsh.model.mesh.getNodes(0, c)[0][0]): i for i, c in enumerate(corners)}\n tags, xyz, _ = gmsh.model.mesh.getNodes(2, -1)\n for t in tags.tolist():\n index[int(t)] = len(index)\n triangles = np.vectorize(index.get)(gmsh.model.mesh.getElementsByType(2)[1].reshape(-1, 3))\n gmsh.model.remove()\n return xyz.reshape(-1, 3)[:, :2], triangles\n\ndef close_slit(ev):\n ell, R = ev.parameters[\"ell\"], ev.parameters[\"radius\"]\n m = ev.mesh()\n points = [m.vertices]\n def add(xyz): # new points; their rows\n points.append(np.atleast_2d(np.asarray(xyz, float)))\n start = sum(len(p) for p in points[:-1])\n return list(range(start, start + len(points[-1])))\n\n # the rim, cut into its six arcs at the corners (vertices on two constraints)\n rim = boundary_loop(m.facets)\n on = np.array([ev.on_constraint(k) for k in range(1, 6)]).T\n rim = np.roll(rim, -np.flatnonzero(on[rim].sum(1) == 2)[0])\n cuts = list(np.flatnonzero(on[rim].sum(1) == 2)) + [len(rim)]\n arcs = {}\n for a, b in zip(cuts[:-1], cuts[1:]):\n arc = list(rim[a:b]) + [rim[b % len(rim)]]\n k = int(np.flatnonzero(on[arc[1]])[0]) + 1 # its constraint\n name = {1: \"bottom\", 2: \"top\", 3: \"bead\", 5: \"far\"}.get(k)\n if k == 4:\n name = \"below\" if m.vertices[arc, 2].mean() < 0 else \"above\"\n arcs[name] = arc\n def corner(a, b):\n return (set(arcs[a]) & set(arcs[b])).pop()\n def run(name, start): # an arc from one of its ends\n return arcs[name] if arcs[name][0] == start else arcs[name][::-1]\n\n # edges that are not on the rim, with points about as far apart as the rim's\n h = np.linalg.norm(np.diff(m.vertices[rim], axis=0), axis=1).mean()\n def line(p, q):\n P = np.vstack(points)\n n = max(1, int(np.ceil(np.linalg.norm(P[q] - P[p])/h)))\n t = np.linspace(0, 1, n + 1)[1:-1, None]\n return [p] + (add(P[p] + t*(P[q] - P[p])) if len(t) else []) + [q]\n def bead_arc(p, q): # a great-circle arc on the bead\n P = np.vstack(points)\n a, b = P[p], P[q]\n angle = np.arccos(np.clip(a @ b/R**2, -1, 1))\n n = max(1, int(np.ceil(R*angle/h)))\n t = np.linspace(0, 1, n + 1)[1:-1]\n inner = [(np.sin((1 - s)*angle)*a + np.sin(s*angle)*b)/np.sin(angle) for s in t]\n return [p] + (add(inner) if inner else []) + [q]\n 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),\n (0, 0, -R), (0, 0, R), (R, 0, 0)])\n c_bb, c_bf = corner(\"bottom\", \"below\"), corner(\"bottom\", \"far\")\n c_tf, c_ta = corner(\"top\", \"far\"), corner(\"top\", \"above\")\n c_ba, c_bb2 = corner(\"bead\", \"above\"), corner(\"bead\", \"below\")\n floor, roof, back = line(b0, b1), line(t1, t0), line(b1, t1)\n axis_b, axis_t = line(pb, b0), line(t0, pt)\n meridian = bead_arc(pt, eq)[:-1] + bead_arc(eq, pb) # bead in y = 0\n gap_b, gap_t = bead_arc(c_bb2, pb), bead_arc(pt, c_ba) # bead in x = 0\n wall_b0, wall_b1 = line(c_bb, b0), line(b1, c_bf)\n wall_t0, wall_t1 = line(c_ta, t0), line(t1, c_tf)\n\n def loop(*chains): # closed polygon from chains that join end to start\n return [i for c in chains for i in c[:-1]]\n flat = lambda d: (lambda p: p[:, d])\n sphere_to_plane = lambda p: R*p[:, [0, 2]]/(p[:, [1]] + R)\n def plane_to_sphere(uv):\n s = 2*R**2/((uv**2).sum(1) + R**2)\n return np.column_stack([s*uv[:, 0], -R + s*R, s*uv[:, 1]])\n faces = { # name: (polygon, to 2D, back to 3D, direction out of the liquid)\n \"bottom wall\": (loop(run(\"bottom\", c_bf), wall_b0, floor, wall_b1), flat([0, 1]), None, (0, 0, -1)),\n \"top wall\": (loop(run(\"top\", c_ta), wall_t1[::-1], roof, wall_t0[::-1]), flat([0, 1]), None, (0, 0, 1)),\n \"mirror x = ell/2\": (loop(run(\"far\", c_tf), wall_b1[::-1], back, wall_t1), flat([1, 2]), None, (1, 0, 0)),\n \"mirror x = 0, below\": (loop(run(\"below\", c_bb), gap_b, axis_b, wall_b0[::-1]), flat([1, 2]), None, (-1, 0, 0)),\n \"mirror x = 0, above\": (loop(run(\"above\", c_ba), wall_t0, axis_t, gap_t), flat([1, 2]), None, (-1, 0, 0)),\n \"mirror y = 0\": (loop(floor, back, roof, axis_t, meridian, axis_b), flat([0, 2]), None, (0, -1, 0)),\n \"bead\": (loop(run(\"bead\", c_bb2), gap_t[::-1], meridian, gap_b[::-1]),\n sphere_to_plane, plane_to_sphere, None),\n }\n out = m.facet_bodies[:, 0] == 1\n cells = [np.where(out[:, None], m.facets, m.facets[:, ::-1])] # free surface\n tags, names = [np.zeros(len(m.facets), int)], {0: \"free surface\"}\n gmsh.initialize()\n gmsh.option.setNumber(\"General.Terminal\", 0)\n gmsh.option.setNumber(\"Mesh.MeshSizeMax\", h)\n for k, (name, (polygon, to_2d, to_3d, outward)) in enumerate(faces.items(), start=1):\n P = np.vstack(points)\n uv, triangles = mesh_polygon(to_2d(P[polygon]), h)\n if to_3d is None: # a flat face: put back its constant coordinate\n d = int(np.flatnonzero(np.abs(outward))[0])\n xyz = np.insert(uv, d, P[polygon[0], d], axis=1)\n else:\n xyz = to_3d(uv)\n rows = np.array(polygon + (add(xyz) if len(xyz) else []))\n T = rows[triangles]\n P = np.vstack(points)\n normal = np.cross(P[T[:, 1]] - P[T[:, 0]], P[T[:, 2]] - P[T[:, 0]])\n want = -P[T].mean(1) if outward is None else np.asarray(outward, float) # bead: inward\n flip = (normal*want).sum(1) < 0\n T[flip] = T[flip][:, ::-1]\n cells.append(T)\n tags.append(np.full(len(T), k))\n names[k] = name\n gmsh.finalize()\n T = np.vstack(cells)\n surface = pyse.BodySurface(1, np.vstack(points), T, \"triangle\", np.zeros(len(T), int),\n pyse.is_watertight(T), 0, cap_ids=np.concatenate(tags))\n return surface, names", "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "s16", "metadata": {}, "source": "A droplet refined twice (fewer elements than above), closed and checked: the\nenclosed volume against the target, and the energy recomputed from the face areas\n($A_\\text{free} - \\cos\\theta_\\text{wall} A_\\text{walls} - \\cos\\theta_\\text{bead} A_\\text{bead}$)\nagainst Evolver's." }, { "cell_type": "code", "id": "s17", "metadata": {}, "source": "ev = load_slit(ell, theta_bead, theta_wall, volume, radius)\nev.relax(levels=2)\nsurface, names = close_slit(ev)\n\narea = surface.areas() # per face: index k is names[k]\nby_face = {name: area[k] for k, name in names.items()}\nenergy = (by_face[\"free surface\"]\n - np.cos(np.radians(theta_wall))*(by_face[\"bottom wall\"] + by_face[\"top wall\"])\n - np.cos(np.radians(theta_bead))*by_face[\"bead\"])\nprint(\"watertight:\", surface.watertight, f\"| {len(surface.cells)} triangles\")\nprint(f\"enclosed volume {surface.volume:.5f} (target {volume})\")\n# the energy is a small difference of larger areas: the faces (flat caps on the walls)\n# and Evolver's line integrals differ at second order, 4e-4 here and 1e-4 one\n# refinement further\nprint(f\"energy from the faces {energy:.5f}, Evolver's {ev.energy:.5f}\")\nfor name, A in by_face.items():\n print(f\" {name:20s} area {A:.4f}\")", "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "s18", "metadata": {}, "source": "The volume mesh, with a physical group per face (named after it) and the body's\nfor the tetrahedra; `.msh` for `dolfinx.io.gmshio`, XDMF for others:" }, { "cell_type": "code", "id": "s19", "metadata": {}, "source": "import os, tempfile\nimport meshio\n\nfolder = tempfile.mkdtemp()\nmesh = surface.volume_mesh(size=0.05, names=names, path=os.path.join(folder, \"slit.msh\"))\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\"]\nfor g, n in zip(*np.unique(groups[\"triangle\"], return_counts=True)):\n print(f\" group {g}: {names[g - 1]:20s} {n} boundary triangles\")\nmeshio.write(os.path.join(folder, \"slit.xdmf\"),\n meshio.Mesh(P, [(\"tetra\", tets)], cell_data={\"region\": [groups[\"tetra\"]]}))\nmeshio.write(os.path.join(folder, \"slit_boundary.xdmf\"),\n meshio.Mesh(P, [(\"triangle\", mesh.cells_dict[\"triangle\"])],\n cell_data={\"boundary\": [groups[\"triangle\"]]}))\n\ngrid = pv.from_meshio(meshio.Mesh(P, [(\"tetra\", tets)]))\nplotter = pv.Plotter()\nplotter.add_mesh(grid.clip(normal=\"z\", invert=True, crinkle=True), color=\"lightskyblue\",\n show_edges=True)\nplotter.camera_position = [(1.4, 1.6, 1.3), (0.5, 0.2, -0.25), (0, 0, 1)]\nplotter.show()", "outputs": [], "execution_count": null }, { "cell_type": "markdown", "id": "s20", "metadata": {}, "source": "## Notes and limits\n\n* **A touching bead.** With $R = 0.5$ the bead touches the walls at single points,\n and a wetting liquid is drawn into the wedge around them, which closes to zero\n thickness; the mesh can't follow it there. A small gap ($R = 0.48$) ends the\n wedge in a thin but finite meniscus.\n* **Pinch-off.** Between walls a distance 1 apart, the meniscus sags inward by\n about $0.5\\,(1 - \\sin\\theta_\\text{wall})/\\cos\\theta_\\text{wall}$ (0.13 at 60°, 0.23\n at 40°). With too little liquid for that, the free surface reaches the $y = 0$\n mirror: physically the droplet splits into a top and a bottom part, which the\n mesh can't represent. Check the closest approach to the mirror; give more volume\n (or a larger wall angle).\n* **Non-wetting beads** ($\\theta_\\text{bead} > 90°$) push the liquid away until it\n pinches in the same way.\n* **Stability.** A relaxed droplet is an equilibrium; the relax result's `stable`\n (or `ev.eigen_counts().negative`) tells whether it is a stable one.\n* **Export.** The liquid's free surface ends on several solids and mirrors at once,\n so `body_surfaces(cap=True)` (one cap per opening, on one constraint) cannot\n close it; `close_slit` above closes it face by face for this geometry.\n\n## Next\n\n[Draining the slit](slit_drainage.ipynb) takes the same geometry through a long\nsequence of equilibria, removing liquid step by step until the liquid snaps." } ], "metadata": { "kernelspec": { "name": "python3", "display_name": "Python 3", "language": "python" }, "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 5 }