solver, model: imperfections from buckling modes in nonlinear steps

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 <[email protected]>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

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");

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");

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;

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<Vec<[f64; 2]>>,
+    /// 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<Imperfection>,
+}
+
+/// 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)?;

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<StaticResult> {
-    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<String>) -> Result<StaticResult> {
+    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<String> = resolved.warnings.iter().map(ToString::to_string).chain(mesh.warnings.iter().cloned()).collect();
+    let mut warnings: Vec<String> =
+        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);

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());
+}

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