diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index ac4e0e1..28a0b42 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -30,7 +30,14 @@ //! - a frequency step is a perturbation of the last general step: its boundary conditions are //! written against that base state and do not carry on to later steps; loads have no effect on //! natural frequencies and are not written in it; -//! - every step restates its output request (`*OUTPUT, FIELD, OP=NEW`). +//! - every step restates its output request (`*OUTPUT, FIELD, OP=NEW`); +//! - a nonlinear static step is `*STEP, NLGEOM=YES` (or `NO` without nonlinear geometry) with +//! `INC=` its increment cap, `*STATIC` with its initial, minimum and maximum increments over a +//! step period of 1, and `*CONTROLS, PARAMETERS=LINE SEARCH` when it searches; its iteration cap +//! 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 //! @@ -96,7 +103,10 @@ //! //! # Materials //! -//! Isotropic elasticity is `*ELASTIC, TYPE=ISOTROPIC`, engineering constants +//! Isotropic elasticity is `*ELASTIC, TYPE=ISOTROPIC` (or, for a neo-Hookean material in a deck +//! with nonlinear steps, `*HYPERELASTIC, NEO HOOKE` with C₁₀ = μ/2 and D₁ = 2/κ, the same strain +//! energy; a linear static step is then refused, since Abaqus's general static step would apply +//! the hyperelastic law), engineering constants //! `TYPE=ENGINEERING CONSTANTS`, density `*DENSITY`. Von Mises plasticity is `*PLASTIC`: a //! hardening table as its (yield stress, plastic strain) lines; bilinear isotropic hardening as //! two lines, the second at a plastic strain of [`BILINEAR_REACH`] (Abaqus holds the last yield @@ -219,6 +229,12 @@ pub fn element_type(kind: ElementKind, section: &Section) -> std::result::Result Ok(if f.hybrid { format!("{plain}H") } else { plain.to_string() }) } +/// Whether the model has a nonlinear static step: hyperelastic materials are then written as +/// `*HYPERELASTIC` (in a deck of linear steps their elastic constants are the whole of it). +fn has_nonlinear_step(model: &Model) -> bool { + model.steps().iter().any(|s| matches!(s.procedure, Procedure::NonlinearStatic(_))) +} + /// Every element, grouped by element type (from its section) and sorted by label. fn elements_by_type<'m>(model: &'m Model, r: &Resolved) -> Result>> { let mut by_type: BTreeMap> = BTreeMap::new(); @@ -287,6 +303,16 @@ fn check_expressible(model: &Model) -> Result<()> { for step in model.steps().iter().filter(|s| matches!(s.procedure, Procedure::ExplicitDynamic(_))) { bad.insert(format!("step `{}`: explicit dynamic steps (not written to Abaqus input yet)", step.name)); } + if has_nonlinear_step(model) + && let Some(material) = model.hyperelastic_in_use().into_iter().next() + { + for step in model.steps().iter().filter(|s| matches!(s.procedure, Procedure::Static {})) { + bad.insert(format!( + "step `{}`: a linear static step beside nonlinear steps of hyperelastic material `{material}` (Abaqus's general static step applies the hyperelastic law; write the linear step to its own deck)", + step.name + )); + } + } for (name, m) in model.materials() { if m.damage.is_some() || m.erosion.is_some() { bad.insert(format!( @@ -759,6 +785,12 @@ impl Deck { for (name, m) in model.materials() { self.line(&format!("*MATERIAL, NAME={name}")); match m.elastic { + Some(Elastic::Isotropic { young, poisson }) if m.hyperelastic.is_some() && has_nonlinear_step(model) => { + // W = C10 (Ī1 − 3) + (J − 1)²/D1 with C10 = μ/2 and D1 = 2/κ. + let (shear, bulk) = fieldcraft_model::Hyperelastic::moduli(young, poisson); + self.line("*HYPERELASTIC, NEO HOOKE"); + self.numbers_line("", &[shear / 2.0, 2.0 / bulk])?; + } Some(Elastic::Isotropic { young, poisson }) => { self.line("*ELASTIC, TYPE=ISOTROPIC"); self.numbers_line("", &[young, poisson])?; @@ -1038,10 +1070,22 @@ impl Deck { self.comment(&format!("Step {}", step.name)); match step.procedure { Procedure::Fluid { .. } => return Err(crate::fluid_step(FORMAT, &step.name)), - Procedure::Static {} => { - self.line(&format!("*STEP, NAME={}, NLGEOM=NO", step.name)); - self.line("*STATIC"); - self.line("1.0, 1.0"); + Procedure::Static {} | Procedure::NonlinearStatic(_) => { + if let Procedure::NonlinearStatic(n) = &step.procedure { + let nlgeom = if n.nonlinear_geometry { "YES" } else { "NO" }; + self.line(&format!("*STEP, NAME={}, NLGEOM={nlgeom}, INC={}", step.name, n.max_increments)); + self.line("*STATIC"); + self.numbers_line("", &[n.initial_increment, 1.0, n.min_increment, n.max_increment])?; + if n.line_search { + self.line("*CONTROLS, PARAMETERS=LINE SEARCH"); + self.line("5,"); + } + self.comment("iterations and convergence: Abaqus's own controls (its residual norms differ from FieldCraft's)"); + } else { + self.line(&format!("*STEP, NAME={}, NLGEOM=NO", step.name)); + self.line("*STATIC"); + self.line("1.0, 1.0"); + } self.step_boundaries(model, r, &base_bcs, &step.boundary_conditions)?; if step.loads == base_loads { self.comment("loads: as in the previous general step"); diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index 5741c9c..50db6d9 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -24,7 +24,13 @@ //! the set's translations, equal weights); //! - 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`); boundary conditions +//! 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 +//! 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 +//! LS-DYNA's defaults; a step without nonlinear geometry and hyperelastic materials in nonlinear +//! steps are refused); boundary conditions //! as `*BOUNDARY_SPC_*_ID` (named) and `*BOUNDARY_PRESCRIBED_MOTION_*_ID` for non-zero //! displacements; forces, moments and tractions as `*LOAD_NODE_SET`/`_POINT` from the //! resolved nodal values; pressure as `*LOAD_SEGMENT_SET`; gravity as `*LOAD_BODY_X/Y/Z` @@ -267,7 +273,21 @@ fn plan<'a>(model: &'a Model, r: &Resolved) -> Result> { if matches!(step.procedure, Procedure::ExplicitDynamic(_)) { refused.insert(format!("step `{}`: explicit dynamic steps (this writer writes implicit analyses so far)", step.name)); } - let is_static = matches!(step.procedure, Procedure::Static {}); + if let Procedure::NonlinearStatic(n) = &step.procedure { + if !n.nonlinear_geometry { + refused.insert(format!( + "step `{}`: a nonlinear step without nonlinear geometry (LS-DYNA's nonlinear implicit analysis always follows the deformed shape)", + step.name + )); + } + for material in model.hyperelastic_in_use() { + refused.insert(format!( + "step `{}`: hyperelastic material `{material}` in a nonlinear step (LS-DYNA's rubber materials define their volumetric energy differently; not written yet)", + step.name + )); + } + } + let is_static = matches!(step.procedure, Procedure::Static {} | Procedure::NonlinearStatic(_)); for b in &step.boundary_conditions { let Some(bc) = model.boundary_conditions().get(b) else { continue }; match bc.target() { @@ -793,18 +813,29 @@ fn dof_code(dof: u8) -> i32 { fn write_step(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan, step: &Step, n: u32, constraint_id: &mut u32) -> Result<()> { d.line("*TITLE"); d.line(&step.name); - let is_static = matches!(step.procedure, Procedure::Static {}); + let is_static = matches!(step.procedure, Procedure::Static {} | Procedure::NonlinearStatic(_)); + let first_increment = if let Procedure::NonlinearStatic(n) = &step.procedure { n.initial_increment } else { 1.0 }; d.line("*CONTROL_IMPLICIT_GENERAL"); d.header(&["imflag", "dt0"]); - d.card([Field::from(1), Field::from(1.0)]); - match step.procedure { + d.card([Field::from(1), Field::from(first_increment)]); + match &step.procedure { Procedure::Fluid { .. } => return Err(crate::fluid_step(FORMAT, &step.name)), Procedure::Static {} => { d.line("*CONTROL_IMPLICIT_SOLUTION"); d.header(&["nsolvr"]); d.card([Field::from(1)]); } - Procedure::Frequency { modes, min_hz, max_hz } => { + Procedure::NonlinearStatic(n) => { + // Nonlinear with a stiffness reformation every iteration (full Newton), at most + // `max_iterations` per increment; LS-DYNA's own convergence norms and line search. + d.line("*CONTROL_IMPLICIT_SOLUTION"); + d.header(&["nsolvr", "ilimit", "maxref"]); + d.card([Field::from(2), Field::from(1), Field::from(n.max_iterations)]); + d.line("*CONTROL_IMPLICIT_AUTO"); + d.header(&["iauto", "iteopt", "itewin", "dtmin", "dtmax"]); + d.card([Field::from(1), Field::Blank, Field::Blank, Field::from(n.min_increment), Field::from(n.max_increment)]); + } + &Procedure::Frequency { modes, min_hz, max_hz } => { d.line("*CONTROL_IMPLICIT_EIGENVALUE"); d.header(&["neig", "center", "lflag", "lftend", "rflag", "rhtend"]); let bound = |b: Option| match b { diff --git a/crates/decks/src/nastran/mod.rs b/crates/decks/src/nastran/mod.rs index a1ecbe1..28d171c 100644 --- a/crates/decks/src/nastran/mod.rs +++ b/crates/decks/src/nastran/mod.rs @@ -280,6 +280,9 @@ fn solution(steps: &[&Step]) -> std::result::Result { if let Some(s) = steps.iter().find(|s| matches!(s.procedure, Procedure::ExplicitDynamic(_))) { return Err(format!("step `{}`: explicit dynamic steps (Nastran's linear solutions have none)", s.name)); } + if let Some(s) = steps.iter().find(|s| matches!(s.procedure, Procedure::NonlinearStatic(_))) { + return Err(format!("step `{}`: nonlinear static steps (Nastran's SOL 400; this writer writes SOL 101 and 103 so far)", s.name)); + } let statics = steps.iter().find(|s| matches!(s.procedure, Procedure::Static {})); let modes = steps.iter().find(|s| matches!(s.procedure, Procedure::Frequency { .. })); match (statics, modes) { diff --git a/crates/decks/tests/abaqus.rs b/crates/decks/tests/abaqus.rs index e21f552..33638d4 100644 --- a/crates/decks/tests/abaqus.rs +++ b/crates/decks/tests/abaqus.rs @@ -753,3 +753,45 @@ fn plasticity_is_written_and_failure_reported() { let Err(Error::Unsupported { what, .. }) = abaqus::write(&failing) else { panic!("damage is not written") }; assert!(what.contains("material `Steel`: damage and erosion"), "{what}"); } + +/// A nonlinear static step is a general step with NLGEOM, its increments and its line search; a +/// neo-Hookean material is then `*HYPERELASTIC` and a linear static step beside it is refused. +#[test] +fn nonlinear_steps_and_hyperelastic_materials_are_written() { + use fieldcraft_model::{Hyperelastic, NonlinearStatic}; + let mut m = cantilever(); + let controls = NonlinearStatic { + initial_increment: 0.05, + min_increment: 1e-4, + max_increment: 0.5, + max_increments: 40, + line_search: true, + ..NonlinearStatic::default() + }; + m.replace_step(step("Bend", Procedure::NonlinearStatic(controls.clone()), &["clamp"], &["end"])).unwrap(); + let deck = abaqus::write(&m).unwrap(); + let bend = block(&deck, "Bend"); + assert_eq!(bend.first(), Some(&"*STEP, NAME=Bend, NLGEOM=YES, INC=40"), "{bend:?}"); + let at = bend.iter().position(|l| *l == "*STATIC").unwrap(); + assert_eq!(bend.get(at + 1..at + 4), Some(&["0.05, 1.0, 0.0001, 0.5", "*CONTROLS, PARAMETERS=LINE SEARCH", "5,"][..]), "{bend:?}"); + assert!(deck.contains("*ELASTIC, TYPE=ISOTROPIC\n210000.0, 0.3\n"), "St Venant–Kirchhoff is the elastic material"); + + let small = NonlinearStatic { nonlinear_geometry: false, line_search: false, ..controls }; + m.replace_step(step("Bend", Procedure::NonlinearStatic(small), &["clamp"], &["end"])).unwrap(); + m.replace_material("Steel", Material { hyperelastic: Some(Hyperelastic::NeoHookean {}), ..steel() }).unwrap(); + let deck = abaqus::write(&m).unwrap(); + let bend = block(&deck, "Bend"); + assert_eq!(bend.first(), Some(&"*STEP, NAME=Bend, NLGEOM=NO, INC=40"), "{bend:?}"); + assert!(!bend.iter().any(|l| l.starts_with("*CONTROLS")), "{bend:?}"); + // μ = E/(2(1 + ν)) = 80769.23…, κ = E/(3(1 − 2ν)) = 175000: C10 = μ/2, D1 = 2/κ. + let (c10, d1) = (210_000.0 / 2.6 / 2.0, 2.0 / 175_000.0); + let line = deck.split("*HYPERELASTIC, NEO HOOKE\n").nth(1).and_then(|s| s.lines().next()).unwrap(); + let values: Vec = line.split(", ").map(|v| v.parse().unwrap()).collect(); + assert!((values[0] - c10).abs() < 1e-9 * c10 && (values[1] - d1).abs() < 1e-9 * d1, "{line}"); + assert!(!deck.contains("*ELASTIC"), "{deck}"); + + m.add_step(step("Linear", Procedure::Static {}, &["clamp"], &["end"]), None).unwrap(); + assert!(unsupported(&m).contains("step `Linear`: a linear static step beside nonlinear steps of hyperelastic material `Steel`")); + 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"); +} diff --git a/crates/decks/tests/lsdyna.rs b/crates/decks/tests/lsdyna.rs index 99eb830..302cc34 100644 --- a/crates/decks/tests/lsdyna.rs +++ b/crates/decks/tests/lsdyna.rs @@ -737,3 +737,37 @@ fn plasticity_is_written_as_plastic_materials() { let Err(Error::Unsupported { what, .. }) = lsdyna::write(&with(jc, None)) else { panic!() }; assert!(what.contains("Johnson–Cook"), "{what}"); } + +/// A nonlinear static step: the ramped implicit run with full Newton and automatic increments; +/// what LS-DYNA cannot take is named. +#[test] +fn nonlinear_static_steps() { + let mut m = cantilever(); + let controls = + NonlinearStatic { initial_increment: 0.05, min_increment: 1e-4, max_increment: 0.5, max_iterations: 12, ..NonlinearStatic::default() }; + m.replace_step(step("Bend", Procedure::NonlinearStatic(controls.clone()), &["clamp"], &["tip_load"])).unwrap(); + let deck = lsdyna::write(&m).unwrap().text; + let after = |keyword: &str| -> Vec { + let mut lines = deck.lines().skip_while(|l| *l != keyword).skip(1).filter(|l| !l.starts_with('$')); + lines.next().unwrap().split_whitespace().map(str::to_string).collect() + }; + let num = |s: &String| s.parse::().unwrap(); + let general = after("*CONTROL_IMPLICIT_GENERAL"); + assert_eq!((num(&general[0]), num(&general[1])), (1.0, 0.05), "{general:?}"); + let solution = after("*CONTROL_IMPLICIT_SOLUTION"); + assert_eq!(solution.iter().map(num).collect::>(), [2.0, 1.0, 12.0], "{solution:?}"); + let auto = after("*CONTROL_IMPLICIT_AUTO"); + 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}"); + + m.replace_step(step("Bend", Procedure::NonlinearStatic(NonlinearStatic { nonlinear_geometry: false, ..controls }), &["clamp"], &["tip_load"])) + .unwrap(); + unsupported(&m, "a nonlinear step without nonlinear geometry"); + m.replace_step(step("Bend", Procedure::NonlinearStatic(NonlinearStatic::default()), &["clamp"], &["tip_load"])).unwrap(); + m.replace_material("Steel", Material { hyperelastic: Some(Hyperelastic::NeoHookean {}), ..steel() }).unwrap(); + unsupported(&m, "hyperelastic material `Steel` in a nonlinear step"); + // In a linear step the elastic constants are the whole of a hyperelastic material. + m.replace_step(step("Bend", Procedure::Static {}, &["clamp"], &["tip_load"])).unwrap(); + assert!(lsdyna::write(&m).unwrap().text.contains("*MAT_ELASTIC_TITLE")); +} diff --git a/crates/decks/tests/nastran.rs b/crates/decks/tests/nastran.rs index 31785a2..39fdb6f 100644 --- a/crates/decks/tests/nastran.rs +++ b/crates/decks/tests/nastran.rs @@ -782,3 +782,13 @@ fn plasticity_is_written_as_mats1() { let Err(Error::Unsupported { what, .. }) = nastran::write(&m) else { panic!() }; assert!(what.contains("Johnson–Cook plasticity") && what.contains("damage and erosion"), "{what}"); } + +#[test] +fn nonlinear_static_steps_are_refused_by_name() { + let mut m = cantilever(); + let name = m.steps()[0].name.clone(); + let mut s = m.steps()[0].clone(); + s.set_procedure(Procedure::NonlinearStatic(fieldcraft_model::NonlinearStatic::default())); + m.replace_step(s).unwrap(); + assert!(unsupported(&m).contains(&format!("step `{name}`: nonlinear static steps (Nastran's SOL 400")), "{}", unsupported(&m)); +} diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs index ff86bee..5cce949 100644 --- a/crates/engine/src/cmd/define.rs +++ b/crates/engine/src/cmd/define.rs @@ -81,7 +81,7 @@ definition_commands!( Material, "material", "Material", - r#"{"name": "Steel", "elastic": {"type": "isotropic", "young": 210000, "poisson": 0.3} | {"type": "engineering_constants", "young": [E1, E2, E3], "poisson": [ν12, ν13, ν23], "shear": [G12, G13, G23]}, "density": 7.85e-9, "fluid": {"type": "newtonian", "viscosity": 1.0e-9} (dynamic viscosity; a fluid needs a density), "plastic": {…} (see model.material.plasticity), "damage": {…}, "erosion": {…} (see model.material.damage)} — in the model's unit system (model.units), never converted"#, + r#"{"name": "Steel", "elastic": {"type": "isotropic", "young": 210000, "poisson": 0.3} | {"type": "engineering_constants", "young": [E1, E2, E3], "poisson": [ν12, ν13, ν23], "shear": [G12, G13, G23]}, "density": 7.85e-9, "fluid": {"type": "newtonian", "viscosity": 1.0e-9} (dynamic viscosity; a fluid needs a density), "plastic": {…} (see model.material.plasticity), "damage": {…}, "erosion": {…} (see model.material.damage), "hyperelastic": {"type": "neo_hookean"} (large strain in nonlinear steps, from the isotropic constants: μ = E/(2(1 + ν)), κ = E/(3(1 − 2ν)); without it an elastic material is St Venant–Kirchhoff there)} — in the model's unit system (model.units), never converted"#, add_material, replace_material, remove_material, @@ -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": "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} (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)}"#, kind: CommandKind::Model, run: step_create, }; diff --git a/crates/engine/src/cmd/solve.rs b/crates/engine/src/cmd/solve.rs index 2e40549..a84f39c 100644 --- a/crates/engine/src/cmd/solve.rs +++ b/crates/engine/src/cmd/solve.rs @@ -142,6 +142,7 @@ fn run_all(s: &mut Session, _: &Value) -> Result { Procedure::Frequency { .. } => run_frequency, Procedure::ExplicitDynamic(_) => run_explicit, Procedure::Fluid { .. } => run_fluid, + Procedure::NonlinearStatic(_) => run_nonlinear, }; (st.name.clone(), st.procedure.name(), solve) }) @@ -193,6 +194,10 @@ fn run_static(s: &mut Session, p: &Value) -> Result { Ok(summary) } +fn run_nonlinear(_: &mut Session, _: &Value) -> Result { + Err(EngineError::Params("not supported yet: nonlinear static steps (the nonlinear solver arrives with solve.nonlinear)".into())) +} + fn run_heat(s: &mut Session, p: &Value) -> Result { let limits = limits(p)?; let step = opt::(p, "step")?; diff --git a/crates/explicit/src/system.rs b/crates/explicit/src/system.rs index 2cdc773..e70677c 100644 --- a/crates/explicit/src/system.rs +++ b/crates/explicit/src/system.rs @@ -139,6 +139,9 @@ fn invalid(m: impl Into) -> Error { fn material(model: &Model, name: &str) -> Result { let m = model.material(name).ok_or_else(|| invalid(format!("material `{name}` is missing")))?; let elastic = m.elastic.as_ref().ok_or_else(|| invalid(format!("material `{name}` has no elastic constants")))?; + if let Some(h) = m.hyperelastic { + return Err(Error::NotSupported(format!("{} hyperelastic material `{name}` in explicit dynamics", h.name()))); + } let density = m .density .filter(|d| *d > 0.0 && d.is_finite()) diff --git a/crates/model/src/check.rs b/crates/model/src/check.rs index 5121442..e279e33 100644 --- a/crates/model/src/check.rs +++ b/crates/model/src/check.rs @@ -322,6 +322,7 @@ impl Model { Procedure::ExplicitDynamic(_) => "is dynamic", Procedure::HeatTransfer { transient: Some(_) } => "is transient heat transfer", Procedure::Static {} + | Procedure::NonlinearStatic(_) | Procedure::Fluid { .. } | Procedure::HeatTransfer { transient: None } | Procedure::CoupledTemperatureDisplacement {} => { @@ -469,7 +470,7 @@ impl Model { } } match &s.procedure { - Procedure::Static {} => { + Procedure::Static {} | Procedure::NonlinearStatic(_) => { if s.loads.is_empty() && !prescribes_motion { d.warning(&subject, "nothing loads the model in this step"); } diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs index 0a1b232..9726c46 100644 --- a/crates/model/src/lib.rs +++ b/crates/model/src/lib.rs @@ -57,12 +57,13 @@ pub use discrete::{Direction, Mass, PointMass, ResolvedDirection, ResolvedMass, pub use fluids::{BoxFace, BoxFaces, FaceCondition, FluidBoundary, FluidDomain, Profile, SurfaceBoundary, SurfaceCondition}; pub use geometry::{Classification, Entity, GeometryRef}; pub use loads::{Distribution, Load, centrifugal_acceleration}; -pub use materials::{Damage, Directional, Elastic, Erosion, Fluid, Hardening, HardeningPoint, Material, Plastic, Thermal}; +pub use materials::{Damage, Directional, Elastic, Erosion, Fluid, Hardening, HardeningPoint, Hyperelastic, Material, Plastic, Thermal}; pub use mesh::{Element, ElementEdge, ElementFace, ElementKind, Node, check_coordinates}; pub use sections::{BeamProfile, Extent, Formulation, Integration, Section, SectionKind}; pub use sets::Set; pub use steps::{ - BulkViscosity, ExplicitDynamic, Hourglass, HourglassKind, InitialVelocity, OutputVariable, Procedure, RigidWall, Step, StressRate, Transient, + BulkViscosity, ExplicitDynamic, Hourglass, HourglassKind, InitialVelocity, MAX_NONLINEAR_INCREMENTS, MAX_NONLINEAR_ITERATIONS, NonlinearStatic, + OutputVariable, Procedure, RigidWall, Step, StressRate, Tolerances, Transient, }; pub use units::UnitSystem; diff --git a/crates/model/src/materials.rs b/crates/model/src/materials.rs index 0571bf1..61dc819 100644 --- a/crates/model/src/materials.rs +++ b/crates/model/src/materials.rs @@ -74,6 +74,37 @@ impl Elastic { } } +/// Hyperelastic behaviour of an isotropic material at large strain, in nonlinear steps. +/// +/// The constants are the material's isotropic elastic ones: the shear modulus μ = E / (2(1 + ν)) +/// and the bulk modulus κ = E / (3(1 − 2ν)), so at small strain the material is that linear +/// elastic one (linear steps use it so). Without this option an elastic material in a nonlinear +/// step is St Venant–Kirchhoff (the second Piola–Kirchhoff stress linear in the Green–Lagrange +/// strain), which is right for large rotations with small strains. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] +pub enum Hyperelastic { + /// Compressible neo-Hookean with the isochoric–volumetric split: the strain energy per + /// reference volume W = μ/2 (Ī₁ − 3) + κ/2 (J − 1)², with Ī₁ = J^(−2/3) tr C (C₁₀ = μ/2, + /// D₁ = 2/κ in the usual polynomial notation). + NeoHookean {}, +} + +impl Hyperelastic { + /// The shear and bulk moduli μ = E / (2(1 + ν)) and κ = E / (3(1 − 2ν)) the model takes from + /// the isotropic constants. + pub fn moduli(young: f64, poisson: f64) -> (f64, f64) { + (young / (2.0 * (1.0 + poisson)), young / (3.0 * (1.0 - 2.0 * poisson))) + } + + /// The name of the model, for messages. + pub fn name(&self) -> &'static str { + match self { + Hyperelastic::NeoHookean {} => "neo-Hookean", + } + } +} + /// Fluid behaviour. #[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] #[serde(tag = "type", rename_all = "snake_case", deny_unknown_fields)] @@ -332,6 +363,9 @@ pub struct Erosion { pub struct Material { #[serde(default, skip_serializing_if = "Option::is_none")] pub elastic: Option, + /// Large-strain behaviour in nonlinear steps, built on the isotropic elastic constants. + #[serde(default, skip_serializing_if = "Option::is_none")] + pub hyperelastic: Option, /// Mass per volume, needed by frequency steps, gravity and transient heat transfer. #[serde(default, skip_serializing_if = "Option::is_none")] pub density: Option, @@ -366,6 +400,14 @@ impl Material { if let Some(e) = &self.elastic { e.validate()?; } + if let Some(h) = &self.hyperelastic { + if !matches!(self.elastic, Some(Elastic::Isotropic { .. })) { + return Err(Error::Invalid(format!("{} hyperelasticity needs isotropic elastic constants", h.name()))); + } + if self.plastic.is_some() { + return Err(Error::Invalid(format!("plasticity of a {} hyperelastic material is not supported yet", h.name()))); + } + } if let Some(d) = self.density { positive("density", d)?; } @@ -412,6 +454,11 @@ impl Model { self.materials.get(name) } + /// The hyperelastic materials some section uses, by name. + pub fn hyperelastic_in_use(&self) -> std::collections::BTreeSet<&str> { + self.sections.values().map(|s| s.material.as_str()).filter(|m| self.material(m).is_some_and(|m| m.hyperelastic.is_some())).collect() + } + pub fn add_material(&mut self, name: &str, material: Material) -> Result<()> { material.validate().map_err(|e| context("material", name, e))?; insert_new(&mut self.materials, "material", name, material) @@ -544,6 +591,25 @@ mod tests { assert_eq!(serde_json::from_str::(&text).unwrap(), Material { plastic: Some(jc()), damage: Some(damage), ..steel() }); } + #[test] + fn hyperelasticity_builds_on_isotropic_constants() { + let m: Material = serde_json::from_value(serde_json::json!({ + "elastic": {"type": "isotropic", "young": 3.0, "poisson": 0.49}, "hyperelastic": {"type": "neo_hookean"} + })) + .unwrap(); + m.validate().unwrap(); + assert_eq!(m.hyperelastic, Some(Hyperelastic::NeoHookean {})); + let back: Material = serde_json::from_value(serde_json::to_value(&m).unwrap()).unwrap(); + assert_eq!(back, m); + assert!(serde_json::from_value::(serde_json::json!({"hyperelastic": {"type": "neo_hookean", "c10": 1}})).is_err()); + let alone = Material { hyperelastic: Some(Hyperelastic::NeoHookean {}), density: Some(1.0), ..Material::default() }; + assert!(alone.validate().is_err(), "needs isotropic constants"); + let ortho = Elastic::EngineeringConstants { young: [1.0; 3], poisson: [0.0; 3], shear: [1.0; 3] }; + assert!(Material { elastic: Some(ortho), ..alone }.validate().is_err()); + let plastic = Plastic::Bilinear { yield_stress: 1.0, hardening_modulus: 0.0, hardening: Hardening::Isotropic }; + assert!(Material { plastic: Some(plastic), ..m }.validate().is_err()); + } + #[test] fn thermal_properties() { let m: Material = serde_json::from_value(serde_json::json!({ diff --git a/crates/model/src/steps.rs b/crates/model/src/steps.rs index 198d70e..ec52171 100644 --- a/crates/model/src/steps.rs +++ b/crates/model/src/steps.rs @@ -86,6 +86,114 @@ pub enum Procedure { /// Steady heat conduction and the static equilibrium of the thermal strains it causes, /// together: every boundary condition and load of the step acts. CoupledTemperatureDisplacement {}, + /// Static equilibrium followed in load increments with Newton–Raphson iterations: large + /// displacements, rotations and strains (`nonlinear_geometry`), hyperelastic materials and + /// follower pressure. Loads and prescribed displacements grow in proportion to the step's + /// load factor, from 0 to 1. + NonlinearStatic(NonlinearStatic), +} + +/// Most load increments a nonlinear static step may take. +pub const MAX_NONLINEAR_INCREMENTS: u32 = 100_000; + +/// Most Newton iterations a nonlinear static step may take in one increment. +pub const MAX_NONLINEAR_ITERATIONS: u32 = 1000; + +/// The settings of a nonlinear static step. +/// +/// The step's load factor λ runs from 0 to 1: every load is λ times its value and every +/// prescribed displacement λ times its value (a linear ramp). Increments are fractions of the +/// step: the first is `initial_increment`; an increment that converges quickly lets the next +/// grow (up to `max_increment`), one that diverges, fails to converge in `max_iterations` or +/// inverts an element is retried with a quarter of its size, and the step fails, naming the +/// load factor and the reason, when an increment would fall below `min_increment` or the step +/// has taken `max_increments` increments. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(default, deny_unknown_fields)] +pub struct NonlinearStatic { + /// Large displacements, rotations and strains (total Lagrangian elements, follower + /// pressure); `false` keeps the kinematics of small displacements (only the materials are + /// nonlinear). + pub nonlinear_geometry: bool, + /// The first increment, as a fraction of the step. + pub initial_increment: f64, + /// The smallest increment a cutback may reach. + pub min_increment: f64, + /// The largest increment. + pub max_increment: f64, + /// Cap on the increments (cutbacks included). + pub max_increments: u32, + /// Cap on the Newton iterations of one increment. + pub max_iterations: u32, + /// Convergence tolerances. + pub tolerances: Tolerances, + /// A line search along each Newton correction (off by default: full Newton steps converge + /// best for elastic large-rotation problems; a search helps where the material softens). + pub line_search: bool, +} + +/// When an increment's Newton iterations have converged: all three hold. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +#[serde(default, deny_unknown_fields)] +pub struct Tolerances { + /// The largest residual force component over the largest force (applied, internal or + /// reaction) the step has seen. + pub force: f64, + /// The largest displacement correction over the largest displacement of the increment. + pub displacement: f64, + /// The work of the last correction against the residual over that of the first. + pub energy: f64, +} + +impl Default for Tolerances { + fn default() -> Tolerances { + Tolerances { force: 1e-6, displacement: 1e-5, energy: 1e-10 } + } +} + +impl Default for NonlinearStatic { + fn default() -> NonlinearStatic { + NonlinearStatic { + nonlinear_geometry: true, + initial_increment: 0.1, + min_increment: 1e-5, + max_increment: 1.0, + max_increments: 100, + max_iterations: 16, + tolerances: Tolerances::default(), + line_search: false, + } + } +} + +impl NonlinearStatic { + pub fn validate(&self) -> Result<()> { + let fraction = |what: &str, v: f64| { + if v > 0.0 && v <= 1.0 { Ok(v) } else { Err(Error::Invalid(format!("the {what} is a fraction of the step in (0, 1], got {v}"))) } + }; + let initial = fraction("initial increment", self.initial_increment)?; + let min = fraction("minimum increment", self.min_increment)?; + let max = fraction("maximum increment", self.max_increment)?; + if !(min <= initial && initial <= max) { + return Err(Error::Invalid(format!("the increments must satisfy minimum ≤ initial ≤ maximum, got {min}, {initial} and {max}"))); + } + if !(1..=MAX_NONLINEAR_INCREMENTS).contains(&self.max_increments) { + return Err(Error::Invalid(format!("a nonlinear step takes 1 to {MAX_NONLINEAR_INCREMENTS} increments, got {}", self.max_increments))); + } + if !(1..=MAX_NONLINEAR_ITERATIONS).contains(&self.max_iterations) { + return Err(Error::Invalid(format!( + "a nonlinear step takes 1 to {MAX_NONLINEAR_ITERATIONS} iterations per increment, got {}", + self.max_iterations + ))); + } + let t = self.tolerances; + for (what, v) in [("force", t.force), ("displacement", t.displacement), ("energy", t.energy)] { + if !(v > 0.0 && v < 1.0) { + return Err(Error::Invalid(format!("the {what} tolerance is in (0, 1), got {v}"))); + } + } + Ok(()) + } } /// Most time increments of a transient heat transfer step. @@ -386,7 +494,7 @@ impl BulkViscosity { } impl Procedure { - /// `static` or `frequency`. + /// The procedure's type, as the model file names it. pub fn name(&self) -> &'static str { match self { Procedure::Static {} => "static", @@ -395,6 +503,7 @@ impl Procedure { Procedure::HeatTransfer { .. } => "heat_transfer", Procedure::CoupledTemperatureDisplacement {} => "coupled_temperature_displacement", Procedure::Fluid { .. } => "fluid", + Procedure::NonlinearStatic(_) => "nonlinear_static", } } @@ -435,6 +544,9 @@ impl Procedure { if let Procedure::ExplicitDynamic(e) = self { return e.validate(); } + if let Procedure::NonlinearStatic(n) = self { + return n.validate(); + } if let Procedure::Fluid { end_time, steps, time_step, frames, magic, .. } = *self { match (end_time, steps) { (Some(t), None) => { @@ -480,7 +592,7 @@ impl Procedure { /// The output a step of this procedure requests unless told otherwise. pub fn default_output(&self) -> BTreeSet { match self { - Procedure::Static {} => OutputVariable::ALL.into_iter().collect(), + Procedure::Static {} | Procedure::NonlinearStatic(_) => OutputVariable::ALL.into_iter().collect(), Procedure::Frequency { .. } => [OutputVariable::Displacement].into_iter().collect(), Procedure::ExplicitDynamic(_) => [OutputVariable::Displacement, OutputVariable::Stress].into_iter().collect(), // Velocity, pressure and wall forces, always. @@ -753,6 +865,43 @@ mod tests { assert!(serde_json::from_value::(serde_json::json!({"type": "explicit_dynamic", "end_time": 1, "dt": 2})).is_err()); } + #[test] + fn nonlinear_static_steps_validate_their_controls() { + let p: Procedure = serde_json::from_value(serde_json::json!({"type": "nonlinear_static"})).unwrap(); + assert_eq!(p, Procedure::NonlinearStatic(NonlinearStatic::default())); + assert_eq!(p.name(), "nonlinear_static"); + assert!(p.validate().is_ok() && p.is_structural() && !p.is_dynamic()); + assert_eq!(p.default_output(), Procedure::Static {}.default_output()); + let p: Procedure = serde_json::from_value(serde_json::json!({ + "type": "nonlinear_static", "nonlinear_geometry": false, "initial_increment": 0.5, "min_increment": 0.01, + "max_increments": 20, "max_iterations": 8, "tolerances": {"force": 1e-4}, "line_search": true + })) + .unwrap(); + let Procedure::NonlinearStatic(n) = &p else { panic!("{p:?}") }; + assert_eq!((n.nonlinear_geometry, n.max_increment, n.tolerances.force, n.tolerances.energy), (false, 1.0, 1e-4, 1e-10)); + let back: Procedure = serde_json::from_value(serde_json::to_value(&p).unwrap()).unwrap(); + assert_eq!(back, p); + let base = NonlinearStatic::default(); + let bad = [ + NonlinearStatic { initial_increment: 0.0, ..base.clone() }, + NonlinearStatic { initial_increment: f64::NAN, ..base.clone() }, + NonlinearStatic { initial_increment: 2.0, max_increment: 2.0, ..base.clone() }, + NonlinearStatic { min_increment: 0.5, initial_increment: 0.1, ..base.clone() }, + NonlinearStatic { max_increment: 0.05, ..base.clone() }, + NonlinearStatic { max_increments: 0, ..base.clone() }, + NonlinearStatic { max_increments: MAX_NONLINEAR_INCREMENTS + 1, ..base.clone() }, + NonlinearStatic { max_iterations: 0, ..base.clone() }, + NonlinearStatic { max_iterations: u32::MAX, ..base.clone() }, + NonlinearStatic { tolerances: Tolerances { force: 0.0, ..Tolerances::default() }, ..base.clone() }, + NonlinearStatic { tolerances: Tolerances { energy: f64::INFINITY, ..Tolerances::default() }, ..base.clone() }, + NonlinearStatic { tolerances: Tolerances { displacement: 1.0, ..Tolerances::default() }, ..base.clone() }, + ]; + for b in bad { + assert!(Procedure::NonlinearStatic(b.clone()).validate().is_err(), "{b:?}"); + } + assert!(serde_json::from_value::(serde_json::json!({"type": "nonlinear_static", "nlgeom": true})).is_err()); + } + #[test] fn heat_transfer_procedures() { let steady: Procedure = serde_json::from_value(serde_json::json!({"type": "heat_transfer"})).unwrap();