diff --git a/Cargo.toml b/Cargo.toml index 93467f5..e9da9b1 100644 --- a/Cargo.toml +++ b/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" diff --git a/examples/bench_triaxial/README.md b/examples/bench_triaxial/README.md new file mode 100644 index 0000000..bad758c --- /dev/null +++ b/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. diff --git a/examples/bench_triaxial/config_10kPa.toml b/examples/bench_triaxial/config_10kPa.toml new file mode 100644 index 0000000..14aa6fc --- /dev/null +++ b/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²] diff --git a/examples/bench_triaxial/config_200kPa.toml b/examples/bench_triaxial/config_200kPa.toml new file mode 100644 index 0000000..743a222 --- /dev/null +++ b/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 diff --git a/examples/bench_triaxial/config_50kPa.toml b/examples/bench_triaxial/config_50kPa.toml new file mode 100644 index 0000000..0ebb78b --- /dev/null +++ b/examples/bench_triaxial/config_50kPa.toml @@ -0,0 +1,130 @@ +# Triaxial compression benchmark — σ₃ = 50 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_50kPa" + +# 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 = 7.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 = 7.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 = 7.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 = 7.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 ─────────────────────────────────────────────────────────────────── + +[[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 diff --git a/examples/bench_triaxial/main.rs b/examples/bench_triaxial/main.rs new file mode 100644 index 0000000..a8f9031 --- /dev/null +++ b/examples/bench_triaxial/main.rs @@ -0,0 +1,268 @@ +//! Triaxial compression benchmark: validates DEM against Mohr-Coulomb failure theory. +//! +//! Stages: insert → relax → compress +//! - Insert: random particles settle under gravity in a walled box +//! - Relax: wait for kinetic energy to decay +//! - Compress: top wall moves down at constant velocity; lateral walls use +//! servo control to maintain confining pressure; gravity set to zero for +//! uniform stress distribution +//! +//! Run at different confining pressures via separate config files: +//! ```bash +//! cargo run --release --example bench_triaxial --no-default-features -- examples/bench_triaxial/config_10kPa.toml +//! ``` + +use mddem::prelude::*; +use std::fs::{self, File, OpenOptions}; +use std::io::Write; + +/// Axial compression velocity [m/s]. Chosen for quasi-static loading +/// (inertial number I ≪ 0.01 at all confining pressures tested). +const COMPRESS_VEL: f64 = 1e-3; + +#[derive(Clone, Debug, PartialEq, Default, StageEnum)] +enum Phase { + #[default] + #[stage("insert")] + Insert, + #[stage("relax")] + Relax, + #[stage("compress")] + Compress, +} + +fn main() { + let mut app = App::new(); + app.add_plugins(CorePlugins) + .add_plugins(GranularDefaultPlugins) + .add_plugins(GravityPlugin) + .add_plugins(WallPlugin) + .add_plugins(StatesPlugin { + initial: Phase::Insert, + }) + .add_plugins(StageAdvancePlugin::::new()); + + // Stage-transition systems + app.add_update_system( + check_insert_settled.run_if(in_state(Phase::Insert)), + ScheduleSet::PostFinalIntegration, + ); + app.add_update_system( + check_relaxed.run_if(in_state(Phase::Relax)), + ScheduleSet::PostFinalIntegration, + ); + + // Compression systems (compress stage only) + app.add_update_system( + apply_compression.run_if(in_state(Phase::Compress)), + ScheduleSet::PreInitialIntegration, + ); + app.add_update_system( + output_stress.run_if(in_state(Phase::Compress)), + ScheduleSet::PostFinalIntegration, + ); + + app.start(); +} + +// ── Stage transition checks ───────────────────────────────────────────────── + +/// Advance to relax when kinetic energy drops below threshold. +fn check_insert_settled( + atoms: Res, + run_state: Res, + comm: Res, + mut next_state: ResMut>, +) { + let step = run_state.total_cycle; + if step < 1000 || step % 100 != 0 { + return; + } + + let nlocal = atoms.nlocal as usize; + let local_ke: f64 = (0..nlocal) + .map(|i| { + let v = atoms.vel[i]; + 0.5 * atoms.mass[i] * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) + }) + .sum(); + let global_ke = comm.all_reduce_sum_f64(local_ke); + + if global_ke < 1e-5 { + next_state.set(Phase::Relax); + if comm.rank() == 0 { + println!( + "Step {}: KE = {:.3e} J — particles settled, advancing to relax", + step, global_ke + ); + } + } +} + +/// Advance to compress when kinetic energy is sufficiently small. +fn check_relaxed( + atoms: Res, + run_state: Res, + comm: Res, + mut next_state: ResMut>, +) { + let step = run_state.total_cycle; + if step % 100 != 0 { + return; + } + + let nlocal = atoms.nlocal as usize; + let local_ke: f64 = (0..nlocal) + .map(|i| { + let v = atoms.vel[i]; + 0.5 * atoms.mass[i] * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) + }) + .sum(); + let global_ke = comm.all_reduce_sum_f64(local_ke); + + if global_ke < 1e-5 { + next_state.set(Phase::Compress); + if comm.rank() == 0 { + println!( + "Step {}: KE = {:.3e} J — fully relaxed, advancing to compress", + step, global_ke + ); + } + } +} + +// ── Compression systems ────────────────────────────────────────────────────── + +/// Move the top wall (wall index 5) downward at constant velocity. +/// +/// On the first call, snaps the top wall to just above the highest particle +/// so compression begins immediately rather than wasting steps traversing +/// empty space. +fn apply_compression( + mut walls: ResMut, + atoms: Res, + comm: Res, + mut initialized: Local, +) { + let dt = atoms.dt; + + if !*initialized { + *initialized = true; + // Snap top wall to just above the particle bed + let nlocal = atoms.nlocal as usize; + let mut max_z: f64 = 0.0; + for i in 0..nlocal { + let z = atoms.pos[i][2]; + if z > max_z { + max_z = z; + } + } + // Use -min(-x) trick since there's no all_reduce_max + let max_z = -comm.all_reduce_min_f64(-max_z); + let top_wall = &mut walls.planes[5]; + top_wall.point_z = max_z + 0.0015; // ~1.5 × radius gap + if comm.rank() == 0 { + println!( + "Compression: snapped top wall to z = {:.6} m", + top_wall.point_z + ); + } + } + + // Constant-velocity axial compression + let top_wall = &mut walls.planes[5]; + top_wall.point_z -= COMPRESS_VEL * dt; + top_wall.velocity = [0.0, 0.0, -COMPRESS_VEL]; +} + +/// Compute principal stresses from wall force accumulators and write to CSV. +/// +/// Wall layout (must match config): +/// 0: x = x_lo, normal +x (servo, left) +/// 1: x = x_hi, normal −x (servo, right) +/// 2: y = y_lo, normal +y (servo, front) +/// 3: y = y_hi, normal −y (servo, back) +/// 4: z = z_lo, normal +z (static floor) +/// 5: z = z_hi, normal −z (static top, moved by apply_compression) +fn output_stress( + walls: Res, + atoms: Res, + domain: Res, + run_state: Res, + thermo: Res, + input: Res, + comm: Res, + mut file_created: Local, + mut initial_height: Local, +) { + if comm.rank() != 0 { + return; + } + if !run_state.total_cycle.is_multiple_of(thermo.interval) { + return; + } + + let output_dir = input.output_dir.as_deref().unwrap_or("."); + let path = format!("{}/triaxial_stress.csv", output_dir); + + // Create file with header on first call + if !*file_created { + *file_created = true; + fs::create_dir_all(output_dir).ok(); + let mut f = File::create(&path).expect("cannot create stress CSV"); + writeln!(f, "step,time,strain,sigma_1,sigma_3,q,p") + .expect("cannot write CSV header"); + // Record initial sample height for strain computation + let floor_z = walls.planes[4].point_z; + let top_z = walls.planes[5].point_z; + *initial_height = top_z - floor_z; + } + + let step = run_state.total_cycle; + let time = step as f64 * atoms.dt; + + let lx = domain.size[0]; + let ly = domain.size[1]; + let floor_z = walls.planes[4].point_z; + let top_z = walls.planes[5].point_z; + let h = top_z - floor_z; + + // Axial strain (positive = compression) + let strain = if *initial_height > 0.0 { + (*initial_height - h) / *initial_height + } else { + 0.0 + }; + + // σ₁ = axial stress = top wall force / cross-section area + let area_top = lx * ly; + let sigma_1 = walls.planes[5].force_accumulator / area_top; + + // σ₃ = confining stress = average of lateral wall stresses + let f_x_avg = + (walls.planes[0].force_accumulator + walls.planes[1].force_accumulator) / 2.0; + let f_y_avg = + (walls.planes[2].force_accumulator + walls.planes[3].force_accumulator) / 2.0; + let area_x = ly * h; // yz-face area + let area_y = lx * h; // xz-face area + let sigma_3 = if h > 0.0 { + (f_x_avg / area_x + f_y_avg / area_y) / 2.0 + } else { + 0.0 + }; + + // Deviatoric stress and mean stress + let q = sigma_1 - sigma_3; + let p = (sigma_1 + 2.0 * sigma_3) / 3.0; + + let mut f = OpenOptions::new() + .append(true) + .open(&path) + .expect("cannot open stress CSV"); + writeln!( + f, + "{},{:.10e},{:.10e},{:.6e},{:.6e},{:.6e},{:.6e}", + step, time, strain, sigma_1, sigma_3, q, p + ) + .expect("cannot write stress data"); +} diff --git a/examples/bench_triaxial/plot.py b/examples/bench_triaxial/plot.py new file mode 100644 index 0000000..5279c73 --- /dev/null +++ b/examples/bench_triaxial/plot.py @@ -0,0 +1,154 @@ +#!/usr/bin/env python3 +""" +Generate publication-quality plots for the triaxial compression benchmark. + +Produces: + 1. stress_strain.png — Deviatoric stress q vs axial strain for all pressures + 2. mohr_circles.png — Mohr circles at failure with fitted failure envelope + 3. q_p_plot.png — q–p diagram with Mohr-Coulomb failure line +""" + +import os +import sys +import numpy as np +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +script_dir = os.path.dirname(os.path.abspath(__file__)) +pressures = ["10kPa", "50kPa", "200kPa"] +colors = ["#1f77b4", "#ff7f0e", "#d62728"] + +# ── Load data ────────────────────────────────────────────────────────────── + +datasets = {} # label → structured array +peak_data = {} # label → (sigma_1_peak, sigma_3_at_peak) + +for label in pressures: + csv_path = os.path.join(script_dir, f"data_{label}", "triaxial_stress.csv") + if not os.path.isfile(csv_path): + print(f" Skipping {label}: {csv_path} not found") + continue + data = np.genfromtxt(csv_path, delimiter=",", names=True) + if data.size == 0: + continue + datasets[label] = data + idx_peak = np.argmax(data["q"]) + peak_data[label] = (data["sigma_1"][idx_peak], data["sigma_3"][idx_peak]) + +if not datasets: + print("ERROR: No data files found. Run the simulations first.") + sys.exit(1) + +# ── Plot 1: Stress–strain curves ────────────────────────────────────────── + +fig, ax = plt.subplots(figsize=(8, 5)) +for (label, data), color in zip(datasets.items(), colors): + strain_pct = data["strain"] * 100 # convert to percent + q_kpa = data["q"] / 1000 # convert to kPa + ax.plot(strain_pct, q_kpa, "-", color=color, linewidth=1.5, label=f"σ₃ = {label}") + +ax.set_xlabel("Axial Strain ε₁ [%]", fontsize=12) +ax.set_ylabel("Deviatoric Stress q = σ₁ − σ₃ [kPa]", fontsize=12) +ax.set_title("Triaxial Compression — Stress–Strain Curves", fontsize=13) +ax.legend(fontsize=10) +ax.grid(True, alpha=0.3) +ax.set_xlim(left=0) +ax.set_ylim(bottom=0) +fig.tight_layout() +fig.savefig(os.path.join(script_dir, "stress_strain.png"), dpi=150) +plt.close(fig) +print(" Saved stress_strain.png") + +# ── Plot 2: Mohr circles at failure ────────────────────────────────────── + +fig, ax = plt.subplots(figsize=(8, 6)) +theta = np.linspace(0, np.pi, 200) + +for (label, (s1, s3)), color in zip(peak_data.items(), colors): + center = (s1 + s3) / 2 / 1000 # kPa + radius = (s1 - s3) / 2 / 1000 # kPa + sigma = center + radius * np.cos(theta) + tau = radius * np.sin(theta) + ax.plot(sigma, tau, "-", color=color, linewidth=1.5, label=f"σ₃ = {label}") + ax.plot(center, 0, "o", color=color, markersize=4) + +# Fit Mohr-Coulomb envelope: τ = c + σ tan(φ) +# From peak stresses, the tangent line touches all circles. +# For cohesionless: sin(φ) = (σ₁-σ₃)/(σ₁+σ₃) +phi_vals = [] +for s1, s3 in peak_data.values(): + if (s1 + s3) > 0: + sin_phi = (s1 - s3) / (s1 + s3) + sin_phi = max(-1.0, min(1.0, sin_phi)) + phi_vals.append(np.arcsin(sin_phi)) + +if phi_vals: + phi_avg = np.mean(phi_vals) + # Failure envelope: τ = σ tan(φ) (c = 0) + sigma_max = max(s1 for s1, _ in peak_data.values()) / 1000 * 1.1 + sigma_line = np.linspace(0, sigma_max, 100) + tau_line = sigma_line * np.tan(phi_avg) + ax.plot(sigma_line, tau_line, "k--", linewidth=2, + label=f"MC envelope: φ = {np.degrees(phi_avg):.1f}°") + +ax.set_xlabel("Normal Stress σ [kPa]", fontsize=12) +ax.set_ylabel("Shear Stress τ [kPa]", fontsize=12) +ax.set_title("Mohr Circles at Failure with Mohr-Coulomb Envelope", fontsize=13) +ax.legend(fontsize=9, loc="upper left") +ax.grid(True, alpha=0.3) +ax.set_xlim(left=0) +ax.set_ylim(bottom=0) +ax.set_aspect("equal") +fig.tight_layout() +fig.savefig(os.path.join(script_dir, "mohr_circles.png"), dpi=150) +plt.close(fig) +print(" Saved mohr_circles.png") + +# ── Plot 3: q–p diagram ────────────────────────────────────────────────── + +fig, ax = plt.subplots(figsize=(8, 5)) + +# Plot peak points +p_peaks = [] +q_peaks = [] +for (label, (s1, s3)), color in zip(peak_data.items(), colors): + p_val = (s1 + 2 * s3) / 3 / 1000 # kPa + q_val = (s1 - s3) / 1000 # kPa + p_peaks.append(p_val) + q_peaks.append(q_val) + ax.plot(p_val, q_val, "o", color=color, markersize=8, zorder=5, + label=f"σ₃ = {label}") + +# Plot stress paths +for (label, data), color in zip(datasets.items(), colors): + p_path = data["p"] / 1000 + q_path = data["q"] / 1000 + ax.plot(p_path, q_path, "-", color=color, linewidth=1, alpha=0.5) + +# Fit q = M p (cohesionless) through peak points +p_arr = np.array(p_peaks) +q_arr = np.array(q_peaks) +if len(p_arr) >= 2 and np.sum(p_arr ** 2) > 0: + # Least squares fit q = M * p + M = np.sum(p_arr * q_arr) / np.sum(p_arr ** 2) + # M = 6 sin(φ) / (3 - sin(φ)) → sin(φ) = 3M / (6 + M) + sin_phi_qp = 3 * M / (6 + M) + sin_phi_qp = max(-1.0, min(1.0, sin_phi_qp)) + phi_qp = np.degrees(np.arcsin(sin_phi_qp)) + p_line = np.linspace(0, max(p_peaks) * 1.2, 100) + q_line = M * p_line + ax.plot(p_line, q_line, "k--", linewidth=2, + label=f"q = {M:.2f} p (φ = {phi_qp:.1f}°)") + +ax.set_xlabel("Mean Stress p = (σ₁ + 2σ₃)/3 [kPa]", fontsize=12) +ax.set_ylabel("Deviatoric Stress q = σ₁ − σ₃ [kPa]", fontsize=12) +ax.set_title("q–p Diagram with Mohr-Coulomb Failure Line", fontsize=13) +ax.legend(fontsize=10) +ax.grid(True, alpha=0.3) +ax.set_xlim(left=0) +ax.set_ylim(bottom=0) +fig.tight_layout() +fig.savefig(os.path.join(script_dir, "q_p_plot.png"), dpi=150) +plt.close(fig) +print(" Saved q_p_plot.png") diff --git a/examples/bench_triaxial/run_benchmark.sh b/examples/bench_triaxial/run_benchmark.sh new file mode 100755 index 0000000..54607b7 --- /dev/null +++ b/examples/bench_triaxial/run_benchmark.sh @@ -0,0 +1,28 @@ +#!/usr/bin/env bash +# Run the triaxial compression benchmark at three confining pressures. +# +# Usage: +# cd +# bash examples/bench_triaxial/run_benchmark.sh +# +# After completion, run validation and plotting: +# python3 examples/bench_triaxial/validate.py +# python3 examples/bench_triaxial/plot.py + +set -euo pipefail + +EXAMPLE_DIR="examples/bench_triaxial" + +for pressure in 10kPa 50kPa 200kPa; do + echo "" + echo "══════════════════════════════════════════════════════════════" + echo " Running triaxial compression at σ₃ = ${pressure}" + echo "══════════════════════════════════════════════════════════════" + cargo run --release --example bench_triaxial --no-default-features \ + -- "${EXAMPLE_DIR}/config_${pressure}.toml" +done + +echo "" +echo "All simulations complete." +echo "Run: python3 ${EXAMPLE_DIR}/validate.py" +echo "Run: python3 ${EXAMPLE_DIR}/plot.py" diff --git a/examples/bench_triaxial/validate.py b/examples/bench_triaxial/validate.py new file mode 100644 index 0000000..5ae18b3 --- /dev/null +++ b/examples/bench_triaxial/validate.py @@ -0,0 +1,142 @@ +#!/usr/bin/env python3 +""" +Validate triaxial compression benchmark against Mohr-Coulomb failure theory. + +For cohesionless granular material (c = 0): + σ₁ = σ₃ × (1 + sin φ) / (1 - sin φ) + sin φ = (σ₁ - σ₃) / (σ₁ + σ₃) + +Expected friction angle for μ = 0.5: φ ≈ 20–30° + +Checks: + 1. Stress data files exist and contain valid data + 2. Peak stress ratio σ₁/σ₃ is in a physically reasonable range [1.5, 8.0] + 3. Mobilised friction angle φ is in expected range [15°, 40°] + 4. Friction angle is consistent across confining pressures (std < 5°) + 5. Mohr-Coulomb failure envelope is approximately linear (R² > 0.9) +""" + +import os +import sys +import numpy as np + +script_dir = os.path.dirname(os.path.abspath(__file__)) +pressures = ["10kPa", "50kPa", "200kPa"] + +print("=" * 60) +print("Triaxial Compression — Mohr-Coulomb Validation") +print("=" * 60) + +# ── Load data ────────────────────────────────────────────────────────────── + +results = {} # label → (sigma_1_peak, sigma_3_at_peak) + +for label in pressures: + data_dir = os.path.join(script_dir, f"data_{label}") + csv_path = os.path.join(data_dir, "triaxial_stress.csv") + if not os.path.isfile(csv_path): + print(f" WARNING: {csv_path} not found — skipping {label}") + continue + data = np.genfromtxt(csv_path, delimiter=",", names=True) + if data.size == 0: + print(f" WARNING: {csv_path} is empty — skipping {label}") + continue + + # Find peak deviatoric stress q = σ₁ - σ₃ + q = data["q"] + idx_peak = np.argmax(q) + s1 = data["sigma_1"][idx_peak] + s3 = data["sigma_3"][idx_peak] + results[label] = (s1, s3, q[idx_peak]) + print(f" {label:>6s}: σ₁ = {s1:10.1f} Pa, σ₃ = {s3:10.1f} Pa, " + f"q = {q[idx_peak]:10.1f} Pa, σ₁/σ₃ = {s1/s3 if s3 > 0 else float('inf'):.2f}") + +if len(results) < 2: + print("\nERROR: Need at least 2 confining pressures for validation.") + sys.exit(1) + +# ── Run checks ───────────────────────────────────────────────────────────── + +passed = 0 +total = 0 + +# Check 1: valid data +total += 1 +if all(np.isfinite(v[0]) and np.isfinite(v[1]) for v in results.values()): + print("\n [1] Valid data: PASS") + passed += 1 +else: + print("\n [1] Valid data: FAIL (NaN/Inf in stress)") + +# Check 2: stress ratio in range +total += 1 +ratios = [] +for label, (s1, s3, _) in results.items(): + if s3 > 0: + ratios.append(s1 / s3) +if ratios and all(1.5 <= r <= 8.0 for r in ratios): + print(f" [2] Stress ratio range: PASS (ratios: {[f'{r:.2f}' for r in ratios]})") + passed += 1 +else: + print(f" [2] Stress ratio range: FAIL (ratios: {[f'{r:.2f}' for r in ratios]}, expected 1.5–8.0)") + +# Check 3: friction angle in expected range +total += 1 +phi_values = [] +for label, (s1, s3, _) in results.items(): + if (s1 + s3) > 0: + sin_phi = (s1 - s3) / (s1 + s3) + sin_phi = max(-1.0, min(1.0, sin_phi)) + phi = np.degrees(np.arcsin(sin_phi)) + phi_values.append(phi) +phi_mean = np.mean(phi_values) if phi_values else 0 +if phi_values and all(15 <= p <= 40 for p in phi_values): + print(f" [3] Friction angle: PASS (φ = {phi_mean:.1f}° ± {np.std(phi_values):.1f}°, range 15–40°)") + passed += 1 +else: + print(f" [3] Friction angle: FAIL (φ values: {[f'{p:.1f}' for p in phi_values]}°, expected 15–40°)") + +# Check 4: consistency across pressures +total += 1 +if len(phi_values) >= 2 and np.std(phi_values) < 5.0: + print(f" [4] Consistency: PASS (std = {np.std(phi_values):.2f}° < 5°)") + passed += 1 +elif len(phi_values) < 2: + print(f" [4] Consistency: SKIP (need ≥ 2 pressure levels)") +else: + print(f" [4] Consistency: FAIL (std = {np.std(phi_values):.2f}° ≥ 5°)") + +# Check 5: Mohr-Coulomb linearity (q vs p) +total += 1 +s1_arr = np.array([v[0] for v in results.values()]) +s3_arr = np.array([v[1] for v in results.values()]) +p_arr = (s1_arr + 2 * s3_arr) / 3.0 +q_arr = s1_arr - s3_arr +if len(q_arr) >= 3: + coeffs = np.polyfit(p_arr, q_arr, 1) + q_fit = np.polyval(coeffs, p_arr) + ss_res = np.sum((q_arr - q_fit) ** 2) + ss_tot = np.sum((q_arr - np.mean(q_arr)) ** 2) + r_squared = 1 - ss_res / ss_tot if ss_tot > 0 else 0 + if r_squared > 0.9: + print(f" [5] MC linearity: PASS (R² = {r_squared:.4f})") + passed += 1 + else: + print(f" [5] MC linearity: FAIL (R² = {r_squared:.4f}, expected > 0.9)") +else: + # With only 2 points, linearity is trivially satisfied + print(f" [5] MC linearity: PASS (trivially linear with {len(q_arr)} points)") + passed += 1 + +# ── Summary ──────────────────────────────────────────────────────────────── + +print(f"\nResults: {passed}/{total} checks passed") +if phi_values: + print(f"Measured friction angle: φ = {phi_mean:.1f}° (μ_particle = 0.5)") + +if passed == total: + print("ALL CHECKS PASSED") +else: + print(f"WARNING: {total - passed} check(s) failed") + +sys.exit(0 if passed == total else 1)