Add triaxial compression benchmark (Mohr-Coulomb validation)

Cargo.toml

@@ -102,6 +102,10 @@ path = "examples/polymer_chain/main.rs"
 name = "dem_triaxial"
 path = "examples/dem_triaxial/main.rs"
 
+[[example]]
+name = "bench_triaxial"
+path = "examples/bench_triaxial/main.rs"
+
 [[example]]
 name = "shrink_wrap_demo"
 path = "examples/shrink_wrap_demo/main.rs"

examples/bench_triaxial/README.md

@@ -0,0 +1,81 @@
+# Triaxial Compression Benchmark
+
+Validates DEM triaxial compression against **Mohr-Coulomb failure theory**.
+
+## Physics
+
+A sample of 100 randomly packed spheres is confined laterally by servo-controlled
+walls at a specified confining pressure σ₃, then compressed axially at constant
+velocity. The Mohr-Coulomb criterion predicts a linear failure envelope:
+
+    σ₁ = σ₃ × (1 + sin φ) / (1 - sin φ) + 2c cos φ / (1 - sin φ)
+
+For cohesionless particles (c = 0):
+
+    sin φ = (σ₁ − σ₃) / (σ₁ + σ₃)
+
+where σ₁ is peak axial stress, σ₃ is confining pressure, and φ is the internal
+friction angle. For μ = 0.5 inter-particle friction, theory predicts φ ≈ 20–35°.
+
+## Setup
+
+| Parameter            | Value           | Units  |
+|---------------------|-----------------|--------|
+| Box dimensions      | 10 × 10 × 40   | mm     |
+| Particle radius     | 1               | mm     |
+| Particle count      | 100             | –      |
+| Particle density    | 2500            | kg/m³  |
+| Young's modulus     | 10 MPa          | Pa     |
+| Poisson's ratio     | 0.3             | –      |
+| Restitution         | 0.3             | –      |
+| Friction (μ)        | 0.5             | –      |
+| Confining pressures | 10, 50, 200     | kPa    |
+
+### Simplifications
+
+- **Monodisperse** spheres (no size distribution)
+- **Soft particles** (E = 10 MPa) for fast timestep; contact overlaps remain < 2%
+- Gravity set to zero during compression for **uniform stress distribution**
+- Servo walls target constant **force** (not constant pressure), so σ₃ varies
+  slightly with sample height changes during compression
+- **100 particles** — small sample for fast runtime (~2 s per confining pressure)
+
+## Running
+
+```bash
+# Run all three confining pressures (~4 s total)
+bash examples/bench_triaxial/run_benchmark.sh
+
+# Or run a single pressure
+cargo run --release --example bench_triaxial --no-default-features \
+    -- examples/bench_triaxial/config_10kPa.toml
+
+# Validate results
+python3 examples/bench_triaxial/validate.py
+
+# Generate plots
+python3 examples/bench_triaxial/plot.py
+```
+
+## Validation
+
+`validate.py` checks:
+
+1. **Valid data** — no NaN/Inf in stress measurements
+2. **Stress ratio** — peak σ₁/σ₃ is in range [1.5, 8.0]
+3. **Friction angle** — φ is in range [15°, 40°]
+4. **Consistency** — φ is consistent across confining pressures (σ < 5°)
+5. **MC linearity** — q–p failure line has R² > 0.9
+
+## Plots
+
+- **stress_strain.png** — Deviatoric stress q vs axial strain at all pressures
+- **mohr_circles.png** — Mohr circles at failure with fitted envelope
+- **q_p_plot.png** — q–p diagram with Mohr-Coulomb failure line
+
+## References
+
+- Cundall, P.A. & Strack, O.D.L. (1979). A discrete numerical model for
+  granular assemblies. *Géotechnique*, 29(1), 47–65.
+- Thornton, C. (2000). Numerical simulations of deviatoric shear deformation
+  of granular media. *Géotechnique*, 50(1), 43–53.

examples/bench_triaxial/config_10kPa.toml

