Flow through a randomly packed bed (peclet.dem → peclet.flow)

Generate a random close packing with the DEM engine, characterise its microstructure, then push Stokes flow through the pore space and measure permeability.

dem

flow

porous-media

multi-method

permeability

Author

Peclet

Published

July 3, 2026

Runs on a free Colab CPU runtime — the first cell installs peclet from PyPI.

What you’ll learn

The gallery’s first cross-method example: two peclet codes in one pipeline. We use peclet.dem (discrete-element method) to build a random close packing of spheres, characterise its microstructure (porosity, coordination number, radial distribution), convert it to a signed distance function, and drive creeping flow through the pore space with peclet.flow to measure the Darcy permeability — compared against the classic Carman–Kozeny relation. This is the realistic, disordered counterpart to the periodic-lattice Zick–Homsy example.

import importlib.util, os, subprocess, sys
_local = os.environ.get("PECLET_LOCAL_BUILD")
if _local:
    for p in _local.split(os.pathsep):
        sys.path.insert(0, p)
elif importlib.util.find_spec("peclet") is None:
    subprocess.run([sys.executable, "-m", "pip", "install", "-q", "peclet"], check=True)

import numpy as np
import math, time
import matplotlib.pyplot as plt
from peclet import dem, flow

plt.rcParams.update({"figure.dpi": 130, "font.size": 10, "axes.grid": True,
                     "grid.alpha": 0.3, "axes.axisbelow": True,
                     "figure.facecolor": "white", "savefig.bbox": "tight"})
BLUE, RED, GREEN = "#1f77b4", "#d62728", "#2ca02c"

Step 1 — Build a random close packing with `peclet.dem`

We grow monodisperse spheres in a periodic box under a Lubachevsky–Stillinger protocol: the particles start small, grow while a thermostat keeps the assembly fluid, and the growth is feedback-controlled by the committed overlap so the system jams near random close packing. A final quench settles it into a rigid contact network.

def pack_bed(N=180, phi_ref=0.63, radius=0.5, seed=3, dt=0.002, limit_time=7.0, iters=60,
             temperature=1.0, scale_init=0.05, criterion=5e-3, cooling_time=5.0, quench=1200,
             growth_accel=1.02, growth_decay=0.85, growth_rate_init=0.5):
    # box sized so that full-scale spheres would give phi_ref; the growth stops *before* that,
    # at the jamming point, so the achieved density (and radii) come out below phi_ref.
    volp = (4 / 3) * math.pi * radius ** 3
    side = (N * volp / phi_ref) ** (1 / 3)
    half = side / 2.0
    gr = growth_rate_init
    cool = min(int(cooling_time / dt), int(limit_time / dt))
    rng = np.random.default_rng(seed)
    s = dem.Simulation(N)
    s.initialize(shape_type=1, radius=radius)
    s.set_domain((-half, -half, -half), (half, half, half))
    s.enable_periodicity(True, True, True); s.set_gravity(0, 0, 0)
    s.set_material_params(1.0, 1.0, 0.0); s.set_solver_iterations(iters, iters)   # restitution 1 while fluid
    pos = rng.uniform(-half, half, (N, 4)).astype(np.float32); pos[:, 3] = 1.0
    s.set_positions(pos)
    s.set_velocities(rng.normal(0.0, math.sqrt(temperature), (N, 3)).astype(np.float32))
    s.set_scales(np.full(N, 1.0, np.float32))
    s.set_growth_params(gr, scale_init); s.set_thermostat(temperature, dt)

    def overlap_frac():
        c = 2 * radius * float(s.get_scales().ravel().mean())
        return float(s.compute_overlaps()) / max(c, 1e-9)

    # MONOTONIC growth. When overlaps build, we RELAX them at the current size (dt = 0 steps grow
    # nothing) and only slow the growth rate — we never shrink the grains back. The old scheme
    # un-grew on every overlap spike, and because the GPU contact solve is atomically non-deterministic
    # (float atomic_add ordering) a transient spike could trigger a runaway of shrink-steps that
    # quenched a loose, unjammed state on some seeds. Growing in one direction only makes the final
    # density reproducible: every seed lands at the box-limited phi_ref with an isostatic Z ~ 6.
    for step in range(int(limit_time / dt)):
        if step == cool:                                    # anneal: mildly dissipative
            s.set_material_params(0.5, 1.0, 0.0); s.set_thermostat(0.0, 1e4 * dt)
        s.step(dt); mo = overlap_frac()
        gf = float(s.get_growth_factor())
        if mo > criterion:                                  # overlaps too big -> relax at CURRENT size
            it = 0; prev = mo
            while it < 40:
                s.step(0.0); it += 1; mn = overlap_frac()   # dt = 0: pure overlap removal, no growth
                if mn < criterion or (it > 8 and mn > 0.98 * prev):
                    break                                   # cleared, or plateaued (jammed at this size)
                prev = mn
            if overlap_frac() > criterion:                  # still jammed -> slow growth, never un-grow
                gr = max(gr * growth_decay, 0.02); s.set_growth_params(gr, gf)
        elif gf >= 1.0 - 1e-6:
            break                                           # reached full size (box sized for phi_ref)
        else:                                               # room to grow -> accelerate
            gr = min(gr * growth_accel, growth_rate_init); s.set_growth_params(gr, gf)
    s.set_material_params(0.0, 0.0, 0.0); s.set_thermostat(0.0, 10 * dt)   # quench into the rigid network
    for _ in range(quench):
        s.step(dt)
    # EFFECTIVE radius = per-particle scale * global growth factor * base radius. The growth factor
    # (<= 1) is essential — omitting it overstates the radii and fakes overlaps in g(r).
    r_eff = radius * s.get_scales().ravel() * float(s.get_growth_factor())
    pos = s.get_positions()[:, :3].astype(float)
    phi = float(np.sum(4 / 3 * np.pi * r_eff ** 3) / side ** 3)
    return pos, r_eff.astype(float), float(side), phi

