diff --git a/ROADMAP.md b/ROADMAP.md index 74ee6a2..d4dc3b4 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -36,7 +36,7 @@ pressure) is a degree of freedom of the existing DOF framework, not a separate s | M7 Assemblies, contact and friction | Needs M15. Part instances from SolidCraft occurrences; tie between non-matching meshes (mortar); implicit contact: surface-to-surface mortar and node-to-surface, penalty and augmented Lagrange, self-contact, automatic pair detection (the model's contact definitions, `interaction.contact.*` with `find`, the deck mapping, `fieldcraft_fem::contact`'s closest-point projection and bucket grid, and the explicit solver's segmentation and friction return map are in place to reuse; the implicit solver refuses contact by name until then); friction: Coulomb with stick–slip, static/kinetic, anisotropic, pressure-, rate- and temperature-dependent coefficients; contact output (pressure, slip, status); thermal contact conductance and frictional heating; wear (Archard); bolt pretension, connectors with behaviour (elastic, damped, stops, locks); amplitudes, multiple steps, restart; nonlinear steps in all three decks. Verified by Hertz (sphere and cylinder), the Cattaneo–Mindlin partial-slip problem, the sliding block, NAFEMS contact benchmarks | planned | | M8 Materials, explicit, shells, composites | Plasticity, hyperelasticity, damage data per solver; explicit dynamics decks (Abaqus/Explicit, LS-DYNA); midsurface shells; composite layups and material orientations, continuum shells on swept hex; hybrid (mixed u–p) hex8/hex20/tet10 for nearly incompressible materials (done); heat transfer (solids done: steady and θ-method transient conduction with flux, convection and sources, thermal strain in static steps (solids and 2D) and explicit dynamics (with plasticity), coupled temperature–displacement, `solve.heat`; deck output, 2D/shell/beam conduction planned) | started | | M9 Post-processing and deck readers | Results readers for public formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`); XY and path plots, cuts, display groups, free-body forces; deck readers that keep unknown keywords; mesh convergence studies | planned | -| M10 Explicit dynamics | Shared element and material kernels with the implicit solver (`fem`); central-difference solver with lumped mass, stable time step, one-point hex8 with hourglass control and nodal-pressure tet4, bulk viscosity, energy accounting (X1 core: done — `fieldcraft-explicit`, `explicit_dynamic` steps with initial velocities, prescribed motion and planar rigid walls, `solve.explicit`/`solve.history`, small-strain linear elastic; bar wave, rigid-wall impact, energy balance and cantilever-period tests); large-deformation materials (X2: done — updated Lagrangian elements, Jaumann/Hughes–Winget and Green–Naghdi/Flanagan–Taylor rates, J2 bilinear/table plasticity, Johnson–Cook strength with heating and damage, erosion; materials gain `plastic`, `damage`, `erosion` with `model.material.plasticity`/`damage`; Dienes shear, rigid rotation, J2 uniaxial and Taylor-bar tests); contact (X3: done — `interaction.contact.pair`/`general`/`find`, segment-to-segment penalty contact both ways with polygon clipping (passes the contact patch test on non-matching meshes), pairs, general (all exterior faces) and self contact, bucket-grid search, deterministic on any thread count, a penalty from the masses and the stable increment that never shortens it, Coulomb friction with exponential decay by a return map, friction on rigid walls, contact energy and friction dissipation in the balance, initial penetrations adjusted or refused; Abaqus `*CONTACT PAIR`/`*CONTACT`/`*FRICTION` and LS-DYNA automatic contacts in the decks; verified by colliding bars, sliding blocks, the patch test and Hertz's sphere impact); hyperelastic materials on the shared material point and mass scaling next; runs on hex and tet meshes | started | +| M10 Explicit dynamics | Shared element and material kernels with the implicit solver (`fem`); central-difference solver with lumped mass, stable time step, one-point hex8 with hourglass control and nodal-pressure tet4, bulk viscosity, energy accounting (X1 core: done — `fieldcraft-explicit`, `explicit_dynamic` steps with initial velocities, prescribed motion and planar rigid walls, `solve.explicit`/`solve.history`, small-strain linear elastic; bar wave, rigid-wall impact, energy balance and cantilever-period tests); large-deformation materials (X2: done — updated Lagrangian elements, Jaumann/Hughes–Winget and Green–Naghdi/Flanagan–Taylor rates, J2 bilinear/table plasticity, Johnson–Cook strength with heating and damage, erosion; materials gain `plastic`, `damage`, `erosion` with `model.material.plasticity`/`damage`; Dienes shear, rigid rotation, J2 uniaxial and Taylor-bar tests); contact (X3: done — `interaction.contact.pair`/`general`/`find`, segment-to-segment penalty contact both ways with polygon clipping (passes the contact patch test on non-matching meshes), pairs, general (all exterior faces) and self contact, bucket-grid search, deterministic on any thread count, a penalty from the masses and the stable increment that never shortens it, Coulomb friction with exponential decay by a return map, friction on rigid walls, contact energy and friction dissipation in the balance, initial penetrations adjusted or refused; Abaqus `*CONTACT PAIR`/`*CONTACT`/`*FRICTION` and LS-DYNA automatic contacts in the decks; verified by colliding bars, sliding blocks, the patch test and Hertz's sphere impact); hyperelastic materials (neo-Hookean, Mooney–Rivlin, Yeoh) on the shared material point with average nodal pressure tets, and opt-in mass scaling with the added mass reported: done; runs on hex and tet meshes | started | | M11 Equations of state and shells | Polynomial, ideal-gas, Mie–Grüneisen, JWL with programmed burn, tabulated; energy coupling; shock verification (Rankine–Hugoniot, Sod, Noh, Sedov); Belytschko–Lin–Tsay shells with self-contact on quad meshes | planned | | M12 CFD: lattice Boltzmann | D3Q19/D3Q27 on grids voxelised from the face-tagged surface; TRT, regularised and cumulant collision; interpolated walls, inlets and outlets; LES (Smagorinsky, WALE) with wall models; thermal LBM and conjugate heat transfer with the structure's thermal solver; non-Newtonian viscosity (power-law, Carreau–Yasuda, Cross, Herschel–Bulkley, Bingham) and temperature-dependent viscosity; porous media (Brinkman, Darcy–Forchheimer); multiphase and free surface (M13); validation by the lid-driven cavity, backward-facing step, cylinder wake (Strouhal) and turbulent channel; octree refinement shared with the hex mesher; wgpu backend (deterministic per device, checked against the CPU) (CPU core: done — `fieldcraft-lbm` with D3Q19 TRT and Guo forcing, halfway bounce-back and moving walls, slip and periodic faces, Zou–He velocity and pressure faces, exact voxelisation of obstacle bodies, momentum-exchange forces per face; fluid materials, fluid domains and fluid steps in the model; `model.fluid_domain.*`, `solve.fluid`; verified by Poiseuille channel and pipe, Couette, Taylor–Green decay, sphere-array Stokes drag and mass conservation; Schäfer–Turek 2D-1 as a benchmark) | started | | M13 Fluid–structure interaction | Immersed-boundary coupling of the lattice Boltzmann fluid with the explicit or implicit structure; interface force transfer; sub-cycling; Aitken then interface quasi-Newton coupling; free-surface flow | planned | diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs index 3437d64..59e10cb 100644 --- a/crates/engine/src/cmd/define.rs +++ b/crates/engine/src/cmd/define.rs @@ -292,7 +292,7 @@ fn fluid_boundary(s: &mut Session, p: &Value) -> Result { pub const STEP_CREATE: CommandSpec = CommandSpec { id: "model.step.create", label: "New Step", - params: r#"{"name": "Bend", "procedure": {"type": "static"} | {"type": "nonlinear_static", "nonlinear_geometry": true, "initial_increment": 0.1, "min_increment": 1e-5, "max_increment": 1.0 (fractions of the step), "max_increments": 100, "max_iterations": 16 (per increment), "tolerances": {"force": 1e-6, "displacement": 1e-5, "energy": 1e-10}, "line_search": false, "load_path": [[0, 0], [0.5, 1], [1, 0]] (optional: the load factor against the step fraction, piecewise linear from [0, 0] to fraction 1; increments end at its points)} (loads and prescribed displacements are the load factor times their values; it ramps from 0 to 1 without a load path) | {"type": "frequency", "modes": 10, "min_hz": 0, "max_hz": 500} | {"type": "fluid", "domain": "Flow", "end_time": 1.0 | "steps": 1000, "time_step": (default: chosen for Mach 0.087 and τ ≤ 1), "frames": 10, "magic": 0.1875 (the two-relaxation-time Λ)} | {"type": "explicit_dynamic", "end_time": 0.001, "dt_scale": 0.9, "max_increments": 10000000, "output_interval", "history_interval", "history_nodes": "set", "initial_velocities": [{"set", "velocity": [vx, vy, vz]}], "rigid_walls": [{"set", "point", "normal", "friction": {"static": 0.3, "kinetic": 0.2, "decay": 10} (Coulomb; default frictionless)}], "hourglass": {"type": "viscous|stiffness", "coefficient": 0.1}, "bulk_viscosity": {"linear": 0.06, "quadratic": 1.5}}, "output": ["displacement", "reaction_force", "stress", "strain"] (default by procedure), "boundary_conditions": [names], "loads": [names] (active in this step; nothing carries over between steps), "before": "step" (default: last)}"#, + params: r#"{"name": "Bend", "procedure": {"type": "static"} | {"type": "nonlinear_static", "nonlinear_geometry": true, "initial_increment": 0.1, "min_increment": 1e-5, "max_increment": 1.0 (fractions of the step), "max_increments": 100, "max_iterations": 16 (per increment), "tolerances": {"force": 1e-6, "displacement": 1e-5, "energy": 1e-10}, "line_search": false, "load_path": [[0, 0], [0.5, 1], [1, 0]] (optional: the load factor against the step fraction, piecewise linear from [0, 0] to fraction 1; increments end at its points)} (loads and prescribed displacements are the load factor times their values; it ramps from 0 to 1 without a load path) | {"type": "frequency", "modes": 10, "min_hz": 0, "max_hz": 500} | {"type": "fluid", "domain": "Flow", "end_time": 1.0 | "steps": 1000, "time_step": (default: chosen for Mach 0.087 and τ ≤ 1), "frames": 10, "magic": 0.1875 (the two-relaxation-time Λ)} | {"type": "explicit_dynamic", "end_time": 0.001, "dt_scale": 0.9, "max_increments": 10000000, "output_interval", "history_interval", "history_nodes": "set", "initial_velocities": [{"set", "velocity": [vx, vy, vz]}], "rigid_walls": [{"set", "point", "normal", "friction": {"static": 0.3, "kinetic": 0.2, "decay": 10} (Coulomb; default frictionless)}], "hourglass": {"type": "viscous|stiffness", "coefficient": 0.1}, "bulk_viscosity": {"linear": 0.06, "quadratic": 1.5}, "mass_scaling": {"target_increment": 1e-6, "set": "…" (optional, default every element)} (opt-in: elements below the target increment get the density that reaches it; the added mass is reported)}, "output": ["displacement", "reaction_force", "stress", "strain"] (default by procedure), "boundary_conditions": [names], "loads": [names] (active in this step; nothing carries over between steps), "before": "step" (default: last)}"#, form: Some(super::forms::STEP_CREATE), kind: CommandKind::Model, run: step_create, diff --git a/crates/engine/src/cmd/forms.rs b/crates/engine/src/cmd/forms.rs index 65d0abd..dfeee85 100644 --- a/crates/engine/src/cmd/forms.rs +++ b/crates/engine/src/cmd/forms.rs @@ -505,6 +505,11 @@ const PROCEDURE: Field = variants( group("bulk_viscosity", "Bulk viscosity", &[positive("linear", "Linear", D::None), positive("quadratic", "Quadratic", D::None)]), flag("nonlinear_geometry", "Nonlinear geometry").def("true"), choice("stress_rate", "Stress rate", &[("jaumann", "Jaumann"), ("green_naghdi", "Green–Naghdi")]), + group( + "mass_scaling", + "Mass scaling", + &[positive("target_increment", "Target increment", D::Time).req(), reference("set", "Set", C::Set)], + ), ], ), variant( diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index 929dad0..737f534 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -355,6 +355,7 @@ fn explicit_summary(r: &ExplicitResult) -> Value { "history_samples": r.history.time.len(), "eroded": r.eroded, "contacts": r.contacts, + "mass_scaling": r.mass_scaling, "warnings": r.warnings, }) } diff --git a/crates/explicit/src/elements.rs b/crates/explicit/src/elements.rs index dee47e5..53da037 100644 --- a/crates/explicit/src/elements.rs +++ b/crates/explicit/src/elements.rs @@ -547,6 +547,8 @@ fn hexahedron(sys: &System, h: &crate::system::Hex, m: &Motion, es: &mut Element let (d, _) = rotation::split(&g); let d = mechanical(d, h.thermal_rate * m.dt); let (rho, c) = current(mat, h.kinematics.volume, now.volume); + // Mass scaling raises the density the dynamics see (not the one adiabatic heating uses). + let (rho_s, c_s) = (rho * h.mass_factor, c / h.mass_factor.sqrt()); let stiffer = match &mat.hyperelastic { Some(model) => { let mut f = h.kinematics.gradient(&e.u); @@ -568,7 +570,7 @@ fn hexahedron(sys: &System, h: &crate::system::Hex, m: &Motion, es: &mut Element // geometry). let (gamma, scales) = if sys.nonlinear { let gamma = now.hourglass_vectors(&x_now); - let scales = hex_scales(&now, &gamma, &h.mass, mat, rho, hourglass); + let scales = hex_scales(&now, &gamma, &h.mass, mat, rho_s, hourglass); (gamma, scales) } else { (h.gamma, h.rest) @@ -593,7 +595,7 @@ fn hexahedron(sys: &System, h: &crate::system::Hex, m: &Motion, es: &mut Element } } } - let element = Element { now: &now, mid: &mid, density: rho, wave_speed: c, length: 2.0 * c / (stiffer * scales.omega2).sqrt() }; + let element = Element { now: &now, mid: &mid, density: rho_s, wave_speed: c_s, length: 2.0 * c_s / (stiffer * scales.omega2).sqrt() }; let mut forces = element.finish(&es.point.stress, &e.v, m.bulk); forces.hourglass = Some(f); Ok(Outcome::Forces(forces)) @@ -632,6 +634,7 @@ fn tetrahedron( let thermal = t.thermal_rate * m.dt; let d = mechanical(d, thermal); let (rho, c) = current(mat, t.kinematics.volume, inc.now.volume); + let (rho_s, c_s) = (rho * t.mass_factor, c / t.mass_factor.sqrt()); let stiffer = match &mat.hyperelastic { // The average nodal pressure of a decoupled material (Bonet & Burton 1998, at finite // strain): each node's pressure U′(J_a) from its tetrahedra's volume ratio, the element's @@ -660,7 +663,7 @@ fn tetrahedron( return Ok(Outcome::Eroded(failure)); } let omega2 = stiffer * if sys.nonlinear { frequency_bound(&inc.now, inc.now.gradient_norm(), mat.modulus, &t.mass) } else { t.omega2 }; - let element = Element { now: &inc.now, mid: &inc.mid, density: rho, wave_speed: c, length: 2.0 * c / omega2.sqrt() }; + let element = Element { now: &inc.now, mid: &inc.mid, density: rho_s, wave_speed: c_s, length: 2.0 * c_s / omega2.sqrt() }; Ok(Outcome::Forces(element.finish(&es.point.stress, &inc.v, m.bulk))) } diff --git a/crates/explicit/src/integrate.rs b/crates/explicit/src/integrate.rs index e8e48fc..27f9239 100644 --- a/crates/explicit/src/integrate.rs +++ b/crates/explicit/src/integrate.rs @@ -201,10 +201,14 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res let mut critical: Vec = sys .hexes .iter() - .map(|h| (h.label, h.length, h.material)) - .chain(sys.tets.iter().map(|t| (t.label, t.length, t.material))) - .chain(sys.rods.iter().map(|r| (r.label, r.length, r.material))) - .map(|(label, length, m)| ControllingElement { label, increment: length / sys.materials.get(m).map_or(f64::INFINITY, |m| m.wave_speed) }) + .map(|h| (h.label, h.length, h.material, h.mass_factor)) + .chain(sys.tets.iter().map(|t| (t.label, t.length, t.material, t.mass_factor))) + .chain(sys.rods.iter().map(|r| (r.label, r.length, r.material, 1.0))) + .map(|(label, length, m, f)| ControllingElement { + label, + // Mass scaling slows the element's waves by √f. + increment: length * f.sqrt() / sys.materials.get(m).map_or(f64::INFINITY, |m| m.wave_speed), + }) .collect(); critical.sort_by(|a, b| a.increment.total_cmp(&b.increment).then(a.label.cmp(&b.label))); critical.truncate(CONTROLLING); @@ -411,6 +415,10 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res energy.external += external; energy.constraint += dt_n * constraint; energy.kinetic = 0.5 * kinetic; + if !sys.added_mass.is_empty() { + let per_node = v_now.as_chunks::().0.iter().zip(&sys.added_mass); + energy.added_kinetic = 0.5 * per_node.map(|(v, m)| m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2])).sum::(); + } if !ws.contact_force.is_empty() { let work: f64 = du.iter().zip(ws.contact_force.iter().zip(&previous.contact_force)).map(|(d, (c, c0))| d * (c0 + c)).sum(); energy.contact -= 0.5 * work; @@ -483,6 +491,7 @@ pub fn solve_explicit(model: &Model, step: Option<&str>, limits: &Limits) -> Res initial_kinetic, eroded, contacts: summaries, + mass_scaling: sys.mass_scaling.clone(), warnings, }) } diff --git a/crates/explicit/src/lib.rs b/crates/explicit/src/lib.rs index 42f0cb0..edd8eae 100644 --- a/crates/explicit/src/lib.rs +++ b/crates/explicit/src/lib.rs @@ -46,7 +46,7 @@ pub mod system; pub use elements::Failure; pub use integrate::solve_explicit; -pub use result::{ContactSummary, ControllingElement, Energies, ErodedElement, ExplicitResult, Frame, History}; +pub use result::{ContactSummary, ControllingElement, Energies, ErodedElement, ExplicitResult, Frame, History, MassScalingSummary}; #[derive(Debug, Clone, PartialEq, thiserror::Error)] pub enum Error { diff --git a/crates/explicit/src/result.rs b/crates/explicit/src/result.rs index 9067a43..70179cb 100644 --- a/crates/explicit/src/result.rs +++ b/crates/explicit/src/result.rs @@ -25,6 +25,8 @@ pub struct Energies { /// Energy taken by contact: the work the model does against the contact forces (the penalty /// energy stored in the contact springs, plus what friction dissipates there). pub contact: f64, + /// The kinetic energy of the mass added by mass scaling (part of `kinetic`). + pub added_kinetic: f64, /// Friction dissipation, in contact and on rigid walls (part of `contact` and of /// `constraint`). pub friction: f64, @@ -106,6 +108,19 @@ pub struct History { pub velocity: Vec>, } +/// What mass scaling added to a run. +#[derive(Debug, Clone, Default, PartialEq, Serialize)] +pub struct MassScalingSummary { + pub target_increment: f64, + /// Elements whose density was raised. + pub elements: usize, + pub added_mass: f64, + /// The added mass over the model's mass (with it). + pub fraction: f64, + /// The added mass per section, by name. + pub sections: Vec<(String, f64)>, +} + /// What a contact did over the run. #[derive(Debug, Clone, Default, PartialEq, Serialize)] pub struct ContactSummary { @@ -158,5 +173,7 @@ pub struct ExplicitResult { pub eroded: Vec, /// Per contact of the model, in name order. pub contacts: Vec, + /// What mass scaling added, when the step asks for it. + pub mass_scaling: Option, pub warnings: Vec, } diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs index dd4ffcf..c344484 100644 --- a/crates/explicit/src/system.rs +++ b/crates/explicit/src/system.rs @@ -68,6 +68,8 @@ pub struct Hex { pub length: f64, /// The rate of its thermal strain, α ΔT / end time (every normal direction). pub thermal_rate: f64, + /// Its density's factor from mass scaling (1 without). + pub mass_factor: f64, } /// A nodal-pressure tetrahedron. @@ -84,6 +86,8 @@ pub struct Tet { pub length: f64, /// The rate of its thermal strain, α ΔT / end time (every normal direction). pub thermal_rate: f64, + /// Its density's factor from mass scaling (1 without). + pub mass_factor: f64, } /// A fixed rigid plane. @@ -120,6 +124,9 @@ pub struct System { /// The springs' and dashpots' stable increment on the lumped mass, and a link that sets it, /// when there are any. pub discrete_increment: Option<(f64, String)>, + /// What mass scaling added: per node, and per section (in name order). + pub added_mass: Vec, + pub mass_scaling: Option, /// Prescribed velocity per degree of freedom (`None` = free). pub prescribed: Vec>, pub walls: Vec, @@ -277,6 +284,8 @@ impl System { tets: Vec::new(), rods: Vec::new(), mass: vec![0.0; n], + added_mass: Vec::new(), + mass_scaling: None, discrete: Discrete::default(), discrete_increment: None, prescribed: vec![None; dofs], @@ -322,7 +331,19 @@ impl System { let rest = hex_scales(&k, &gamma, &mass, mat, rho, settings.hourglass); sys.lump(&nodes, &mass); let length = 2.0 * c / rest.omega2.sqrt(); - sys.hexes.push(Hex { label, index, nodes, material: mi, kinematics: k, gamma, mass, rest, length, thermal_rate: 0.0 }); + sys.hexes.push(Hex { + label, + index, + nodes, + material: mi, + kinematics: k, + gamma, + mass, + rest, + length, + thermal_rate: 0.0, + mass_factor: 1.0, + }); } ElementKind::Tet4 => { let nodes: [usize; 4] = idx.try_into().map_err(|_| invalid(format!("element {label} does not have 4 nodes")))?; @@ -341,6 +362,7 @@ impl System { omega2, length: 2.0 * c / omega2.sqrt(), thermal_rate: 0.0, + mass_factor: 1.0, }); } ElementKind::Rod2 | ElementKind::Rod3 => { @@ -369,6 +391,9 @@ impl System { other => return Err(Error::Internal(format!("{} element kind", other.name()))), } } + if let Some(scaling) = &settings.mass_scaling { + sys.scale_mass(model, resolved, scaling, settings)?; + } sys.discrete = Discrete::new(resolved, sys.dofs(), &sys.xyz, |l| sys.node(l))?; sys.discrete_increment = sys.discrete.critical_increment(&sys.dof_mass(), sys.nonlinear); sys.prescribe(model, resolved, step, settings.end_time)?; @@ -379,6 +404,90 @@ impl System { Ok(sys) } + /// Mass scaling to a target increment (see [`fieldcraft_model::MassScaling`]): each solid + /// element (of the set) whose stable increment times the step's scale is below the target + /// gets its density times (target / increment)², its lumped masses and frequency bound with + /// it; the added mass is kept per node and per section. + fn scale_mass(&mut self, model: &Model, resolved: &Resolved, scaling: &fieldcraft_model::MassScaling, settings: &ExplicitDynamic) -> Result<()> { + let chosen: Option<&[u32]> = match &scaling.set { + Some(s) => Some(resolved.sets.get(s).map(|r| r.elements.as_slice()).ok_or_else(|| invalid(format!("set `{s}` is not resolved")))?), + None => None, + }; + let wanted = |label: u32| chosen.is_none_or(|l| l.binary_search(&label).is_ok()); + let factor = |length: f64, material: usize, materials: &[Material]| -> f64 { + let c = materials.get(material).map_or(f64::INFINITY, |m| m.wave_speed); + let dt = settings.dt_scale * length / c; + if dt > 0.0 && dt < scaling.target_increment { (scaling.target_increment / dt).powi(2) } else { 1.0 } + }; + let hourglass = fieldcraft_model::Hourglass { kind: settings.hourglass.kind, coefficient: settings.hourglass.coefficient }; + let mut added = vec![0.0; self.labels.len()]; + let mut by_element: Vec<(u32, f64)> = Vec::new(); + for h in self.hexes.iter_mut().filter(|h| wanted(h.label)) { + let f = factor(h.length, h.material, &self.materials); + if f == 1.0 { + continue; + } + let extra = h.mass.map(|m| m * (f - 1.0)); + for (n, x) in h.nodes.iter().zip(extra) { + if let Some(a) = added.get_mut(*n) { + *a += x; + } + } + h.mass = h.mass.map(|m| m * f); + h.mass_factor = f; + if let Some(mat) = self.materials.get(h.material) { + h.rest = hex_scales(&h.kinematics, &h.gamma, &h.mass, mat, mat.density * f, hourglass); + } + by_element.push((h.label, extra.iter().sum())); + } + for t in self.tets.iter_mut().filter(|t| wanted(t.label)) { + let f = factor(t.length, t.material, &self.materials); + if f == 1.0 { + continue; + } + let extra = t.mass.map(|m| m * (f - 1.0)); + for (n, x) in t.nodes.iter().zip(extra) { + if let Some(a) = added.get_mut(*n) { + *a += x; + } + } + t.mass = t.mass.map(|m| m * f); + t.omega2 /= f; + t.mass_factor = f; + by_element.push((t.label, extra.iter().sum())); + } + for (m, a) in self.mass.iter_mut().zip(&added) { + *m += a; + } + let total_added: f64 = added.iter().sum(); + let total: f64 = self.mass.iter().sum(); + let mut sections: Vec<(String, f64)> = Vec::new(); + for (name, section) in model.sections() { + let elements = resolved.sets.get(§ion.set).map(|r| r.elements.as_slice()).unwrap_or_default(); + let mass: f64 = by_element.iter().filter(|(l, _)| elements.binary_search(l).is_ok()).map(|(_, m)| m).sum(); + if mass > 0.0 { + sections.push((name.clone(), mass)); + } + } + if total_added > 0.0 { + self.warnings.push(format!( + "mass scaling to an increment of {:.3e} added {total_added:.4e} of mass ({:.2}% of the model's) to {} elements", + scaling.target_increment, + 100.0 * total_added / total.max(f64::MIN_POSITIVE), + by_element.len() + )); + } + self.mass_scaling = Some(crate::result::MassScalingSummary { + target_increment: scaling.target_increment, + elements: by_element.len(), + added_mass: total_added, + fraction: total_added / total.max(f64::MIN_POSITIVE), + sections, + }); + self.added_mass = added; + Ok(()) + } + /// Refuse what hyperelastic elements cannot take yet: small-displacement steps (the /// material is a finite-strain law) and thermal strain. fn check_hyperelastic(&self) -> Result<()> { diff --git a/crates/explicit/tests/one_d.rs b/crates/explicit/tests/one_d.rs index 084f6e7..6babf04 100644 --- a/crates/explicit/tests/one_d.rs +++ b/crates/explicit/tests/one_d.rs @@ -305,3 +305,55 @@ fn whirling_mass_on_a_corotational_rod() { // A quarter turn: the mass ends near (0, L). assert!((l + end[0]).abs() < 1e-2 && (end[1] - l).abs() < 1e-2, "{end:?}"); } + +/// Mass scaling to twice a uniform bar's stable increment quadruples its density: the added mass +/// is three times the bar's, the stable increment doubles with the halved wave speed +/// c = √(E/ρ), and the kinetic energy of the added mass is three quarters of the total. +#[test] +fn mass_scaling_slows_the_waves_as_the_added_density_says() { + use fieldcraft_explicit::{Limits, solve_explicit}; + use fieldcraft_model::*; + let build = |scaled: Option| { + let mut m = Model::default(); + let n = 40; + for i in 0..=n { + for (j, (y, z)) in [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)].iter().enumerate() { + m.set_node((4 * i + j + 1) as u32, [i as f64, *y, *z]).unwrap(); + } + } + for i in 0..n { + let a = |j: usize| (4 * i + j + 1) as u32; + let b = |j: usize| (4 * (i + 1) + j + 1) as u32; + m.set_element(i as u32 + 1, ElementKind::Hex8, vec![a(0), b(0), b(1), a(1), a(3), b(3), b(2), a(2)]).unwrap(); + } + m.add_set("all", Set { elements: (1..=n as u32).collect(), ..Set::default() }).unwrap(); + m.add_material("Mat", Material { elastic: Some(Elastic::Isotropic { young: 1.0, poisson: 0.0 }), density: Some(1.0), ..Material::default() }) + .unwrap(); + let formulation = Formulation { integration: Integration::Reduced, ..Formulation::default() }; + m.add_section("bar", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation }).unwrap(); + let mut step = ExplicitDynamic::new(10.0); + step.initial_velocities = vec![InitialVelocity { set: "all".into(), velocity: [0.01, 0.0, 0.0] }]; + step.mass_scaling = scaled.map(|t| MassScaling { target_increment: t, set: None }); + m.add_step(Step::new("Run", Procedure::ExplicitDynamic(step)), None).unwrap(); + solve_explicit(&m, None, &Limits::default()).unwrap() + }; + let plain = build(None); + let scaled = build(Some(2.0 * plain.initial_increment)); + let s = scaled.mass_scaling.as_ref().unwrap(); + eprintln!( + "mass scaling: {s:?}; increments {} vs {}; Δt {:.4e} vs {:.4e}", + scaled.increments, plain.increments, scaled.initial_increment, plain.initial_increment + ); + assert_eq!(s.elements, 40); + assert!((s.added_mass / plain.total_mass - 3.0).abs() < 1e-9 && (s.fraction - 0.75).abs() < 1e-9, "{s:?}"); + assert_eq!(s.sections.len(), 1); + assert!((scaled.initial_increment / plain.initial_increment - 2.0).abs() < 1e-9); + assert!((scaled.total_mass / plain.total_mass - 4.0).abs() < 1e-9); + assert!((scaled.energy.added_kinetic / scaled.energy.kinetic - 0.75).abs() < 1e-9, "{:?}", scaled.energy); + assert!(scaled.warnings.iter().any(|w| w.contains("mass scaling"))); + assert!( + scaled.controlling.first().is_some_and(|c| (0.9 * c.increment / scaled.initial_increment - 1.0).abs() < 1e-6), + "{:?}", + scaled.controlling + ); +} diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs index f5ba1de..4c9c21f 100644 --- a/crates/model/src/lib.rs +++ b/crates/model/src/lib.rs @@ -68,8 +68,8 @@ pub use sections::{BeamProfile, Extent, Formulation, Integration, Section, Secti pub use sets::Set; pub use steps::{ ArcConstraint, ArcLength, BulkViscosity, DisplacementLimit, ExplicitDynamic, Hourglass, HourglassKind, Imperfection, InitialVelocity, - MAX_LOAD_PATH_POINTS, MAX_NONLINEAR_INCREMENTS, MAX_NONLINEAR_ITERATIONS, NonlinearStatic, OutputVariable, Procedure, RigidWall, Step, - StressRate, Tolerances, Transient, + MAX_LOAD_PATH_POINTS, MAX_NONLINEAR_INCREMENTS, MAX_NONLINEAR_ITERATIONS, MassScaling, NonlinearStatic, OutputVariable, Procedure, RigidWall, + Step, StressRate, Tolerances, Transient, }; pub use units::UnitSystem; diff --git a/crates/model/src/steps.rs b/crates/model/src/steps.rs index 10bd7ce..ae456d3 100644 --- a/crates/model/src/steps.rs +++ b/crates/model/src/steps.rs @@ -425,6 +425,22 @@ pub struct ExplicitDynamic { /// The objective stress rate of large deformation. #[serde(default)] pub stress_rate: StressRate, + /// Opt-in mass scaling: the elements whose stable increment is below a target get the + /// density that brings them to it (the added mass is reported). + #[serde(default, skip_serializing_if = "Option::is_none")] + pub mass_scaling: Option, +} + +/// Mass scaling to a target time increment: every element (of `set`, or of the whole model) whose +/// stable increment, times the step's scale, is below `target_increment` has its density raised +/// by (target / increment)², which brings it to the target. It changes the dynamics: the added +/// mass is reported, per section, with its share of the kinetic energy. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct MassScaling { + pub target_increment: f64, + #[serde(default, skip_serializing_if = "Option::is_none")] + pub set: Option, } fn yes() -> bool { @@ -467,6 +483,7 @@ impl ExplicitDynamic { bulk_viscosity: BulkViscosity::default(), nonlinear_geometry: true, stress_rate: StressRate::default(), + mass_scaling: None, } } @@ -487,6 +504,7 @@ impl ExplicitDynamic { .map(String::as_str) .chain(self.initial_velocities.iter().map(|v| v.set.as_str())) .chain(self.rigid_walls.iter().map(|w| w.set.as_str())) + .chain(self.mass_scaling.iter().filter_map(|m| m.set.as_deref())) } pub fn validate(&self) -> Result<()> { @@ -524,6 +542,12 @@ impl ExplicitDynamic { } w.friction.validate().map_err(|e| Error::Invalid(format!("rigid wall on set `{}`: {e}", w.set)))?; } + if let Some(m) = &self.mass_scaling { + positive("mass scaling target increment", m.target_increment)?; + if m.target_increment > end { + return Err(Error::Invalid(format!("the mass scaling target increment {} is longer than the step", m.target_increment))); + } + } self.hourglass.validate()?; self.bulk_viscosity.validate() } @@ -1012,6 +1036,9 @@ mod tests { }, ExplicitDynamic { hourglass: Hourglass { kind: HourglassKind::Viscous, coefficient: 2.0 }, ..base.clone() }, ExplicitDynamic { bulk_viscosity: BulkViscosity { linear: -1.0, quadratic: 1.0 }, ..base.clone() }, + ExplicitDynamic { mass_scaling: Some(MassScaling { target_increment: 0.0, set: None }), ..base.clone() }, + ExplicitDynamic { mass_scaling: Some(MassScaling { target_increment: f64::NAN, set: None }), ..base.clone() }, + ExplicitDynamic { mass_scaling: Some(MassScaling { target_increment: 2.0, set: None }), ..base.clone() }, ]; for b in bad { assert!(Procedure::ExplicitDynamic(b.clone()).validate().is_err(), "{b:?}"); diff --git a/docs/fem.md b/docs/fem.md index 497be23..2b067e0 100644 --- a/docs/fem.md +++ b/docs/fem.md @@ -468,3 +468,15 @@ geometry and take no thermal strain yet (both refused by name). Verified by unia λ = 1.5 of all three laws on hex8 and tet4 against the material point (`crates/explicit/tests/hyperelastic.rs`) and by `tests/oracle/explicit-neo-hookean-implicit`, a quasi-static explicit pull against the implicit nonlinear solver. + +## Mass scaling (explicit dynamics) + +Opt-in per explicit step (`mass_scaling: {target_increment, set}`): every solid element (of the +set, or of the model) whose stable increment times the step's scale is below the target has its +density raised by (target / increment)², which brings its increment to the target (its lumped +masses, frequency bound and hourglass and bulk-viscosity scales follow; adiabatic heating keeps +the physical density). The run reports the elements scaled, the added mass in all and per +section, its fraction of the model's mass, a warning, and the kinetic energy of the added mass +(`energy.added_kinetic`, part of the kinetic energy, so the balance holds). Verified on a uniform +bar scaled to twice its increment: added mass three times the bar's, the increment doubled, the +added mass's share of the kinetic energy three quarters (`crates/explicit/tests/one_d.rs`). diff --git a/docs/parity.md b/docs/parity.md index 9f87d39..4fc2ce8 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -64,7 +64,7 @@ products are run, no documentation text is copied). Status: **done**, **started* | Decks | Readers keeping unknown keywords | M9 | | planned | | Solve | Own linear static and modal solver, verified against analytic solutions | M4 | yes | started (solid elements: tet4, tet10, wedge6, hex8, hex20; 2D plane stress, plane strain and axisymmetric tri3, tri6, quad4 (full, incompatible modes), quad8 (full, reduced); quad4/tri3 MITC shells and Timoshenko beams with profiles, mixed with solids; reference points, kinematic and distributing couplings, moments; centrifugal loads) | | Solve | Own nonlinear implicit solver: Newton–Raphson with line search and automatic increments, total Lagrangian solids and 2D elements, mixed u–p at finite strain, follower pressure | M15 | yes | started (solids, plane strain, axisymmetric; St Venant–Kirchhoff, neo-Hookean and J2 plasticity (Simo's exponential return map at finite strain, radial return at small strain, consistent tangents, the return shared with explicit dynamics); verified by the elastica, tube inflation, patch and objectivity tests, uniaxial tension to large strain, Hill's thick cylinder, a cantilever's collapse and the Bauschinger effect; Mooney–Rivlin and Yeoh; thermal strain at finite and small strain; a symmetric indefinite LDLᵀ factorization with inertia, Newton through indefinite tangents (hybrid elements at ν = 0.4999 under prescribed displacement); linear buckling of solids (Euler columns, Bryan's plate) and imperfections from its modes; arc-length with critical points from the inertia (shallow-arch snap-through); shells, implicit contact next) | -| Solve | Own explicit dynamics solver (central differences, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, energy balance), verified against analytic wave and impact solutions | M10 | yes | started (X1 core and X2: large deformation, objective rates, J2 and Johnson–Cook plasticity, damage, erosion; X3: deformable contact and friction, verified by bar collision, sliding blocks, the contact patch test and Hertz) | +| Solve | Own explicit dynamics solver (central differences, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, energy balance), verified against analytic wave and impact solutions | M10 | yes | started (X1 core and X2: large deformation, objective rates, J2 and Johnson–Cook plasticity, damage, erosion; X3: deformable contact and friction, verified by bar collision, sliding blocks, the contact patch test and Hertz; hyperelastic materials on the shared material point; opt-in mass scaling) | | Results | Deformed shape, contours, mode shapes, probes, reaction sums | M4 | | started (results in the session, `solve.probe`, reaction sums; display planned) | | Results | XY and path plots, cuts, display groups, free-body forces, animation | M9 | | planned | | Results | Readers for public result formats (Nastran OP2, LS-DYNA d3plot, Abaqus `.fil`) | M9 | | planned |