From 0bb1c2462ff4ddba3e6e11df9f37ef7fe81d94b7 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Wed, 14 Oct 2026 13:17:31 -0400 Subject: [PATCH] solver, model: imperfections from buckling modes in nonlinear steps MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit NonlinearStatic.imperfection {step, mode, scale}: the buckling step's mode, scaled to a largest component of `scale`, is added to the node coordinates before the step is solved, with a warning naming it; Abaqus and LS-DYNA decks refuse it by name. A clamped–free column with its first mode as imperfection under half its buckling load amplifies it by 1/(1 − α) within 2.8 %. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/decks/src/abaqus.rs | 9 +++++++ crates/decks/src/lsdyna/mod.rs | 6 +++++ crates/model/src/lib.rs | 2 +- crates/model/src/steps.rs | 21 ++++++++++++++++ crates/solver/src/nonlinear.rs | 46 ++++++++++++++++++++++++++++++---- crates/solver/tests/buckle.rs | 34 +++++++++++++++++++++++++ docs/fem.md | 7 ++++++ 7 files changed, 119 insertions(+), 6 deletions(-) diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index 79a28e3..2c19212 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -1191,6 +1191,15 @@ impl Deck { Procedure::Fluid { .. } => return Err(crate::fluid_step(FORMAT, &step.name)), Procedure::Static {} | Procedure::NonlinearStatic(_) => { if let Procedure::NonlinearStatic(n) = &step.procedure { + if n.imperfection.is_some() { + return Err(Error::Unsupported { + format: FORMAT, + what: format!( + "step `{}`: an imperfection from a buckling mode (FieldCraft computes the mode when it solves)", + step.name + ), + }); + } 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"); diff --git a/crates/decks/src/lsdyna/mod.rs b/crates/decks/src/lsdyna/mod.rs index 43bb79a..8d86150 100644 --- a/crates/decks/src/lsdyna/mod.rs +++ b/crates/decks/src/lsdyna/mod.rs @@ -874,6 +874,12 @@ fn write_step(d: &mut Deck, model: &Model, r: &Resolved, plan: &Plan, step: &Ste d.card([Field::from(1)]); } Procedure::NonlinearStatic(n) => { + if n.imperfection.is_some() { + return Err(Error::Unsupported { + format: FORMAT, + what: format!("step `{}`: an imperfection from a buckling mode (FieldCraft computes the mode when it solves)", step.name), + }); + } // 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"); diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs index cbee8d5..361e868 100644 --- a/crates/model/src/lib.rs +++ b/crates/model/src/lib.rs @@ -67,7 +67,7 @@ pub use mesh::{ pub use sections::{BeamProfile, Extent, Formulation, Integration, Section, SectionKind}; pub use sets::Set; pub use steps::{ - BulkViscosity, ExplicitDynamic, Hourglass, HourglassKind, InitialVelocity, MAX_LOAD_PATH_POINTS, MAX_NONLINEAR_INCREMENTS, + BulkViscosity, ExplicitDynamic, Hourglass, HourglassKind, Imperfection, InitialVelocity, MAX_LOAD_PATH_POINTS, MAX_NONLINEAR_INCREMENTS, MAX_NONLINEAR_ITERATIONS, NonlinearStatic, OutputVariable, Procedure, RigidWall, Step, StressRate, Tolerances, Transient, }; pub use units::UnitSystem; diff --git a/crates/model/src/steps.rs b/crates/model/src/steps.rs index 73ed648..2fb30bd 100644 --- a/crates/model/src/steps.rs +++ b/crates/model/src/steps.rs @@ -145,6 +145,20 @@ pub struct NonlinearStatic { /// them (cyclic loading). Increments end at every point of the path. #[serde(skip_serializing_if = "Option::is_none")] pub load_path: Option>, + /// A geometric imperfection: the mode shape of a buckling step, scaled, added to the node + /// coordinates before this step is solved (the step and its results are on that geometry). + #[serde(skip_serializing_if = "Option::is_none")] + pub imperfection: Option, +} + +/// Mode `mode` (from 1) of buckling step `step`, scaled so its largest component is `scale` +/// (a length), as an initial imperfection. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Imperfection { + pub step: String, + pub mode: u32, + pub scale: f64, } /// Most points of a nonlinear step's load path. @@ -181,6 +195,7 @@ impl Default for NonlinearStatic { tolerances: Tolerances::default(), line_search: false, load_path: None, + imperfection: None, } } } @@ -190,6 +205,12 @@ impl NonlinearStatic { 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}"))) } }; + if let Some(imp) = &self.imperfection { + if !(1..=MAX_MODES).contains(&imp.mode) { + return Err(Error::Invalid(format!("an imperfection takes mode 1 to {MAX_MODES}, got {}", imp.mode))); + } + finite("imperfection scale", imp.scale)?; + } let initial = fraction("initial increment", self.initial_increment)?; let min = fraction("minimum increment", self.min_increment)?; let max = fraction("maximum increment", self.max_increment)?; diff --git a/crates/solver/src/nonlinear.rs b/crates/solver/src/nonlinear.rs index e9e452b..4651d3f 100644 --- a/crates/solver/src/nonlinear.rs +++ b/crates/solver/src/nonlinear.rs @@ -98,7 +98,7 @@ use fieldcraft_fem::finite::{Kinematics, Temperature}; use fieldcraft_fem::plane::Theory; use fieldcraft_fem::quadrature::Rule; use fieldcraft_fem::{DofKind, DofSet, ElementType, PointStress, face, von_mises}; -use fieldcraft_model::{BoundaryCondition, Extent, Hardening, Load, Model, NonlinearStatic, Plastic, Procedure, SectionKind, Step}; +use fieldcraft_model::{BoundaryCondition, Extent, Hardening, Imperfection, Load, Model, NonlinearStatic, Plastic, Procedure, SectionKind, Step}; use fieldcraft_parallel as parallel; use serde::Serialize; @@ -837,14 +837,50 @@ fn followers(model: &Model, resolved: &fieldcraft_model::Resolved, mesh: &Mesh, /// Solve a nonlinear static step: the step named `step`, or the model's only nonlinear static /// step. pub fn solve_nonlinear(model: &Model, step: Option<&str>, limits: &Limits) -> Result { - let resolved = model.resolve()?; - let step = model.select_step(step, |p| matches!(p, Procedure::NonlinearStatic(_)), "nonlinear static")?; - let Procedure::NonlinearStatic(settings) = &step.procedure else { + let chosen = model.select_step(step, |p| matches!(p, Procedure::NonlinearStatic(_)), "nonlinear static")?; + let Procedure::NonlinearStatic(settings) = &chosen.procedure else { return Err(Error::Internal("nonlinear static step".into())); }; settings.validate()?; + let Some(imperfection) = &settings.imperfection else { return solve_on(model, chosen.name.as_str(), limits, Vec::new()) }; + let (perturbed, note) = imperfect(model, imperfection, limits)?; + solve_on(&perturbed, chosen.name.as_str(), limits, vec![note]) +} + +/// The model with a buckling mode, scaled to a largest component of `scale`, added to its node +/// coordinates (the buckling step solved first), and a note saying so. +fn imperfect(model: &Model, imperfection: &Imperfection, limits: &Limits) -> Result<(Model, String)> { + let Imperfection { step, mode, scale } = imperfection; + if !matches!(model.step(step).map(|s| &s.procedure), Some(Procedure::Buckle { .. })) { + return Err(Error::Invalid(format!("the imperfection names step `{step}`, which is not a buckling step"))); + } + let buckled = crate::buckle::solve_buckle(model, Some(step), limits)?; + let shape = + usize::try_from(*mode).ok().and_then(|m| m.checked_sub(1)).and_then(|m| buckled.modes.get(m)).ok_or_else(|| { + Error::Invalid(format!("the imperfection takes mode {mode} of step `{step}`, which found {} modes", buckled.modes.len())) + })?; + let mut out = model.clone(); + for (label, u) in buckled.nodes.iter().zip(&shape.shape) { + let node = model.node(*label).ok_or_else(|| Error::Internal(format!("node {label}")))?; + let xyz: [f64; 3] = std::array::from_fn(|k| node.xyz[k] + scale * u[k]); + out.set_node(*label, xyz)?; + } + let note = format!( + "solved on an imperfect geometry: mode {mode} of buckling step `{step}` (load factor {:.6e}) scaled to {scale} added to the node coordinates", + shape.eigenvalue + ); + Ok((out, note)) +} + +fn solve_on(model: &Model, step: &str, limits: &Limits, notes: Vec) -> Result { + let resolved = model.resolve()?; + let step = model.step(step).ok_or_else(|| Error::Internal("nonlinear static step".into()))?; + let Procedure::NonlinearStatic(settings) = &step.procedure else { + return Err(Error::Internal("nonlinear static step".into())); + }; 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(); + let mut warnings: Vec = + notes.into_iter().chain(resolved.warnings.iter().map(ToString::to_string)).chain(mesh.warnings.iter().cloned()).collect(); check_supported(model, &mesh, step, settings, &mut warnings)?; let (graph, order) = pattern(&mesh); ignored_dashpots(&mesh, &mut warnings); diff --git a/crates/solver/tests/buckle.rs b/crates/solver/tests/buckle.rs index fe94921..1b7fb27 100644 --- a/crates/solver/tests/buckle.rs +++ b/crates/solver/tests/buckle.rs @@ -91,3 +91,37 @@ fn buckling_steps_refuse_what_they_cannot_do() { assert!(Procedure::Buckle { modes: 0, shift: None }.validate().is_err()); assert!(Procedure::Buckle { modes: 1, shift: Some(f64::NAN) }.validate().is_err()); } + +/// A clamped–free column given its first buckling mode as an imperfection δ₀ at the top, under a +/// dead end force of half its buckling load: the top's lateral deflection grows to +/// δ₀ / (1 − P/P_cr) = 2 δ₀ (Timoshenko & Gere 1961, §1.12, for the mode-shaped imperfection). +#[test] +fn an_imperfection_from_the_buckling_mode() { + let mut g = column(L / 2.0, 10, None); + g.fixed("base", |x| near(x[2], 0.0)); + let critical = solve_buckle(&g.model, Some("Step"), &Limits::default()).unwrap().modes[0].eigenvalue * B * B; + let delta0 = 0.01; + g.node_set("tip", |x| near(x[2], 1.0)); + g.model + .add_load("push", Load::Force { target: Target::Set("tip".into()), vector: [0.0, 0.0, -0.5 * critical], distribution: Distribution::Total }) + .unwrap(); + let settings = NonlinearStatic { imperfection: Some(Imperfection { step: "Step".into(), mode: 1, scale: delta0 }), ..NonlinearStatic::default() }; + let mut nl = Step::new("Push", Procedure::NonlinearStatic(settings)); + nl.boundary_conditions.insert("base".into()); + nl.loads.insert("push".into()); + g.model.add_step(nl, None).unwrap(); + let r = fieldcraft_solver::solve_nonlinear(&g.model, Some("Push"), &Limits::default()).unwrap(); + assert!(r.warnings.iter().any(|w| w.contains("imperfect geometry")), "{:?}", r.warnings); + let lateral = r.displacement.iter().map(|u| u[0].hypot(u[1])).fold(0.0, f64::max); + // The deflection added to δ₀: δ₀ α/(1 − α) = δ₀ at α = ½. + eprintln!("imperfection {delta0}: added top deflection {lateral:.6e}"); + assert!((lateral / delta0 - 1.0).abs() < 0.05, "{lateral}"); + // A mode the buckling step does not have, or a step that does not buckle. + let mut bad = g.model.clone(); + let mut s = bad.step("Push").unwrap().clone(); + if let Procedure::NonlinearStatic(n) = &mut s.procedure { + n.imperfection = Some(Imperfection { step: "Step".into(), mode: 9, scale: delta0 }); + } + bad.replace_step(s).unwrap(); + assert!(fieldcraft_solver::solve_nonlinear(&bad, Some("Push"), &Limits::default()).is_err()); +} diff --git a/docs/fem.md b/docs/fem.md index f6d09d8..882720a 100644 --- a/docs/fem.md +++ b/docs/fem.md @@ -193,6 +193,13 @@ buckle from, and with a shift); Nastran SOL 105 and LS-DYNA refuse it by name. V columns of hex20 (pinned and clamped–free, 0.18 % from Engesser's shear-corrected load) and the shifted form against the plain one. +**Imperfections.** A nonlinear static step's `"imperfection": {"step": "Buckle", "mode": 1, +"scale": 0.01}` solves that buckling step first and adds its mode, scaled to a largest component +of `scale`, to the node coordinates; the step is solved and reported on that geometry, with a +warning naming the mode. Abaqus and LS-DYNA decks refuse it by name (their imperfection reads a +results file). Verified: a clamped–free column with its first mode as imperfection under half its +buckling load deflects by δ₀ α/(1 − α) (2.8 %). + ## Nonlinear static analysis (finite strain) A step `{"type": "nonlinear_static", …}` (`Procedure::NonlinearStatic`, solved by