t0 = time.time()
pos, r, side, phi = pack_bed()
print(f"random close packing: N={len(pos)}  φ={phi:.3f}  porosity ε={1-phi:.3f}  "
      f"mean radius={r.mean():.4f}  box={side:.2f}  ({time.time()-t0:.0f}s)")

random close packing: N=180  φ=0.630  porosity ε=0.370  mean radius=0.5000  box=5.31  (21s)

Step 2 — Characterise the microstructure

Porosity $\varepsilon = 1-\phi$, the mean coordination number $Z$ (contacts per sphere — near the isostatic value 6 for a frictionless jammed packing), and the radial distribution function $g(r)$ with its sharp contact peak at one diameter.

def microstructure(pos, r, side):
    box = np.array([side] * 3)
    wp = np.mod(pos + side / 2, box)                        # wrap into [0, side)
    # coordination number (periodic contacts within 1% of touching)
    Z = np.zeros(len(pos))
    for i in range(len(pos)):
        d = wp - wp[i]; d -= box * np.round(d / box)
        dist = np.linalg.norm(d, axis=1)
        Z[i] = ((dist > 1e-6) & (dist < (r + r[i]) * 1.01)).sum()
    # g(r): histogram of periodic pair distances, normalised by the ideal-gas count
    dd = []
    for i in range(len(pos)):
        d = wp[i + 1:] - wp[i]; d -= box * np.round(d / box)
        dd.append(np.linalg.norm(d, axis=1))
    dd = np.concatenate(dd)
    rmax = side / 2
    h, edges = np.histogram(dd, bins=120, range=(0, rmax))
    rc = 0.5 * (edges[1:] + edges[:-1])
    shell = 4 * np.pi * rc ** 2 * (edges[1] - edges[0])
    rho = len(pos) / side ** 3
    g = h / (shell * rho * len(pos) / 2)
    return Z, rc, g

Z, rc, g = microstructure(pos, r, side)
print(f"mean coordination Z = {Z[Z>0].mean():.2f} (isostatic 6)   rattlers = {int((Z==0).sum())}")

mean coordination Z = 6.27 (isostatic 6)   rattlers = 0

d = 2 * r.mean()                                  # sphere diameter (effective)
fig, (a0, a1) = plt.subplots(1, 2, figsize=(9, 3.9))
zc = 0.0
for k in range(len(pos)):
    dz = pos[k, 2] - zc; dz -= side * round(dz / side)
    if abs(dz) < r[k]:
        rr = math.sqrt(max(r[k] ** 2 - dz ** 2, 0))
        a0.add_patch(plt.Circle((pos[k, 0], pos[k, 1]), rr, color=BLUE, alpha=0.55))
