Add Beverloo hopper discharge benchmark

Cargo.toml

@@ -126,6 +126,10 @@ path = "examples/example_template/main.rs"
 name = "cone_hopper"
 path = "examples/cone_hopper/main.rs"
 
+[[example]]
+name = "bench_beverloo_hopper"
+path = "examples/bench_beverloo_hopper/main.rs"
+
 [profile.release]
 opt-level = 3
 lto = "fat"

examples/bench_beverloo_hopper/README.md

@@ -0,0 +1,75 @@
+# Beverloo Hopper Discharge Benchmark
+
+Validates MDDEM hopper discharge mass flow rate against the **Beverloo correlation** for quasi-2D granular flow.
+
+## Physics
+
+The Beverloo equation predicts the steady-state mass flow rate through a hopper orifice:
+
+**3D (circular):** W = C · rho_bulk · sqrt(g) · (D - k·d)^(5/2)
+
+**2D (slot):** W = C · rho_bulk · sqrt(g) · (D - k·d)^(3/2) · depth
+
+where:
+- `D` = orifice width [m]
+- `d` = particle diameter [m]
+- `g` = gravitational acceleration [m/s^2]
+- `rho_bulk` = bulk density = rho_particle x phi (packing fraction ~0.58)
+- `C ~ 0.58` = empirical discharge coefficient
+- `k ~ 1.4` = empty annulus correction (particles can't flow within ~k·d/2 of the edge)
+- `depth` = slot length (periodic y extent for quasi-2D)
+
+## Setup
+
+- **Geometry**: Flat-bottom rectangular hopper with a central orifice, periodic in y (quasi-2D slab, 2d thick)
+- **Particles**: 800 monodisperse glass spheres, d = 2 mm, rho = 2500 kg/m^3
+- **Container**: 80 mm wide x 80 mm tall
+- **Orifice widths**: 6d, 10d, 15d (12 mm, 20 mm, 30 mm)
+- **Stages**: (1) Fill and settle under gravity, (2) Remove blocker wall and discharge
+- **Settling criterion**: Per-particle KE < 1e-7 J (avoids scaling issues with total KE)
+
+## Running
+
+### Single orifice width (default D = 10d):
+```bash
+cargo run --release --example bench_beverloo_hopper --no-default-features \
+    -- examples/bench_beverloo_hopper/config.toml
+```
+
+### Full sweep (3 orifice widths):
+```bash
+bash examples/bench_beverloo_hopper/run_sweep.sh
+```
+
+### Analysis:
+```bash
+python3 examples/bench_beverloo_hopper/validate.py   # PASS/FAIL checks
+python3 examples/bench_beverloo_hopper/plot.py        # Generate plots
+```
+
+## Expected Results
+
+The measured mass flow rate should agree with the Beverloo prediction within ~50% for each orifice width (small DEM systems with ~800 particles have significant statistical noise). The key physics test is the correct **scaling exponent of 3/2**: on a log-log plot of W vs (D - k·d), the data points should follow the theory line.
+
+![Beverloo comparison](beverloo_comparison.png)
+
+![Mass vs time](mass_vs_time.png)
+
+## Assumptions and Simplifications
+
+- **Quasi-2D**: Periodic boundary in y with 2d slab thickness. The 2D Beverloo formula is used with slab depth as the slot length.
+- **Monodisperse**: All particles have the same diameter.
+- **Flat bottom**: No converging hopper walls (funnel angle = 90 deg). Beverloo is valid for flat-bottom or steep hoppers.
+- **Hertz-Mindlin contacts**: Standard DEM contact model with friction = 0.3, restitution = 0.3.
+- **Reduced Young's modulus**: E = 5 MPa (vs. 70 GPa for real glass) for computational efficiency. This does not significantly affect steady-state flow rates.
+- **Small system**: 800 particles provides adequate statistics for Beverloo validation but with ~50% noise.
+
+## Files
+
+| File | Description |
+|------|-------------|
+| `main.rs` | Simulation: fill, settle (per-particle KE), discharge with tracking |
+| `config.toml` | Default config (D = 10d), fully documented |
+| `run_sweep.sh` | Shell script to run 3 orifice widths |
+| `validate.py` | Quantitative validation: PASS/FAIL per orifice |
+| `plot.py` | Generates beverloo_comparison.png and mass_vs_time.png |

examples/bench_beverloo_hopper/config.toml

@@ -0,0 +1,166 @@
+# Beverloo hopper discharge benchmark
+# Validates mass flow rate against W = C * rho_bulk * sqrt(g) * (D - k*d)^(3/2) * depth (2D)
+#
+# Flat-bottom rectangular hopper, periodic in y (quasi-2D slice).
+# Central orifice at bottom. Particles fill, settle, then discharge.
+#
+# Key parameters (change wall bounds to test different D values):
+#   Particle diameter d = 0.002 m
+#   Default orifice width D = 0.020 m (= 10d)
+#   Container width = 0.08 m (= 40d)
+#   Container height = 0.08 m (= 40d)
+#   Y depth = 0.004 m (= 2d, periodic -- quasi-2D)
+#
+# To sweep orifice widths, use run_sweep.sh
+
+[comm]
+processors_x = 1
+processors_y = 1
+processors_z = 1
+
+[domain]
+# Simulation box extents [m]
+x_low = -0.04
+x_high = 0.04
+y_low = 0.0
+y_high = 0.004               # Thin slab (2d) for quasi-2D
+z_low = -0.03
+z_high = 0.08
+periodic_x = false
+periodic_y = true             # Periodic in y for quasi-2D
+periodic_z = false
+
+[neighbor]
+skin_fraction = 1.1           # Neighbor search radius multiplier [dimensionless]
+bin_size = 0.005              # Neighbor-list bin width [m], >= largest diameter
+
+[gravity]
+gz = -9.81                    # Gravitational acceleration [m/s^2]
+
+[[dem.materials]]
+name = "glass"
+youngs_mod = 5.0e6            # Young's modulus [Pa] -- reduced for faster contacts
+poisson_ratio = 0.3           # Poisson's ratio [dimensionless]
+restitution = 0.3             # Coefficient of restitution [0-1] -- low for faster settling
+friction = 0.3                # Sliding friction coefficient [dimensionless]
+
+[[particles.insert]]
+material = "glass"
+count = 800                   # Number of particles
+radius = 0.001                # Particle radius [m] -> diameter d = 0.002 m
+density = 2500.0              # Particle density [kg/m^3]
+velocity_z = -0.5             # Initial downward velocity [m/s]
+# Insertion region: full width, thin y, upper part of box
+region = { type = "block", min = [-0.035, 0.0, 0.01], max = [0.035, 0.004, 0.075] }
+
+# -- Walls --
+
+# Left side wall (x = -0.04, normal +x)
+[[wall]]
+point_x = -0.04
+point_y = 0.0
+point_z = 0.0
+normal_x = 1.0
+normal_y = 0.0
+normal_z = 0.0
+material = "glass"
+
+# Right side wall (x = 0.04, normal -x)
+[[wall]]
+point_x = 0.04
+point_y = 0.0
+point_z = 0.0
+normal_x = -1.0
+normal_y = 0.0
+normal_z = 0.0
+material = "glass"
+
+# Ceiling (z = 0.08, normal -z)
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.08
+normal_x = 0.0
+normal_y = 0.0
+normal_z = -1.0
+material = "glass"
+
+# Floor LEFT (z = 0, normal +z) -- from x = -0.04 to x = -0.010 (orifice at center)
+# For D = 10d = 0.020 m, half-orifice = 0.010 m
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.0
+normal_x = 0.0
+normal_y = 0.0
+normal_z = 1.0
+material = "glass"
+bound_x_high = -0.010
+
+# Floor RIGHT (z = 0, normal +z) -- from x = +0.010 to x = +0.04
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.0
+normal_x = 0.0
+normal_y = 0.0
+normal_z = 1.0
+material = "glass"
+bound_x_low = 0.010
+
+# Blocker wall -- covers the orifice during filling (removed when settled)
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.0
+normal_x = 0.0
+normal_y = 0.0
+normal_z = 1.0
+material = "glass"
+name = "blocker"
+bound_x_low = -0.010
+bound_x_high = 0.010
+
+# Bottom catch wall (z = -0.03, normal +z) -- prevents particles from falling forever
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = -0.03
+normal_x = 0.0
+normal_y = 0.0
+normal_z = 1.0
+material = "glass"
+
+# -- Measure plane at orifice --
+
+[[measure_plane]]
+name = "orifice"
+point = [0.0, 0.0, -0.002]   # Just below the orifice floor
+normal = [0.0, 0.0, -1.0]    # Downward (positive = exiting hopper)
+report_interval = 1000        # Average over 1000 steps
+
+# -- Output --
+
+[output]
+dir = "examples/bench_beverloo_hopper"
+
+[vtp]
+interval = 10000              # VTP visualization every N steps
+
+# -- Run stages --
+
+# Stage 1: fill and settle
+[[run]]
+name = "filling"
+dt = 5.0e-6                   # Timestep [s]
+steps = 200000                # Allow up to 200k steps for settling
+thermo = 2000
+thermo_keys = ["step", "atoms", "ke", "time"]
+
+# Stage 2: discharge through orifice
+[[run]]
+name = "discharge"
+dt = 5.0e-6                   # Timestep [s]
+steps = 400000                # Discharge period (2.0 s simulation time)
+thermo = 2000
+thermo_keys = ["step", "atoms", "ke", "time", "flow_rate_orifice", "crossings_orifice"]