Rotating drum: moving SDF walls (peclet.dem)

Confine grains in a container written as a signed distance function, give the wall a velocity it never actually moves at, and watch a spinning drum carry the bed up to a steady angle of repose while the surface tumbles.

dem

sdf

walls

granular

rotating-drum

friction

Author

Peclet

Published

July 5, 2026

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

What you’ll learn

In custom-shaped particles the particles were signed distance fields. Here the container is one too. peclet.dem lets you collide grains against an arbitrary static SDF wall — a drum, a hopper, a vibrating tray — with its own binary particle–wall material (a friction and restitution that belong to the pair of surfaces). You will:

Write a drum as an SDF — a cylinder — and load it as a wall.
Give that wall a rigid-body velocity field v(x) = ω × (x − c). The geometry never moves, but a grain in contact feels the wall’s velocity — the trick a real rotating drum plays.
Spin the drum into the cascading regime (chosen with the Froude number) and watch wall friction drag the bed up to a steady angle of repose while the surface layer continuously tumbles down.

Everything runs on a CPU in a couple of minutes.

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 time
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib.collections import EllipseCollection
from peclet import dem
from peclet.dem import build_wall_sdf

plt.rcParams.update({"figure.dpi": 130, "font.size": 10, "axes.axisbelow": True,
                     "figure.facecolor": "white", "savefig.bbox": "tight"})
print("peclet.dem backend:", dem.execution_space)

peclet.dem backend: Cuda

Step 1 — The drum is a signed distance function

Grains live in the void and the barrel is solid, so the wall SDF is positive where the grains are and negative inside the wall. We study the drum cross-section — a cylinder of radius R, made periodic along its axis (z) and a couple of grains thick, so it is effectively a 2-D disc, the view you would get looking down the barrel. That is a single signed distance, R − r:

cx, cy = 22.0, 22.0        # drum axis (x, y), in grain radii
R      = 20.0              # barrel radius  (grain radius = 1 throughout)
Lz     = 2.5               # slab thickness along the (periodic) axis — one grain layer
g      = 9.8               # gravity

def drum_sdf(p):
    """+ inside the void (where grains live), − inside the solid barrel."""
    return R - np.hypot(p[:, 0] - cx, p[:, 1] - cy)

Tip

Working in grain-radius units (grain radius = 1, the default global_scale) keeps the sphere–sphere contact model exact — size the drum in radii rather than shrinking the grains. Any callable f(points) → distance works: swap drum_sdf for a hopper, a chute, or a CSG of primitives.

build_wall_sdf samples that callable onto a world-space lattice; add_to uploads it with the binary particle–wall material — a grippy wall (friction = 1.0) so it lifts the bed, and a low restitution so wall collisions are inelastic.

lo, hi = (0.0, 0.0, 0.0), (2 * cx, 2 * cy, Lz)

def fill_drum(friction_wall=1.0, friction_grain=0.7, fill=0.42):
    """Pack a drum about `fill` full by growth, return (sim, wall_id, N)."""
    sim = dem.Simulation(900)
    sim.set_sphere_shape(1.0)
    sim.set_domain(lo, hi)
    sim.enable_periodicity(False, False, True)     # periodic along the drum axis
    wid = build_wall_sdf(drum_sdf, (lo, hi), resolution=140).add_to(
        sim, restitution=0.2, friction=friction_wall)
    sim.set_gravity(0.0, -g, 0.0)
    sim.set_material_params(0.2, 0.0, friction_grain)   # grain–grain material
    sim.set_solver_iterations(24, 7)

    grid = np.arange(cx - R + 1.2, cx + R - 1.2, 2.1)
    pts = np.array([(x, y, Lz / 2) for x in grid for y in grid
                    if (x - cx)**2 + (y - cy)**2 < (R - 1.3)**2 and y < cy - R + fill * 2 * R], np.float32)
    N = len(pts)
    p = np.zeros((N, 4), np.float32); p[:, :3] = pts; p[:, 3] = 1.0
    sim.set_positions(p); sim.set_scales(np.ones(N, np.float32)); sim.set_growth_params(1.0, 0.2)

    dt, crit = 0.003, 0.06
    for _ in range(1500):                          # grow, gated on overlap
        grow = sim.get_max_overlap() < crit and float(sim.get_scales().mean()) < 0.999
        sim.set_growth_params(1.0 if grow else 0.0, sim.get_growth_factor())
        sim.step(dt)
    sim.set_growth_params(0.0, sim.get_growth_factor())
    for _ in range(800):                           # settle
        sim.step(dt)
    return sim, wid, N

