forked from 1b-agent/fieldcraft-upstream
jyvikoandClaude Opus 5.5
a30d0b5c5b
Segment-to-segment penalty contact between element faces: each secondary face and each main face facing it are projected onto the secondary face's plane, split into triangles, clipped (Sutherland–Hodgman) and integrated with a three-point rule per sub-triangle (Puso–Laursen segmentation, the penalty at the integration points as Fischer–Wriggers), so a uniform pressure gives the consistent nodal forces on both meshes and the contact patch test passes on non-matching meshes. Every contact runs both ways at half the penalty: pairs with each surface as the secondary, general and self contact face against face. A bucket grid of face boxes every increment, faces integrated in parallel and their forces added in face order: the same numbers on any thread count. The penalty per area ε = s γ min(m/A)/Δt² comes from the lumped masses, each node's directional area and the stable increment, bounded by Gershgorin's row sums so contact never shortens the increment. Coulomb friction with exponential decay by a regularised return map at the secondary nodes; rigid walls take Coulomb impulses. Contact work and friction dissipation join the energy balance; each contact reports its peak force, first and last contact, deepest penetration, friction and the nodes adjusted at the start (initial penetrations adjusted strain-free or refused with the node, face and depth); a cap on contact faces. Verified (crates/explicit/tests/contact.rs, tests/oracle/explicit-*): colliding bars (2L/c +1.2%, ρcAv₀/2, KE kept to 0.005%, contact energy back to 8e-5 of it), blocks sliding to a stop on a rigid wall and on a plate under constant and decaying friction (0.005–0.2%), the patch test (σ_zz uniform to 0.06%), two spheres colliding against Hertz (peak +2.3%, duration −0.4%; the oracle case waits for a graded sphere mesher), a free body untouched by general contact, determinism across threads. Oracle measures gain `over` (ratios). Existing fixtures do not move. Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm
483 lines
24 KiB
Rust
483 lines
24 KiB
Rust
//! Large deformation and plasticity in explicit dynamics against closed-form solutions
|
||
//! (docs/solvers/literature.md, A4, A6, A10):
|
||
//! - simple shear of a hypoelastic material (J. K. Dienes, *On the analysis of rotation and
|
||
//! stress rate in deforming bodies*, Acta Mech. 32:217–232, 1979): the Jaumann rate gives the
|
||
//! oscillating σ12 = G sin γ, σ11 = −σ22 = G(1 − cos γ), the Green–Naghdi rate the monotonic
|
||
//! solution in β = arctan(γ/2);
|
||
//! - a stressed element turned rigidly by 90° in many increments: the stress turns with it and
|
||
//! keeps its magnitude (objectivity), where the small-strain update does not;
|
||
//! - uniaxial stress under large stretch with J2 plasticity, loaded and unloaded: Cauchy stress
|
||
//! against the logarithmic strain follows the bilinear flow curve and unloads elastically;
|
||
//! - the Taylor bar (G. I. Taylor, Proc. R. Soc. A 194:289–299, 1948): the yield stress Taylor's
|
||
//! one-dimensional estimate recovers from the computed final shape.
|
||
//!
|
||
//! Run with `--nocapture` to see the errors.
|
||
|
||
use fieldcraft_explicit::elements::{Motion, State, Workspace, internal_forces};
|
||
use fieldcraft_explicit::system::System;
|
||
use fieldcraft_explicit::{ExplicitResult, Failure, Limits, solve_explicit};
|
||
use fieldcraft_model::*;
|
||
|
||
/// A block [0, size] of `n` cells: one hex8 per cell, or six tet4 (Kuhn's split).
|
||
fn block(n: [usize; 3], size: [f64; 3], kind: ElementKind) -> Model {
|
||
let mut m = Model::default();
|
||
let label = |i: usize, j: usize, k: usize| (1 + i + (n[0] + 1) * (j + (n[1] + 1) * k)) as u32;
|
||
for k in 0..=n[2] {
|
||
for j in 0..=n[1] {
|
||
for i in 0..=n[0] {
|
||
let x = [i as f64 * size[0] / n[0] as f64, j as f64 * size[1] / n[1] as f64, k as f64 * size[2] / n[2] as f64];
|
||
m.set_node(label(i, j, k), x).unwrap();
|
||
}
|
||
}
|
||
}
|
||
let corners = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]];
|
||
let mut next = 1;
|
||
for k in 0..n[2] {
|
||
for j in 0..n[1] {
|
||
for i in 0..n[0] {
|
||
let at = |c: [usize; 3]| label(i + c[0], j + c[1], k + c[2]);
|
||
if kind == ElementKind::Hex8 {
|
||
m.set_element(next, kind, corners.iter().map(|c| at(*c)).collect()).unwrap();
|
||
next += 1;
|
||
continue;
|
||
}
|
||
for perm in [[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]] {
|
||
let mut v = [[0usize; 3]; 4];
|
||
for s in 0..3 {
|
||
v[s + 1] = v[s];
|
||
v[s + 1][perm[s]] = 1;
|
||
}
|
||
let p = v.map(|c| [0, 1, 2].map(|d| c[d] as f64 * size[d] / n[d] as f64));
|
||
let e = |a: usize| [0, 1, 2].map(|d| p[a][d] - p[0][d]);
|
||
let (a, b, c) = (e(1), e(2), e(3));
|
||
let vol = a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0]);
|
||
if vol < 0.0 {
|
||
v.swap(1, 2);
|
||
}
|
||
m.set_element(next, kind, v.iter().map(|c| at(*c)).collect()).unwrap();
|
||
next += 1;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
let all: Vec<u32> = m.elements().keys().copied().collect();
|
||
m.add_set("all", Set { elements: all, ..Set::default() }).unwrap();
|
||
m
|
||
}
|
||
|
||
fn node_set(m: &mut Model, name: &str, f: impl Fn([f64; 3]) -> bool) {
|
||
let nodes: Vec<u32> = m.nodes().iter().filter(|(_, n)| f(n.xyz)).map(|(l, _)| *l).collect();
|
||
m.add_set(name, Set { nodes, ..Set::default() }).unwrap();
|
||
}
|
||
|
||
fn define(m: &mut Model, kind: ElementKind, material: Material, step: ExplicitDynamic) {
|
||
m.add_material("Mat", material).unwrap();
|
||
let integration = if kind == ElementKind::Hex8 { Integration::Reduced } else { Integration::Full };
|
||
let formulation = Formulation { integration, ..Formulation::default() };
|
||
m.add_section("solid", Section { set: "all".into(), material: "Mat".into(), kind: SectionKind::Solid {}, formulation }).unwrap();
|
||
m.add_step(Step::new("Dynamic", Procedure::ExplicitDynamic(step)), None).unwrap();
|
||
}
|
||
|
||
fn bc(m: &mut Model, name: &str, bc: BoundaryCondition) {
|
||
m.add_boundary_condition(name, bc).unwrap();
|
||
m.set_boundary_condition_active("Dynamic", name, true).unwrap();
|
||
}
|
||
|
||
fn elastic(young: f64, poisson: f64, density: f64) -> Material {
|
||
Material { elastic: Some(Elastic::Isotropic { young, poisson }), density: Some(density), ..Material::default() }
|
||
}
|
||
|
||
/// Dienes' closed forms: (σ11, σ12) / G at shear γ.
|
||
fn dienes(rate: StressRate, gamma: f64) -> (f64, f64) {
|
||
match rate {
|
||
StressRate::Jaumann => (1.0 - gamma.cos(), gamma.sin()),
|
||
StressRate::GreenNaghdi => {
|
||
let beta = (gamma / 2.0).atan();
|
||
let lc = beta.cos().ln();
|
||
let s11 = 4.0 * ((2.0 * beta).cos() * lc + beta * (2.0 * beta).sin() - beta.sin().powi(2));
|
||
let s12 = 2.0 * (2.0 * beta).cos() * (2.0 * beta - 2.0 * (2.0 * beta).tan() * lc - beta.tan());
|
||
(s11, s12)
|
||
}
|
||
}
|
||
}
|
||
|
||
/// A unit cube sheared to γ = 6 by prescribing every node: x = X + γ(t) Y e_x.
|
||
#[test]
|
||
fn dienes_simple_shear_with_both_rates() {
|
||
let (young, poisson) = (1.0, 0.3);
|
||
let g = young / (2.0 * (1.0 + poisson));
|
||
let gamma_end = 6.0;
|
||
for kind in [ElementKind::Hex8, ElementKind::Tet4] {
|
||
for rate in [StressRate::Jaumann, StressRate::GreenNaghdi] {
|
||
let mut m = block([1, 1, 1], [1.0; 3], kind);
|
||
let end = 4000.0;
|
||
let mut step = ExplicitDynamic::new(end);
|
||
step.stress_rate = rate;
|
||
step.output_interval = Some(end / 60.0);
|
||
step.max_increments = 100_000;
|
||
define(&mut m, kind, elastic(young, poisson, 1.0), step);
|
||
node_set(&mut m, "bottom", |x| x[1] == 0.0);
|
||
node_set(&mut m, "top", |x| x[1] == 1.0);
|
||
bc(&mut m, "hold", BoundaryCondition::Fixed { target: Target::Set("bottom".into()) });
|
||
bc(
|
||
&mut m,
|
||
"shear",
|
||
BoundaryCondition::Displacement {
|
||
target: Target::Set("top".into()),
|
||
values: [Some(gamma_end), Some(0.0), Some(0.0), None, None, None],
|
||
},
|
||
);
|
||
let r = solve_explicit(&m, None, &Limits::default()).unwrap();
|
||
assert!(r.completed);
|
||
let mut worst: f64 = 0.0;
|
||
for f in &r.frames {
|
||
let gamma = gamma_end * f.time / end;
|
||
let (s11, s12) = dienes(rate, gamma);
|
||
for s in &f.stress {
|
||
let e = [(s[0] / g - s11).abs(), (s[1] / g + s11).abs(), (s[3] / g - s12).abs(), s[2].abs() / g, s[4].abs() / g, s[5].abs() / g];
|
||
worst = worst.max(e.into_iter().fold(0.0, f64::max));
|
||
}
|
||
}
|
||
let last = r.frames.last().unwrap();
|
||
eprintln!(
|
||
"Dienes {} {rate:?}: σ/G at γ = 6: ({:.5}, {:.5}) vs ({:.5}, {:.5}); largest error over the run {worst:.2e} G; {} increments",
|
||
kind.name(),
|
||
last.stress[0][0] / g,
|
||
last.stress[0][3] / g,
|
||
dienes(rate, gamma_end).0,
|
||
dienes(rate, gamma_end).1,
|
||
r.increments
|
||
);
|
||
assert!(worst < 2e-3, "{} {rate:?}: {worst}", kind.name());
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Drive one element through a sequence of nodal displacements with the solver's own element
|
||
/// update, returning its stress after each.
|
||
struct Driver {
|
||
sys: System,
|
||
state: State,
|
||
ws: Workspace,
|
||
u: Vec<f64>,
|
||
time: f64,
|
||
}
|
||
|
||
impl Driver {
|
||
fn new(model: &Model) -> Driver {
|
||
let resolved = model.resolve().unwrap();
|
||
let step = model.step("Dynamic").unwrap();
|
||
let Procedure::ExplicitDynamic(settings) = &step.procedure else { panic!() };
|
||
let sys = System::new(model, &resolved, step, settings, &Limits { max_dofs: 1000, ..Limits::default() }).unwrap();
|
||
let (state, ws) = (State::new(&sys), Workspace::new(&sys));
|
||
let u = vec![0.0; sys.dofs()];
|
||
Driver { sys, state, ws, u, time: 0.0 }
|
||
}
|
||
|
||
/// Move the nodes to `x(X)` over an increment of length `dt`.
|
||
fn go(&mut self, to: impl Fn([f64; 3]) -> [f64; 3], dt: f64) -> [f64; 6] {
|
||
self.try_go(to, dt).unwrap()
|
||
}
|
||
|
||
fn try_go(&mut self, to: impl Fn([f64; 3]) -> [f64; 3], dt: f64) -> fieldcraft_explicit::Result<[f64; 6]> {
|
||
let mut du = vec![0.0; self.u.len()];
|
||
for (i, x0) in self.sys.xyz.iter().enumerate() {
|
||
let x = to(*x0);
|
||
for k in 0..3 {
|
||
du[3 * i + k] = x[k] - x0[k] - self.u[3 * i + k];
|
||
self.u[3 * i + k] = x[k] - x0[k];
|
||
}
|
||
}
|
||
self.time += dt;
|
||
let v: Vec<f64> = du.iter().map(|d| d / dt).collect();
|
||
let bulk = BulkViscosity { linear: 0.0, quadratic: 0.0 };
|
||
let m = Motion { u: &self.u, du: &du, v: &v, dt, time: self.time, bulk: &bulk };
|
||
internal_forces(&self.sys, &m, &mut self.state, &mut self.ws)?;
|
||
Ok(self.state.elements[0].point.stress)
|
||
}
|
||
}
|
||
|
||
fn rotation_z(angle: f64, x: [f64; 3]) -> [f64; 3] {
|
||
let (s, c) = angle.sin_cos();
|
||
[c * x[0] - s * x[1], s * x[0] + c * x[1], x[2]]
|
||
}
|
||
|
||
/// Stretch an element 1% along x, then turn it rigidly by 90° about z in 50 increments: the
|
||
/// stress turns with it (σ_yy takes σ_xx's value), exactly for both rates, and the internal
|
||
/// forces stay those of the turned stress. The small-strain update keeps σ_xx and adds a
|
||
/// spurious stress from the rotation's "strain".
|
||
#[test]
|
||
fn rigid_rotation_turns_the_stress_with_the_body() {
|
||
for (rate, nonlinear) in [(StressRate::Jaumann, true), (StressRate::GreenNaghdi, true), (StressRate::Jaumann, false)] {
|
||
let mut m = block([1, 1, 1], [1.0; 3], ElementKind::Hex8);
|
||
let mut step = ExplicitDynamic::new(1.0);
|
||
step.stress_rate = rate;
|
||
step.nonlinear_geometry = nonlinear;
|
||
define(&mut m, ElementKind::Hex8, elastic(1000.0, 0.3, 1.0), step);
|
||
let mut d = Driver::new(&m);
|
||
let stretch = |x: [f64; 3]| [1.01 * x[0], x[1], x[2]];
|
||
let before = d.go(stretch, 1.0);
|
||
let n = 50;
|
||
let mut after = before;
|
||
for k in 1..=n {
|
||
let angle = std::f64::consts::FRAC_PI_2 * k as f64 / n as f64;
|
||
after = d.go(|x| rotation_z(angle, stretch(x)), 1.0);
|
||
}
|
||
// R σ Rᵀ for 90° about z: xx ↔ yy, xy → −xy, xz → yz, yz → −xz.
|
||
let want = [before[1], before[0], before[2], -before[3], -before[5], before[4]];
|
||
let error = after.iter().zip(want).map(|(a, b)| (a - b).abs()).fold(0.0, f64::max) / before[0].abs();
|
||
eprintln!("rigid rotation {rate:?} nonlinear {nonlinear}: σ before {before:.4?}, after {after:.4?}; error {error:.2e} of σ_xx");
|
||
if nonlinear {
|
||
assert!(error < 1e-12, "{rate:?}: {error}");
|
||
// Forces of the turned stress balance on the turned element: ∑ f = 0, ∑ x × f = 0.
|
||
let f = &d.ws.stress_force;
|
||
let total: f64 = (0..3).map(|k| f.iter().skip(k).step_by(3).sum::<f64>().abs()).sum();
|
||
assert!(total < 1e-10 * before[0].abs());
|
||
} else {
|
||
assert!(error > 0.5, "the small-strain update is not objective: {error}");
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Uniaxial stress to a logarithmic strain of 0.3 and back by 0.006 (inside the elastic range) in a single hexahedron with
|
||
/// bilinear isotropic hardening: the lateral stretch is iterated each increment so the lateral
|
||
/// stresses vanish. The Cauchy stress follows σ = σ_y + H ε_p with ln λ = σ/E + ε_p (for a
|
||
/// hypoelastic–plastic law under coaxial stretching the integral of D is the logarithmic
|
||
/// strain), then unloads with slope E.
|
||
#[test]
|
||
fn j2_uniaxial_stress_strain_curve_and_unloading_at_large_strain() {
|
||
let (young, poisson, y0, h) = (200_000.0, 0.3, 250.0, 2000.0);
|
||
for kind in [ElementKind::Hex8, ElementKind::Tet4] {
|
||
let mut m = block([1, 1, 1], [1.0; 3], kind);
|
||
let plastic = Plastic::Bilinear { yield_stress: y0, hardening_modulus: h, hardening: Hardening::Isotropic };
|
||
define(&mut m, kind, Material { plastic: Some(plastic), ..elastic(young, poisson, 7.8e-9) }, ExplicitDynamic::new(1.0));
|
||
let mut d = Driver::new(&m);
|
||
let mut path: Vec<f64> = (1..=600).map(|k| 0.3 * k as f64 / 600.0).collect();
|
||
path.extend((1..=200).map(|k| 0.3 - 0.006 * k as f64 / 200.0));
|
||
let (mut lateral, mut worst, mut peak) = (0.0f64, 0.0f64, (0.0, 0.0, 0.0));
|
||
let mut last = [0.0; 6];
|
||
for (step, axial) in path.iter().enumerate() {
|
||
// Secant iteration on the lateral log strain for σ_yy = 0 (σ_zz follows by symmetry).
|
||
let saved = (d.state.clone(), d.u.clone(), d.time);
|
||
let trial = |lat: f64, d: &mut Driver| {
|
||
(d.state, d.u, d.time) = saved.clone();
|
||
let (a, l) = (axial.exp(), lat.exp());
|
||
d.go(|x| [a * x[0], l * x[1], l * x[2]], 1e-6)
|
||
};
|
||
let (mut a, mut b) = (lateral, lateral - 1e-4);
|
||
let (mut fa, mut fb) = (trial(a, &mut d)[1], trial(b, &mut d)[1]);
|
||
for _ in 0..60 {
|
||
if fb.abs() < 1e-9 || fb == fa {
|
||
break;
|
||
}
|
||
let c = b - fb * (b - a) / (fb - fa);
|
||
(a, fa) = (b, fb);
|
||
b = c;
|
||
fb = trial(b, &mut d)[1];
|
||
}
|
||
lateral = b;
|
||
let s = trial(b, &mut d);
|
||
let p = d.state.elements.iter().map(|e| e.point.plastic_strain).fold(0.0, f64::max);
|
||
if step < 600 {
|
||
// Loading: σ = y0 + H ε_p and ε = σ/E + ε_p, so σ = y0 + E_t (ε − y0/E).
|
||
let et = young * h / (young + h);
|
||
let exact = if *axial * young <= y0 { young * axial } else { y0 + et * (axial - y0 / young) };
|
||
worst = worst.max((s[0] - exact).abs() / exact);
|
||
peak = (*axial, s[0], p);
|
||
} else {
|
||
// Unloading: elastic from the peak.
|
||
let exact = peak.1 - young * (peak.0 - axial);
|
||
worst = worst.max((s[0] - exact).abs() / peak.1);
|
||
assert!((p - peak.2).abs() < 1e-12, "no plastic flow while unloading");
|
||
}
|
||
last = s;
|
||
}
|
||
eprintln!(
|
||
"J2 uniaxial {}: peak σ {:.4} at ln λ = 0.3 (exact {:.4}); after unloading to 0.294 σ = {:.4}; largest relative error {worst:.2e}",
|
||
kind.name(),
|
||
peak.1,
|
||
y0 + young * h / (young + h) * (0.3 - y0 / young),
|
||
last[0]
|
||
);
|
||
assert!(worst < 1e-6, "{}: {worst}", kind.name());
|
||
}
|
||
}
|
||
|
||
/// The Taylor bar: a copper-like bar of square section hits a rigid wall at 190 m/s (units mm,
|
||
/// t, s, MPa) and mushrooms. Taylor's one-dimensional analysis of a rigid–perfectly plastic bar
|
||
/// (Taylor 1948) relates the final length L₁ and the length X of the undeformed rear part to the
|
||
/// yield stress: σ_y = ρ U² (L − X) / (2 (L − L₁) ln(L / X)). Applied to the computed shape it
|
||
/// should give back the yield stress the material has, within the accuracy of the estimate
|
||
/// itself (which neglects elastic waves and radial inertia): 25%.
|
||
#[test]
|
||
fn taylor_bar_recovers_its_yield_stress() {
|
||
let (rho, young, poisson, yield_stress) = (8.93e-9, 117_000.0, 0.35, 400.0);
|
||
let (length, width, speed) = (25.4, 7.6, 190_000.0);
|
||
let end = 80e-6;
|
||
let mut step = ExplicitDynamic::new(end);
|
||
step.output_interval = Some(end / 4.0);
|
||
step.history_interval = Some(end / 400.0);
|
||
step.initial_velocities = vec![InitialVelocity { set: "nodes".into(), velocity: [0.0, 0.0, -speed] }];
|
||
step.rigid_walls = vec![RigidWall { set: "nodes".into(), point: [0.0; 3], normal: [0.0, 0.0, 1.0], friction: Friction::NONE }];
|
||
step.hourglass = Hourglass { kind: HourglassKind::Stiffness, coefficient: 0.1 };
|
||
let mut m = block([4, 4, 30], [width, width, length], ElementKind::Hex8);
|
||
node_set(&mut m, "nodes", |_| true);
|
||
let plastic = Plastic::Bilinear { yield_stress, hardening_modulus: 0.0, hardening: Hardening::Isotropic };
|
||
define(&mut m, ElementKind::Hex8, Material { plastic: Some(plastic), ..elastic(young, poisson, rho) }, step);
|
||
let r = solve_explicit(&m, None, &Limits::default()).unwrap();
|
||
let f = r.frames.last().unwrap();
|
||
// Final length: the top face's mean height; undeformed part: the axis above the last
|
||
// element with plastic strain, along the centre column of elements.
|
||
let z = |label: u32| m.node(label).unwrap().xyz[2] + f.displacement[r.nodes.binary_search(&label).unwrap()][2];
|
||
let top: Vec<u32> = m.nodes().iter().filter(|(_, n)| n.xyz[2] == length).map(|(l, _)| *l).collect();
|
||
let l1 = top.iter().map(|l| z(*l)).sum::<f64>() / top.len() as f64;
|
||
let column: Vec<usize> = (0..30).map(|k| 4 * 4 * k + 4 + 1).collect();
|
||
// The rear part counts as undeformed from the first element (from the impact end) whose
|
||
// plastic strain is below 1%, small against the mushroom's.
|
||
let deformed = column.iter().position(|&i| f.plastic_strain[i] < 0.01).unwrap_or(30);
|
||
let x = length * (30 - deformed) as f64 / 30.0;
|
||
let estimate = rho * speed * speed * (length - x) / (2.0 * (length - l1) * (length / x).ln());
|
||
let kinetic = r.energy.kinetic / r.initial_kinetic;
|
||
eprintln!(
|
||
"Taylor bar: L₁/L = {:.4}, X/L = {:.4}; Taylor's estimate σ_y = {estimate:.1} vs {yield_stress} ({:+.1}%); final KE/KE₀ {kinetic:.4}; plastic work {:.3e} of {:.3e}; balance {:.2e}; {} increments",
|
||
l1 / length,
|
||
x / length,
|
||
100.0 * (estimate / yield_stress - 1.0),
|
||
r.energy.plastic,
|
||
r.initial_kinetic,
|
||
r.energy.relative_error(r.initial_kinetic),
|
||
r.increments
|
||
);
|
||
assert!(r.completed && l1 < 0.95 * length && l1 > 0.5 * length);
|
||
assert!((estimate / yield_stress - 1.0).abs() < 0.25, "{estimate}");
|
||
assert!(r.energy.relative_error(r.initial_kinetic) < 0.01, "{:?}", r.energy);
|
||
assert!(r.energy.plastic > 0.5 * r.initial_kinetic, "most of the impact energy goes into plastic work: {:?}", r.energy);
|
||
}
|
||
|
||
fn taylor_with(material: Material) -> (Model, ExplicitResult) {
|
||
let end = 40e-6;
|
||
let mut step = ExplicitDynamic::new(end);
|
||
step.output_interval = Some(end / 2.0);
|
||
step.initial_velocities = vec![InitialVelocity { set: "nodes".into(), velocity: [0.0, 0.0, -300_000.0] }];
|
||
step.rigid_walls = vec![RigidWall { set: "nodes".into(), point: [0.0; 3], normal: [0.0, 0.0, 1.0], friction: Friction::NONE }];
|
||
step.hourglass = Hourglass { kind: HourglassKind::Stiffness, coefficient: 0.1 };
|
||
let mut m = block([2, 2, 12], [7.6, 7.6, 25.4], ElementKind::Hex8);
|
||
node_set(&mut m, "nodes", |_| true);
|
||
define(&mut m, ElementKind::Hex8, material, step);
|
||
let r = solve_explicit(&m, None, &Limits::default()).unwrap();
|
||
(m, r)
|
||
}
|
||
|
||
/// Johnson–Cook strength, adiabatic heating, damage and erosion in a fast impact: elements at
|
||
/// the impact face reach the erosion strain (compression keeps their damage low, as the
|
||
/// fracture strain grows with confinement), are removed and reported; the run stays balanced
|
||
/// and repeatable bit for bit.
|
||
#[test]
|
||
fn johnson_cook_impact_erodes_elements() {
|
||
// OFHC copper constants as published by Johnson & Cook (1983, 1985), units mm, t, s, MPa, K.
|
||
let plastic = Plastic::JohnsonCook {
|
||
a: 90.0,
|
||
b: 292.0,
|
||
n: 0.31,
|
||
c: 0.025,
|
||
m: 1.09,
|
||
reference_strain_rate: 1.0,
|
||
room_temperature: 294.0,
|
||
melt_temperature: 1356.0,
|
||
specific_heat: Some(3.83e8),
|
||
taylor_quinney: 0.9,
|
||
};
|
||
let damage = Damage::JohnsonCook { d1: 0.54, d2: 4.89, d3: -3.03, d4: 0.014, d5: 1.12 };
|
||
let material = Material {
|
||
plastic: Some(plastic),
|
||
damage: Some(damage),
|
||
erosion: Some(Erosion { plastic_strain: 0.8 }),
|
||
..elastic(124_000.0, 0.34, 8.96e-9)
|
||
};
|
||
let (_, r) = taylor_with(material.clone());
|
||
let f = r.frames.last().unwrap();
|
||
let hottest = f.temperature.iter().copied().fold(0.0, f64::max);
|
||
eprintln!(
|
||
"Johnson–Cook impact: {} elements eroded ({:?}…), eroded energy {:.3e}; hottest element {hottest:.1} K; largest damage {:.3}; balance {:.2e}",
|
||
r.eroded.len(),
|
||
r.eroded.first(),
|
||
r.energy.eroded,
|
||
f.damage.iter().copied().fold(0.0, f64::max),
|
||
r.energy.relative_error(r.initial_kinetic)
|
||
);
|
||
assert!(!r.eroded.is_empty() && r.eroded.iter().all(|e| matches!(e.failure, Failure::Damage | Failure::PlasticStrain)));
|
||
assert_eq!(f.failed.iter().filter(|x| **x).count(), r.eroded.len());
|
||
assert!(hottest > 294.0 + 10.0);
|
||
assert!(r.energy.eroded > 0.0 && r.warnings.iter().any(|w| w.contains("eroded")));
|
||
assert!(r.energy.relative_error(r.initial_kinetic) < 0.02, "{:?}", r.energy);
|
||
let (_, again) = taylor_with(material);
|
||
assert!(r == again, "bit-identical runs");
|
||
}
|
||
|
||
/// Damage with a constant fracture strain D1 (D2 = D4 = D5 = 0): an element stretched in
|
||
/// tension fails when its plastic strain reaches D1, and from then on carries no stress.
|
||
#[test]
|
||
fn damage_fails_an_element_at_its_fracture_strain() {
|
||
let mut m = block([1, 1, 1], [1.0; 3], ElementKind::Hex8);
|
||
let plastic = Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 1000.0, hardening: Hardening::Isotropic };
|
||
let damage = Damage::JohnsonCook { d1: 0.2, d2: 0.0, d3: 0.0, d4: 0.0, d5: 0.0 };
|
||
define(
|
||
&mut m,
|
||
ElementKind::Hex8,
|
||
Material { plastic: Some(plastic), damage: Some(damage), ..elastic(200_000.0, 0.3, 7.8e-9) },
|
||
ExplicitDynamic::new(1.0),
|
||
);
|
||
let mut d = Driver::new(&m);
|
||
let mut failed_at = None;
|
||
for k in 1..=1000 {
|
||
let stretch = 1.0 + 0.5 * k as f64 / 1000.0;
|
||
let before = d.state.elements[0].point.plastic_strain;
|
||
d.go(|x| [stretch * x[0], x[1] / stretch.sqrt(), x[2] / stretch.sqrt()], 1e-6);
|
||
let e = &d.state.elements[0];
|
||
if e.failed.is_some() {
|
||
failed_at = Some((before, e.point.plastic_strain, e.failed));
|
||
break;
|
||
}
|
||
}
|
||
let (before, at, why) = failed_at.expect("the element fails");
|
||
eprintln!("damage: failed between plastic strain {before:.5} and {at:.5} (D1 = 0.2)");
|
||
assert_eq!(why, Some(Failure::Damage));
|
||
assert!(before < 0.2 && at >= 0.2 - 1e-12, "{before} {at}");
|
||
// Failed: no forces from it.
|
||
d.go(|x| [1.6 * x[0], x[1], x[2]], 1e-6);
|
||
assert!(d.ws.stress_force.iter().all(|f| *f == 0.0));
|
||
}
|
||
|
||
/// Inelastic models that cannot run are refused by name, never panics.
|
||
#[test]
|
||
fn inelastic_hostile_models_are_errors() {
|
||
// Plasticity with orthotropic elasticity is refused by the model itself.
|
||
let mut m = block([1, 1, 1], [1.0; 3], ElementKind::Hex8);
|
||
let ortho = Elastic::EngineeringConstants { young: [1.0; 3], poisson: [0.0; 3], shear: [0.5; 3] };
|
||
let plastic = Plastic::Bilinear { yield_stress: 1.0, hardening_modulus: 0.0, hardening: Hardening::Isotropic };
|
||
assert!(m.add_material("Bad", Material { elastic: Some(ortho), density: Some(1.0), plastic: Some(plastic), ..Material::default() }).is_err());
|
||
// An element turned inside out stops the run with an error naming it.
|
||
define(&mut m, ElementKind::Hex8, elastic(1.0, 0.3, 1.0), ExplicitDynamic::new(1.0));
|
||
let mut d = Driver::new(&m);
|
||
d.go(|x| [x[0], x[1], 0.5 * x[2]], 1.0);
|
||
let e = d.try_go(|x| [x[0], x[1], 0.5 - x[2]], 1.0).unwrap_err().to_string();
|
||
assert!(e.contains("element 1") && e.contains("inside out"), "{e}");
|
||
// A crushed element shrinks the increment until the cap stops the run, with a warning.
|
||
let mut m = block([1, 1, 1], [1.0; 3], ElementKind::Hex8);
|
||
let mut step = ExplicitDynamic::new(1.0);
|
||
step.max_increments = 20_000;
|
||
define(&mut m, ElementKind::Hex8, elastic(1.0, 0.3, 1.0), step);
|
||
node_set(&mut m, "top", |x| x[2] == 1.0);
|
||
node_set(&mut m, "bottom", |x| x[2] == 0.0);
|
||
bc(&mut m, "hold", BoundaryCondition::Fixed { target: Target::Set("bottom".into()) });
|
||
bc(
|
||
&mut m,
|
||
"crush",
|
||
BoundaryCondition::Displacement { target: Target::Set("top".into()), values: [Some(0.0), Some(0.0), Some(-2.0), None, None, None] },
|
||
);
|
||
match solve_explicit(&m, None, &Limits::default()) {
|
||
Ok(r) => assert!(!r.completed && r.warnings.iter().any(|w| w.contains("cap"))),
|
||
Err(e) => assert!(e.to_string().contains("element 1"), "{e}"),
|
||
}
|
||
}
|