a0.set(xlim=(-side / 2, side / 2), ylim=(-side / 2, side / 2), aspect="equal",
       title="packing slice (z=0)", xlabel="x", ylabel="y"); a0.grid(False)
a1.plot(rc / d, g, "-", color=RED)
a1.axhline(1, ls=":", color="0.6")
a1.set(xlabel="r / d", ylabel="g(r)", title="radial distribution function", xlim=(0, rc.max() / d))
fig.tight_layout()
plt.show()

Figure 1: Left: a slice through the packing (spheres intersecting the mid-plane). Right: the radial distribution function g(r/d), with the contact peak at one diameter and the split second peak characteristic of random close packing.

Step 3 — Pore-space flow and permeability with `peclet.flow`

Sampling the packing onto a grid gives the SDF (negative inside a sphere). A body force drives Stokes flow through the pores; the Darcy permeability is $k = \mu\langle u\rangle/F$ (here in physical length units via the cell size).

def sdf_from_pack(Ng, pos, r, side):
    g = (np.arange(Ng) + 0.5) / Ng * side
    X, Y, Z = np.meshgrid(g, g, g, indexing="ij")
    best = np.full((Ng, Ng, Ng), 1e30)
    for k in range(len(pos)):
        dx = X - (pos[k, 0] + side / 2); dx -= side * np.round(dx / side)
        dy = Y - (pos[k, 1] + side / 2); dy -= side * np.round(dy / side)
        dz = Z - (pos[k, 2] + side / 2); dz -= side * np.round(dz / side)
        best = np.minimum(best, np.sqrt(dx * dx + dy * dy + dz * dz) - r[k])
    return best