@@ -0,0 +1,143 @@
+# Triaxial compression benchmark — σ₃ = 10 kPa confining pressure
+#
+# Box: 0.01 × 0.01 × 0.04 m (walled, non-periodic)
+# Particles: 100 spheres, r = 1 mm, ρ = 2500 kg/m³
+# Material: E = 1e7 Pa (soft for fast timestep), ν = 0.3, e = 0.3, μ = 0.5
+# Servo walls maintain lateral confining pressure; top wall compressed axially.
+
+[comm]
+processors_x = 1
+processors_y = 1
+processors_z = 1
+
+[domain]
+x_low  = 0.0
+x_high = 0.01        # 10 mm width [m]
+y_low  = 0.0
+y_high = 0.01        # 10 mm depth [m]
+z_low  = 0.0
+z_high = 0.04        # 40 mm tall for insertion [m]
+periodic_x = false
+periodic_y = false
+periodic_z = false
+
+[neighbor]
+skin_fraction = 1.1  # neighbor search radius multiplier [–]
+bin_size = 0.003     # bin width [m], ≥ largest particle diameter
+
+[gravity]
+gz = -9.81           # gravitational acceleration [m/s²]
+
+[dem]
+contact_model = "hertz"
+
+[[dem.materials]]
+name = "granular"
+youngs_mod = 1e7     # Young's modulus [Pa] — soft for fast auto-dt (~1.3e-5 s)
+poisson_ratio = 0.3  # Poisson's ratio [–]
+restitution = 0.3    # coefficient of restitution [–] — high damping for fast settling
+friction = 0.5       # sliding friction coefficient [–]
+
+[[particles.insert]]
+material = "granular"
+count = 100          # number of particles
+radius = 0.001       # particle radius [m]
+density = 2500.0     # particle density [kg/m³]
+
+[output]
+dir = "examples/bench_triaxial/data_10kPa"
+
+# ── Walls ────────────────────────────────────────────────────────────────────
+# Walls 0–3: lateral servo walls for confining pressure
+# Wall  4:   static floor
+# Wall  5:   top wall (moved by custom system during compress stage)
+#
+# Servo target_force = σ₃ × wall_area
+#   σ₃ = 10 kPa = 10000 Pa
+#   Wall area ≈ Ly × H_bed ≈ 0.01 × 0.015 = 1.5e-4 m²
+#   target_force = 10000 × 1.5e-4 = 1.5 N
+
+# Wall 0: x = 0, normal +x (left, servo)
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.0
+normal_x = 1.0
+normal_y = 0.0
+normal_z = 0.0
+material = "granular"
+servo = { target_force = 1.5, max_velocity = 0.005, gain = 0.005 }
+
+# Wall 1: x = 0.01, normal −x (right, servo)
+[[wall]]
+point_x = 0.01
+point_y = 0.0
+point_z = 0.0
+normal_x = -1.0
+normal_y = 0.0
+normal_z = 0.0
+material = "granular"
+servo = { target_force = 1.5, max_velocity = 0.005, gain = 0.005 }
+
+# Wall 2: y = 0, normal +y (front, servo)
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.0
+normal_x = 0.0
+normal_y = 1.0
+normal_z = 0.0
+material = "granular"
+servo = { target_force = 1.5, max_velocity = 0.005, gain = 0.005 }
+
+# Wall 3: y = 0.01, normal −y (back, servo)
+[[wall]]
+point_x = 0.0
+point_y = 0.01
+point_z = 0.0
+normal_x = 0.0
+normal_y = -1.0
+normal_z = 0.0
+material = "granular"
+servo = { target_force = 1.5, max_velocity = 0.005, gain = 0.005 }
+
+# Wall 4: z = 0, normal +z (floor, static)
+[[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 = "granular"
+
+# Wall 5: z = 0.04, normal −z (top, static — moved by custom system)
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.04
+normal_x = 0.0
+normal_y = 0.0
+normal_z = -1.0
+material = "granular"
+
+# ── Stages ───────────────────────────────────────────────────────────────────
+
+# Stage 1: insert and settle under gravity
+[[run]]
+name = "insert"
+steps = 50000
+thermo = 5000
+
+# Stage 2: relax to low KE
+[[run]]
+name = "relax"
+steps = 30000
+thermo = 2000
+
+# Stage 3: axial compression with zero gravity for uniform stress
+[[run]]
+name = "compress"
+steps = 100000
+thermo = 500
+gravity.gz = 0.0     # zero gravity during compression for uniform stress [m/s²]

examples/bench_triaxial/config_200kPa.toml

@@ -0,0 +1,130 @@
+# Triaxial compression benchmark — σ₃ = 200 kPa confining pressure
+#
+# Box: 0.01 × 0.01 × 0.04 m (walled, non-periodic)
+# Particles: 100 spheres, r = 1 mm, ρ = 2500 kg/m³
+# Material: E = 1e7 Pa (soft for fast timestep), ν = 0.3, e = 0.3, μ = 0.5
+# Servo walls maintain lateral confining pressure; top wall compressed axially.
+
+[comm]
+processors_x = 1
+processors_y = 1
+processors_z = 1
+
+[domain]
+x_low  = 0.0
+x_high = 0.01
+y_low  = 0.0
+y_high = 0.01
+z_low  = 0.0
+z_high = 0.04
+periodic_x = false
+periodic_y = false
+periodic_z = false
+
+[neighbor]
+skin_fraction = 1.1
+bin_size = 0.003
+
+[gravity]
+gz = -9.81
+
+[dem]
+contact_model = "hertz"
+
+[[dem.materials]]
+name = "granular"
+youngs_mod = 1e7
+poisson_ratio = 0.3
+restitution = 0.3
+friction = 0.5
+
+[[particles.insert]]
+material = "granular"
+count = 100
+radius = 0.001
+density = 2500.0
+
+[output]
+dir = "examples/bench_triaxial/data_200kPa"
+
+# Wall 0: x = 0, normal +x (left, servo)
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.0
+normal_x = 1.0
+normal_y = 0.0
+normal_z = 0.0
+material = "granular"
+servo = { target_force = 30.0, max_velocity = 0.005, gain = 0.005 }
+
+# Wall 1: x = 0.01, normal −x (right, servo)
+[[wall]]
+point_x = 0.01
+point_y = 0.0
+point_z = 0.0
+normal_x = -1.0
+normal_y = 0.0
+normal_z = 0.0
+material = "granular"
+servo = { target_force = 30.0, max_velocity = 0.005, gain = 0.005 }
+
+# Wall 2: y = 0, normal +y (front, servo)
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.0
+normal_x = 0.0
+normal_y = 1.0
+normal_z = 0.0
+material = "granular"
+servo = { target_force = 30.0, max_velocity = 0.005, gain = 0.005 }
+
+# Wall 3: y = 0.01, normal −y (back, servo)
+[[wall]]
+point_x = 0.0
+point_y = 0.01
+point_z = 0.0
+normal_x = 0.0
+normal_y = -1.0
+normal_z = 0.0
+material = "granular"
+servo = { target_force = 30.0, max_velocity = 0.005, gain = 0.005 }
+
+# Wall 4: z = 0, normal +z (floor, static)
+[[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 = "granular"
+
+# Wall 5: z = 0.04, normal −z (top, static — moved by custom system)
+[[wall]]
+point_x = 0.0
+point_y = 0.0
+point_z = 0.04
+normal_x = 0.0
+normal_y = 0.0
+normal_z = -1.0
+material = "granular"
+
+# ── Stages ───────────────────────────────────────────────────────────────────
+
+[[run]]
+name = "insert"
+steps = 50000
+thermo = 5000
+
+[[run]]
+name = "relax"
+steps = 30000
+thermo = 2000
+
+[[run]]
+name = "compress"
+steps = 100000
+thermo = 500
+gravity.gz = 0.0