diff --git a/crates/decks/src/abaqus.rs b/crates/decks/src/abaqus.rs index 89f922a..99e3051 100644 --- a/crates/decks/src/abaqus.rs +++ b/crates/decks/src/abaqus.rs @@ -108,7 +108,8 @@ //! //! 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 +//! energy, and likewise `MOONEY-RIVLIN` (C₁₀, C₀₁, D₁) and `YEOH` (C₁₀, C₂₀, C₃₀, D₁, the higher +//! volumetric terms zero); 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 @@ -129,8 +130,8 @@ use std::collections::{BTreeMap, BTreeSet}; use std::fmt::Write as _; use fieldcraft_model::{ - BeamProfile, BoundaryCondition, CouplingKind, Elastic, Element, ElementKind, Hardening, Integration, Load, Model, OutputVariable, Plastic, - Procedure, Resolved, Section, SectionKind, Step, Symmetry, Target, + BeamProfile, BoundaryCondition, CouplingKind, Elastic, Element, ElementKind, Hardening, Hyperelastic, Integration, Load, Model, OutputVariable, + Plastic, Procedure, Resolved, Section, SectionKind, Step, Symmetry, Target, }; use crate::{Error, Result, unsupported}; @@ -809,11 +810,20 @@ 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 }) + if has_nonlinear_step(model) + && let Some(h) = m.hyperelastic => + { + // The polynomial constants of the option and U = (J − 1)²/D1 with D1 = 2/κ. + let (shear, bulk) = Hyperelastic::moduli(young, poisson); + let [c10, c01, c20, c30] = h.constants(shear); + let (keyword, values) = match h { + Hyperelastic::NeoHookean {} => ("NEO HOOKE", vec![c10, 2.0 / bulk]), + Hyperelastic::MooneyRivlin { .. } => ("MOONEY-RIVLIN", vec![c10, c01, 2.0 / bulk]), + Hyperelastic::Yeoh { .. } => ("YEOH", vec![c10, c20, c30, 2.0 / bulk]), + }; + self.line(&format!("*HYPERELASTIC, {keyword}")); + self.numbers_line("", &values)?; } Some(Elastic::Isotropic { young, poisson }) => { self.line("*ELASTIC, TYPE=ISOTROPIC"); diff --git a/crates/decks/tests/abaqus.rs b/crates/decks/tests/abaqus.rs index 543d601..54c694f 100644 --- a/crates/decks/tests/abaqus.rs +++ b/crates/decks/tests/abaqus.rs @@ -790,6 +790,19 @@ fn nonlinear_steps_and_hyperelastic_materials_are_written() { assert!((values[0] - c10).abs() < 1e-9 * c10 && (values[1] - d1).abs() < 1e-9 * d1, "{line}"); assert!(!deck.contains("*ELASTIC"), "{deck}"); + // Mooney–Rivlin and Yeoh: C10 = μ/2 − C01, C01, D1 and C10, C20, C30, D1. + m.replace_material("Steel", Material { hyperelastic: Some(Hyperelastic::MooneyRivlin { c01: 10000.0 }), ..steel() }).unwrap(); + let deck = abaqus::write(&m).unwrap(); + let line = deck.split("*HYPERELASTIC, MOONEY-RIVLIN\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 - 10000.0)).abs() < 1e-9 * c10 && values[1] == 10000.0 && (values[2] - d1).abs() < 1e-9 * d1, "{line}"); + m.replace_material("Steel", Material { hyperelastic: Some(Hyperelastic::Yeoh { c20: -50.0, c30: 2.0 }), ..steel() }).unwrap(); + let deck = abaqus::write(&m).unwrap(); + let line = deck.split("*HYPERELASTIC, YEOH\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..3] == [-50.0, 2.0] && (values[3] - d1).abs() < 1e-9 * d1, "{line}"); + m.replace_material("Steel", Material { hyperelastic: Some(Hyperelastic::NeoHookean {}), ..steel() }).unwrap(); + 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(); diff --git a/crates/engine/src/cmd/define.rs b/crates/engine/src/cmd/define.rs index f4bf3de..b1ebe58 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), "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"#, + 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"} | {"type": "mooney_rivlin", "c01": 0.1} (C10 = μ/2 − C01) | {"type": "yeoh", "c20": -0.01, "c30": 0.001} (C10 = μ/2) (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, diff --git a/crates/fem/src/constitutive.rs b/crates/fem/src/constitutive.rs index 529bd2e..6cba2bf 100644 --- a/crates/fem/src/constitutive.rs +++ b/crates/fem/src/constitutive.rs @@ -36,6 +36,9 @@ //! §6.6, eqs. 6.165–6.169 and 6.193, in Kirchhoff form). The split is what the mixed (hybrid) //! elements need: they take the isochoric part from [`MaterialModel::isochoric`] and the //! volumetric one from [`MaterialModel::volumetric`] with the element's own dilatation. +//! - **Mooney–Rivlin and Yeoh**: decoupled polynomials +//! W̄ = C₁₀ (Ī₁ − 3) + C₀₁ (Ī₂ − 3) + C₂₀ (Ī₁ − 3)² + C₃₀ (Ī₁ − 3)³ with the neo-Hookean volume; their fictitious stress and +//! elasticity from Holzapfel 2000, eqs. 6.193 and 6.168 ([`MaterialModel::Polynomial`]). //! - **J2 plasticity** ([`crate::elastoplastic`]): Simo's multiplicative plasticity with the //! exponential return map in the principal logarithmic strains (finite strain), or the //! additive radial return (small strain, [`MaterialModel::small_strain`]); decoupled too, the @@ -73,6 +76,14 @@ pub enum MaterialModel { shear: f64, bulk: f64, }, + /// A decoupled polynomial in the isochoric invariants, + /// W̄ = C₁₀ (Ī₁ − 3) + C₀₁ (Ī₂ − 3) + C₂₀ (Ī₁ − 3)² + C₃₀ (Ī₁ − 3)³ with `c` = (C₁₀, C₀₁, C₂₀, C₃₀), and U = κ/2 (J − 1)²: + /// Mooney–Rivlin and Yeoh. + Polynomial { + name: &'static str, + c: [f64; 4], + bulk: f64, + }, /// Von Mises plasticity (shared: one per material, whatever the element count). Plastic(Arc), } @@ -123,27 +134,35 @@ fn jacobian(f: &Mat3) -> Result { } impl MaterialModel { - /// The model of a material's elastic behaviour in a nonlinear step: neo-Hookean when the - /// material asks for it (its isotropic constants give μ and κ), St Venant–Kirchhoff + /// The model of a material's elastic behaviour in a nonlinear step: the hyperelastic option + /// when the material has one (its isotropic constants give μ and κ), St Venant–Kirchhoff /// otherwise. pub fn from_model(elastic: &Elastic, hyperelastic: Option<&Hyperelastic>) -> Result { match hyperelastic { None => Ok(MaterialModel::StVenantKirchhoff(Box::new(elasticity(elastic)?))), - Some(Hyperelastic::NeoHookean {}) => { - let (young, poisson) = crate::material::isotropic_constants(elastic, "neo-Hookean")?; + Some(h) => { + let (young, poisson) = crate::material::isotropic_constants(elastic, h.name())?; let (shear, bulk) = Hyperelastic::moduli(young, poisson); - Ok(MaterialModel::NeoHookean { shear, bulk }) + Ok(Self::hyperelastic(h, shear, bulk)) } } } + /// The hyperelastic model of an option with small-strain moduli μ and κ. + fn hyperelastic(h: &Hyperelastic, shear: f64, bulk: f64) -> MaterialModel { + match h { + Hyperelastic::NeoHookean {} => MaterialModel::NeoHookean { shear, bulk }, + _ => MaterialModel::Polynomial { name: h.name(), c: h.constants(shear), bulk }, + } + } + /// The model of an element whose elasticity is `d` (the small-strain one) and whose material - /// has the large-strain option `hyperelastic`: neo-Hookean takes μ = D₄₄ and κ = D₁₂ + ⅔ μ of - /// an isotropic `d`, St Venant–Kirchhoff `d` itself. + /// has the large-strain option `hyperelastic`: the hyperelastic models take μ = D₄₄ and + /// κ = D₁₂ + ⅔ μ of an isotropic `d`, St Venant–Kirchhoff `d` itself. pub fn from_elasticity(d: &Elasticity, hyperelastic: Option<&Hyperelastic>) -> MaterialModel { match hyperelastic { None => MaterialModel::StVenantKirchhoff(Box::new(*d)), - Some(Hyperelastic::NeoHookean {}) => MaterialModel::NeoHookean { shear: d[3][3], bulk: d[0][1] + 2.0 / 3.0 * d[3][3] }, + Some(h) => Self::hyperelastic(h, d[3][3], d[0][1] + 2.0 / 3.0 * d[3][3]), } } @@ -151,6 +170,7 @@ impl MaterialModel { match self { MaterialModel::StVenantKirchhoff(_) => "St Venant–Kirchhoff", MaterialModel::NeoHookean { .. } => "neo-Hookean", + MaterialModel::Polynomial { name, .. } => name, MaterialModel::Plastic(_) => "J2 plastic", } } @@ -180,7 +200,7 @@ impl MaterialModel { /// Whether the strain energy splits into isochoric and volumetric parts (what mixed /// elements need). pub fn is_decoupled(&self) -> bool { - matches!(self, MaterialModel::NeoHookean { .. } | MaterialModel::Plastic(_)) + !matches!(self, MaterialModel::StVenantKirchhoff(_)) } /// The Kirchhoff stress and spatial tangent at deformation gradient `f` and (for @@ -197,7 +217,7 @@ impl MaterialModel { let t = push_forward(f); Ok(Response { stress: std::array::from_fn(|a| (0..6).map(|b| t[a][b] * s[b]).sum()), tangent: sandwich(&t, d) }) } - MaterialModel::NeoHookean { .. } | MaterialModel::Plastic(_) => { + MaterialModel::NeoHookean { .. } | MaterialModel::Polynomial { .. } | MaterialModel::Plastic(_) => { let iso = self.isochoric(f, temperature, committed, updated)?; let (p, dp) = self.volumetric(j, false).ok_or_else(|| Error::Invalid("volumetric response".into()))?; Ok(add_pressure(&iso, j, p, j * dp)) @@ -212,6 +232,7 @@ impl MaterialModel { let mu = match self { MaterialModel::NeoHookean { shear, .. } => *shear, MaterialModel::Plastic(j2) => return elastoplastic::finite_isochoric(j2, f, temperature, committed, updated), + MaterialModel::Polynomial { c, .. } => return polynomial_isochoric(c, f), MaterialModel::StVenantKirchhoff(_) => { return Err(Error::NotSupported(format!("the isochoric–volumetric split of a {} material", self.name()))); } @@ -241,7 +262,7 @@ impl MaterialModel { /// `small`, J = 1 + tr ε of small strain, where every material's is K/2 (J − 1)²). pub fn volumetric(&self, j: f64, small: bool) -> Option<(f64, f64)> { match self { - MaterialModel::NeoHookean { bulk, .. } => Some((bulk * (j - 1.0), *bulk)), + MaterialModel::NeoHookean { bulk, .. } | MaterialModel::Polynomial { bulk, .. } => Some((bulk * (j - 1.0), *bulk)), MaterialModel::Plastic(j2) => elastoplastic::volumetric(j2, j, small), MaterialModel::StVenantKirchhoff(_) => None, } @@ -262,6 +283,51 @@ impl MaterialModel { } } +/// The isochoric response of the polynomial energy W̄(Ī₁, Ī₂) with `c` = (C₁₀, C₀₁, C₂₀, C₃₀) +/// (Holzapfel 2000, §6.5–6.6, in Kirchhoff form): the fictitious stress +/// τ̄ = 2 (Ψ₁ + Ī₁ Ψ₂) b̄ − 2 Ψ₂ b̄², τ_iso = dev τ̄, and with the fictitious elasticity +/// c̄ = 4 (Ψ₁₁ + Ψ₂) b̄ ⊗ b̄ − 4 Ψ₂ b̄ ⊙ b̄ (eq. 6.193 pushed forward, Ψ₁₂ = Ψ₂₂ = 0 here), +/// c_iso = ℙ : c̄ : ℙ + ⅔ tr τ̄ ℙ − ⅔ (τ_iso ⊗ I + I ⊗ τ_iso) (eq. 6.168). With C₁₀ alone it is +/// the neo-Hookean response. +fn polynomial_isochoric(c: &[f64; 4], f: &Mat3) -> Result { + let [c10, c01, c20, c30] = *c; + let j = jacobian(f)?; + let scale = j.powf(-2.0 / 3.0); + let b = mul(f, &transpose(f)).map(|row| row.map(|v| scale * v)); + let b2 = mul(&b, &b); + let i1 = b[0][0] + b[1][1] + b[2][2]; + let x = i1 - 3.0; + let (psi1, psi2, psi11) = (c10 + 2.0 * c20 * x + 3.0 * c30 * x * x, c01, 2.0 * c20 + 6.0 * c30 * x); + let (bv, b2v) = (voigt(&b), voigt(&b2)); + let fict: [f64; 6] = std::array::from_fn(|a| 2.0 * (psi1 + i1 * psi2) * bv[a] - 2.0 * psi2 * b2v[a]); + let tr = fict[0] + fict[1] + fict[2]; + let mut tau = fict; + for v in tau.iter_mut().take(3) { + *v -= tr / 3.0; + } + let mut cbar = [[0.0; 6]; 6]; + for (a, &(i, jj)) in PAIRS.iter().enumerate() { + for (bb, &(k, l)) in PAIRS.iter().enumerate() { + cbar[a][bb] = 4.0 * (psi11 + psi2) * bv[a] * bv[bb] - 2.0 * psi2 * (b[i][k] * b[jj][l] + b[i][l] * b[jj][k]); + } + } + // ℙ : c̄ : ℙ = c̄ − ⅓ I ⊗ (I : c̄) − ⅓ (c̄ : I) ⊗ I + ⅑ (I : c̄ : I) I ⊗ I. + let row: [f64; 6] = std::array::from_fn(|a| (0..3).map(|bb| cbar[a][bb]).sum()); + let col: [f64; 6] = std::array::from_fn(|bb| (0..3).map(|a| cbar[a][bb]).sum()); + let all: f64 = (0..3).map(|a| row[a]).sum(); + let i4 = identity4(); + let mut out = [[0.0; 6]; 6]; + for a in 0..6 { + for bb in 0..6 { + let (ia, ib) = (if a < 3 { 1.0 } else { 0.0 }, if bb < 3 { 1.0 } else { 0.0 }); + let projected = cbar[a][bb] - ia * col[bb] / 3.0 - row[a] * ib / 3.0 + ia * ib * all / 9.0; + let projector = i4[a][bb] - ia * ib / 3.0; + out[a][bb] = projected + 2.0 / 3.0 * tr * projector - 2.0 / 3.0 * (tau[a] * ib + ia * tau[bb]); + } + } + Ok(Response { stress: tau, tangent: out }) +} + /// An isochoric response with a pressure p acting on dilatation J: τ + J p I, and the tangent /// c + J (p̃ I ⊗ I − 2p 𝕀) with p̃ = p + `j_dp` (J dp/dJ, zero when p is held). pub fn add_pressure(iso: &Response, j: f64, p: f64, j_dp: f64) -> Response { @@ -361,6 +427,21 @@ mod tests { } } } + // Mooney–Rivlin and Yeoh: consistent and symmetric; with C₁₀ alone the neo-Hookean model. + for m in [ + MaterialModel::Polynomial { name: "Mooney–Rivlin", c: [120.0, 30.0, 0.0, 0.0], bulk: 5000.0 }, + MaterialModel::Polynomial { name: "Yeoh", c: [150.0, 0.0, -12.0, 1.5], bulk: 5000.0 }, + ] { + check_tangent(&m, &deformation()); + check_tangent(&m, &IDENTITY); + let c = m.response(&deformation(), None, &[], &mut []).unwrap().tangent; + assert!((0..36).all(|k| (c[k / 6][k % 6] - c[k % 6][k / 6]).abs() < 1e-9 * 5000.0), "{} symmetric", m.name()); + } + let neo = MaterialModel::NeoHookean { shear: 300.0, bulk: 5000.0 }; + let poly = MaterialModel::Polynomial { name: "neo", c: [150.0, 0.0, 0.0, 0.0], bulk: 5000.0 }; + let (a, b) = (neo.response(&deformation(), None, &[], &mut []).unwrap(), poly.response(&deformation(), None, &[], &mut []).unwrap()); + assert!(a.stress.iter().zip(&b.stress).all(|(x, y)| (x - y).abs() < 1e-11 * 300.0)); + assert!(a.tangent.iter().flatten().zip(b.tangent.iter().flatten()).all(|(x, y)| (x - y).abs() < 1e-10 * 5000.0)); let ortho = crate::material::orthotropic([1000.0, 300.0, 500.0], [0.3, 0.2, 0.25], [100.0, 120.0, 90.0]).unwrap(); check_tangent(&MaterialModel::StVenantKirchhoff(Box::new(ortho)), &deformation()); } diff --git a/crates/fem/src/finite.rs b/crates/fem/src/finite.rs index 2ba5d08..d65fe10 100644 --- a/crates/fem/src/finite.rs +++ b/crates/fem/src/finite.rs @@ -597,6 +597,10 @@ mod tests { for (kind, x, np) in [(ElementKind::Hex8, &x8, 1), (ElementKind::Hex20, &x20, 4), (ElementKind::Tet10, &x10, 1)] { let hybrid = Solid::new(kind, x, Formulation::Hybrid { pressures: np }).unwrap(); check_tangent(&hybrid, &neo, &displacement(x.len(), 0.08)); + let mr = MaterialModel::Polynomial { name: "Mooney–Rivlin", c: [100.0, 50.0, 0.0, 0.0], bulk: 3000.0 }; + let yeoh = MaterialModel::Polynomial { name: "Yeoh", c: [150.0, 0.0, -10.0, 2.0], bulk: 3000.0 }; + check_tangent(&hybrid, &mr, &displacement(x.len(), 0.08)); + check_tangent(&hybrid, &yeoh, &displacement(x.len(), 0.08)); assert!(hybrid.finite_strain(&svk, &displacement(x.len(), 0.0), &[], &mut [], false).is_err(), "SVK has no split"); } let im = Solid::new(ElementKind::Hex8, &x8, Formulation::IncompatibleModes).unwrap(); diff --git a/crates/model/src/materials.rs b/crates/model/src/materials.rs index ecbf320..2cdc840 100644 --- a/crates/model/src/materials.rs +++ b/crates/model/src/materials.rs @@ -88,6 +88,16 @@ pub enum Hyperelastic { /// reference volume W = μ/2 (Ī₁ − 3) + κ/2 (J − 1)², with Ī₁ = J^(−2/3) tr C (C₁₀ = μ/2, /// D₁ = 2/κ in the usual polynomial notation). NeoHookean {}, + /// Mooney–Rivlin (M. Mooney, J. Appl. Phys. 11:582–592, 1940; R. S. Rivlin, Phil. Trans. R. + /// Soc. A 241:379–397, 1948), decoupled: W = C₁₀ (Ī₁ − 3) + C₀₁ (Ī₂ − 3) + κ/2 (J − 1)². The + /// small-strain shear modulus 2 (C₁₀ + C₀₁) is the elastic μ, so `c01` sets the second + /// invariant's share and C₁₀ = μ/2 − C₀₁ (0 ≤ C₀₁ ≤ μ/2, both constants non-negative). + MooneyRivlin { c01: f64 }, + /// Yeoh (O. H. Yeoh, Rubber Chem. Technol. 66(5):754–771, 1993), decoupled: + /// W = C₁₀ (Ī₁ − 3) + C₂₀ (Ī₁ − 3)² + C₃₀ (Ī₁ − 3)³ + κ/2 (J − 1)², with C₁₀ = μ/2 from the + /// elastic constants (the small-strain shear modulus); `c20` (often negative) and `c30` + /// shape the large-strain response. + Yeoh { c20: f64, c30: f64 }, } impl Hyperelastic { @@ -101,8 +111,39 @@ impl Hyperelastic { pub fn name(&self) -> &'static str { match self { Hyperelastic::NeoHookean {} => "neo-Hookean", + Hyperelastic::MooneyRivlin { .. } => "Mooney–Rivlin", + Hyperelastic::Yeoh { .. } => "Yeoh", } } + + /// The polynomial constants (C₁₀, C₀₁, C₂₀, C₃₀) for the shear modulus μ. + pub fn constants(&self, shear: f64) -> [f64; 4] { + match *self { + Hyperelastic::NeoHookean {} => [shear / 2.0, 0.0, 0.0, 0.0], + Hyperelastic::MooneyRivlin { c01 } => [shear / 2.0 - c01, c01, 0.0, 0.0], + Hyperelastic::Yeoh { c20, c30 } => [shear / 2.0, 0.0, c20, c30], + } + } + + /// Check the option's constants against the elastic shear modulus μ. + pub fn validate(&self, shear: f64) -> Result<()> { + match *self { + Hyperelastic::NeoHookean {} => {} + Hyperelastic::MooneyRivlin { c01 } => { + if !(c01 >= 0.0 && c01 <= shear / 2.0) { + return Err(Error::Invalid(format!( + "the Mooney–Rivlin C01 = {c01} must be within [0, μ/2 = {}] (C10 = μ/2 − C01 and C01 both non-negative)", + shear / 2.0 + ))); + } + } + Hyperelastic::Yeoh { c20, c30 } => { + finite("Yeoh C20", c20)?; + finite("Yeoh C30", c30)?; + } + } + Ok(()) + } } /// Fluid behaviour. @@ -407,6 +448,9 @@ impl Material { if self.plastic.is_some() { return Err(Error::Invalid(format!("plasticity of a {} hyperelastic material is not supported yet", h.name()))); } + if let Some(Elastic::Isotropic { young, poisson }) = self.elastic { + h.validate(Hyperelastic::moduli(young, poisson).0)?; + } } if let Some(d) = self.density { positive("density", d)?; @@ -607,7 +651,28 @@ mod tests { 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()); + assert!(Material { plastic: Some(plastic), ..m.clone() }.validate().is_err()); + // Mooney–Rivlin and Yeoh: μ = E/(2(1 + ν)) = 1.00671 here, so C01 is within [0, μ/2]. + let mr: Material = serde_json::from_value(serde_json::json!({ + "elastic": {"type": "isotropic", "young": 3.0, "poisson": 0.49}, "hyperelastic": {"type": "mooney_rivlin", "c01": 0.1} + })) + .unwrap(); + mr.validate().unwrap(); + let mu = Hyperelastic::moduli(3.0, 0.49).0; + let c = Hyperelastic::MooneyRivlin { c01: 0.1 }.constants(mu); + assert!((2.0 * (c[0] + c[1]) - mu).abs() < 1e-15 && c[1] == 0.1); + let yeoh = Hyperelastic::Yeoh { c20: -0.01, c30: 0.001 }; + assert_eq!(yeoh.constants(mu), [mu / 2.0, 0.0, -0.01, 0.001]); + assert!(Material { hyperelastic: Some(yeoh), ..m.clone() }.validate().is_ok()); + for bad in [ + Hyperelastic::MooneyRivlin { c01: -0.1 }, + Hyperelastic::MooneyRivlin { c01: 0.6 }, + Hyperelastic::MooneyRivlin { c01: f64::NAN }, + Hyperelastic::Yeoh { c20: f64::INFINITY, c30: 0.0 }, + Hyperelastic::Yeoh { c20: 0.0, c30: f64::NAN }, + ] { + assert!(Material { hyperelastic: Some(bad), ..m.clone() }.validate().is_err(), "{bad:?}"); + } } #[test] diff --git a/crates/solver/tests/nonlinear.rs b/crates/solver/tests/nonlinear.rs index 69d5468..0946da9 100644 --- a/crates/solver/tests/nonlinear.rs +++ b/crates/solver/tests/nonlinear.rs @@ -433,3 +433,81 @@ fn nonlinear_results_are_deterministic() { assert_eq!(fieldcraft_parallel::with_threads(threads, || solve(&g)), a, "{threads} threads"); } } + +/// The principal Cauchy stresses of a decoupled polynomial energy at principal stretches `l`: +/// σᵢ = (λᵢ/J) ∂W/∂λᵢ by central differences of W itself (independent of the solver's +/// stress and tangent), W = C₁₀ (Ī₁ − 3) + C₀₁ (Ī₂ − 3) + C₂₀ (Ī₁ − 3)² + C₃₀ (Ī₁ − 3)³ + κ/2 (J − 1)². +fn polynomial_stress(c: [f64; 4], kappa: f64, l: [f64; 3]) -> [f64; 3] { + let w = |l: [f64; 3]| { + let j = l[0] * l[1] * l[2]; + let b: Vec = l.iter().map(|x| x * x * j.powf(-2.0 / 3.0)).collect(); + let i1 = b[0] + b[1] + b[2]; + let i2 = b[0] * b[1] + b[1] * b[2] + b[2] * b[0]; + let x = i1 - 3.0; + c[0] * x + c[1] * (i2 - 3.0) + c[2] * x * x + c[3] * x * x * x + 0.5 * kappa * (j - 1.0) * (j - 1.0) + }; + let j = l[0] * l[1] * l[2]; + std::array::from_fn(|i| { + let h = 1e-6 * l[i]; + let (mut p, mut m) = (l, l); + p[i] += h; + m[i] -= h; + l[i] / j * (w(p) - w(m)) / (2.0 * h) + }) +} + +/// Mooney–Rivlin and Yeoh blocks in the homogeneous modes of rubber testing (L. R. G. Treloar, +/// *The Physics of Rubber Elasticity*, 3rd ed., Oxford 1975, ch. 10): uniaxial and equibiaxial +/// tension and pure shear, against the closed form of the compressible decoupled energy (the +/// free stretches found by bisection so the free faces carry no stress). +#[test] +fn mooney_rivlin_and_yeoh_in_the_homogeneous_modes() { + let (e, nu) = (3.0, 0.49); + let (mu, kappa) = (e / (2.0 * (1.0 + nu)), e / (3.0 * (1.0 - 2.0 * nu))); + let hybrid = Formulation { hybrid: true, ..Formulation::default() }; + let (c01, c20, c30) = (0.2, -0.02, 0.004); + for option in [Hyperelastic::MooneyRivlin { c01 }, Hyperelastic::Yeoh { c20, c30 }] { + let c = option.constants(mu); + for (mode, lam) in [("uniaxial", 2.0), ("equibiaxial", 1.5), ("pure shear", 2.0)] { + // The stretches with the free one(s) l from σ_free(l) = 0. + let stretches = |l: f64| match mode { + "uniaxial" => [lam, l, l], + "equibiaxial" => [lam, lam, l], + _ => [lam, 1.0, l], + }; + let (mut lo, mut hi) = (0.05, 1.5); + for _ in 0..200 { + let mid = 0.5 * (lo + hi); + if polynomial_stress(c, kappa, stretches(mid))[2] > 0.0 { hi = mid } else { lo = mid } + } + let want = polynomial_stress(c, kappa, stretches(0.5 * (lo + hi))); + let mut g = grid([1, 1, 1], ElementKind::Hex20, |x| x); + g.define(iso(e, nu), None, hybrid, nonlinear()); + g.model.replace_material("Mat", Material { elastic: Some(iso(e, nu)), hyperelastic: Some(option), ..Material::default() }).unwrap(); + g.symmetry("x0", Symmetry::X, |x| near(x[0], 0.0)); + g.symmetry("y0", Symmetry::Y, |x| near(x[1], 0.0)); + g.symmetry("z0", Symmetry::Z, |x| near(x[2], 0.0)); + g.node_set("x1", |x| near(x[0], 1.0)); + g.bc( + "pull", + BoundaryCondition::Displacement { target: Target::Set("x1".into()), values: [Some(lam - 1.0), None, None, None, None, None] }, + ); + g.node_set("y1", |x| near(x[1], 1.0)); + let y = match mode { + "equibiaxial" => Some(lam - 1.0), + "pure shear" => Some(0.0), + _ => None, + }; + if let Some(v) = y { + g.bc("y1", BoundaryCondition::Displacement { target: Target::Set("y1".into()), values: [None, Some(v), None, None, None, None] }); + } + let r = solve(&g); + for s in &r.stress { + for k in 0..3 { + assert!((s[k] - want[k]).abs() < 1e-7 * want[0], "{} {mode}: {s:?} vs {want:?}", option.name()); + } + } + eprintln!("{} {mode} λ = {lam}: σ = {:.6?} (closed form {want:.6?})", option.name(), &r.stress[0][..3]); + } + } +}