def permeability(Ng, pos, r, side, mu=0.1, F=1e-3, dt=80.0, max_steps=500, tol=1e-4, want_field=False):
    sdf = sdf_from_pack(Ng, pos, r, side)
    lv = max(2, int(np.log2(Ng)) - 1)
    s = flow.Solver(Ng, Ng, Ng)
    s.set_rho(1.0); s.set_mu(mu); s.set_dt(dt); s.set_body_force(F, 0, 0); s.set_advection(False)
    s.set_velocity_solver_params(150)
    s.set_pressure_multigrid(True, levels=lv)
    s.set_pressure_pcg(True, 200, 1e-8)
    s.set_solid(np.asfortranarray(sdf), cutcell_pressure=True, pressure_coarse="rediscretized")
    prev = 0.0
    for it in range(max_steps):
        s.step()
        if it % 5 == 4:
            m = float(s.get_u().mean())
            if it > 10 and abs(m - prev) < tol * (abs(m) + 1e-30):
                break
            prev = m
    umean = float(s.get_u().mean())
    k = mu * umean / F * (side / Ng) ** 2
    out = dict(Ng=Ng, k=k, umean=umean, steps=it + 1)
    if want_field:
        out["u"] = s.get_u()[:, :, Ng // 2].copy()
        out["sdf"] = sdf[:, :, Ng // 2].copy()
    return out

main = permeability(56, pos, r, side, want_field=True)
eps = 1 - phi
k_CK = eps ** 3 * (2 * r.mean()) ** 2 / (180 * (1 - eps) ** 2)    # Carman–Kozeny estimate (d = diameter)
print(f"Ng=56:  permeability k = {main['k']:.3e}   Carman–Kozeny ≈ {k_CK:.3e}   "
      f"ε={eps:.3f}  ({main['steps']} steps)")

Ng=56:  permeability k = 7.619e-04   Carman–Kozeny ≈ 7.090e-04   ε=0.370  (170 steps)

u = np.where(main["sdf"] < 0, np.nan, main["u"])
fig, ax = plt.subplots(figsize=(4.6, 4.2))
im = ax.imshow(u.T, origin="lower", cmap="magma")
ax.set(title="pore-space velocity (mid-plane)", xlabel="x", ylabel="y"); ax.grid(False)
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label="u")
plt.show()

Figure 2: Streamwise velocity on a mid-plane slice through the packed bed: flow threads through the connected pore space and stagnates against the grains (masked white).

Step 4 — Grid convergence (with a caveat)

At this porosity the permeability is already fairly grid-insensitive across these resolutions — the connected pores are well enough resolved that $k$ barely moves from $N=32$ to $56$. That is not guaranteed in general: denser beds have near-touching grains whose limiting throats are only a few cells wide, and there finer grids (a GPU) are needed. We show the trend against the Carman–Kozeny estimate.

conv = [permeability(Ng, pos, r, side) for Ng in (32, 44, 56)]
Ngs = [c["Ng"] for c in conv]; ks = [c["k"] for c in conv]

fig, ax = plt.subplots(figsize=(5.0, 3.8))
ax.plot(Ngs, ks, "o-", color=BLUE, label="peclet")
ax.axhline(k_CK, ls="--", color="0.4", label="Carman–Kozeny")
ax.set(xlabel="grid resolution N", ylabel="permeability k", title="permeability vs resolution")
ax.legend()
plt.show()
for c in conv:
    print(f"  N={c['Ng']:3d}  k={c['k']:.3e}  ({c['steps']} steps)")

Figure 3: Permeability vs grid resolution for one packing, with the Carman–Kozeny estimate (dashed). The tight throats of a random packing make convergence demanding — contrast the smooth Zick–Homsy lattice.

  N= 32  k=7.262e-04  (315 steps)
  N= 44  k=7.649e-04  (215 steps)
  N= 56  k=7.619e-04  (170 steps)

Step 5 — Statistics over independent realizations

Permeability is a property of the random microstructure, so it fluctuates between independent packings. We generate several realizations (different random seeds), characterise each, and report the spread.

real = [dict(seed=3, phi=phi, k=conv[1]["k"], Z=Z[Z > 0].mean())]   # reuse the main packing (seed 3, Ng=44)
for sd in (4, 5):
    p2, r2, s2, ph2 = pack_bed(seed=sd)
    Z2, _, _ = microstructure(p2, r2, s2)
    k2 = permeability(44, p2, r2, s2)["k"]
    real.append(dict(seed=sd, phi=ph2, k=k2, Z=Z2[Z2 > 0].mean()))

ks = np.array([x["k"] for x in real]); phis = np.array([x["phi"] for x in real])
fig, (a0, a1) = plt.subplots(1, 2, figsize=(9, 3.6))
a0.bar([x["seed"] for x in real], ks, color=BLUE, alpha=0.8)
a0.axhline(ks.mean(), ls="--", color=RED, label=f"mean {ks.mean():.2e}")
a0.set(xlabel="realization (seed)", ylabel="permeability k", title="permeability spread"); a0.legend(fontsize=8)
a1.scatter(phis, ks, color=GREEN, s=40)
a1.set(xlabel="solid fraction φ", ylabel="permeability k", title="k vs φ (per realization)")
fig.tight_layout()
plt.show()
print(f"permeability over {len(real)} realizations: mean={ks.mean():.3e}  std={ks.std():.2e}  "
      f"({100*ks.std()/ks.mean():.1f}%)   mean φ={phis.mean():.3f}  mean Z={np.mean([x['Z'] for x in real]):.2f}")

Figure 4: Permeability and porosity across independent random packings. The permeability scatters by a few percent around the mean at fixed target density.

permeability over 3 realizations: mean=7.481e-04  std=1.70e-05  (2.3%)   mean φ=0.630  mean Z=6.20

The takeaway

Two peclet codes compose cleanly: dem builds the geometry, flow solves the transport, joined by an SDF — the same signed-distance representation used throughout the suite.
The generated packing shows the hallmarks of a random close packing: $\phi\approx0.62$, coordination $Z\approx6$–$7$ (at the frictionless isostatic value, no rattlers), and a $g(r)$ that is zero inside one diameter with a sharp contact peak and the split second peak of dense random packing.
Permeability at this porosity is well-behaved (grid-insensitive across the resolutions shown) and a few percent above the Carman–Kozeny estimate. Denser beds, with near-touching grains and few-cell throats, are the genuinely hard case where fine grids (a GPU) matter. For the validated second-order convergence methodology on a smooth geometry, see Zick–Homsy.

Adapt this yourself

Polydisperse beds. Give the spheres a size distribution (per-particle radius) and watch $\phi$ and permeability shift.
Wall-bounded columns. Replace the periodic box with confining walls to study wall channelling.
Bigger, on a GPU. The identical code runs on the CUDA/HIP peclet build; there the grid can be fine enough for a grid-converged permeability across many realizations.

Reproduce this

pip install peclet
quarto render examples/random-packed-bed/index.qmd --execute

--- title: "Flow through a randomly packed bed (peclet.dem → peclet.flow)" subtitle: "Generate a random close packing with the DEM engine, characterise its microstructure, then push Stokes flow through the pore space and measure permeability." author: "Peclet" date: "2026-07-03" categories: [dem, flow, porous-media, multi-method, permeability] jupyter: python3 --- [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/computational-chemical-engineering/peclet-examples/blob/main/examples/random-packed-bed/index.ipynb){target="_blank"}  Runs on a free Colab CPU runtime — the first cell installs `peclet` from PyPI. ## What you'll learn The gallery's first **cross-method** example: two `peclet` codes in one pipeline. We use **`peclet.dem`** (discrete-element method) to build a random close packing of spheres, characterise its microstructure (porosity, coordination number, radial distribution), convert it to a **signed distance function**, and drive creeping flow through the pore space with **`peclet.flow`** to measure the **Darcy permeability** — compared against the classic Carman–Kozeny relation. This is the realistic, disordered counterpart to the periodic-lattice [Zick–Homsy](../zick-homsy/index.qmd) example. ```{python} #| label: bootstrap #| code-summary: "Environment bootstrap (installs peclet from PyPI on Colab/Binder)" import importlib.util, os, subprocess, sys _local = os.environ.get("PECLET_LOCAL_BUILD") if _local: for p in _local.split(os.pathsep): sys.path.insert(0, p) elif importlib.util.find_spec("peclet") is None: subprocess.run([sys.executable, "-m", "pip", "install", "-q", "peclet"], check=True) ``` ```{python} import numpy as np import math, time import matplotlib.pyplot as plt from peclet import dem, flow plt.rcParams.update({"figure.dpi": 130, "font.size": 10, "axes.grid": True, "grid.alpha": 0.3, "axes.axisbelow": True, "figure.facecolor": "white", "savefig.bbox": "tight"}) BLUE, RED, GREEN = "#1f77b4", "#d62728", "#2ca02c" ``` ## Step 1 — Build a random close packing with `peclet.dem` We grow monodisperse spheres in a periodic box under a **Lubachevsky–Stillinger** protocol: the particles start small, grow while a thermostat keeps the assembly fluid, and the growth is feedback-controlled by the committed overlap so the system jams near random close packing. A final quench settles it into a rigid contact network. ```{python} def pack_bed(N=180, phi_ref=0.63, radius=0.5, seed=3, dt=0.002, limit_time=7.0, iters=60, temperature=1.0, scale_init=0.05, criterion=5e-3, cooling_time=5.0, quench=1200, growth_accel=1.02, growth_decay=0.85, growth_rate_init=0.5): # box sized so that full-scale spheres would give phi_ref; the growth stops *before* that, # at the jamming point, so the achieved density (and radii) come out below phi_ref. volp = (4 / 3) * math.pi * radius ** 3 side = (N * volp / phi_ref) ** (1 / 3) half = side / 2.0 gr = growth_rate_init cool = min(int(cooling_time / dt), int(limit_time / dt)) rng = np.random.default_rng(seed) s = dem.Simulation(N) s.initialize(shape_type=1, radius=radius) s.set_domain((-half, -half, -half), (half, half, half)) s.enable_periodicity(True, True, True); s.set_gravity(0, 0, 0) s.set_material_params(1.0, 1.0, 0.0); s.set_solver_iterations(iters, iters) # restitution 1 while fluid pos = rng.uniform(-half, half, (N, 4)).astype(np.float32); pos[:, 3] = 1.0 s.set_positions(pos) s.set_velocities(rng.normal(0.0, math.sqrt(temperature), (N, 3)).astype(np.float32)) s.set_scales(np.full(N, 1.0, np.float32)) s.set_growth_params(gr, scale_init); s.set_thermostat(temperature, dt) def overlap_frac(): c = 2 * radius * float(s.get_scales().ravel().mean()) return float(s.compute_overlaps()) / max(c, 1e-9) # MONOTONIC growth. When overlaps build, we RELAX them at the current size (dt = 0 steps grow # nothing) and only slow the growth rate — we never shrink the grains back. The old scheme # un-grew on every overlap spike, and because the GPU contact solve is atomically non-deterministic # (float atomic_add ordering) a transient spike could trigger a runaway of shrink-steps that # quenched a loose, unjammed state on some seeds. Growing in one direction only makes the final # density reproducible: every seed lands at the box-limited phi_ref with an isostatic Z ~ 6. for step in range(int(limit_time / dt)): if step == cool: # anneal: mildly dissipative s.set_material_params(0.5, 1.0, 0.0); s.set_thermostat(0.0, 1e4 * dt) s.step(dt); mo = overlap_frac() gf = float(s.get_growth_factor()) if mo > criterion: # overlaps too big -> relax at CURRENT size it = 0; prev = mo while it < 40: s.step(0.0); it += 1; mn = overlap_frac() # dt = 0: pure overlap removal, no growth if mn < criterion or (it > 8 and mn > 0.98 * prev): break # cleared, or plateaued (jammed at this size) prev = mn if overlap_frac() > criterion: # still jammed -> slow growth, never un-grow gr = max(gr * growth_decay, 0.02); s.set_growth_params(gr, gf) elif gf >= 1.0 - 1e-6: break # reached full size (box sized for phi_ref) else: # room to grow -> accelerate gr = min(gr * growth_accel, growth_rate_init); s.set_growth_params(gr, gf) s.set_material_params(0.0, 0.0, 0.0); s.set_thermostat(0.0, 10 * dt) # quench into the rigid network for _ in range(quench): s.step(dt) # EFFECTIVE radius = per-particle scale * global growth factor * base radius. The growth factor # (<= 1) is essential — omitting it overstates the radii and fakes overlaps in g(r). r_eff = radius * s.get_scales().ravel() * float(s.get_growth_factor()) pos = s.get_positions()[:, :3].astype(float) phi = float(np.sum(4 / 3 * np.pi * r_eff ** 3) / side ** 3) return pos, r_eff.astype(float), float(side), phi t0 = time.time() pos, r, side, phi = pack_bed() print(f"random close packing: N={len(pos)} φ={phi:.3f} porosity ε={1-phi:.3f} " f"mean radius={r.mean():.4f} box={side:.2f} ({time.time()-t0:.0f}s)") ``` ## Step 2 — Characterise the microstructure Porosity $\varepsilon = 1-\phi$, the mean **coordination number** $Z$ (contacts per sphere — near the isostatic value 6 for a frictionless jammed packing), and the **radial distribution function** $g(r)$ with its sharp contact peak at one diameter. ```{python} def microstructure(pos, r, side): box = np.array([side] * 3) wp = np.mod(pos + side / 2, box) # wrap into [0, side) # coordination number (periodic contacts within 1% of touching) Z = np.zeros(len(pos)) for i in range(len(pos)): d = wp - wp[i]; d -= box * np.round(d / box) dist = np.linalg.norm(d, axis=1) Z[i] = ((dist > 1e-6) & (dist < (r + r[i]) * 1.01)).sum() # g(r): histogram of periodic pair distances, normalised by the ideal-gas count dd = [] for i in range(len(pos)): d = wp[i + 1:] - wp[i]; d -= box * np.round(d / box) dd.append(np.linalg.norm(d, axis=1)) dd = np.concatenate(dd) rmax = side / 2 h, edges = np.histogram(dd, bins=120, range=(0, rmax)) rc = 0.5 * (edges[1:] + edges[:-1]) shell = 4 * np.pi * rc ** 2 * (edges[1] - edges[0]) rho = len(pos) / side ** 3 g = h / (shell * rho * len(pos) / 2) return Z, rc, g Z, rc, g = microstructure(pos, r, side) print(f"mean coordination Z = {Z[Z>0].mean():.2f} (isostatic 6) rattlers = {int((Z==0).sum())}") ``` ```{python} #| label: fig-microstructure #| fig-cap: "Left: a slice through the packing (spheres intersecting the mid-plane). Right: the radial distribution function g(r/d), with the contact peak at one diameter and the split second peak characteristic of random close packing." d = 2 * r.mean() # sphere diameter (effective) fig, (a0, a1) = plt.subplots(1, 2, figsize=(9, 3.9)) zc = 0.0 for k in range(len(pos)): dz = pos[k, 2] - zc; dz -= side * round(dz / side) if abs(dz) < r[k]: rr = math.sqrt(max(r[k] ** 2 - dz ** 2, 0)) a0.add_patch(plt.Circle((pos[k, 0], pos[k, 1]), rr, color=BLUE, alpha=0.55)) a0.set(xlim=(-side / 2, side / 2), ylim=(-side / 2, side / 2), aspect="equal", title="packing slice (z=0)", xlabel="x", ylabel="y"); a0.grid(False) a1.plot(rc / d, g, "-", color=RED) a1.axhline(1, ls=":", color="0.6") a1.set(xlabel="r / d", ylabel="g(r)", title="radial distribution function", xlim=(0, rc.max() / d)) fig.tight_layout() plt.show() ``` ## Step 3 — Pore-space flow and permeability with `peclet.flow` Sampling the packing onto a grid gives the SDF (negative inside a sphere). A body force drives Stokes flow through the pores; the Darcy permeability is $k = \mu\langle u\rangle/F$ (here in physical length units via the cell size). ```{python} def sdf_from_pack(Ng, pos, r, side): g = (np.arange(Ng) + 0.5) / Ng * side X, Y, Z = np.meshgrid(g, g, g, indexing="ij") best = np.full((Ng, Ng, Ng), 1e30) for k in range(len(pos)): dx = X - (pos[k, 0] + side / 2); dx -= side * np.round(dx / side) dy = Y - (pos[k, 1] + side / 2); dy -= side * np.round(dy / side) dz = Z - (pos[k, 2] + side / 2); dz -= side * np.round(dz / side) best = np.minimum(best, np.sqrt(dx * dx + dy * dy + dz * dz) - r[k]) return best def permeability(Ng, pos, r, side, mu=0.1, F=1e-3, dt=80.0, max_steps=500, tol=1e-4, want_field=False): sdf = sdf_from_pack(Ng, pos, r, side) lv = max(2, int(np.log2(Ng)) - 1) s = flow.Solver(Ng, Ng, Ng) s.set_rho(1.0); s.set_mu(mu); s.set_dt(dt); s.set_body_force(F, 0, 0); s.set_advection(False) s.set_velocity_solver_params(150) s.set_pressure_multigrid(True, levels=lv) s.set_pressure_pcg(True, 200, 1e-8) s.set_solid(np.asfortranarray(sdf), cutcell_pressure=True, pressure_coarse="rediscretized") prev = 0.0 for it in range(max_steps): s.step() if it % 5 == 4: m = float(s.get_u().mean()) if it > 10 and abs(m - prev) < tol * (abs(m) + 1e-30): break prev = m umean = float(s.get_u().mean()) k = mu * umean / F * (side / Ng) ** 2 out = dict(Ng=Ng, k=k, umean=umean, steps=it + 1) if want_field: out["u"] = s.get_u()[:, :, Ng // 2].copy() out["sdf"] = sdf[:, :, Ng // 2].copy() return out main = permeability(56, pos, r, side, want_field=True) eps = 1 - phi k_CK = eps ** 3 * (2 * r.mean()) ** 2 / (180 * (1 - eps) ** 2) # Carman–Kozeny estimate (d = diameter) print(f"Ng=56: permeability k = {main['k']:.3e} Carman–Kozeny ≈ {k_CK:.3e} " f"ε={eps:.3f} ({main['steps']} steps)") ``` ```{python} #| label: fig-flowfield #| fig-cap: "Streamwise velocity on a mid-plane slice through the packed bed: flow threads through the connected pore space and stagnates against the grains (masked white)." u = np.where(main["sdf"] < 0, np.nan, main["u"]) fig, ax = plt.subplots(figsize=(4.6, 4.2)) im = ax.imshow(u.T, origin="lower", cmap="magma") ax.set(title="pore-space velocity (mid-plane)", xlabel="x", ylabel="y"); ax.grid(False) fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label="u") plt.show() ``` ## Step 4 — Grid convergence (with a caveat) At this porosity the permeability is already fairly **grid-insensitive** across these resolutions — the connected pores are well enough resolved that $k$ barely moves from $N=32$ to $56$. That is *not* guaranteed in general: denser beds have near-touching grains whose limiting throats are only a few cells wide, and there finer grids (a GPU) are needed. We show the trend against the Carman–Kozeny estimate. ```{python} #| label: fig-convergence #| fig-cap: "Permeability vs grid resolution for one packing, with the Carman–Kozeny estimate (dashed). The tight throats of a random packing make convergence demanding — contrast the smooth Zick–Homsy lattice." conv = [permeability(Ng, pos, r, side) for Ng in (32, 44, 56)] Ngs = [c["Ng"] for c in conv]; ks = [c["k"] for c in conv] fig, ax = plt.subplots(figsize=(5.0, 3.8)) ax.plot(Ngs, ks, "o-", color=BLUE, label="peclet") ax.axhline(k_CK, ls="--", color="0.4", label="Carman–Kozeny") ax.set(xlabel="grid resolution N", ylabel="permeability k", title="permeability vs resolution") ax.legend() plt.show() for c in conv: print(f" N={c['Ng']:3d} k={c['k']:.3e} ({c['steps']} steps)") ``` ## Step 5 — Statistics over independent realizations Permeability is a property of the *random* microstructure, so it fluctuates between independent packings. We generate several realizations (different random seeds), characterise each, and report the spread. ```{python} #| label: fig-statistics #| fig-cap: "Permeability and porosity across independent random packings. The permeability scatters by a few percent around the mean at fixed target density." real = [dict(seed=3, phi=phi, k=conv[1]["k"], Z=Z[Z > 0].mean())] # reuse the main packing (seed 3, Ng=44) for sd in (4, 5): p2, r2, s2, ph2 = pack_bed(seed=sd) Z2, _, _ = microstructure(p2, r2, s2) k2 = permeability(44, p2, r2, s2)["k"] real.append(dict(seed=sd, phi=ph2, k=k2, Z=Z2[Z2 > 0].mean())) ks = np.array([x["k"] for x in real]); phis = np.array([x["phi"] for x in real]) fig, (a0, a1) = plt.subplots(1, 2, figsize=(9, 3.6)) a0.bar([x["seed"] for x in real], ks, color=BLUE, alpha=0.8) a0.axhline(ks.mean(), ls="--", color=RED, label=f"mean {ks.mean():.2e}") a0.set(xlabel="realization (seed)", ylabel="permeability k", title="permeability spread"); a0.legend(fontsize=8) a1.scatter(phis, ks, color=GREEN, s=40) a1.set(xlabel="solid fraction φ", ylabel="permeability k", title="k vs φ (per realization)") fig.tight_layout() plt.show() print(f"permeability over {len(real)} realizations: mean={ks.mean():.3e} std={ks.std():.2e} " f"({100*ks.std()/ks.mean():.1f}%) mean φ={phis.mean():.3f} mean Z={np.mean([x['Z'] for x in real]):.2f}") ``` ## The takeaway - Two `peclet` codes compose cleanly: `dem` builds the geometry, `flow` solves the transport, joined by an SDF — the same signed-distance representation used throughout the suite. - The generated packing shows the hallmarks of a **random close packing**: $\phi\approx0.62$, coordination $Z\approx6$–$7$ (at the frictionless isostatic value, no rattlers), and a $g(r)$ that is zero inside one diameter with a sharp contact peak and the split second peak of dense random packing. - Permeability at this porosity is well-behaved (grid-insensitive across the resolutions shown) and a few percent above the Carman–Kozeny estimate. *Denser* beds, with near-touching grains and few-cell throats, are the genuinely hard case where fine grids (a GPU) matter. For the *validated* second-order convergence methodology on a smooth geometry, see [Zick–Homsy](../zick-homsy/index.qmd). ## Adapt this yourself - **Polydisperse beds.** Give the spheres a size distribution (per-particle radius) and watch $\phi$ and permeability shift. - **Wall-bounded columns.** Replace the periodic box with confining walls to study wall channelling. - **Bigger, on a GPU.** The identical code runs on the CUDA/HIP `peclet` build; there the grid can be fine enough for a grid-converged permeability across many realizations. ## Reproduce this ```bash pip install peclet quarto render examples/random-packed-bed/index.qmd --execute ```