From 8e41a104b1becf0d88bdf7d435f0257d9a0de380 Mon Sep 17 00:00:00 2001 From: elizabeth-suehr Date: Mon, 16 Mar 2026 00:35:51 -0700 Subject: [PATCH] Add hopper discharge, rotating drum examples and measurement plane enhancements MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Measurement plane enhancements (dem_measure_plane): - Circular measurement planes (shape="circular" + radius) for cylindrical geometries - Rolling average statistics (rolling_windows) for smoothed flow rate tracking - Per-type measurement (per_type=true) to track different particle types separately - New tests for circular bounds checking and rolling history Cylinder wall angular velocity (dem_wall): - Add angular_velocity field to WallCylinder for rotating cylinder walls - Compute wall surface velocity at contact point via omega × r - Apply Coulomb-limited tangential friction force opposing relative sliding - Compute torque on particles from wall friction Hopper discharge example: - Slot hopper with funnel walls, blocker wall removal on settling - Measurement plane at outlet tracks mass flow rate (Beverloo validation case) - Rolling average flow rate for smoother statistics - Two-stage simulation: filling → flowing Rotating drum example: - Z-axis cylinder with angular_velocity = 1 rad/s (~10 RPM) - Particles inside drum, gravity in -y creates avalanching flow - Demonstrates dynamic angle of repose in rolling/cascading regime - Froude number ≈ 0.005 (well below centrifuging transition) Co-Authored-By: Claude Opus 4.6 --- Cargo.lock | 14 ++ Cargo.toml | 8 + crates/dem_measure_plane/src/lib.rs | 266 +++++++++++++++++++++++++- crates/dem_wall/src/lib.rs | 60 ++++++ examples/hopper_discharge/config.toml | 190 ++++++++++++++++++ examples/hopper_discharge/main.rs | 82 ++++++++ examples/rotating_drum/config.toml | 132 +++++++++++++ examples/rotating_drum/main.rs | 25 +++ 8 files changed, 774 insertions(+), 3 deletions(-) create mode 100644 examples/hopper_discharge/config.toml create mode 100644 examples/hopper_discharge/main.rs create mode 100644 examples/rotating_drum/config.toml create mode 100644 examples/rotating_drum/main.rs diff --git a/Cargo.lock b/Cargo.lock index f47fb00..487d8a3 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -155,6 +155,17 @@ dependencies = [ "mddem_verlet", ] +[[package]] +name = "dem_measure_plane" +version = "0.1.3" +dependencies = [ + "mddem_app", + "mddem_core", + "mddem_print", + "mddem_scheduler", + "serde", +] + [[package]] name = "dem_thermal" version = "0.1.3" @@ -423,6 +434,7 @@ dependencies = [ "dem_atom", "dem_bond", "dem_granular", + "dem_measure_plane", "dem_thermal", "dem_wall", "md_bond", @@ -447,7 +459,9 @@ dependencies = [ name = "mddem-bin" version = "0.1.3" dependencies = [ + "dem_measure_plane", "libffi-sys", + "md_type_rdf", "mddem", "mddem_derive", "mddem_fire", diff --git a/Cargo.toml b/Cargo.toml index c667b19..261850e 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -94,6 +94,14 @@ path = "examples/granular_gas_vdist/main.rs" name = "polymer_chain" path = "examples/polymer_chain/main.rs" +[[example]] +name = "hopper_discharge" +path = "examples/hopper_discharge/main.rs" + +[[example]] +name = "rotating_drum" +path = "examples/rotating_drum/main.rs" + [profile.release] opt-level = 3 diff --git a/crates/dem_measure_plane/src/lib.rs b/crates/dem_measure_plane/src/lib.rs index d4dd980..06b76ed 100644 --- a/crates/dem_measure_plane/src/lib.rs +++ b/crates/dem_measure_plane/src/lib.rs @@ -11,6 +11,10 @@ //! point = [0.1, 0.0, 0.0] //! normal = [1.0, 0.0, 0.0] //! report_interval = 1000 +//! shape = "circular" # optional: "infinite" (default) or "circular" +//! radius = 0.01 # required when shape = "circular" +//! rolling_windows = 5 # optional: number of windows for rolling average +//! per_type = true # optional: report per-atom-type statistics //! ``` //! //! Multiple `[[measure_plane]]` blocks can be defined. Each plane tracks @@ -18,8 +22,11 @@ //! - `crossings_` — total crossing count (positive direction) //! - `flow_rate_` — mass flow rate (kg/s) averaged over `report_interval` //! - `cross_rate_` — particle crossing rate (1/s) averaged over `report_interval` +//! - `avg_flow_rate_` — rolling average mass flow rate (when `rolling_windows > 1`) +//! - `crossings__type` — per-type crossing count (when `per_type = true`) +//! - `flow_rate__type` — per-type mass flow rate (when `per_type = true`) -use std::collections::HashMap; +use std::collections::{HashMap, VecDeque}; use mddem_app::prelude::*; use mddem_core::{Atom, CommResource, Config, RunState}; @@ -33,6 +40,14 @@ fn default_report_interval() -> usize { 1000 } +fn default_shape() -> String { + "infinite".to_string() +} + +fn default_rolling_windows() -> usize { + 1 +} + #[derive(Deserialize, Clone)] #[serde(deny_unknown_fields)] /// TOML `[[measure_plane]]` — defines a measurement plane for particle crossing detection. @@ -46,6 +61,19 @@ pub struct MeasurePlaneDef { /// Report interval in timesteps. Crossing rates are averaged over this window. #[serde(default = "default_report_interval")] pub report_interval: usize, + /// Shape of the measurement plane: "infinite" (default) or "circular". + #[serde(default = "default_shape")] + pub shape: String, + /// Radius of the circular measurement plane [m]. Required when shape = "circular". + #[serde(default)] + pub radius: Option, + /// Number of reporting windows to average over for rolling statistics. + /// 1 = no averaging (default). >1 = rolling average over the last N windows. + #[serde(default = "default_rolling_windows")] + pub rolling_windows: usize, + /// If true, report crossings broken down by atom_type. + #[serde(default)] + pub per_type: bool, } // ── Runtime state ─────────────────────────────────────────────────────────── @@ -58,6 +86,15 @@ struct MeasurePlaneState { normal: [f64; 3], report_interval: usize, + /// Whether this plane is circular (true) or infinite (false). + circular: bool, + /// Squared radius for circular planes (avoids sqrt in hot path). + radius_sq: f64, + /// Number of windows for rolling average. + rolling_windows: usize, + /// Whether to track per-type statistics. + per_type: bool, + /// Signed distance of each tracked particle at the previous step. /// Key: atom tag, Value: signed distance from plane. prev_signed_dist: HashMap, @@ -70,6 +107,14 @@ struct MeasurePlaneState { total_crossings: u64, /// Step at which the last report window started. window_start_step: usize, + + /// Ring buffer of past window values for rolling average: (crossings, mass, window_time). + history: VecDeque<(f64, f64, f64)>, + + /// Per-type crossing count within current window. Key: atom_type. + crossings_by_type: HashMap, + /// Per-type mass within current window. Key: atom_type. + mass_by_type: HashMap, } impl MeasurePlaneState { @@ -81,16 +126,38 @@ impl MeasurePlaneState { } else { [1.0, 0.0, 0.0] // fallback }; + + let circular = def.shape == "circular"; + let radius_sq = if circular { + let r = def.radius.unwrap_or_else(|| { + eprintln!( + "ERROR: measure_plane '{}' has shape='circular' but no 'radius' specified", + def.name + ); + std::process::exit(1); + }); + r * r + } else { + 0.0 + }; + MeasurePlaneState { name: def.name.clone(), point: def.point, normal, report_interval: def.report_interval, + circular, + radius_sq, + rolling_windows: def.rolling_windows.max(1), + per_type: def.per_type, prev_signed_dist: HashMap::new(), crossings_window: 0, mass_window: 0.0, total_crossings: 0, window_start_step: 0, + history: VecDeque::new(), + crossings_by_type: HashMap::new(), + mass_by_type: HashMap::new(), } } @@ -102,6 +169,22 @@ impl MeasurePlaneState { let dz = pos[2] - self.point[2]; dx * self.normal[0] + dy * self.normal[1] + dz * self.normal[2] } + + /// Check if a position is within the circular region of this plane. + /// Always returns true for infinite planes. + #[inline] + fn in_bounds(&self, pos: &[f64; 3]) -> bool { + if !self.circular { + return true; + } + let dx = pos[0] - self.point[0]; + let dy = pos[1] - self.point[1]; + let dz = pos[2] - self.point[2]; + // Project out the normal component to get in-plane displacement + let proj_n = dx * self.normal[0] + dy * self.normal[1] + dz * self.normal[2]; + let in_plane_sq = dx * dx + dy * dy + dz * dz - proj_n * proj_n; + in_plane_sq <= self.radius_sq + } } /// Resource holding all measurement plane states. @@ -122,7 +205,11 @@ impl Plugin for MeasurePlanePlugin { # name = "outlet" # point = [0.1, 0.0, 0.0] # normal = [1.0, 0.0, 0.0] -# report_interval = 1000"#, +# report_interval = 1000 +# shape = "infinite" # or "circular" (requires 'radius') +# radius = 0.01 +# rolling_windows = 1 # >1 for rolling average statistics +# per_type = false # true to report per-atom-type statistics"#, ) } @@ -162,7 +249,14 @@ fn measure_plane_detect_crossings(atoms: Res, mut planes: ResMut, mut planes: ResMut 1 { + plane + .history + .push_back((global_crossings, global_mass, window_time)); + while plane.history.len() > plane.rolling_windows { + plane.history.pop_front(); + } + + let (total_cross, total_mass, total_time) = + plane + .history + .iter() + .fold((0.0, 0.0, 0.0), |(c, m, t), &(wc, wm, wt)| { + (c + wc, m + wm, t + wt) + }); + + let avg_flow_rate = if total_time > 0.0 { + total_mass / total_time + } else { + 0.0 + }; + let avg_cross_rate = if total_time > 0.0 { + total_cross / total_time + } else { + 0.0 + }; + + thermo.set( + &format!("avg_flow_rate_{}", plane.name), + avg_flow_rate, + ); + thermo.set( + &format!("avg_cross_rate_{}", plane.name), + avg_cross_rate, + ); + } + + // Per-type statistics + if plane.per_type { + // Collect all type keys, reduce each + let type_keys: Vec = plane.crossings_by_type.keys().copied().collect(); + for atype in &type_keys { + let local_tc = *plane.crossings_by_type.get(atype).unwrap_or(&0) as f64; + let local_tm = *plane.mass_by_type.get(atype).unwrap_or(&0.0); + let global_tc = comm.all_reduce_sum_f64(local_tc); + let global_tm = comm.all_reduce_sum_f64(local_tm); + + let type_flow_rate = if window_time > 0.0 { + global_tm / window_time + } else { + 0.0 + }; + + thermo.set( + &format!("crossings_{}_type{}", plane.name, atype), + global_tc, + ); + thermo.set( + &format!("flow_rate_{}_type{}", plane.name, atype), + type_flow_rate, + ); + } + + plane.crossings_by_type.clear(); + plane.mass_by_type.clear(); + } + if comm.rank() == 0 { println!( " [{}] crossings={}, mass_flow_rate={:.6e} kg/s, crossing_rate={:.1} /s", @@ -251,6 +420,10 @@ mod tests { point: [0.5, 0.0, 0.0], normal: [1.0, 0.0, 0.0], report_interval: 100, + shape: "infinite".to_string(), + radius: None, + rolling_windows: 1, + per_type: false, }; let state = MeasurePlaneState::new(&def); @@ -269,6 +442,10 @@ mod tests { point: [0.0, 0.0, 0.0], normal: [3.0, 4.0, 0.0], report_interval: 100, + shape: "infinite".to_string(), + radius: None, + rolling_windows: 1, + per_type: false, }; let state = MeasurePlaneState::new(&def); let mag = (state.normal[0].powi(2) + state.normal[1].powi(2) + state.normal[2].powi(2)).sqrt(); @@ -295,4 +472,87 @@ mod tests { let curr_dist3 = -0.1_f64; assert!(!(prev_dist3 <= 0.0 && curr_dist3 > 0.0)); } + + #[test] + fn test_circular_plane_in_bounds() { + let def = MeasurePlaneDef { + name: "circ".to_string(), + point: [0.0, 0.0, 0.5], + normal: [0.0, 0.0, 1.0], + report_interval: 100, + shape: "circular".to_string(), + radius: Some(0.01), + rolling_windows: 1, + per_type: false, + }; + let state = MeasurePlaneState::new(&def); + + // Point within circular region (at center) + assert!(state.in_bounds(&[0.0, 0.0, 0.5])); + // Point within circular region (near edge) + assert!(state.in_bounds(&[0.005, 0.005, 0.5])); + // Point outside circular region + assert!(!state.in_bounds(&[0.02, 0.0, 0.5])); + // Point at different z but within radial bounds + assert!(state.in_bounds(&[0.005, 0.0, 0.6])); + } + + #[test] + fn test_infinite_plane_always_in_bounds() { + let def = MeasurePlaneDef { + name: "inf".to_string(), + point: [0.0, 0.0, 0.0], + normal: [1.0, 0.0, 0.0], + report_interval: 100, + shape: "infinite".to_string(), + radius: None, + rolling_windows: 1, + per_type: false, + }; + let state = MeasurePlaneState::new(&def); + + assert!(state.in_bounds(&[100.0, 200.0, 300.0])); + assert!(state.in_bounds(&[-50.0, -50.0, -50.0])); + } + + #[test] + fn test_rolling_history() { + let def = MeasurePlaneDef { + name: "test".to_string(), + point: [0.0; 3], + normal: [1.0, 0.0, 0.0], + report_interval: 100, + shape: "infinite".to_string(), + radius: None, + rolling_windows: 3, + per_type: false, + }; + let mut state = MeasurePlaneState::new(&def); + assert_eq!(state.rolling_windows, 3); + + // Simulate pushing window data + state.history.push_back((10.0, 0.5, 1.0)); + state.history.push_back((20.0, 1.0, 1.0)); + state.history.push_back((30.0, 1.5, 1.0)); + assert_eq!(state.history.len(), 3); + + // Push one more — oldest should be evicted when we check + state.history.push_back((40.0, 2.0, 1.0)); + while state.history.len() > state.rolling_windows { + state.history.pop_front(); + } + assert_eq!(state.history.len(), 3); + + // Rolling average should be over windows 2,3,4 + let (total_cross, total_mass, total_time) = + state + .history + .iter() + .fold((0.0, 0.0, 0.0), |(c, m, t), &(wc, wm, wt)| { + (c + wc, m + wm, t + wt) + }); + assert!((total_cross - 90.0).abs() < 1e-10); + assert!((total_mass - 4.5).abs() < 1e-10); + assert!((total_time - 3.0).abs() < 1e-10); + } } diff --git a/crates/dem_wall/src/lib.rs b/crates/dem_wall/src/lib.rs index 1b20267..a897830 100644 --- a/crates/dem_wall/src/lib.rs +++ b/crates/dem_wall/src/lib.rs @@ -101,6 +101,10 @@ pub struct WallDef { /// Servo control: adjust velocity to reach target force. #[serde(default)] pub servo: Option, + /// Angular velocity for cylinder walls [rad/s]. Positive = counterclockwise + /// when looking along the +axis direction. + #[serde(default)] + pub angular_velocity: Option, } // ── Runtime types ─────────────────────────────────────────────────────────── @@ -169,6 +173,8 @@ pub struct WallCylinder { pub material_index: usize, pub name: Option, pub force_accumulator: f64, + /// Angular velocity around the cylinder axis [rad/s]. + pub angular_velocity: f64, } /// Runtime sphere wall. @@ -321,6 +327,7 @@ impl Plugin for WallPlugin { let lo = w.lo.unwrap_or(f64::NEG_INFINITY); let hi = w.hi.unwrap_or(f64::INFINITY); let inside = w.inside.unwrap_or(false); + let angular_velocity = w.angular_velocity.unwrap_or(0.0); cylinders.push(WallCylinder { axis, center, @@ -331,6 +338,7 @@ impl Plugin for WallPlugin { material_index: mat_idx, name: w.name.clone(), force_accumulator: 0.0, + angular_velocity, }); } "sphere" => { @@ -733,6 +741,54 @@ pub fn wall_contact_force( atoms.force[i][1] += f_net * ny; atoms.force[i][2] += f_net * nz; cyl_forces[cyl_idx] += f_net; + + // ── Tangential friction from cylinder angular velocity ── + // Compute wall surface velocity at contact point: v_wall = omega × r + // where r is the radial vector from cylinder axis to contact point. + let omega = cyl.angular_velocity; + let mu = material_table.friction_ij[mat_i][wall_mat]; + if omega.abs() > 1e-30 && mu > 0.0 { + let (vw_x, vw_y, vw_z) = match cyl.axis { + // axis=X: omega=(omega,0,0), r=(0,d0,d1) → v=(0,-omega*d1,omega*d0) + 0 => (0.0, -omega * d1, omega * d0), + // axis=Y: omega=(0,omega,0), r=(d0,0,d1) → v=(-omega*d1,0,omega*d0) + 1 => (-omega * d1, 0.0, omega * d0), + // axis=Z: omega=(0,0,omega), r=(d0,d1,0) → v=(-omega*d1,omega*d0,0) + _ => (-omega * d1, omega * d0, 0.0), + }; + + // Relative tangential velocity + let vr_x = atoms.vel[i][0] - vw_x; + let vr_y = atoms.vel[i][1] - vw_y; + let vr_z = atoms.vel[i][2] - vw_z; + let v_n_comp = vr_x * nx + vr_y * ny + vr_z * nz; + let vt_x = vr_x - v_n_comp * nx; + let vt_y = vr_y - v_n_comp * ny; + let vt_z = vr_z - v_n_comp * nz; + let vt_mag = (vt_x * vt_x + vt_y * vt_y + vt_z * vt_z).sqrt(); + + if vt_mag > 1e-30 { + // Coulomb friction: F_t = mu * |F_n|, opposing relative sliding + let f_t = mu * f_net.abs(); + let inv_vt = 1.0 / vt_mag; + let ft_x = -f_t * vt_x * inv_vt; + let ft_y = -f_t * vt_y * inv_vt; + let ft_z = -f_t * vt_z * inv_vt; + + atoms.force[i][0] += ft_x; + atoms.force[i][1] += ft_y; + atoms.force[i][2] += ft_z; + + // Torque on particle: tau = r_contact × F_t + // r_contact from particle center to contact point = -radius * n_hat + let rc_x = -radius * nx; + let rc_y = -radius * ny; + let rc_z = -radius * nz; + dem.torque[i][0] += rc_y * ft_z - rc_z * ft_y; + dem.torque[i][1] += rc_z * ft_x - rc_x * ft_z; + dem.torque[i][2] += rc_x * ft_y - rc_y * ft_x; + } + } } } for (idx, &f) in cyl_forces.iter().enumerate() { @@ -1261,6 +1317,7 @@ mod tests { material_index: 0, name: None, force_accumulator: 0.0, + angular_velocity: 0.0, }); let mut app = App::new(); @@ -1307,6 +1364,7 @@ mod tests { material_index: 0, name: None, force_accumulator: 0.0, + angular_velocity: 0.0, }); let mut app = App::new(); @@ -1442,6 +1500,7 @@ mod tests { material_index: 0, name: None, force_accumulator: 0.0, + angular_velocity: 0.0, }); let mut app = App::new(); @@ -1518,6 +1577,7 @@ mod tests { material_index: 0, name: None, force_accumulator: 0.0, + angular_velocity: 0.0, }); let mut app = App::new(); diff --git a/examples/hopper_discharge/config.toml b/examples/hopper_discharge/config.toml new file mode 100644 index 0000000..6a4d767 --- /dev/null +++ b/examples/hopper_discharge/config.toml @@ -0,0 +1,190 @@ +# ═══════════════════════════════════════════════════════════════════════════════ +# Hopper Discharge with Mass Flow Rate Measurement +# ═══════════════════════════════════════════════════════════════════════════════ +# +# Classic DEM validation case: slot hopper discharge through an orifice. +# Particles fill a funnel under gravity, settle, then the blocker wall is +# removed and particles flow out. A measurement plane at the outlet tracks +# mass flow rate for comparison with the Beverloo correlation. +# +# Geometry (side view, XZ plane): +# +# z=0.08 ┌────────────────────────────┐ ceiling +# │ │ +# │ ◊ ◊ ◊ ◊ ◊ ◊ ◊ ◊ ◊ ◊ ◊ │ particles inserted here +# z=0.06 │ ◊ ◊ ◊ ◊ ◊ ◊ ◊ ◊ ◊ ◊ ◊ │ +# │ │ +# z=0.05 │\ /│ funnel walls (angled) +# │ \ / │ +# │ \ / │ +# z=0.015 │ \ ┌──blocker──┐ / │ ← removed when settled +# │ \ │ │ / │ +# z=0.01 │─────measure plane──────────│ ← tracks mass flow +# │ │ +# z=0.0 └────────────────────────────┘ floor +# x=0 x=0.04 +# +# Beverloo correlation for slot hoppers (2D): +# W = C * ρ_bulk * sqrt(g) * (D - k*d)^(3/2) * L +# where D = orifice width ≈ 0.010 m, d = 0.002 m, L = y-depth = 0.02 m +# +# ═══════════════════════════════════════════════════════════════════════════════ + +[comm] +processors = [1, 1, 1] + +[domain] +x_lo = 0.0 +x_hi = 0.04 +y_lo = 0.0 +y_hi = 0.02 +z_lo = 0.0 +z_hi = 0.08 +x_periodic = false +y_periodic = true # Periodic in y to simulate a 2D-like slot +z_periodic = false + +[neighbor] +skin_fraction = 1.1 +bin_size = 0.005 + +[gravity] +gz = -9.81 + +# ── Material ────────────────────────────────────────────────────────────────── + +[dem] +contact_model = "hertz" + +[[dem.materials]] +name = "glass" +youngs_mod = 8.7e9 +poisson_ratio = 0.3 +restitution = 0.3 # Low restitution for heavy damping (faster settling) +friction = 0.5 + +# ── Particle insertion ──────────────────────────────────────────────────────── + +[[particles.insert]] +material = "glass" +count = 300 +radius = 0.001 # d = 2mm particles +density = 2500.0 +velocity_z = -1.0 # Initial downward velocity for faster filling +region = { type = "block", min = [0.005, 0.0, 0.050], max = [0.035, 0.02, 0.075] } + +# ── Walls ───────────────────────────────────────────────────────────────────── + +# Left side wall (x = 0, normal +x) +[[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 = "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 (z = 0, normal +z) +[[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" + +# Left funnel wall: from (0, 0, 0.05) angled inward to the orifice +# Orifice at z=0.015, centered at x=0.02, half-width = 0.005 (full width D=0.010) +# Left edge at (0.015, 0, 0.015), wall goes from x=0 z=0.05 to x=0.015 z=0.015 +# Direction along wall: (0.015, 0, -0.035), inward normal: (0.035, 0, 0.015) +[[wall]] +point_x = 0.0 +point_y = 0.0 +point_z = 0.05 +normal_x = 0.035 +normal_y = 0.0 +normal_z = 0.015 +material = "glass" +bound_z_low = 0.014 +bound_z_high = 0.052 + +# Right funnel wall: mirror of left +# Right edge at (0.025, 0, 0.015) +# Inward normal: (-0.035, 0, 0.015) +[[wall]] +point_x = 0.04 +point_y = 0.0 +point_z = 0.05 +normal_x = -0.035 +normal_y = 0.0 +normal_z = 0.015 +material = "glass" +bound_z_low = 0.014 +bound_z_high = 0.052 + +# Blocker wall — horizontal at z = 0.015 covering the orifice (removed when settled) +[[wall]] +point_x = 0.0 +point_y = 0.0 +point_z = 0.015 +normal_x = 0.0 +normal_y = 0.0 +normal_z = 1.0 +material = "glass" +name = "blocker" + +# ── Measurement plane ──────────────────────────────────────────────────────── +# +# Placed at z = 0.010, below the orifice, counting particles flowing downward. +# Normal points -z so that particles falling through are counted as positive crossings. + +[[measure_plane]] +name = "outlet" +point = [0.02, 0.0, 0.010] +normal = [0.0, 0.0, -1.0] +report_interval = 5000 +rolling_windows = 3 # Rolling average over 3 reporting windows + +# ── Output ──────────────────────────────────────────────────────────────────── + +[vtp] +interval = 2000 + +[thermo] +columns = ["step", "atoms", "ke", "crossings_outlet", "flow_rate_outlet", "avg_flow_rate_outlet", "walltime"] + +# ── Stages ──────────────────────────────────────────────────────────────────── + +# Stage 1: Fill the hopper — particles settle under gravity onto the blocker +[[run]] +name = "filling" +steps = 500000 +thermo = 5000 + +# Stage 2: Blocker removed — particles discharge through the funnel opening +[[run]] +name = "flowing" +steps = 1000000 +thermo = 5000 diff --git a/examples/hopper_discharge/main.rs b/examples/hopper_discharge/main.rs new file mode 100644 index 0000000..a530ff5 --- /dev/null +++ b/examples/hopper_discharge/main.rs @@ -0,0 +1,82 @@ +//! Hopper discharge with mass flow rate measurement. +//! +//! A classic DEM validation case: particles fill a slot hopper under gravity, +//! then a blocker wall is removed and particles discharge through the orifice. +//! A measurement plane at the outlet tracks the mass flow rate, which can be +//! compared against the Beverloo correlation: +//! +//! W = C * ρ_bulk * sqrt(g) * (D - k*d)^(5/2) +//! +//! where D = orifice width, d = particle diameter, C ≈ 0.58, k ≈ 1.4. +//! +//! ```bash +//! cargo run --example hopper_discharge --no-default-features -- examples/hopper_discharge/config.toml +//! ``` + +use mddem::prelude::*; + +#[derive(Clone, Debug, PartialEq, Default, StageEnum)] +enum Phase { + #[default] + #[stage("filling")] + Filling, + #[stage("flowing")] + Flowing, +} + +fn main() { + let mut app = App::new(); + app.add_plugins(CorePlugins) + .add_plugins(GranularDefaultPlugins) + .add_plugins(GravityPlugin) + .add_plugins(WallPlugin) + .add_plugins(MeasurePlanePlugin) + .add_plugins(StatesPlugin { + initial: Phase::Filling, + }) + .add_plugins(StageAdvancePlugin::::new()); + + app.add_update_system( + check_settled.run_if(in_state(Phase::Filling)), + ScheduleSet::PostFinalIntegration, + ); + + app.start(); +} + +/// Check if particles have settled (KE near zero) and remove the blocker wall. +fn check_settled( + atoms: Res, + run_state: Res, + comm: Res, + mut walls: ResMut, + mut next_state: ResMut>, +) { + let step = run_state.total_cycle; + // Wait at least 1000 steps for particles to start moving, then check every 100 steps + if step < 1000 || step % 100 != 0 { + return; + } + + let nlocal = atoms.nlocal as usize; + let local_ke: f64 = (0..nlocal) + .map(|i| { + let vx = atoms.vel[i][0]; + let vy = atoms.vel[i][1]; + let vz = atoms.vel[i][2]; + 0.5 * atoms.mass[i] * (vx * vx + vy * vy + vz * vz) + }) + .sum(); + let global_ke = comm.all_reduce_sum_f64(local_ke); + + if global_ke < 1e-5 { + walls.deactivate_by_name("blocker"); + next_state.set(Phase::Flowing); + if comm.rank() == 0 { + println!( + "Step {}: KE = {:.3e} J — particles settled, removing blocker wall", + step, global_ke + ); + } + } +} diff --git a/examples/rotating_drum/config.toml b/examples/rotating_drum/config.toml new file mode 100644 index 0000000..c069e3e --- /dev/null +++ b/examples/rotating_drum/config.toml @@ -0,0 +1,132 @@ +# ═══════════════════════════════════════════════════════════════════════════════ +# Rotating Drum DEM Simulation +# ═══════════════════════════════════════════════════════════════════════════════ +# +# A classic industrial DEM test case. A cylindrical drum rotates at constant +# angular velocity, dragging particles upward via friction. Gravity pulls +# particles back down, creating avalanching flow and a dynamic angle of repose. +# +# Geometry (cross-section, XY plane): +# +# ╭───────────────╮ +# ╱ ╲ +# │ ◊ ◊ ◊ │ +# │ ◊ ◊ ◊ ◊ ◊ │ ← drum rotates CCW +# │ ◊ ◊ ◊ ◊ ◊ ◊ ◊ │ at 1 rad/s (~10 RPM) +# │ ◊ ◊ ◊ ◊ ◊ ◊ ◊ │ +# │ ◊ ◊ ◊ ◊ ◊ │ +# ╲ ╱ +# ╰───────────────╯ +# +# The drum is a Z-axis cylinder centered at (0, 0) with radius 0.05 m. +# Flat end caps close the drum at z = 0 and z = 0.02. +# Particles are 2 mm diameter glass beads inserted in the lower half. +# +# At ~10 RPM (rolling/cascading regime), particles exhibit: +# - A free surface with dynamic angle of repose (~25-35° for glass beads) +# - Continuous avalanching from the raised side +# - Minimal centrifugal segregation at low speeds +# +# ═══════════════════════════════════════════════════════════════════════════════ + +[comm] +processors = [1, 1, 1] + +[domain] +# Box must fully enclose the drum cylinder +x_lo = -0.055 +x_hi = 0.055 +y_lo = -0.055 +y_hi = 0.055 +z_lo = 0.0 +z_hi = 0.02 +x_periodic = false +y_periodic = false +z_periodic = false + +[neighbor] +skin_fraction = 1.1 +bin_size = 0.005 + +[gravity] +gx = 0.0 +gy = -9.81 # Gravity in -y direction (downward in cross-section view) +gz = 0.0 + +# ── Material ────────────────────────────────────────────────────────────────── + +[dem] +contact_model = "hertz" + +[[dem.materials]] +name = "glass" +youngs_mod = 5.0e8 # Reduced Young's modulus for faster simulation +poisson_ratio = 0.3 +restitution = 0.4 +friction = 0.5 # Sliding friction — important for drum dynamics + +# ── Particle insertion ──────────────────────────────────────────────────────── +# +# Insert particles in the lower half of the drum (y < 0). +# They will settle under gravity before the drum rotation takes effect. + +[[particles.insert]] +material = "glass" +count = 400 +radius = 0.001 # d = 2 mm particles +density = 2500.0 # Glass density +velocity = 0.1 # Small random velocity for initial spreading +region = { type = "block", min = [-0.035, -0.045, 0.002], max = [0.035, -0.005, 0.018] } + +# ── Walls ───────────────────────────────────────────────────────────────────── + +# Rotating cylinder wall (Z-axis, particles inside) +# angular_velocity = 1.0 rad/s ≈ 9.5 RPM (rolling/cascading regime) +# Froude number = ω²R/g = 1.0² × 0.05 / 9.81 ≈ 0.005 (well below centrifuging) +[[wall]] +type = "cylinder" +axis = "z" +center = [0.0, 0.0] +radius = 0.05 +lo = 0.0 +hi = 0.02 +inside = true +material = "glass" +name = "drum" +angular_velocity = 1.0 + +# Bottom end cap (z = 0, normal +z) +[[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 = "bottom_cap" + +# Top end cap (z = 0.02, normal -z) +[[wall]] +point_x = 0.0 +point_y = 0.0 +point_z = 0.02 +normal_x = 0.0 +normal_y = 0.0 +normal_z = -1.0 +material = "glass" +name = "top_cap" + +# ── Output ──────────────────────────────────────────────────────────────────── + +[vtp] +interval = 2000 + +[thermo] +columns = ["step", "atoms", "ke", "walltime", "stepps"] + +# ── Run ─────────────────────────────────────────────────────────────────────── + +[run] +steps = 2000000 +thermo = 5000 diff --git a/examples/rotating_drum/main.rs b/examples/rotating_drum/main.rs new file mode 100644 index 0000000..a6c2781 --- /dev/null +++ b/examples/rotating_drum/main.rs @@ -0,0 +1,25 @@ +//! Rotating drum DEM simulation — a classic industrial test case. +//! +//! A cylinder wall rotates at constant angular velocity with particles inside. +//! Friction between the drum wall and particles drags particles upward, +//! creating avalanching behavior and a dynamic angle of repose. +//! +//! This demonstrates: +//! - Cylinder walls with prescribed angular velocity +//! - Tangential friction from rotating walls +//! - Gravity-driven granular flow inside a drum +//! +//! ```bash +//! cargo run --example rotating_drum --no-default-features -- examples/rotating_drum/config.toml +//! ``` + +use mddem::prelude::*; + +fn main() { + let mut app = App::new(); + app.add_plugins(CorePlugins) + .add_plugins(GranularDefaultPlugins) + .add_plugins(GravityPlugin) + .add_plugins(WallPlugin); + app.start(); +}