forked from 1b-agent/fieldcraft-upstream
jyvikoandClaude Opus 5.5
0bb1c2462f
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
128 lines
6.7 KiB
Rust
128 lines
6.7 KiB
Rust
//! Linear buckling against closed forms: Euler columns of hex20 solids (pinned–pinned by a
|
||
//! symmetry plane at midspan, and clamped–free), with Engesser's shear correction
|
||
//! P = P_e / (1 + P_e/(κ G A)) (S. P. Timoshenko, J. M. Gere, *Theory of Elastic Stability*,
|
||
//! 2nd ed., McGraw-Hill 1961, §2.1, §2.17), the shifted (LDLᵀ) form against the plain one, and
|
||
//! what buckling steps refuse.
|
||
|
||
mod common;
|
||
|
||
use common::*;
|
||
use fieldcraft_model::*;
|
||
use fieldcraft_solver::{Limits, solve_buckle};
|
||
|
||
const E: f64 = 1000.0;
|
||
const B: f64 = 1.0;
|
||
const L: f64 = 20.0;
|
||
|
||
fn buckle(modes: u32, shift: Option<f64>) -> Procedure {
|
||
Procedure::Buckle { modes, shift }
|
||
}
|
||
|
||
/// Engesser's critical load of a square column of side B and effective length `le`, ν = 0.
|
||
fn engesser(le: f64) -> f64 {
|
||
let (a, i, g) = (B * B, B.powi(4) / 12.0, E / 2.0);
|
||
let pe = std::f64::consts::PI.powi(2) * E * i / (le * le);
|
||
pe / (1.0 + pe / (5.0 / 6.0 * g * a))
|
||
}
|
||
|
||
/// A column of `length` along z, 2 × 2 × `cells` hex20, pressed by a unit pressure on its top.
|
||
fn column(length: f64, cells: usize, shift: Option<f64>) -> Grid {
|
||
let mut g = grid([2, 2, cells], ElementKind::Hex20, |x| [B * x[0], B * x[1], length * x[2]]);
|
||
g.define(iso(E, 0.0), None, Formulation::default(), buckle(3, shift));
|
||
g.face_set("top", |x| near(x[2], 1.0));
|
||
g.load("press", Load::Pressure { set: "top".into(), magnitude: 1.0 });
|
||
g
|
||
}
|
||
|
||
#[test]
|
||
fn euler_columns() {
|
||
// Pinned–pinned: half the column, the pinned end (z = 0 here is the loaded top) held
|
||
// laterally over its face (which leaves its rotation free), midspan a symmetry plane.
|
||
let mut g = grid([2, 2, 10], ElementKind::Hex20, |x| [B * x[0], B * x[1], L / 2.0 * x[2]]);
|
||
g.define(iso(E, 0.0), None, Formulation::default(), buckle(2, None));
|
||
g.face_set("end", |x| near(x[2], 0.0));
|
||
g.load("press", Load::Pressure { set: "end".into(), magnitude: 1.0 });
|
||
g.node_set("pin", |x| near(x[2], 0.0));
|
||
g.bc("pin", BoundaryCondition::Displacement { target: Target::Set("pin".into()), values: [Some(0.0), Some(0.0), None, None, None, None] });
|
||
g.symmetry("mid", Symmetry::Z, |x| near(x[2], 1.0));
|
||
let r = solve_buckle(&g.model, Some("Step"), &Limits::default()).unwrap();
|
||
let want = engesser(L) / (B * B);
|
||
let got = r.modes[0].eigenvalue;
|
||
eprintln!("pinned: λ = {got:.6} (Engesser {want:.6}, {:+.3} %), second {:.6}", 100.0 * (got / want - 1.0), r.modes[1].eigenvalue);
|
||
assert!((got / want - 1.0).abs() < 5e-3);
|
||
assert!((r.modes[1].eigenvalue / got - 1.0).abs() < 1e-6, "the square section's two planes");
|
||
assert!(r.buckling.as_ref().unwrap().negative.is_empty() || r.buckling.as_ref().unwrap().negative[0] < -got);
|
||
|
||
// Clamped–free: the base held, the top free: le = 2L.
|
||
let mut g = column(L / 2.0, 10, None);
|
||
g.fixed("base", |x| near(x[2], 0.0));
|
||
let r = solve_buckle(&g.model, Some("Step"), &Limits::default()).unwrap();
|
||
let want = engesser(L) / (B * B);
|
||
let got = r.modes[0].eigenvalue;
|
||
eprintln!("clamped–free: λ = {got:.6} (Engesser {want:.6}, {:+.3} %)", 100.0 * (got / want - 1.0));
|
||
assert!((got / want - 1.0).abs() < 5e-3);
|
||
// Scaled to a largest component of 1, which is lateral at the free top.
|
||
let big = r.modes[0].shape.iter().flatten().fold(0.0f64, |m, v| m.max(v.abs()));
|
||
assert!((big - 1.0).abs() < 1e-12);
|
||
|
||
// Shifted between the first and the higher modes (LDLᵀ of K + σ K_σ): its inertia counts
|
||
// the two load factors below the shift, and the modes above it come out.
|
||
let mut s = column(L / 2.0, 10, Some(2.0 * got));
|
||
s.fixed("base", |x| near(x[2], 0.0));
|
||
let shifted = solve_buckle(&s.model, Some("Step"), &Limits::default()).unwrap();
|
||
let info = shifted.buckling.as_ref().unwrap();
|
||
assert_eq!(info.below_shift, Some(2));
|
||
let third = r.modes[2].eigenvalue;
|
||
assert!((shifted.modes[0].eigenvalue / third - 1.0).abs() < 1e-7, "{} vs {third}", shifted.modes[0].eigenvalue);
|
||
}
|
||
|
||
#[test]
|
||
fn buckling_steps_refuse_what_they_cannot_do() {
|
||
// No load: nothing to buckle.
|
||
let mut g = column(L, 4, None);
|
||
g.model.set_load_active("Step", "press", false).unwrap();
|
||
g.fixed("base", |x| near(x[2], 0.0));
|
||
assert!(solve_buckle(&g.model, Some("Step"), &Limits::default()).is_err());
|
||
// Unsupported: a located singular stiffness.
|
||
let g = column(L, 4, None);
|
||
let e = solve_buckle(&g.model, Some("Step"), &Limits::default()).unwrap_err();
|
||
assert!(e.to_string().contains("singular"), "{e}");
|
||
// A shift that is not finite, and a step of another kind.
|
||
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());
|
||
}
|