From 17b2c837f2c71c98c9a8deea74996202d414b792 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Wed, 14 Oct 2026 05:07:52 -0400 Subject: [PATCH] solver: J2 plasticity in nonlinear static steps, finite and small strain Elements of a plastic material take the J2 material point: total Lagrangian with nonlinear geometry, the small-strain loops without (elastic elements keep their linear stiffness, so elastic results do not move). The plastic history is committed per converged increment. Increments follow the step's load path, ending at its points, the extrapolated predictor scaled by the change of the load factor. Steps with path-dependent materials start an increment without a predictor from the tangent predictor (the prescribed motion spread through the converged tangent), and a first tangent that is not positive definite is a cutback when the undeformed model's is. Refused by name: the Johnson-Cook rate term (a static step has no time) and kinematic hardening at finite strain; damage, erosion and adiabatic heating are left out with a warning. Results carry the equivalent plastic strain at the integration points (solve.probe, the step summary). LS-DYNA writes a load path as the step's curve; Abaqus refuses it and a plastic deck of several general steps (Abaqus carries the plastic state between them). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm --- crates/decks/src/abaqus.rs | 26 +++- crates/decks/src/lsdyna/mod.rs | 12 +- crates/decks/tests/abaqus.rs | 20 +++ crates/decks/tests/lsdyna.rs | 13 ++ crates/engine/src/cmd/define.rs | 2 +- crates/engine/src/cmd/solve.rs | 11 +- crates/model/src/materials.rs | 10 +- crates/solver/src/analysis.rs | 7 +- crates/solver/src/nonlinear.rs | 231 ++++++++++++++++++++++----- crates/solver/src/system.rs | 71 +++++++-- crates/solver/tests/plasticity.rs | 249 ++++++++++++++++++++++++++++++ 11 files changed, 588 insertions(+), 64 deletions(-) create mode 100644 crates/solver/tests/plasticity.rs diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index 28a0b42..89f922a 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -37,7 +37,10 @@ //! and tolerances are left to Abaqus's controls, whose residual norms are defined differently. //! Abaqus starts each general step from the end of the previous one, while every model step //! starts from the undeformed state; for elastic materials the end state depends on the final -//! loads alone, so the two agree. +//! loads alone, so the two agree. A plastic material's state depends on the path, so a deck +//! with nonlinear steps, a plastic material and more than one general step is refused, as is +//! a nonlinear step's load path (an amplitude on every load and boundary condition, not +//! written yet). //! //! # Loads //! @@ -313,6 +316,27 @@ fn check_expressible(model: &Model) -> Result<()> { )); } } + for step in model.steps() { + if let Procedure::NonlinearStatic(n) = &step.procedure + && n.load_path.is_some() + { + bad.insert(format!( + "step `{}`: a load path (it would be an amplitude on every load and boundary condition of the step; not written yet)", + step.name + )); + } + } + let general = model.steps().iter().filter(|s| matches!(s.procedure, Procedure::Static {} | Procedure::NonlinearStatic(_))).count(); + let plastic_in_use = model.sections().values().find(|s| model.material(&s.material).is_some_and(|m| m.plastic.is_some())); + if has_nonlinear_step(model) + && general > 1 + && let Some(section) = plastic_in_use + { + bad.insert(format!( + "plastic material `{}` in a deck of {general} general steps (Abaqus carries the plastic state from one general step to the next, while every model step starts from the undeformed state; write each step to its own deck)", + section.material + )); + } for (name, m) in model.materials() { if m.damage.is_some() || m.erosion.is_some() { bad.insert(format!( diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index 50db6d9..40d26c4 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -25,7 +25,8 @@ //! - per step: `*CONTROL_IMPLICIT_GENERAL`, `*CONTROL_IMPLICIT_SOLUTION` (linear) or //! `*CONTROL_IMPLICIT_EIGENVALUE`, `*CONTROL_TERMINATION` (a static step is one implicit step //! from time 0 to 1 with every load ramped by one curve, `*DEFINE_CURVE`); a nonlinear static -//! step runs over the same time with the same ramp, its first increment as `dt0`, the nonlinear +//! step runs over the same time with the same ramp (or its load path as the curve), its first +//! increment as `dt0`, the nonlinear //! solver with a stiffness reformation every iteration (`ilimit` 1, full Newton) and its //! iteration cap as `maxref`, and `*CONTROL_IMPLICIT_AUTO` with its minimum and maximum //! increments (its increment cap and tolerances have no LS-DYNA counterpart and are left to @@ -872,8 +873,13 @@ fn write_step(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan, step: &Ste d.header(&["lcid"]); d.card([Field::from(n)]); d.header(&["a1", "o1"]); - d.card([Field::from(0.0), Field::from(0.0)]); - d.card([Field::from(1.0), Field::from(1.0)]); + let path = match &step.procedure { + Procedure::NonlinearStatic(s) => s.load_path.clone(), + _ => None, + }; + for [t, f] in path.unwrap_or_else(|| vec![[0.0, 0.0], [1.0, 1.0]]) { + d.card([Field::from(t), Field::from(f)]); + } } for b in &step.boundary_conditions { diff --git a/crates/decks/tests/abaqus.rs b/crates/decks/tests/abaqus.rs index 33638d4..543d601 100644 --- a/crates/decks/tests/abaqus.rs +++ b/crates/decks/tests/abaqus.rs @@ -795,3 +795,23 @@ fn nonlinear_steps_and_hyperelastic_materials_are_written() { m.remove_step("Bend").unwrap(); assert!(abaqus::write(&m).unwrap().contains("*ELASTIC, TYPE=ISOTROPIC\n"), "with linear steps only, the elastic constants are the whole of it"); } + +/// A plastic material in nonlinear steps is `*PLASTIC` in a general step; what Abaqus would +/// read differently is refused: a load path, and several general steps (Abaqus carries the +/// plastic state between them). +#[test] +fn nonlinear_plastic_steps_and_what_is_refused() { + use fieldcraft_model::{Hardening, NonlinearStatic, Plastic}; + let mut m = cantilever(); + let plastic = Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 2000.0, hardening: Hardening::Isotropic }; + m.replace_material("Steel", Material { plastic: Some(plastic), ..steel() }).unwrap(); + m.replace_step(step("Bend", Procedure::NonlinearStatic(NonlinearStatic::default()), &["clamp"], &["end"])).unwrap(); + let deck = abaqus::write(&m).unwrap(); + assert!(deck.contains("*PLASTIC\n250.0, 0.0\n") && deck.contains("NLGEOM=YES"), "{deck}"); + let path = NonlinearStatic { load_path: Some(vec![[0.0, 0.0], [0.5, 1.0], [1.0, 0.0]]), ..NonlinearStatic::default() }; + m.replace_step(step("Bend", Procedure::NonlinearStatic(path), &["clamp"], &["end"])).unwrap(); + assert!(unsupported(&m).contains("step `Bend`: a load path")); + m.replace_step(step("Bend", Procedure::NonlinearStatic(NonlinearStatic::default()), &["clamp"], &["end"])).unwrap(); + m.add_step(step("Again", Procedure::NonlinearStatic(NonlinearStatic::default()), &["clamp"], &["end"]), None).unwrap(); + assert!(unsupported(&m).contains("plastic material `Steel` in a deck of 2 general steps")); +} diff --git a/crates/decks/tests/lsdyna.rs b/crates/decks/tests/lsdyna.rs index 302cc34..2a86c63 100644 --- a/crates/decks/tests/lsdyna.rs +++ b/crates/decks/tests/lsdyna.rs @@ -760,6 +760,19 @@ fn nonlinear_static_steps() { assert_eq!(auto.iter().map(num).collect::>(), [1.0, 1e-4, 0.5], "blank fields are skipped: {auto:?}"); assert!(deck.contains("*DEFINE_CURVE_TITLE\nBend ramp\n"), "the loads follow the ramp"); assert!(deck.contains("*LOAD_NODE"), "{deck}"); + // A load path is the curve. + let path = NonlinearStatic { load_path: Some(vec![[0.0, 0.0], [0.5, 1.0], [1.0, -1.0]]), ..controls.clone() }; + m.replace_step(step("Bend", Procedure::NonlinearStatic(path), &["clamp"], &["tip_load"])).unwrap(); + let deck = lsdyna::write(&m).unwrap().text; + let curve: Vec> = deck + .lines() + .skip_while(|l| *l != "Bend ramp") + .filter(|l| !l.starts_with('$')) + .skip(2) + .take(3) + .map(|l| l.split_whitespace().map(|v| v.parse().unwrap()).collect()) + .collect(); + assert_eq!(curve, [vec![0.0, 0.0], vec![0.5, 1.0], vec![1.0, -1.0]], "{deck}"); m.replace_step(step("Bend", Procedure::NonlinearStatic(NonlinearStatic { nonlinear_geometry: false, ..controls }), &["clamp"], &["tip_load"])) .unwrap(); diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs index 5cce949..f4bf3de 100644 --- a/crates/engine/src/cmd/define.rs +++ b/crates/engine/src/cmd/define.rs @@ -266,7 +266,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} (loads and prescribed displacements ramp with the load factor from 0 to 1) | {"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"}], "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"}], "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)}"#, kind: CommandKind::Model, run: step_create, }; diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index e642569..bfe52b9 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -34,7 +34,7 @@ pub const STATIC: CommandSpec = CommandSpec { pub const NONLINEAR: CommandSpec = CommandSpec { id: "solve.nonlinear", label: "Solve Nonlinear Static Step", - params: r#"{"step": "Inflate", "max_dofs": 6000000, "max_factor_entries": 400000000, "max_newton_iterations": 100000, "threads": 8} — solves a nonlinear static step (the only one when `step` is left out) in load increments with Newton–Raphson iterations and keeps the result in the session: the final displacements, reactions, Cauchy stresses and logarithmic strains, every increment's load factor, iterations and residuals (force, displacement and energy, relative to their tolerances' scales), and a displacement frame per converged increment; the caps may only be lowered (`max_newton_iterations` over the whole step); `threads` as for solve.static"#, + params: r#"{"step": "Inflate", "max_dofs": 6000000, "max_factor_entries": 400000000, "max_newton_iterations": 100000, "threads": 8} — solves a nonlinear static step (the only one when `step` is left out) in load increments with Newton–Raphson iterations and keeps the result in the session: the final displacements, reactions, Cauchy stresses and logarithmic strains (small strains without nonlinear geometry), the equivalent plastic strain of plastic materials (the largest in the summary), every increment's load factor, iterations and residuals (force, displacement and energy, relative to their tolerances' scales), and a displacement frame per converged increment; the caps may only be lowered (`max_newton_iterations` over the whole step); `threads` as for solve.static"#, kind: CommandKind::Model, run: run_nonlinear, }; @@ -82,7 +82,7 @@ pub const FLUID: CommandSpec = CommandSpec { pub const PROBE: CommandSpec = CommandSpec { id: "solve.probe", label: "Probe Results", - params: r#"{"step": "Load", "nodes": [1, 2], "mode": 1, "frame": 0} — the stored result at nodes: displacement, rotation, reaction, reaction moment, stress and von Mises (static, with `elements` each element's results: solid, 2D and truss nodal stresses, shell section forces (N, M, Q per length in the element's axes) and bottom, middle and top stresses, beam end forces (N, V_y, V_z, T, M_y, M_z in its axes)), the shape of mode `mode` with its rotations (frequency, from 1), or displacement and velocity in field frame `frame` (explicit, from 0; the last by default) at `nodes` (or at the nodes nearest to points: "at": [[x, y, z], ...], lowest label on a tie, each entry then giving its `xyz` and `distance`) and, with `elements` (explicit), each element's stress, von Mises stress, strain, equivalent plastic strain, damage, temperature and whether it has eroded; a coupled reference point's label gives its six displacements (translations, rotations) and reactions (forces, moments). Fluid steps: {"step": "Flow", "points": [[x, y, z]], "frame": 10 (default the last)} gives velocity and pressure in the cells at the points; {"step": "Flow", "forces": ["x_min", "F2:top", "Body1/F2:top"]} gives those walls' force histories"#, + params: r#"{"step": "Load", "nodes": [1, 2], "mode": 1, "frame": 0} — the stored result at nodes: displacement, rotation, reaction, reaction moment, stress and von Mises (static, with `elements` each element's results: solid, 2D and truss nodal stresses (and, in a nonlinear step, the equivalent plastic strain at the integration points of a plastic element), shell section forces (N, M, Q per length in the element's axes) and bottom, middle and top stresses, beam end forces (N, V_y, V_z, T, M_y, M_z in its axes)), the shape of mode `mode` with its rotations (frequency, from 1), or displacement and velocity in field frame `frame` (explicit, from 0; the last by default) at `nodes` (or at the nodes nearest to points: "at": [[x, y, z], ...], lowest label on a tie, each entry then giving its `xyz` and `distance`) and, with `elements` (explicit), each element's stress, von Mises stress, strain, equivalent plastic strain, damage, temperature and whether it has eroded; a coupled reference point's label gives its six displacements (translations, rotations) and reactions (forces, moments). Fluid steps: {"step": "Flow", "points": [[x, y, z]], "frame": 10 (default the last)} gives velocity and pressure in the cells at the points; {"step": "Flow", "forces": ["x_min", "F2:top", "Body1/F2:top"]} gives those walls' force histories"#, kind: CommandKind::Model, run: probe, }; @@ -220,6 +220,7 @@ fn run_nonlinear(s: &mut Session, p: &Value) -> Result { "frames": history.map(|h| h.load_factors()), "max_displacement": largest(r.nodes.iter().zip(&r.displacement).map(|(l, u)| (*l, crate::results::magnitude(*u)))), "max_von_mises": largest(r.nodes.iter().zip(&r.von_mises).map(|(l, v)| (*l, *v))), + "max_plastic_strain": largest(r.elements.iter().flat_map(|e| e.plastic_strain.iter().map(|v| (e.label, *v)))), "applied_total": r.applied_total, "reaction_total": r.reaction_total, "reference_points": r.reference_points.iter().map(|p| json!({ @@ -613,7 +614,11 @@ fn by_label(list: &[T], label: u32, key: impl Fn(&T) -> u32) -> Option<&T> { /// The stored results of one element of a static step. fn element_result(r: &StaticResult, label: u32) -> Option { if let Some(e) = by_label(&r.elements, label, |e| e.label) { - return Some(json!({"element": label, "kind": e.kind.name(), "stress": e.nodal})); + let mut v = json!({"element": label, "kind": e.kind.name(), "stress": e.nodal}); + if !e.plastic_strain.is_empty() { + v["plastic_strain"] = json!(e.plastic_strain); + } + return Some(v); } if let Some(e) = by_label(&r.shells, label, |e| e.label) { let s = &e.stress; diff --git a/crates/model/src/materials.rs b/crates/model/src/materials.rs index 61dc819..ecbf320 100644 --- a/crates/model/src/materials.rs +++ b/crates/model/src/materials.rs @@ -1,7 +1,7 @@ //! Materials: linear elasticity, mass density, fluid behaviour, thermal properties (conductivity, -//! specific heat, expansion), and the inelastic behaviours of -//! explicit dynamics: von Mises (J2) plasticity with a bilinear, tabulated or Johnson–Cook flow stress, -//! Johnson–Cook damage, and element erosion criteria. +//! specific heat, expansion), and the inelastic behaviours: von Mises (J2) plasticity with a +//! bilinear, tabulated or Johnson–Cook flow stress (explicit dynamics and nonlinear static +//! steps), Johnson–Cook damage and element erosion criteria (explicit dynamics). //! //! The records are solver-neutral and named for what they describe, never after a product's //! keyword or material number; the deck writers map them to each format. @@ -372,8 +372,8 @@ pub struct Material { /// How it flows, for fluid domains (needs a density). #[serde(default, skip_serializing_if = "Option::is_none")] pub fluid: Option, - /// Von Mises plasticity (explicit dynamics; the linear implicit solver uses the elastic - /// constants alone). + /// Von Mises plasticity (explicit dynamics and nonlinear static steps; the linear solvers use + /// the elastic constants alone). #[serde(default, skip_serializing_if = "Option::is_none")] pub plastic: Option, /// Damage to failure (explicit dynamics: a failed element erodes). diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 1cf6337..9aa1bb3 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -61,6 +61,9 @@ pub struct ElementStress { pub points: Vec, /// Extrapolated to the element's nodes, in its node order. pub nodal: Vec<[f64; 6]>, + /// The equivalent plastic strain at the integration points (nonlinear steps with a plastic + /// material; empty otherwise). + pub plastic_strain: Vec, } /// The motion and reaction of a coupled reference point: translations along X, Y, Z, then @@ -431,7 +434,7 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul *c += 1; } } - elements.push(ElementStress { label: e.label, kind: e.kind, points, nodal }); + elements.push(ElementStress { label: e.label, kind: e.kind, points, nodal, plastic_strain: Vec::new() }); } for t in &mesh.trusses { let truss = t.truss()?; @@ -448,7 +451,7 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul *c += 1; } } - elements.push(ElementStress { label: t.label, kind: t.kind, points, nodal }); + elements.push(ElementStress { label: t.label, kind: t.kind, points, nodal, plastic_strain: Vec::new() }); } elements.sort_by_key(|e| e.label); let (mut shells, mut beams) = (Vec::new(), Vec::new()); diff --git a/crates/solver/src/nonlinear.rs b/crates/solver/src/nonlinear.rs index 690a281..d2ac241 100644 --- a/crates/solver/src/nonlinear.rs +++ b/crates/solver/src/nonlinear.rs @@ -25,7 +25,9 @@ //! for the solution of nonlinear finite element equations*, CMAME 22:59–85, 1980); //! //! or when the force residual is a thousandth of its tolerance (round-off: no correction can -//! improve it). +//! improve it). With path-dependent materials, an increment that starts within the force +//! tolerance (a predictor that already solved it) waives the energy criterion, whose reference +//! would be round-off. //! //! **Automatic incrementation.** Each increment starts from the linear extrapolation of the //! last converged one (the free motion scaled by the ratio of their sizes), which keeps the @@ -50,10 +52,21 @@ //! //! **Elements.** Solids (tet4, tet10, wedge6, hex8, hex20 full and reduced, hybrid hex8, hex20 //! and tet10) and 2D plane strain and axisymmetric elements, total Lagrangian -//! ([`fieldcraft_fem::finite`]); without nonlinear geometry the elements are the linear ones -//! (every material is elastic today, so such a step is the linear solution reached by Newton). +//! ([`fieldcraft_fem::finite`]); without nonlinear geometry an element of a plastic material +//! takes the same loops with small strain, and an elastic one its linear stiffness (K u). //! Springs act along their fixed directions. Trusses, shells and beams (large rotations), plane -//! stress, incompatible modes, temperature and centrifugal loads are refused by name. +//! stress, incompatible modes (with nonlinear geometry, or with an inelastic material), +//! temperature and centrifugal loads are refused by name. +//! +//! **Plasticity** (`Material::plastic`, J2 with the radial return explicit dynamics uses, +//! [`fieldcraft_fem::elastoplastic`]): Simo's exponential return map at finite strain, the +//! additive radial return at small strain, both with their consistent tangents, so Newton keeps +//! its quadratic rate. The plastic history is committed when an increment converges and +//! discarded when it is cut back. A static step has no time, so the Johnson–Cook strain-rate term +//! is refused (C = 0 is the rate-independent law), and kinematic hardening at finite strain is +//! refused (its back stress breaks the principal-frame return); damage, erosion and adiabatic +//! heating act in explicit dynamics only, with a warning. Results add the equivalent plastic +//! strain at the integration points. //! Couplings stay the linear (small-rotation) constraints of the linear solver, with a warning. //! //! Element forces and tangents are computed in parallel and added in label order, so results are @@ -63,9 +76,10 @@ use fieldcraft_fem::constitutive::MaterialModel; use fieldcraft_fem::dense::Matrix; use fieldcraft_fem::element::{Extrapolation, extrapolation}; use fieldcraft_fem::finite::Kinematics; +use fieldcraft_fem::plane::Theory; use fieldcraft_fem::quadrature::Rule; use fieldcraft_fem::{DofKind, DofSet, ElementType, PointStress, face, von_mises}; -use fieldcraft_model::{Extent, Load, Model, NonlinearStatic, Procedure, Step}; +use fieldcraft_model::{Extent, Hardening, Load, Model, NonlinearStatic, Plastic, Procedure, SectionKind, Step}; use fieldcraft_parallel as parallel; use serde::Serialize; @@ -76,7 +90,7 @@ use crate::analysis::{ use crate::cholesky::FactorError; use crate::constraints::Reduction; use crate::sparse::{Csr, add_block}; -use crate::system::{Mesh, inelastic_warnings, loads, translations}; +use crate::system::{Mesh, loads, translations}; use crate::{Error, Limits, Result}; /// One Newton iteration of an increment. @@ -237,6 +251,12 @@ impl Problem<'_> { self.element_dofs.get(e).map(|dofs| dofs.iter().map(|d| u.get(*d).copied().unwrap_or(0.0)).collect()).unwrap_or_default() } + /// How element `e` is evaluated: total Lagrangian with nonlinear geometry, with small + /// strain when its material is inelastic, and by its linear stiffness otherwise (`None`). + fn kinematics(&self, e: usize) -> Option { + if self.finite { Some(Kinematics::Finite) } else { self.materials.get(e).filter(|m| m.is_path_dependent()).map(|_| Kinematics::Small) } + } + /// Internal forces and the tangent of element `e` at displacements `u` (all degrees of /// freedom). fn element(&self, e: usize, u: &[f64], tangent: bool, committed: &[f64]) -> Result<(Vec, Option, Vec)> { @@ -244,14 +264,14 @@ impl Problem<'_> { let at = |err| Error::element(el.label, err); let ue = self.element_vector(e, u); let kernel = el.kernel()?; - if !self.finite { + let Some(kin) = self.kinematics(e) else { let k = kernel.stiffness(&el.d).map_err(at)?; let f = k.mul_vec(&ue); return Ok((f, tangent.then_some(k), committed.to_vec())); - } + }; let material = self.materials.get(e).ok_or_else(|| Error::Internal("element material".into()))?; let mut updated = committed.to_vec(); - let r = kernel.nonlinear(material, &ue, Kinematics::Finite, committed, &mut updated, tangent).map_err(at)?; + let r = kernel.nonlinear(material, &ue, kin, committed, &mut updated, tangent).map_err(at)?; Ok((r.force, r.tangent, updated)) } @@ -382,6 +402,10 @@ struct Run<'a> { /// Every element's material history at the end of the last converged increment /// (integration points × the material's history length). committed: Vec>, + /// Whether a material keeps history (the tangent predictor and its convergence rules). + path_dependent: bool, + /// Whether the undeformed model's tangent factors, once asked. + undeformed_stable: Option, } /// The outcome of a converged increment. @@ -405,12 +429,53 @@ impl Run<'_> { Ok(self.problem.evaluate(&self.problem.red.expand_scaled(q, lambda), lambda, tangent, &self.committed)?) } + /// Whether the tangent of the undeformed, unloaded model is positive definite: a first + /// tangent that is not, at a state the first increment has already strained, is then the + /// material or the geometry (a cutback), not a model that can move freely (an error). + fn undeformed_is_stable(&mut self) -> std::result::Result { + if let Some(stable) = self.undeformed_stable { + return Ok(stable); + } + let mut e = self.problem.evaluate(&vec![0.0; self.problem.mesh.dofs()], 0.0, true, &self.committed)?; + let k = e.tangent.take().ok_or_else(|| Failure::Fatal(Error::Internal("tangent".into())))?; + let [kff] = self.problem.red.reduce([&k])?; + let stable = factor_raw(&kff, self.limits).is_ok(); + self.undeformed_stable = Some(stable); + Ok(stable) + } + + /// The tangent predictor of an increment from the converged `q0` at load factor `from` to + /// `to`, and its work δq · r (the increment's reference for the energy criterion, since the + /// Newton iterations after it only correct it): one linear solve with the tangent of the converged state for the change of the + /// loads and of the prescribed displacements, K δq = T_fᵀ (f_ext(to) − f_int(u₀) − K T_p Δλ u_p). + /// It spreads a prescribed motion into the interior before the material sees it, which a + /// path-dependent material needs: started from the boundary's jump alone, the elements next + /// to it would flow plastically far beyond the solution. `None` when the tangent cannot be + /// factored (the increment then starts from `q0`). + fn tangent_predictor(&mut self, q0: &[f64], from: f64, to: f64) -> std::result::Result, f64)>, Failure> { + let red = self.problem.red; + let mut e = self.problem.evaluate(&red.expand_scaled(q0, from), to, true, &self.committed)?; + let k = e.tangent.take().ok_or_else(|| Failure::Fatal(Error::Internal("tangent".into())))?; + let jump = red.expand_scaled(&vec![0.0; red.neq()], to - from); + let mut kj = vec![0.0; jump.len()]; + k.mul(&jump, &mut kj); + let rf: Vec = free_residual(red, &e).iter().zip(red.restrict(&kj)).map(|(a, b)| a - b).collect(); + let mut dq = rf.clone(); + let [kff] = red.reduce([&k])?; + drop(k); + let Ok(chol) = factor_raw(&kff, self.limits) else { return Ok(None) }; + chol.solve(&mut dq)?; + let energy = dot(&dq, &rf).abs(); + Ok(dq.iter().all(|v| v.is_finite()).then(|| (q0.iter().zip(&dq).map(|(a, d)| a + d).collect(), energy))) + } + /// Newton iterations from the converged `q0` (at the previous load factor) to load factor /// `lambda`. fn increment( &mut self, q0: &[f64], predictor: Option<&[f64]>, + predicted_energy: Option, lambda: f64, size: f64, iterations: &mut Vec, @@ -425,8 +490,14 @@ impl Run<'_> { Some(Err(Failure::Cutback(_))) | None => (q0.to_vec(), self.state(q0, lambda, true)?), }; let mut r = self.residual(&eval); + // With path-dependent materials, an increment that starts within the force tolerance (a + // predictor that already solved it) has only round-off left for the energy criterion to + // measure against: it is waived. + let started_converged = self.path_dependent && self.force_floor.is_some() && norm_inf(&r) <= tol.force * self.scale(); let prescribed = size * red.largest_prescribed(); - let mut first_energy: Option = None; + // The energy criterion's reference: the first correction's work, or the tangent + // predictor's when the increment started from one that is still in use. + let mut first_energy: Option = predicted_energy.filter(|_| q.as_slice() != q0); let mut history: Vec = Vec::new(); let max_iterations = usize::try_from(self.settings.max_iterations).unwrap_or(usize::MAX); for _ in 0..max_iterations { @@ -453,7 +524,7 @@ impl Run<'_> { drop(k); let chol = match factor_raw(&kff, self.limits) { Ok(c) => c, - Err(FactorError::NotPositiveDefinite { .. }) if !first => { + Err(FactorError::NotPositiveDefinite { .. }) if !first || self.undeformed_is_stable()? => { return Err(Failure::Cutback( "the tangent stiffness is not positive definite (a limit point, a bifurcation or an unstable material)".into(), )); @@ -488,7 +559,8 @@ impl Run<'_> { if !it.force.is_finite() { return Err(Failure::Cutback("the residual is not finite".into())); } - if (it.force <= tol.force && it.displacement <= tol.displacement && it.energy <= tol.energy) || it.force <= ROUND_OFF * tol.force { + let energy_met = it.energy <= tol.energy || started_converged; + if (it.force <= tol.force && it.displacement <= tol.displacement && energy_met) || it.force <= ROUND_OFF * tol.force { return Ok(Converged { q, evaluation: eval }); } history.push(it.force); @@ -582,6 +654,7 @@ fn check_supported(model: &Model, mesh: &Mesh, step: &Step, settings: &Nonlinear _ => {} } } + plasticity(model, step, settings, warnings)?; if settings.nonlinear_geometry && let Some(name) = model.couplings().keys().next() { @@ -592,6 +665,45 @@ fn check_supported(model: &Model, mesh: &Mesh, step: &Step, settings: &Nonlinear Ok(()) } +/// Refuse the plasticity a static step cannot take, and warn about what it leaves out, for the +/// materials of the solid and 2D sections. +fn plasticity(model: &Model, step: &Step, settings: &NonlinearStatic, warnings: &mut Vec) -> Result<()> { + let at = |what: String| Error::NotSupported(format!("{what} in nonlinear steps (step `{}`)", step.name)); + let used: std::collections::BTreeSet<&str> = model + .sections() + .values() + .filter(|s| matches!(s.kind, SectionKind::Solid {}) || Theory::from_section(&s.kind).is_some()) + .map(|s| s.material.as_str()) + .collect(); + for name in used { + let Some(m) = model.material(name) else { continue }; + match &m.plastic { + Some(Plastic::JohnsonCook { c, specific_heat, .. }) => { + if *c != 0.0 { + return Err(at(format!( + "the Johnson–Cook strain-rate term of material `{name}` (C = {c}): a static step has no time, so no strain rate; set C to 0 or use explicit dynamics" + ))); + } + if specific_heat.is_some() { + warnings.push(format!( + "material `{name}`: a static step is isothermal, so plastic work does not heat it (adiabatic heating acts in explicit dynamics)" + )); + } + } + Some(Plastic::Bilinear { hardening: Hardening::Kinematic, .. }) if settings.nonlinear_geometry => { + return Err(at(format!( + "the kinematic hardening of material `{name}` at finite strain (its back stress is not coaxial with the elastic stretch, which the exponential return map needs); use nonlinear_geometry false or isotropic hardening" + ))); + } + _ => {} + } + if m.plastic.is_some() && (m.damage.is_some() || m.erosion.is_some()) { + warnings.push(format!("material `{name}`: damage and erosion act in explicit dynamics; this static step leaves them out")); + } + } + Ok(()) +} + /// The follower pressures of a step (with nonlinear geometry), and the step without them for /// the dead loads. fn followers(model: &Model, resolved: &fieldcraft_model::Resolved, mesh: &Mesh, step: &Step) -> Result<(Vec, Step)> { @@ -646,7 +758,6 @@ pub fn solve_nonlinear(model: &Model, step: Option<&str>, limits: &Limits) -> Re settings.validate()?; let mesh = Mesh::new(model, &resolved, limits.max_dofs)?; let mut warnings: Vec = resolved.warnings.iter().map(ToString::to_string).chain(mesh.warnings.iter().cloned()).collect(); - warnings.extend(inelastic_warnings(model)); check_supported(model, &mesh, step, settings, &mut warnings)?; let (graph, order) = pattern(&mesh); ignored_dashpots(&mesh, &mut warnings); @@ -683,7 +794,14 @@ pub fn solve_nonlinear(model: &Model, step: Option<&str>, limits: &Limits) -> Re finite: settings.nonlinear_geometry, element_dofs, element_kinds, - materials: mesh.elements.iter().map(|e| MaterialModel::from_elasticity(&e.d, e.hyperelastic.as_ref())).collect(), + materials: mesh + .elements + .iter() + .map(|e| match &e.plastic { + Some(j2) => MaterialModel::Plastic(j2.clone()), + None => MaterialModel::from_elasticity(&e.d, e.hyperelastic.as_ref()), + }) + .collect(), pattern: graph.zero_matrix(&mesh.layout), }; let (lo, hi) = mesh.xyz.iter().fold(([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]), |(lo, hi), p| { @@ -699,8 +817,20 @@ pub fn solve_nonlinear(model: &Model, step: Option<&str>, limits: &Limits) -> Re // The undeformed state first: an element that is bad before anything moves is an error in // the model, not a reason to cut increments back. problem.evaluate(&vec![0.0; mesh.dofs()], 0.0, false, &committed)?; - let mut run = - Run { problem: &problem, settings, limits, force_scale: 0.0, force_floor: None, length, iterations: 0, factor_entries: 0, committed }; + let path_dependent = problem.materials.iter().any(MaterialModel::is_path_dependent); + let mut run = Run { + problem: &problem, + settings, + limits, + force_scale: 0.0, + force_floor: None, + length, + iterations: 0, + factor_entries: 0, + committed, + path_dependent, + undeformed_stable: None, + }; let mut history = History::default(); let (mut frame_values, mut frames_full) = (0usize, false); @@ -719,32 +849,50 @@ pub fn solve_nonlinear(model: &Model, step: Option<&str>, limits: &Limits) -> Re }; keep_frame(&mut history, 0.0, &vec![0.0; mesh.dofs()], &mut warnings); let mut q = vec![0.0; red.neq()]; - let mut lambda = 0.0f64; + // The step fraction t and the load factor λ(t) of the step's load path (λ = t without one). + let (mut t, mut lambda) = (0.0f64, 0.0f64); let mut size = settings.initial_increment; let mut fast = 0usize; let mut last: Option = None; + // The last converged increment's free motion and the change of the load factor over it. let mut previous: Option<(Vec, f64)> = None; let max_increments = usize::try_from(settings.max_increments).unwrap_or(usize::MAX); - while lambda < 1.0 { + while t < 1.0 { if history.increments.len() >= max_increments { return Err(Error::Limit(format!( "step `{}`: {} increments reached the step's cap at load factor {lambda:.6} (raise max_increments)", step.name, settings.max_increments ))); } - let remaining = 1.0 - lambda; - let (target, this) = if size >= remaining * (1.0 - 1e-9) { (1.0, remaining) } else { (lambda + size, size) }; + // An increment ends at the next point of the load path at the latest. + let turn = settings.next_turn(t); + let remaining = turn - t; + let (target, this) = if size >= remaining * (1.0 - 1e-9) { (turn, remaining) } else { (t + size, size) }; + let target_lambda = settings.load_factor(target); + let change = if settings.load_path.is_some() { target_lambda - lambda } else { this }; let mut iterations = Vec::new(); let number = history.increments.len() + 1; - // Linear extrapolation of the last converged increment's free motion. - let predictor: Option> = - previous.as_ref().map(|(dq, prev): &(Vec, f64)| q.iter().zip(dq).map(|(a, d)| a + this / prev * d).collect()); - match run.increment(&q, predictor.as_deref(), target, this, &mut iterations) { + // Linear extrapolation of the last converged increment's free motion, in proportion to + // the change of the load factor (so it turns with the load path). + let mut predictor: Option> = previous + .as_ref() + .filter(|(_, prev)| *prev != 0.0) + .map(|(dq, prev): &(Vec, f64)| q.iter().zip(dq).map(|(a, d)| a + change / prev * d).collect()); + // Path-dependent materials start without one from the tangent predictor. + let mut predicted_energy = None; + if predictor.is_none() && path_dependent { + match run.tangent_predictor(&q, lambda, target_lambda) { + Ok(Some((p, energy))) => (predictor, predicted_energy) = (Some(p), Some(energy)), + Ok(None) | Err(Failure::Cutback(_)) => {} + Err(Failure::Fatal(e)) => return Err(e), + } + } + match run.increment(&q, predictor.as_deref(), predicted_energy, target_lambda, change.abs(), &mut iterations) { Ok(c) => { - previous = Some((c.q.iter().zip(&q).map(|(a, b)| a - b).collect(), this)); + previous = Some((c.q.iter().zip(&q).map(|(a, b)| a - b).collect(), change)); let count = iterations.len(); - history.increments.push(Increment { number, load_factor: target, size: this, iterations, converged: true, cutback: None }); - lambda = target; + history.increments.push(Increment { number, load_factor: target_lambda, size: this, iterations, converged: true, cutback: None }); + (t, lambda) = (target, target_lambda); q = c.q; keep_frame(&mut history, lambda, &red.expand_scaled(&q, lambda), &mut warnings); run.committed.clone_from(&c.evaluation.history); @@ -761,7 +909,7 @@ pub fn solve_nonlinear(model: &Model, step: Option<&str>, limits: &Limits) -> Re Err(Failure::Cutback(reason)) => { history.increments.push(Increment { number, - load_factor: target, + load_factor: target_lambda, size: this, iterations, converged: false, @@ -771,7 +919,7 @@ pub fn solve_nonlinear(model: &Model, step: Option<&str>, limits: &Limits) -> Re size = this * CUTBACK; if size < settings.min_increment * (1.0 - 1e-9) { return Err(Error::NotConverged(format!( - "step `{}` stopped at load factor {lambda:.6}: the increment to {target:.6} failed ({reason}) and a smaller one would be below the minimum increment {}", + "step `{}` stopped at load factor {lambda:.6} (step fraction {t:.6}): the increment to {target_lambda:.6} failed ({reason}) and a smaller one would be below the minimum increment {}", step.name, settings.min_increment ))); } @@ -827,27 +975,32 @@ fn finish( .zip(point_motion.iter().zip(&point_reaction)) .map(|(p, (d, r))| PointResult { name: p.name.clone(), label: p.label, displacement: *d, reaction: *r }) .collect(); - let all_points = parallel::map_grained(mesh.elements.len(), GRAIN, |i| -> Result<(Vec, &'static ElementType, Rule)> { + #[allow(clippy::type_complexity)] + let all_points = parallel::map_grained(mesh.elements.len(), GRAIN, |i| -> Result<(Vec, Vec, &'static ElementType, Rule)> { let el = mesh.elements.get(i).ok_or_else(|| Error::Internal("element index".into()))?; let at = |err| Error::element(el.label, err); let ue = problem.element_vector(i, &u); let kernel = el.kernel()?; - let points = if problem.finite { - let material = problem.materials.get(i).ok_or_else(|| Error::Internal("element material".into()))?; - kernel.nonlinear_stresses(material, &ue, Kinematics::Finite, committed.get(i).map(Vec::as_slice).unwrap_or_default()) - } else { - let per_node: Vec<[f64; 3]> = el.nodes.iter().map(|n| displacement.get(*n).copied().unwrap_or_default()).collect(); - kernel.stresses_with(&el.d, &per_node, None) + let history = committed.get(i).map(Vec::as_slice).unwrap_or_default(); + let material = problem.materials.get(i).ok_or_else(|| Error::Internal("element material".into()))?; + let points = match problem.kinematics(i) { + Some(kin) => kernel.nonlinear_stresses(material, &ue, kin, history), + None => { + let per_node: Vec<[f64; 3]> = el.nodes.iter().map(|n| displacement.get(*n).copied().unwrap_or_default()).collect(); + kernel.stresses_with(&el.d, &per_node, None) + } } .map_err(at)?; - Ok((points, kernel.ty(), kernel.rule())) + let len = material.history_len(); + let plastic: Vec = if len > 0 { history.chunks(len).filter_map(|h| material.plastic_strain(h)).collect() } else { Vec::new() }; + Ok((points, plastic, kernel.ty(), kernel.rule())) }); let mut extrapolations: Vec<(&'static ElementType, Rule, Extrapolation)> = Vec::new(); let mut elements = Vec::with_capacity(mesh.elements.len()); let mut sum = vec![[0.0; 6]; mesh.labels.len()]; let mut count = vec![0u32; mesh.labels.len()]; for (el, computed) in mesh.elements.iter().zip(all_points) { - let (points, ty, rule) = computed?; + let (points, plastic_strain, ty, rule) = computed?; let values: Vec<[f64; 6]> = points.iter().map(|p| p.stress).collect(); let index = match extrapolations.iter().position(|(t, r, _)| std::ptr::eq(*t, ty) && *r == rule) { Some(i) => i, @@ -865,7 +1018,7 @@ fn finish( *c += 1; } } - elements.push(ElementStress { label: el.label, kind: el.kind, points, nodal }); + elements.push(ElementStress { label: el.label, kind: el.kind, points, nodal, plastic_strain }); } let stress: Vec<[f64; 6]> = sum.iter().zip(&count).map(|(s, c)| s.map(|v| if *c > 0 { v / f64::from(*c) } else { 0.0 })).collect(); let von_mises = stress.iter().map(von_mises).collect(); diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index 7d624eb..0b34469 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -13,11 +13,13 @@ use fieldcraft_fem::dense::Matrix; use fieldcraft_fem::element::{Formulation, Solid, ThermalStrain}; use fieldcraft_fem::finite::Kinematics; use fieldcraft_fem::plane::{Plane, Theory}; +use fieldcraft_fem::plasticity::J2; use fieldcraft_fem::shell::{Shell, ShellSection}; use fieldcraft_fem::truss::{Truss, young_modulus}; use fieldcraft_fem::{DofKind, DofLayout, DofSet, Elasticity, ElementDofs, PointStress, elasticity, face, plane}; use fieldcraft_model::nodal::edge_pressure; use std::collections::BTreeSet; +use std::sync::Arc; use fieldcraft_model::{Elastic, ElementKind, Hyperelastic, Integration, Load, Model, Resolved, SectionKind, Step, Target, centrifugal_acceleration}; @@ -42,6 +44,9 @@ pub struct Elem { pub expansion: Option<([f64; 3], f64)>, /// The large-strain option of the element's material (nonlinear steps). pub hyperelastic: Option, + /// The element's material's plasticity (nonlinear steps), shared by the elements of a + /// section. + pub plastic: Option>, } /// An element's kernel: a solid or a 2D element. @@ -208,6 +213,21 @@ fn thermal_strain(expansion: Option<([f64; 3], f64)>, nodes: &[usize], field: &[ delta.iter().any(|d| *d != 0.0).then_some(ThermalStrain { alpha, delta }) } +/// A solid or 2D element of a section, by model labels, before the mesh numbers its nodes. +struct SolidEntry { + label: u32, + kind: ElementKind, + nodes: Vec, + formulation: Formulation, + d: Elasticity, + density: Option, + theory: Option, + dofs: ElementDofs, + expansion: Option<([f64; 3], f64)>, + hyperelastic: Option, + plastic: Option>, +} + /// What a shell or beam element is made of. #[derive(Debug, Clone, Copy)] pub enum Section { @@ -472,6 +492,13 @@ impl Mesh { continue; } let d = elasticity(elastic).map_err(|e| Error::Invalid(format!("material `{}`: {e}", section.material)))?; + let plastic = match &material.plastic { + Some(p) => Some(Arc::new( + J2::from_model(elastic, Some(p), material.damage.as_ref()) + .map_err(|e| Error::Invalid(format!("material `{}`: {e}", section.material)))?, + )), + None => None, + }; for &label in &set.elements { let el = model.element(label).ok_or_else(|| Error::Invalid(format!("element {label} is missing")))?; let (formulation, dofs) = match theory { @@ -480,20 +507,32 @@ impl Mesh { } .map_err(|e| Error::element(label, e))?; let expansion = material.thermal.and_then(|t| t.expansion.map(|a| (a.axes(), t.reference_temperature))); - elements.push((label, el.kind, el.nodes.clone(), formulation, d, material.density, theory, dofs, expansion, material.hyperelastic)); + elements.push(SolidEntry { + label, + kind: el.kind, + nodes: el.nodes.clone(), + formulation, + d, + density: material.density, + theory, + dofs, + expansion, + hyperelastic: material.hyperelastic, + plastic: plastic.clone(), + }); } } if let (Some(_), Some(name)) = (planar, model.couplings().keys().next()) { return Err(Error::NotSupported(format!("couplings in 2D models (coupling `{name}`)"))); } - elements.sort_by_key(|e| e.0); + elements.sort_by_key(|e| e.label); trusses.sort_by_key(|e| e.0); structural.sort_by_key(|e| e.0); let translations = translation_kinds(planar.is_some()); // What discrete elements need: at mesh nodes (some only they use), and at reference // points, which then get their six degrees of freedom like coupled ones. let needs: Vec<(u32, DofSet)> = crate::discrete::needs(resolved, translations).collect(); - let mut labels: Vec = elements.iter().flat_map(|e| e.2.iter().copied()).collect(); + let mut labels: Vec = elements.iter().flat_map(|e| e.nodes.iter().copied()).collect(); labels.extend(trusses.iter().flat_map(|e| e.2.iter().copied())); labels.extend(structural.iter().flat_map(|e| e.2.iter().copied())); labels.extend(needs.iter().map(|n| n.0).filter(|l| model.node(*l).is_some())); @@ -511,7 +550,7 @@ impl Mesh { .collect::>()?; let mut layout = DofLayout::builder(); for e in &elements { - layout.element(&e.2, e.7); + layout.element(&e.nodes, e.dofs); } for e in &trusses { layout.element(&e.2, fieldcraft_fem::ElementDofs::new(translations)); @@ -546,10 +585,22 @@ impl Mesh { warnings, discrete: crate::discrete::Discrete::default(), }; - for (label, kind, nodes, formulation, d, density, theory, _, expansion, hyperelastic) in elements { - let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); + for e in elements { + let idx: Vec = e.nodes.iter().filter_map(|n| mesh.node(*n)).collect(); let exyz: Vec<[f64; 3]> = idx.iter().filter_map(|i| mesh.xyz.get(*i).copied()).collect(); - mesh.elements.push(Elem { label, kind, nodes: idx, xyz: exyz, formulation, d, density, theory, expansion, hyperelastic }); + mesh.elements.push(Elem { + label: e.label, + kind: e.kind, + nodes: idx, + xyz: exyz, + formulation: e.formulation, + d: e.d, + density: e.density, + theory: e.theory, + expansion: e.expansion, + hyperelastic: e.hyperelastic, + plastic: e.plastic, + }); } for (label, kind, nodes, young, area, density, expansion) in trusses { let idx: Vec = nodes.iter().filter_map(|n| mesh.node(*n)).collect(); @@ -940,15 +991,15 @@ fn body(f: &mut [f64], mesh: &Mesh, e: &Elem, b: [f64; 3]) -> Result<()> { Ok(()) } -/// The linear solver uses the elastic constants alone: a warning for each material in use whose -/// plasticity, damage or erosion it leaves out. +/// The linear solvers use the elastic constants alone: a warning for each material in use whose +/// plasticity, damage or erosion they leave out. pub fn inelastic_warnings(model: &Model) -> Vec { let used: BTreeSet<&str> = model.sections().values().map(|s| s.material.as_str()).collect(); used.into_iter() .filter(|name| model.material(name).is_some_and(|m| m.plastic.is_some() || m.damage.is_some() || m.erosion.is_some())) .map(|name| { format!( - "material `{name}`: the linear solver uses its elastic constants only (its plasticity, damage and erosion act in explicit dynamics)" + "material `{name}`: the linear solver uses its elastic constants only (its plasticity acts in nonlinear static steps and explicit dynamics, its damage and erosion in explicit dynamics)" ) }) .collect() diff --git a/crates/solver/tests/plasticity.rs b/crates/solver/tests/plasticity.rs new file mode 100644 index 0000000..8fd200e --- /dev/null +++ b/crates/solver/tests/plasticity.rs @@ -0,0 +1,249 @@ +//! Implicit J2 plasticity in nonlinear static steps against closed forms: uniaxial tension at +//! finite strain (Kirchhoff stress against the logarithmic strain) and at small strain in every +//! element kind, the Bauschinger effect of linear kinematic hardening on a load path that comes +//! back, Newton's quadratic rate with the consistent tangent, and what plastic steps refuse or +//! stop at. Run with `--nocapture` to see the residuals. + +mod common; + +use common::*; +use fieldcraft_model::*; +use fieldcraft_solver::{Limits, StaticResult, solve_nonlinear}; + +const E: f64 = 200_000.0; +const NU: f64 = 0.3; +const Y0: f64 = 250.0; +const H: f64 = 2000.0; + +fn bilinear(hardening: Hardening) -> Plastic { + Plastic::Bilinear { yield_stress: Y0, hardening_modulus: H, hardening } +} + +fn plastic(g: &mut Grid, p: Plastic) { + let m = Material { elastic: Some(iso(E, NU)), plastic: Some(p), ..Material::default() }; + g.model.replace_material("Mat", m).unwrap(); +} + +/// A unit cube on three symmetry planes, its x = 1 face pulled by `stretch` (uniaxial stress). +fn cube(kind: ElementKind, formulation: Formulation, step: NonlinearStatic, p: Plastic, stretch: f64) -> Grid { + let mut g = grid([2, 2, 2], kind, |x| x); + g.define(iso(E, NU), None, formulation, Procedure::NonlinearStatic(step)); + plastic(&mut g, p); + g.symmetry("x0", Symmetry::X, |x| near(x[0], 0.0)); + g.symmetry("y0", Symmetry::Y, |x| near(x[1], 0.0)); + g.symmetry("z0", Symmetry::Z, |x| near(x[2], 0.0)); + g.node_set("x1", |x| near(x[0], 1.0)); + g.bc("pull", BoundaryCondition::Displacement { target: Target::Set("x1".into()), values: [Some(stretch), None, None, None, None, None] }); + g +} + +fn solve(g: &Grid) -> StaticResult { + solve_nonlinear(&g.model, Some("Step"), &Limits::default()).unwrap() +} + +/// The uniaxial stress of linear hardening at strain ε past yield: σ_y + E_t (ε − σ_y/E). +fn uniaxial(strain: f64) -> f64 { + let et = E * H / (E + H); + Y0 + et * (strain - Y0 / E) +} + +fn hybrid() -> Formulation { + Formulation { hybrid: true, ..Formulation::default() } +} + +/// Every element kind (and hybrid hex8, hex20, tet10) stretched to a logarithmic strain of 0.2: +/// the Kirchhoff stress τ = J σ follows the small-strain curve in the logarithmic strain +/// exactly (Simo's exponential map integrates a proportional path exactly), and the +/// equivalent plastic strain is (τ − σ_y)/H. +#[test] +fn finite_strain_uniaxial_tension_follows_the_logarithmic_closed_form() { + let strain: f64 = 0.2; + let tau = uniaxial(strain); + let mut cases: Vec<(ElementKind, Formulation, &str)> = variants().into_iter().filter(|v| v.1 != incompatible()).collect(); + cases.extend([ + (ElementKind::Hex8, hybrid(), "hex8 hybrid"), + (ElementKind::Hex20, hybrid(), "hex20 hybrid"), + (ElementKind::Tet10, hybrid(), "tet10 hybrid"), + ]); + for (kind, formulation, name) in cases { + let g = cube(kind, formulation, NonlinearStatic::default(), bilinear(Hardening::Isotropic), strain.exp() - 1.0); + let r = solve(&g); + for e in &r.elements { + for (p, ep) in e.points.iter().zip(&e.plastic_strain) { + // J from the logarithmic strains, σ the Cauchy stress. + let j = (p.strain[0] + p.strain[1] + p.strain[2]).exp(); + assert!((p.stress[0] * j - tau).abs() < 1e-7 * tau, "{name}: τ = {} vs {tau}", p.stress[0] * j); + assert!(p.stress[1].abs() < 1e-7 * tau && p.stress[3].abs() < 1e-7 * tau, "{name}: {:?}", p.stress); + assert!((ep - (tau - Y0) / H).abs() < 1e-9, "{name}: ε_p = {ep}"); + assert!((p.strain[0] - strain).abs() < 1e-12, "{name}: {:?}", p.strain); + } + assert_eq!(e.plastic_strain.len(), e.points.len(), "{name}"); + } + let h = r.nonlinear.as_ref().unwrap(); + eprintln!("{name}: {} increments, {} iterations", h.increments.len(), h.iterations); + } +} + +/// Without nonlinear geometry: σ = σ_y + E_t (ε − σ_y/E) in the small strain, isotropic and +/// kinematic hardening alike under monotonic loading, hybrid elements too. +#[test] +fn small_strain_uniaxial_tension_follows_the_closed_form() { + let strain = 0.01; + let sigma = uniaxial(strain); + let small = NonlinearStatic { nonlinear_geometry: false, ..NonlinearStatic::default() }; + for hardening in [Hardening::Isotropic, Hardening::Kinematic] { + for (kind, formulation, name) in [ + (ElementKind::Hex8, full(), "hex8"), + (ElementKind::Tet10, full(), "tet10"), + (ElementKind::Hex20, reduced(), "hex20 reduced"), + (ElementKind::Hex20, hybrid(), "hex20 hybrid"), + ] { + let r = solve(&cube(kind, formulation, small.clone(), bilinear(hardening), strain)); + for e in &r.elements { + for (p, ep) in e.points.iter().zip(&e.plastic_strain) { + assert!((p.stress[0] - sigma).abs() < 1e-8 * sigma && p.stress[1].abs() < 1e-8 * sigma, "{name} {hardening:?}: {:?}", p.stress); + assert!((ep - (sigma - Y0) / H).abs() < 1e-10, "{name}: {ep}"); + } + } + } + } +} + +/// The Bauschinger effect: pulled to ε_a and brought back to zero strain on a load path, +/// linear kinematic hardening yields again in reverse at σ_a − 2σ_y and ends at the residual +/// stress σ_a − 2σ_y − E_t (ε_a − 2σ_y/E); isotropic hardening yields again only at −σ_a and +/// ends at −σ_a − E_t (ε_a − 2σ_a/E). +#[test] +fn kinematic_hardening_shows_the_bauschinger_effect_on_unloading() { + let strain = 0.01; + let sa = uniaxial(strain); + let et = E * H / (E + H); + let path = NonlinearStatic { + nonlinear_geometry: false, + initial_increment: 0.05, + max_increment: 0.05, + load_path: Some(vec![[0.0, 0.0], [0.5, 1.0], [1.0, 0.0]]), + ..NonlinearStatic::default() + }; + for (hardening, want) in + [(Hardening::Kinematic, sa - 2.0 * Y0 - et * (strain - 2.0 * Y0 / E)), (Hardening::Isotropic, -sa - et * (strain - 2.0 * sa / E))] + { + let r = solve(&cube(ElementKind::Hex8, full(), path.clone(), bilinear(hardening), strain)); + let s = r.elements[0].points[0].stress[0]; + assert!((s - want).abs() < 1e-8 * sa, "{hardening:?}: residual stress {s} vs {want}"); + assert!(r.displacement.iter().all(|u| u[0].abs() < 1e-12), "back at zero strain"); + let h = r.nonlinear.unwrap(); + assert!(h.increments.iter().any(|i| i.converged && i.load_factor == 1.0), "the path's turn is an increment's end"); + assert_eq!(h.frames.last().map(|f| f.load_factor), Some(0.0)); + eprintln!("{hardening:?}: residual stress {s:.6} (closed form {want:.6})"); + } +} + +/// Newton with the consistent tangent converges quadratically in a plastic structure: a notched +/// bar pulled past general yield, and a cantilever bent into plastic flow at finite strain. +#[test] +fn plastic_newton_converges_quadratically() { + for nonlinear_geometry in [false, true] { + let mut step = NonlinearStatic { initial_increment: 0.5, nonlinear_geometry, ..NonlinearStatic::default() }; + step.tolerances.force = 1e-8; + step.tolerances.displacement = 1e-8; + step.tolerances.energy = 1e-16; + let mut g = grid([8, 1, 2], ElementKind::Hex20, |q| [10.0 * q[0], q[1], q[2]]); + g.define(iso(E, NU), None, full(), Procedure::NonlinearStatic(step)); + plastic(&mut g, bilinear(Hardening::Isotropic)); + g.fixed("root", |q| near(q[0], 0.0)); + g.node_set("tip", |q| near(q[0], 1.0)); + // A tip deflection of ten times the first yield's. + let first_yield = 2.0 * Y0 * 100.0 / (3.0 * E * 2.0); + g.bc( + "push", + BoundaryCondition::Displacement { target: Target::Set("tip".into()), values: [None, None, Some(-10.0 * first_yield), None, None, None] }, + ); + let r = solve(&g); + assert!(r.elements.iter().flat_map(|e| &e.plastic_strain).any(|p| *p > 0.0), "plastic"); + let h = r.nonlinear.unwrap(); + let mut best = 0.0f64; + for inc in h.increments.iter().filter(|i| i.converged) { + let res: Vec = inc.iterations.iter().map(|i| i.force).collect(); + let orders: Vec = res.windows(3).filter(|w| w[0] < 1e-2 && w[2] > 1e-11).map(|w| (w[2] / w[1]).ln() / (w[1] / w[0]).ln()).collect(); + eprintln!("geometry {nonlinear_geometry}: increment {} (λ = {}): residuals {res:?}, orders {orders:?}", inc.number, inc.load_factor); + best = orders.iter().fold(best, |m, o| m.max(*o)); + } + assert!(best > 1.8, "nonlinear geometry {nonlinear_geometry}: no quadratic convergence (best order {best})"); + } +} + +/// What a static step cannot take is refused by name; a material that softens faster than it +/// unloads makes the step cut back and stop cleanly; the results do not depend on the threads. +#[test] +fn plastic_steps_refuse_and_fail_cleanly() { + let jc = |c: f64| Plastic::JohnsonCook { + a: 350.0, + b: 275.0, + n: 0.36, + c, + m: 1.0, + reference_strain_rate: 1.0, + room_temperature: 293.0, + melt_temperature: 1356.0, + specific_heat: Some(4.5e8), + taylor_quinney: 0.9, + }; + let err = solve_nonlinear(&cube(ElementKind::Hex8, full(), NonlinearStatic::default(), jc(0.02), 0.01).model, None, &Limits::default()) + .unwrap_err() + .to_string(); + assert!(err.contains("strain-rate term") && err.contains("C = 0.02"), "{err}"); + // C = 0 is the rate-independent law: it solves, warning that a static step does not heat. + let r = solve(&cube(ElementKind::Hex8, full(), NonlinearStatic::default(), jc(0.0), 0.01)); + assert!(r.warnings.iter().any(|w| w.contains("isothermal")), "{:?}", r.warnings); + let err = solve_nonlinear( + &cube(ElementKind::Hex8, full(), NonlinearStatic::default(), bilinear(Hardening::Kinematic), 0.01).model, + None, + &Limits::default(), + ) + .unwrap_err() + .to_string(); + assert!(err.contains("kinematic hardening") && err.contains("finite strain"), "{err}"); + let small = NonlinearStatic { nonlinear_geometry: false, ..NonlinearStatic::default() }; + let err = solve_nonlinear( + &cube(ElementKind::Hex8, incompatible(), small.clone(), bilinear(Hardening::Isotropic), 0.01).model, + None, + &Limits::default(), + ) + .unwrap_err() + .to_string(); + assert!(err.contains("incompatible modes with an inelastic material"), "{err}"); + // Softening steeper than −3G past a small plastic strain: the point localises, the step cuts + // back to its minimum increment and stops, naming the element. + let g_mod = E / (2.0 * (1.0 + NU)); + let softening = + Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: Y0 }, HardeningPoint { strain: 1e-4, stress: Y0 - 4.0 * g_mod * 1e-4 }] }; + let soft = NonlinearStatic { min_increment: 1e-3, ..small.clone() }; + let err = solve_nonlinear(&cube(ElementKind::Hex8, full(), soft, softening, 0.01).model, None, &Limits::default()).unwrap_err().to_string(); + assert!(err.contains("stopped at load factor") && err.contains("element"), "{err}"); + // Damage is left out with a warning. + let mut g = cube(ElementKind::Hex8, full(), small, bilinear(Hardening::Isotropic), 0.01); + let mut m = g.model.material("Mat").unwrap().clone(); + m.damage = Some(Damage::JohnsonCook { d1: 0.5, d2: 0.0, d3: 0.0, d4: 0.0, d5: 0.0 }); + g.model.replace_material("Mat", m).unwrap(); + let r = solve(&g); + assert!(r.warnings.iter().any(|w| w.contains("damage and erosion")), "{:?}", r.warnings); + // Hostile plasticity never reaches the solver: non-finite or negative tables and yield. + let mut g = cube(ElementKind::Hex8, full(), NonlinearStatic::default(), bilinear(Hardening::Isotropic), 0.01); + for bad in [ + Plastic::Bilinear { yield_stress: -250.0, hardening_modulus: H, hardening: Hardening::Isotropic }, + Plastic::Bilinear { yield_stress: Y0, hardening_modulus: f64::NAN, hardening: Hardening::Isotropic }, + Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: f64::INFINITY }] }, + Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: Y0 }, HardeningPoint { strain: f64::NAN, stress: Y0 }] }, + Plastic::Table { points: vec![HardeningPoint { strain: 0.0, stress: Y0 }, HardeningPoint { strain: 0.1, stress: -1.0 }] }, + ] { + let m = Material { elastic: Some(iso(E, NU)), plastic: Some(bad.clone()), ..Material::default() }; + assert!(g.model.replace_material("Mat", m).is_err(), "{bad:?}"); + } + // Bit for bit on any number of threads. + let g = cube(ElementKind::Tet10, full(), NonlinearStatic::default(), bilinear(Hardening::Isotropic), 0.1); + let a = fieldcraft_parallel::with_threads(1, || solve(&g)); + for threads in [2, 7] { + assert_eq!(fieldcraft_parallel::with_threads(threads, || solve(&g)), a, "{threads} threads"); + } +}