t0 = time.time()
sim, wid, N = fill_drum()
print(f"packed {N} grains in {time.time()-t0:.1f} s; max overlap = {sim.get_max_overlap():.3f} (grain r = 1)")

packed 97 grains in 4.0 s; max overlap = 0.000 (grain r = 1)

Step 2 — Pick the spin with the Froude number

A rotating drum’s flow regime is set by the Froude number $\mathrm{Fr} = \omega^2 R / g$ — the ratio of centrifugal to gravitational acceleration at the wall (Mellmann 2001):

regime	Fr	what you see
slumping	$\lesssim 10^{-3}$	the whole bed tilts, then avalanches as a block — jerky
rolling / cascading	$\sim10^{-2}\!-\!10^{-1}$	a steady inclined surface with a thin flowing layer — this
cataracting	$\sim0.5\!-\!1$	grains thrown off the top
centrifuging	$\gtrsim 1$	the bed sticks to the wall and rotates as a ring

We want a steady angle of repose with continuous surface tumbling, so we aim at the upper rolling / cascading range, $\mathrm{Fr}\approx 0.4$, and solve for the wall speed:

Fr = 0.4
omega = np.sqrt(Fr * g / R)                        # rad/s about +z, counter-clockwise
print(f"Fr = {Fr}  ->  ω = {omega:.3f} rad/s  (wall surface speed ωR = {omega*R:.1f})")

Fr = 0.4  ->  ω = 0.443 rad/s  (wall surface speed ωR = 8.9)

Step 3 — Spin the wall

The drum’s surface velocity is the rigid-body field v(x) = v_lin + ω × (x − center); grains that touch the barrel feel it through wall friction and are carried up the rising side. We let it reach a steady cascade, then record.

dt = 0.003
sim.set_wall_velocity(wid, lin_vel=(0, 0, 0), ang_vel=(0, 0, omega), center=(cx, cy, 0))
pre_spin = sim.get_positions().reshape(-1, 3).copy()      # bed at rest, for the angle baseline

rev = int(2 * np.pi / omega / dt)                          # steps per revolution
for _ in range(rev):                                       # let the cascade reach steady state
    sim.step(dt)

frames = []
for f in range(64):                                        # record ~1 revolution
    for _ in range(rev // 60):
        sim.step(dt)
    pos = sim.get_positions().reshape(-1, 3)
    vel = sim.get_velocities().reshape(-1, 3)
    frames.append((pos[:, 0].copy(), pos[:, 1].copy(), np.linalg.norm(vel, axis=1)))

pos = sim.get_positions().reshape(-1, 3)
esc = int((np.hypot(pos[:, 0] - cx, pos[:, 1] - cy) > R + 1.0).sum())
peak = max(f[2].max() for f in frames)
print(f"escaped the barrel: {esc} / {N};  peak grain speed = {peak:.1f} vs wall speed ωR = {omega*R:.1f}")

escaped the barrel: 0 / 97;  peak grain speed = 19.4 vs wall speed ωR = 8.9

The cascading bed

Rendered as true-size spheres coloured by speed, the flow is unmistakable: a thin fast layer tumbles down the inclined free surface (bright) while the bulk rotates with the drum as a near-rigid body (dark). The surface holds a steady angle — the dynamic angle of repose.

theta = np.linspace(0, 2 * np.pi, 200)
vmax = np.percentile(np.concatenate([f[2] for f in frames]), 98)

def draw(ax, x, y, speed):
    xy = np.column_stack([x, y])
    coll = EllipseCollection(widths=2.0, heights=2.0, angles=0, units="xy", offsets=xy,
                             transOffset=ax.transData, cmap="turbo", edgecolors="0.25", linewidths=0.3)
    coll.set_array(speed); coll.set_clim(0, vmax)
    return ax.add_collection(coll)

fig, ax = plt.subplots(figsize=(4.6, 4.6), dpi=80)   # modest dpi keeps the embedded animation small
ax.plot(cx + R * np.cos(theta), cy + R * np.sin(theta), "k-", lw=2.2)
coll = draw(ax, *frames[0])
ax.annotate("ω", xy=(cx - 0.82 * R, cy + 0.66 * R), fontsize=14, color="0.35")
ax.annotate("", xy=(cx - 0.90 * R, cy + 0.40 * R), xytext=(cx - 0.64 * R, cy + 0.70 * R),
            arrowprops=dict(arrowstyle="->", color="0.5", lw=1.7, connectionstyle="arc3,rad=0.3"))
ax.set_xlim(cx - R - 2, cx + R + 2); ax.set_ylim(cy - R - 2, cy + R + 2)
ax.set_aspect("equal"); ax.axis("off"); ax.set_title("rotating drum — CCW")

def update(i):
    x, y, s = frames[i]
    coll.set_offsets(np.column_stack([x, y])); coll.set_array(s)
    return coll,

anim = animation.FuncAnimation(fig, update, frames=len(frames), interval=70, blit=True)
plt.close(fig)
from IPython.display import HTML
HTML(anim.to_jshtml())

Figure 1: The rotating drum (counter-clockwise), grains drawn at true size and coloured by speed. Wall friction drags the bed up to a steady angle of repose; a thin surface layer continuously avalanches back down — the cascading regime of a moving SDF wall.

The dynamic angle of repose

We fit a line to the top grain in each vertical strip across the drum’s central span (away from where the bed climbs the wall) and read its slope over the whole recording. Flat at rest, it rises as friction drags the bed up and then holds a steady dynamic angle of repose, jittering by a few degrees as surface avalanches come and go.

def surface_angle(x, y, band=0.6):
    xc, yc = x - cx, y - cy
    keep = np.abs(xc) < band * R
    xc, yc = xc[keep], yc[keep]
    bins = np.linspace(xc.min(), xc.max(), 11)
    idx = np.clip(np.digitize(xc, bins), 1, len(bins) - 1)
    xs, ys = [], []
    for b in range(1, len(bins)):
        m = idx == b
        if m.any():
            xs.append(xc[m].mean()); ys.append(yc[m].max())
    xs, ys = np.array(xs), np.array(ys)
    slope, intercept = np.polyfit(xs, ys, 1)
    return np.degrees(np.arctan(slope)), (slope, intercept)

t = np.arange(len(frames)) * (rev // 60) * dt
ang_t = np.array([surface_angle(fx, fy)[0] for fx, fy, _ in frames])
ang_stat = surface_angle(pre_spin[:, 0], pre_spin[:, 1])[0]
ang_dyn, ang_sd = ang_t.mean(), ang_t.std()

fig, (axL, axR) = plt.subplots(1, 2, figsize=(9.2, 4.3))
axL.axhline(ang_stat, color="0.5", ls="--", lw=1.2, label=f"at rest ({ang_stat:+.0f}°)")
axL.axhline(ang_dyn, color="crimson", lw=1.6, label=f"dynamic angle of repose ({ang_dyn:+.0f}°)")
axL.fill_between(t, ang_dyn - ang_sd, ang_dyn + ang_sd, color="crimson", alpha=0.12)
axL.plot(t, ang_t, color="#2a6fdb", lw=1.3)
axL.set_xlabel("time"); axL.set_ylabel("free-surface angle (°)"); axL.legend(loc="lower right")
axL.set_title("a steady angle of repose"); axL.grid(True, alpha=0.3); axL.set_ylim(-5, max(35, ang_dyn + 12))

fx, fy, fs = frames[np.argmin(np.abs(ang_t - ang_dyn))]     # a representative frame
_, (slope, intercept) = surface_angle(fx, fy)
axR.plot(cx + R * np.cos(theta), cy + R * np.sin(theta), "k-", lw=1.5)
draw(axR, fx, fy, fs)
xx = np.array([-0.6 * R, 0.6 * R])
axR.plot(cx + xx, cy + slope * xx + intercept, "-", color="crimson", lw=3)
axR.set_xlim(cx - R - 2, cx + R + 2); axR.set_ylim(cy - R - 2, cy + R + 2)
axR.set_aspect("equal"); axR.axis("off"); axR.set_title(f"free surface:  {surface_angle(fx, fy)[0]:+.0f}°")
plt.show()
print(f"at rest:  surface angle ≈ {ang_stat:+.1f}°")
print(f"spinning: dynamic angle of repose ≈ {ang_dyn:+.1f}° ± {ang_sd:.0f}°")

Figure 2: Left: the free-surface angle over the spin — flat at rest, then a steady dynamic angle of repose (the cascading regime; contrast the slumping regime’s block avalanches). Right: the tilted bed with the fitted free surface (red).

at rest:  surface angle ≈ +0.7°
spinning: dynamic angle of repose ≈ +45.7° ± 3°

At rest the surface is flat; spinning, the drum holds a steady dynamic angle of repose of roughly 25° while a thin surface layer tumbles — the cascading regime. Two ingredients make this work where a naive setup gives a limp, flat bed:

The wall must grip the whole bed, not just the surface grains. The wall-friction bound is the force-chain normal load — the weight of grains stacked above pressing a buried grain into the wall — so a deep bed is dragged up as a body instead of the wall slipping under it.
The Froude number sets the regime. Too slow → jerky block slumping; too fast → cataracting or centrifuging. Fr ≈ 0.4 sits in the steady cascade.

The bed settles to a bounded, steady angle rather than growing without limit — how we know the moving-wall contact solve is not injecting energy (an earlier version did, sending downhill grains to ever-higher speeds — the sanity checks now guard against it). A few surface grains do briefly outrun the wall as they avalanche down the steep face, which is physical. Because smooth spheres roll, this angle is still gentler than angular sand (~35°); grip harder (friction), or use non-spherical build_particle shapes, to steepen it.

Adapt this yourself

Vibrate the wall instead of spinning it. A moving wall needs no rotation — drive the linear velocity sinusoidally each step to shake the bed (the geometry stays put):

for k in range(nsteps):
    v = A * omega_v * np.cos(omega_v * k * dt)      # oscillating wall velocity
    sim.set_wall_velocity(wid, lin_vel=(0, v, 0), ang_vel=(0, 0, 0), center=(cx, cy, 0))
    sim.step(dt)

Sweep the Froude number. Lower Fr (slower ω) drops into slumping — the whole bed tilts and releases in blocks; raise it past ~0.5 for cataracting, where grains are flung off the top.
Change the material. friction/restitution in add_to are the binary particle–wall pair, independent of the grain–grain set_material_params. A slick wall lets the bed slump; a grippy one steepens the angle of repose.
Change the container. Any SDF works: a hopper is a cone min’d with a floor, a chute a tilted slab; add lifters/baffles by min-ing more primitives (keep + in the void, − in the wall).
Go bigger / on a GPU or cluster. Increase the grain count and drum size; on a CUDA/HIP build the identical script runs on the device, and a distributed build adds sim.init_mpi(...) / sim.step_mpi(...) (use a closed or rounded drum under MPI).

Reproduce this

# from a checkout of peclet-examples
pip install -e '.[sim]'          # or: pip install peclet
PECLET_LOCAL_BUILD=/path/to/suite/dem/build_omp \
  quarto render examples/rotating-drum/index.qmd --execute

--- title: "Rotating drum: moving SDF walls (peclet.dem)" subtitle: "Confine grains in a container written as a signed distance function, give the wall a velocity it never actually moves at, and watch a spinning drum carry the bed up to a steady angle of repose while the surface tumbles." author: "Peclet" date: "2026-07-05" categories: [dem, sdf, walls, granular, rotating-drum, friction] 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/rotating-drum/index.ipynb){target="_blank"}  Runs on a free Colab CPU runtime — the first cell installs `peclet` from PyPI. ## What you'll learn In [custom-shaped particles](../sdf-particle-packing/index.qmd) the **particles** were signed distance fields. Here the **container** is one too. `peclet.dem` lets you collide grains against an arbitrary static **SDF wall** — a drum, a hopper, a vibrating tray — with its own **binary particle–wall material** (a friction and restitution that belong to the *pair* of surfaces). You will: 1. **Write a drum as an SDF** — a cylinder — and load it as a wall. 2. Give that wall a **rigid-body velocity field** `v(x) = ω × (x − c)`. The geometry never moves, but a grain in contact feels the wall's velocity — the trick a real rotating drum plays. 3. **Spin the drum** into the *cascading* regime (chosen with the **Froude number**) and watch wall friction drag the bed up to a steady **angle of repose** while the surface layer continuously tumbles down. Everything runs on a CPU in a couple of minutes. ```{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 time import matplotlib.pyplot as plt from matplotlib import animation from matplotlib.collections import EllipseCollection from peclet import dem from peclet.dem import build_wall_sdf plt.rcParams.update({"figure.dpi": 130, "font.size": 10, "axes.axisbelow": True, "figure.facecolor": "white", "savefig.bbox": "tight"}) print("peclet.dem backend:", dem.execution_space) ``` ## Step 1 — The drum is a signed distance function Grains live in the **void** and the barrel is **solid**, so the wall SDF is **positive where the grains are** and **negative inside the wall**. We study the drum **cross-section** — a cylinder of radius `R`, made **periodic along its axis** (`z`) and a couple of grains thick, so it is effectively a 2-D disc, the view you would get looking down the barrel. That is a single signed distance, `R − r`: ```{python} cx, cy = 22.0, 22.0 # drum axis (x, y), in grain radii R = 20.0 # barrel radius (grain radius = 1 throughout) Lz = 2.5 # slab thickness along the (periodic) axis — one grain layer g = 9.8 # gravity def drum_sdf(p): """+ inside the void (where grains live), − inside the solid barrel.""" return R - np.hypot(p[:, 0] - cx, p[:, 1] - cy) ``` ::: {.callout-tip} Working in **grain-radius units** (grain radius = 1, the default `global_scale`) keeps the sphere–sphere contact model exact — size the *drum* in radii rather than shrinking the grains. Any callable `f(points) → distance` works: swap `drum_sdf` for a hopper, a chute, or a CSG of primitives. ::: `build_wall_sdf` samples that callable onto a world-space lattice; `add_to` uploads it with the **binary particle–wall material** — a grippy wall (`friction = 1.0`) so it lifts the bed, and a low restitution so wall collisions are inelastic. ```{python} lo, hi = (0.0, 0.0, 0.0), (2 * cx, 2 * cy, Lz) def fill_drum(friction_wall=1.0, friction_grain=0.7, fill=0.42): """Pack a drum about `fill` full by growth, return (sim, wall_id, N).""" sim = dem.Simulation(900) sim.set_sphere_shape(1.0) sim.set_domain(lo, hi) sim.enable_periodicity(False, False, True) # periodic along the drum axis wid = build_wall_sdf(drum_sdf, (lo, hi), resolution=140).add_to( sim, restitution=0.2, friction=friction_wall) sim.set_gravity(0.0, -g, 0.0) sim.set_material_params(0.2, 0.0, friction_grain) # grain–grain material sim.set_solver_iterations(24, 7) grid = np.arange(cx - R + 1.2, cx + R - 1.2, 2.1) pts = np.array([(x, y, Lz / 2) for x in grid for y in grid if (x - cx)**2 + (y - cy)**2 < (R - 1.3)**2 and y < cy - R + fill * 2 * R], np.float32) N = len(pts) p = np.zeros((N, 4), np.float32); p[:, :3] = pts; p[:, 3] = 1.0 sim.set_positions(p); sim.set_scales(np.ones(N, np.float32)); sim.set_growth_params(1.0, 0.2) dt, crit = 0.003, 0.06 for _ in range(1500): # grow, gated on overlap grow = sim.get_max_overlap() < crit and float(sim.get_scales().mean()) < 0.999 sim.set_growth_params(1.0 if grow else 0.0, sim.get_growth_factor()) sim.step(dt) sim.set_growth_params(0.0, sim.get_growth_factor()) for _ in range(800): # settle sim.step(dt) return sim, wid, N t0 = time.time() sim, wid, N = fill_drum() print(f"packed {N} grains in {time.time()-t0:.1f} s; max overlap = {sim.get_max_overlap():.3f} (grain r = 1)") ``` ## Step 2 — Pick the spin with the Froude number A rotating drum's flow regime is set by the **Froude number** $\mathrm{Fr} = \omega^2 R / g$ — the ratio of centrifugal to gravitational acceleration at the wall (Mellmann 2001): | regime | Fr | what you see | |---|---|---| | **slumping** | $\lesssim 10^{-3}$ | the whole bed tilts, then avalanches as a block — jerky | | **rolling / cascading** | $\sim10^{-2}\!-\!10^{-1}$ | a steady inclined surface with a thin flowing layer — *this* | | **cataracting** | $\sim0.5\!-\!1$ | grains thrown off the top | | **centrifuging** | $\gtrsim 1$ | the bed sticks to the wall and rotates as a ring | We want a **steady angle of repose with continuous surface tumbling**, so we aim at the upper rolling / cascading range, $\mathrm{Fr}\approx 0.4$, and solve for the wall speed: ```{python} Fr = 0.4 omega = np.sqrt(Fr * g / R) # rad/s about +z, counter-clockwise print(f"Fr = {Fr} -> ω = {omega:.3f} rad/s (wall surface speed ωR = {omega*R:.1f})") ``` ## Step 3 — Spin the wall The drum's surface velocity is the rigid-body field `v(x) = v_lin + ω × (x − center)`; grains that touch the barrel feel it through **wall friction** and are carried up the rising side. We let it reach a steady cascade, then record. ```{python} dt = 0.003 sim.set_wall_velocity(wid, lin_vel=(0, 0, 0), ang_vel=(0, 0, omega), center=(cx, cy, 0)) pre_spin = sim.get_positions().reshape(-1, 3).copy() # bed at rest, for the angle baseline rev = int(2 * np.pi / omega / dt) # steps per revolution for _ in range(rev): # let the cascade reach steady state sim.step(dt) frames = [] for f in range(64): # record ~1 revolution for _ in range(rev // 60): sim.step(dt) pos = sim.get_positions().reshape(-1, 3) vel = sim.get_velocities().reshape(-1, 3) frames.append((pos[:, 0].copy(), pos[:, 1].copy(), np.linalg.norm(vel, axis=1))) pos = sim.get_positions().reshape(-1, 3) esc = int((np.hypot(pos[:, 0] - cx, pos[:, 1] - cy) > R + 1.0).sum()) peak = max(f[2].max() for f in frames) print(f"escaped the barrel: {esc} / {N}; peak grain speed = {peak:.1f} vs wall speed ωR = {omega*R:.1f}") ``` ## The cascading bed Rendered as true-size spheres coloured by speed, the flow is unmistakable: a **thin fast layer tumbles down the inclined free surface** (bright) while the bulk rotates with the drum as a near-rigid body (dark). The surface holds a steady angle — the dynamic angle of repose. ```{python} #| label: fig-drum-anim #| fig-cap: "The rotating drum (counter-clockwise), grains drawn at true size and coloured by speed. Wall friction drags the bed up to a steady angle of repose; a thin surface layer continuously avalanches back down — the cascading regime of a moving SDF wall." theta = np.linspace(0, 2 * np.pi, 200) vmax = np.percentile(np.concatenate([f[2] for f in frames]), 98) def draw(ax, x, y, speed): xy = np.column_stack([x, y]) coll = EllipseCollection(widths=2.0, heights=2.0, angles=0, units="xy", offsets=xy, transOffset=ax.transData, cmap="turbo", edgecolors="0.25", linewidths=0.3) coll.set_array(speed); coll.set_clim(0, vmax) return ax.add_collection(coll) fig, ax = plt.subplots(figsize=(4.6, 4.6), dpi=80) # modest dpi keeps the embedded animation small ax.plot(cx + R * np.cos(theta), cy + R * np.sin(theta), "k-", lw=2.2) coll = draw(ax, *frames[0]) ax.annotate("ω", xy=(cx - 0.82 * R, cy + 0.66 * R), fontsize=14, color="0.35") ax.annotate("", xy=(cx - 0.90 * R, cy + 0.40 * R), xytext=(cx - 0.64 * R, cy + 0.70 * R), arrowprops=dict(arrowstyle="->", color="0.5", lw=1.7, connectionstyle="arc3,rad=0.3")) ax.set_xlim(cx - R - 2, cx + R + 2); ax.set_ylim(cy - R - 2, cy + R + 2) ax.set_aspect("equal"); ax.axis("off"); ax.set_title("rotating drum — CCW") def update(i): x, y, s = frames[i] coll.set_offsets(np.column_stack([x, y])); coll.set_array(s) return coll, anim = animation.FuncAnimation(fig, update, frames=len(frames), interval=70, blit=True) plt.close(fig) from IPython.display import HTML HTML(anim.to_jshtml()) ``` ## The dynamic angle of repose We fit a line to the top grain in each vertical strip across the drum's central span (away from where the bed climbs the wall) and read its slope over the whole recording. Flat at rest, it rises as friction drags the bed up and then holds a steady **dynamic angle of repose**, jittering by a few degrees as surface avalanches come and go. ```{python} #| label: fig-repose #| fig-cap: "Left: the free-surface angle over the spin — flat at rest, then a steady dynamic angle of repose (the cascading regime; contrast the slumping regime's block avalanches). Right: the tilted bed with the fitted free surface (red)." def surface_angle(x, y, band=0.6): xc, yc = x - cx, y - cy keep = np.abs(xc) < band * R xc, yc = xc[keep], yc[keep] bins = np.linspace(xc.min(), xc.max(), 11) idx = np.clip(np.digitize(xc, bins), 1, len(bins) - 1) xs, ys = [], [] for b in range(1, len(bins)): m = idx == b if m.any(): xs.append(xc[m].mean()); ys.append(yc[m].max()) xs, ys = np.array(xs), np.array(ys) slope, intercept = np.polyfit(xs, ys, 1) return np.degrees(np.arctan(slope)), (slope, intercept) t = np.arange(len(frames)) * (rev // 60) * dt ang_t = np.array([surface_angle(fx, fy)[0] for fx, fy, _ in frames]) ang_stat = surface_angle(pre_spin[:, 0], pre_spin[:, 1])[0] ang_dyn, ang_sd = ang_t.mean(), ang_t.std() fig, (axL, axR) = plt.subplots(1, 2, figsize=(9.2, 4.3)) axL.axhline(ang_stat, color="0.5", ls="--", lw=1.2, label=f"at rest ({ang_stat:+.0f}°)") axL.axhline(ang_dyn, color="crimson", lw=1.6, label=f"dynamic angle of repose ({ang_dyn:+.0f}°)") axL.fill_between(t, ang_dyn - ang_sd, ang_dyn + ang_sd, color="crimson", alpha=0.12) axL.plot(t, ang_t, color="#2a6fdb", lw=1.3) axL.set_xlabel("time"); axL.set_ylabel("free-surface angle (°)"); axL.legend(loc="lower right") axL.set_title("a steady angle of repose"); axL.grid(True, alpha=0.3); axL.set_ylim(-5, max(35, ang_dyn + 12)) fx, fy, fs = frames[np.argmin(np.abs(ang_t - ang_dyn))] # a representative frame _, (slope, intercept) = surface_angle(fx, fy) axR.plot(cx + R * np.cos(theta), cy + R * np.sin(theta), "k-", lw=1.5) draw(axR, fx, fy, fs) xx = np.array([-0.6 * R, 0.6 * R]) axR.plot(cx + xx, cy + slope * xx + intercept, "-", color="crimson", lw=3) axR.set_xlim(cx - R - 2, cx + R + 2); axR.set_ylim(cy - R - 2, cy + R + 2) axR.set_aspect("equal"); axR.axis("off"); axR.set_title(f"free surface: {surface_angle(fx, fy)[0]:+.0f}°") plt.show() print(f"at rest: surface angle ≈ {ang_stat:+.1f}°") print(f"spinning: dynamic angle of repose ≈ {ang_dyn:+.1f}° ± {ang_sd:.0f}°") ``` At rest the surface is flat; spinning, the drum holds a **steady dynamic angle of repose** of roughly 25° while a thin surface layer tumbles — the cascading regime. Two ingredients make this work where a naive setup gives a limp, flat bed: - **The wall must grip the *whole* bed, not just the surface grains.** The wall-friction bound is the *force-chain* normal load — the weight of grains stacked above pressing a buried grain into the wall — so a deep bed is dragged up as a body instead of the wall slipping under it. - **The Froude number sets the regime.** Too slow → jerky block slumping; too fast → cataracting or centrifuging. `Fr ≈ 0.4` sits in the steady cascade. The bed settles to a **bounded, steady angle** rather than growing without limit — how we know the moving-wall contact solve is not injecting energy (an earlier version did, sending downhill grains to ever-higher speeds — the [sanity checks](../../sanity-checks/index.qmd) now guard against it). A few surface grains *do* briefly outrun the wall as they avalanche down the steep face, which is physical. Because smooth spheres *roll*, this angle is still gentler than angular sand (~35°); grip harder (`friction`), or use non-spherical [`build_particle`](../sdf-particle-packing/index.qmd) shapes, to steepen it. ## Adapt this yourself - **Vibrate the wall instead of spinning it.** A moving wall needs no rotation — drive the *linear* velocity sinusoidally each step to shake the bed (the geometry stays put): ```python for k in range(nsteps): v = A * omega_v * np.cos(omega_v * k * dt) # oscillating wall velocity sim.set_wall_velocity(wid, lin_vel=(0, v, 0), ang_vel=(0, 0, 0), center=(cx, cy, 0)) sim.step(dt) ``` - **Sweep the Froude number.** Lower `Fr` (slower `ω`) drops into slumping — the whole bed tilts and releases in blocks; raise it past ~0.5 for cataracting, where grains are flung off the top. - **Change the material.** `friction`/`restitution` in `add_to` are the *binary* particle–wall pair, independent of the grain–grain `set_material_params`. A slick wall lets the bed slump; a grippy one steepens the angle of repose. - **Change the container.** Any SDF works: a **hopper** is a cone `min`'d with a floor, a **chute** a tilted slab; add lifters/baffles by `min`-ing more primitives (keep + in the void, − in the wall). - **Go bigger / on a GPU or cluster.** Increase the grain count and drum size; on a CUDA/HIP build the identical script runs on the device, and a distributed build adds `sim.init_mpi(...)` / `sim.step_mpi(...)` (use a closed or rounded drum under MPI). ## Reproduce this ```bash # from a checkout of peclet-examples pip install -e '.[sim]' # or: pip install peclet PECLET_LOCAL_BUILD=/path/to/suite/dem/build_omp \ quarto render examples/rotating-drum/index.qmd --execute ```