diff --git a/crates/fem/src/finite.rs b/crates/fem/src/finite.rs index 698abd5..e340de2 100644 --- a/crates/fem/src/finite.rs +++ b/crates/fem/src/finite.rs @@ -515,6 +515,21 @@ impl Plane<'_> { } } +/// The kinematics of one point of a 2D element in a nonlinear step: F, det F_m, the volume share, the strain columns (current at finite strain), the small strain, +/// the temperature and the reference position. +struct PlaneKinematics { + f: Mat3, + /// det F_m (small strain: 1), and the dilatation the mixed element averages (small strain: + /// 1 + tr ε). + j: f64, + dilatation: f64, + dv: f64, + current: Vec<[f64; 3]>, + strain: [f64; 6], + heat: Heat, + x: [f64; 3], +} + /// One integration point of a 2D element in a nonlinear step: the total deformation gradient, /// det F_m, the reference position and the response. struct PlanePoint { @@ -525,8 +540,14 @@ struct PlanePoint { } impl Plane<'_> { - fn check_nonlinear(&self, u: &[[f64; 2]], kin: Kinematics) -> Result<()> { + fn check_nonlinear(&self, m: &MaterialModel, u: &[[f64; 2]], kin: Kinematics) -> Result<()> { check_displacements(u, self.ty.nodes)?; + if matches!(self.formulation, Formulation::Hybrid { .. }) && !m.is_decoupled() { + return Err(Error::NotSupported(format!( + "hybrid elements of a {} material in nonlinear steps (the mixed formulation needs a volumetric–isochoric split of the strain energy; give the material the neo-Hookean option)", + m.name() + ))); + } if let Theory::PlaneStress { .. } = self.theory { let why = match kin { Kinematics::Finite => "the thickness change needs a local iteration", @@ -558,7 +579,7 @@ impl Plane<'_> { updated: &mut [f64], tangent: bool, ) -> Result<(Finite, Vec)> { - self.check_nonlinear(u, kin)?; + self.check_nonlinear(m, u, kin)?; let n = self.ty.nodes; if let Some(t) = thermal { t.check(n)?; @@ -570,8 +591,9 @@ impl Plane<'_> { let mut k = tangent.then(|| Matrix::zeros(ndof)); let axisymmetric = self.theory == Theory::Axisymmetric; let small = kin == Kinematics::Small; - let mut out = Vec::with_capacity(rule.len()); - for (i, p) in rule.iter().enumerate() { + // The kinematics of every point first: the mixed element's pressure needs them all. + let mut kinematics = Vec::with_capacity(rule.len()); + for p in &rule { let at = self.at(p)?; let cols = at.g.get(..n).unwrap_or_default(); let mut f = IDENTITY; @@ -585,12 +607,11 @@ impl Plane<'_> { f[2][2] += ua[0] * g[2]; } } - let updated_point = updated.get_mut(i * len..(i + 1) * len).unwrap_or_default(); let heat = Heat::at(thermal, at.shape.values())?; - let (j, dv, current, r) = if small { + let point = if small { let strain = heat.mechanical_strain(&small_strain(&f)); - let r = m.small_strain(&strain, heat.temperature, point_history(committed, len, i), updated_point, true)?; - (1.0, at.dv, cols.to_vec(), r) + let volume = 1.0 + strain[0] + strain[1] + strain[2]; + PlaneKinematics { f, j: 1.0, dilatation: volume, dv: at.dv, current: cols.to_vec(), strain, heat, x: at.x } } else { let det2 = f[0][0] * f[1][1] - f[0][1] * f[1][0]; let j = det2 * f[2][2]; @@ -600,9 +621,29 @@ impl Plane<'_> { let inv = [[f[1][1] / det2, -f[0][1] / det2], [-f[1][0] / det2, f[0][0] / det2]]; let current: Vec<[f64; 3]> = cols.iter().map(|g| [g[0] * inv[0][0] + g[1] * inv[1][0], g[0] * inv[0][1] + g[1] * inv[1][1], g[2] / f[2][2]]).collect(); - let r = m.response(&heat.mechanical(&f), heat.temperature, point_history(committed, len, i), updated_point)?; - (j / heat.volume(), at.dv * heat.volume(), current, r) + let j = j / heat.volume(); + PlaneKinematics { f, j, dilatation: j, dv: at.dv * heat.volume(), current, strain: [0.0; 6], heat, x: at.x } }; + kinematics.push(point); + } + let responses: Vec = if let Formulation::Hybrid { .. } = self.formulation { + self.mean_dilatation(m, &kinematics, small, committed, updated, k.as_mut())? + } else { + let mut out = Vec::with_capacity(kinematics.len()); + for (i, p) in kinematics.iter().enumerate() { + let updated_point = updated.get_mut(i * len..(i + 1) * len).unwrap_or_default(); + let history = point_history(committed, len, i); + out.push(if small { + m.small_strain(&p.strain, p.heat.temperature, history, updated_point, true)? + } else { + m.response(&p.heat.mechanical(&p.f), p.heat.temperature, history, updated_point)? + }); + } + out + }; + let mut out = Vec::with_capacity(rule.len()); + for (p, r) in kinematics.into_iter().zip(responses) { + let (dv, current) = (p.dv, &p.current); let tau = [r.stress[0], r.stress[1], r.stress[2], r.stress[3]]; for (a, g) in current.iter().enumerate() { if let Some(x) = force.get_mut(2 * a) { @@ -614,7 +655,7 @@ impl Plane<'_> { } if let Some(k) = k.as_mut() { let c4: [[f64; 4]; 4] = std::array::from_fn(|a| std::array::from_fn(|b| r.tangent[a][b])); - crate::plane::accumulate(&mut k.data, ¤t, &c4, dv); + crate::plane::accumulate(&mut k.data, current, &c4, dv); if !small { for (a, ga) in current.iter().enumerate() { for (b, gb) in current.iter().enumerate() { @@ -625,11 +666,56 @@ impl Plane<'_> { } } } - out.push(PlanePoint { f, j, x: at.x, response: r }); + out.push(PlanePoint { f: p.f, j: p.j, x: p.x, response: r }); } Ok((Finite { force, tangent: k }, out)) } + /// The mean dilatation responses of a hybrid 2D element (one constant pressure): θ = ∫ J dV₀ + /// / V₀ (hoop stretch included), p = U′(θ), τ = τ_iso + J p I (small strain: s + p I), and + /// the dilatation stiffness U″(θ) g gᵀ / V₀ added to `k`, g = ∫ J ∇ₓ·N dV₀ (small strain + /// ∫ ∇₀·N dV₀), the divergence taking the hoop column N/r too (the constant-pressure case of + /// the solids' three-field form, Simo, Taylor & Pister 1985). + fn mean_dilatation( + &self, + m: &MaterialModel, + points: &[PlaneKinematics], + small: bool, + committed: &[f64], + updated: &mut [f64], + k: Option<&mut Matrix>, + ) -> Result> { + let len = m.history_len(); + let volume: f64 = points.iter().map(|p| p.dv).sum(); + if !(volume > 0.0 && volume.is_finite()) { + return Err(inverted(&format!("2D {}", self.ty.kind.name()), volume)); + } + let theta = points.iter().map(|p| p.dv * p.dilatation).sum::() / volume; + if !(theta > 0.0 && theta.is_finite()) { + return Err(inverted(&format!("2D {}", self.ty.kind.name()), theta)); + } + let (pressure, stiffness) = m.volumetric(theta, small).ok_or_else(|| Error::Invalid("volumetric response".into()))?; + let mut out = Vec::with_capacity(points.len()); + for (i, p) in points.iter().enumerate() { + let (c, u) = (point_history(committed, len, i), updated.get_mut(i * len..(i + 1) * len).unwrap_or_default()); + out.push(if small { + let mut r = m.small_strain(&p.strain, p.heat.temperature, c, u, false)?; + r.stress.iter_mut().take(3).for_each(|s| *s += pressure); + r + } else { + add_pressure(&m.isochoric(&p.heat.mechanical(&p.f), p.heat.temperature, c, u)?, p.j, pressure, 0.0) + }); + } + if let Some(k) = k { + let mut g = vec![0.0; k.n]; + for p in points { + crate::plane::divergence(&mut g, &p.current, p.dv * p.j); + } + crate::plane::add_outer(k, &g, stiffness / volume); + } + Ok(out) + } + /// Internal forces (u, v per node) and (with `tangent`) the tangent stiffness at nodal /// displacements `u`, as for [`Solid::finite_strain`]. pub fn finite_strain(&self, m: &MaterialModel, u: &[[f64; 2]], committed: &[f64], updated: &mut [f64], tangent: bool) -> Result { @@ -1016,6 +1102,63 @@ mod tests { assert!(matches!(stress.finite_strain(&svk, &[[0.0; 2]; 4], &[], &mut [], true), Err(Error::NotSupported(_)))); } + #[test] + fn hybrid_2d_tangents_are_the_derivative_of_the_forces() { + let quad8 = [ + [1.0, 0.0, 0.0], + [2.0, 0.1, 0.0], + [2.2, 1.0, 0.0], + [0.9, 1.1, 0.0], + [1.5, 0.0, 0.0], + [2.1, 0.55, 0.0], + [1.5, 1.1, 0.0], + [0.95, 0.5, 0.0], + ]; + let (shear, bulk) = (300.0, 30000.0); + let neo = MaterialModel::NeoHookean { shear, bulk }; + let mr = MaterialModel::Polynomial { name: "Mooney–Rivlin", c: [100.0, 50.0, 0.0, 0.0], bulk }; + let hybrid = Formulation::Hybrid { pressures: 1 }; + for theory in [Theory::PlaneStrain { thickness: 2.0 }, Theory::Axisymmetric] { + for (kind, x) in [(ElementKind::Quad8, &quad8[..]), (ElementKind::Quad4, &quad8[..4])] { + let plane = Plane::new(kind, x, hybrid, theory).unwrap(); + for m in [&neo, &mr] { + let u: Vec<[f64; 2]> = displacement(x.len(), 0.06).iter().map(|v| [v[0], v[1]]).collect(); + let k = plane.finite_strain(m, &u, &[], &mut [], true).unwrap().tangent.unwrap(); + let scale = k.data.iter().fold(0.0f64, |a, v| a.max(v.abs())); + let h = 1e-6; + for j in 0..2 * x.len() { + let shifted = |s: f64| { + let mut v = u.clone(); + v[j / 2][j % 2] += s; + plane.finite_strain(m, &v, &[], &mut [], false).unwrap().force + }; + let (p, q) = (shifted(h), shifted(-h)); + for i in 0..2 * x.len() { + let numeric = (p[i] - q[i]) / (2.0 * h); + assert!((numeric - k.get(i, j)).abs() < 1e-6 * scale, "{theory:?} {kind:?} {i} {j}: {numeric} vs {}", k.get(i, j)); + } + } + } + // At rest the tangent is the linear hybrid stiffness of the same moduli. + let young = 9.0 * bulk * shear / (3.0 * bulk + shear); + let poisson = (3.0 * bulk - 2.0 * shear) / (2.0 * (3.0 * bulk + shear)); + let d = isotropic(young, poisson).unwrap(); + let k = plane.finite_strain(&neo, &vec![[0.0; 2]; x.len()], &[], &mut [], true).unwrap().tangent.unwrap(); + let linear = plane.stiffness(&d).unwrap(); + let scale = linear.data.iter().fold(0.0f64, |a, v| a.max(v.abs())); + for (a, b) in k.data.iter().zip(&linear.data) { + assert!((a - b).abs() < 1e-9 * scale, "{a} {b}"); + } + } + } + // The linear triangle and plane stress have none; St Venant–Kirchhoff has no split. + assert!(Plane::new(ElementKind::Tri3, &quad8[..3], hybrid, Theory::PlaneStrain { thickness: 1.0 }).is_err()); + assert!(Plane::new(ElementKind::Quad4, &quad8[..4], hybrid, Theory::PlaneStress { thickness: 1.0 }).is_err()); + let svk = MaterialModel::StVenantKirchhoff(Box::new(isotropic(1000.0, 0.3).unwrap())); + let plane = Plane::new(ElementKind::Quad4, &quad8[..4], hybrid, Theory::Axisymmetric).unwrap(); + assert!(matches!(plane.finite_strain(&svk, &[[0.0; 2]; 4], &[], &mut [], true), Err(Error::NotSupported(_)))); + } + /// uᵀ K_σ u = ∫ σ_xx (∂u_y/∂x)² dV = σ_xx V for u_y = x under a uniform σ_xx, in a distorted /// hex20 and a plane strain quad8, and K_σ annihilates rigid translations. #[test] diff --git a/crates/fem/src/plane.rs b/crates/fem/src/plane.rs index 7c1a2c5..5c94b00 100644 --- a/crates/fem/src/plane.rs +++ b/crates/fem/src/plane.rs @@ -35,6 +35,15 @@ //! leaves the plane case unchanged (its mean is zero) and makes the axisymmetric element pass the //! patch test. //! +//! **Hybrid** (mixed displacement–pressure, plane strain and axisymmetric quad4, quad8 and tri6): +//! one constant pressure per element, the mean dilatation method (J. C. Nagtegaal, D. M. Parks, +//! J. R. Rice, *On numerically accurate finite element solutions in the fully plastic range*, +//! CMAME 4(2):153–177, 1974; T. J. R. Hughes, *The Finite Element Method*, Dover 2000, §4.5 for +//! the axisymmetric B-bar): the volumetric strain, hoop part included, is replaced by its volume +//! mean, so K = ∫ Bᵀ D_dev B dV + κ g gᵀ / V with g = ∫ Bᵀ m dV, m = (1, 1, 1, 0). The linear +//! triangle is refused (one pressure per element still locks it), plane stress too (σ_zz = 0 +//! leaves no volume constraint to relax). The finite-strain form is [`crate::finite`]'s. +//! //! On the axis (r = 0) the hoop strain u/r is evaluated only at integration points, which lie //! inside the element (r > 0) for every element with some node off the axis; the solver holds //! the radial displacement of nodes on the axis at zero, as symmetry requires. @@ -165,11 +174,24 @@ pub fn plane_type(kind: ElementKind) -> Result<&'static ElementType> { .ok_or_else(|| Error::NotSupported(format!("{} elements in a 2D section (2D sections take tri3, tri6, quad4 and quad8)", kind.name()))) } +/// The 2D element kinds with a hybrid (constant pressure) formulation. +const HYBRID_PRESSURES: [ElementKind; 3] = [ElementKind::Quad4, ElementKind::Quad8, ElementKind::Tri6]; + /// The formulation a 2D section asks for, if this element type has it. pub fn formulation(kind: ElementKind, f: &fieldcraft_model::Formulation) -> Result { let ty = plane_type(kind)?; if f.hybrid { - return Err(Error::NotSupported("hybrid (mixed displacement–pressure) 2D elements".into())); + if f.incompatible_modes || f.integration == fieldcraft_model::Integration::Reduced { + return Err(Error::NotSupported("hybrid 2D elements with incompatible modes or reduced integration".into())); + } + return if HYBRID_PRESSURES.contains(&kind) { + Ok(Formulation::Hybrid { pressures: 1 }) + } else { + Err(Error::NotSupported(format!( + "hybrid 2D {} elements (one pressure per element still locks the linear triangle; use quad4, quad8 or tri6)", + kind.name() + ))) + }; } if f.incompatible_modes { return if ty.incompatible_modes { @@ -270,11 +292,14 @@ impl<'a> Plane<'a> { Formulation::Full => true, Formulation::Reduced => ty.reduced.is_some(), Formulation::IncompatibleModes => ty.incompatible_modes, - Formulation::Hybrid { .. } => false, + Formulation::Hybrid { pressures } => pressures == 1 && HYBRID_PRESSURES.contains(&kind), }; if !ok { return Err(Error::NotSupported(format!("{formulation:?} formulation of 2D {} elements", kind.name()))); } + if matches!(formulation, Formulation::Hybrid { .. }) && matches!(theory, Theory::PlaneStress { .. }) { + return Err(Error::NotSupported("hybrid plane stress elements (σ_zz = 0 leaves no volume constraint; use full integration)".into())); + } Ok(Plane { ty, xyz, formulation, theory }) } @@ -282,7 +307,8 @@ impl<'a> Plane<'a> { /// modes it condenses out. pub fn dofs(&self) -> ElementDofs { let incompatible_modes = if self.formulation == Formulation::IncompatibleModes { INCOMPATIBLE_DOFS } else { 0 }; - ElementDofs { node: DofSet::PLANE, internal: Internal { incompatible_modes, pressure: 0 } } + let pressure = if let Formulation::Hybrid { pressures } = self.formulation { pressures } else { 0 }; + ElementDofs { node: DofSet::PLANE, internal: Internal { incompatible_modes, pressure } } } /// The stiffness integration rule. @@ -411,6 +437,13 @@ impl<'a> Plane<'a> { /// The stiffness matrix, 2n × 2n, degrees of freedom ordered node by node (u, v). pub fn stiffness(&self, d: &Elasticity) -> Result { + if let Some(dev) = self.deviatoric(d)? { + // K_dev + κ g gᵀ / V. + let (mut k, points, _) = self.full_stiffness(&dev)?; + let mixed = self.mixed(d, &points, None)?; + add_outer(&mut k, &mixed.g, mixed.bulk / mixed.volume); + return Ok(k); + } let d2 = self.theory.elasticity(d)?; let (full, _, modes) = self.full_stiffness(&d2)?; if modes.is_none() { @@ -433,6 +466,34 @@ impl<'a> Plane<'a> { Ok(k) } + /// A hybrid element's deviatoric elasticity in the theory's form; `None` for the other + /// formulations. + fn deviatoric(&self, d: &Elasticity) -> Result> { + if !matches!(self.formulation, Formulation::Hybrid { .. }) { + return Ok(None); + } + Ok(Some(self.theory.elasticity(&crate::hybrid::split(d)?.0)?)) + } + + /// A hybrid element's parts over its stiffness `points`: the bulk modulus, g = ∫ Bᵀ m dV + /// (2n), the volume and, with a thermal strain, ∫ tr ε_th dV. + pub(crate) fn mixed(&self, d: &Elasticity, points: &[(Point, At)], thermal: Option<&ThermalStrain>) -> Result { + let bulk = crate::hybrid::split(d)?.1; + let (mut g, mut volume, mut thermal_volume) = (vec![0.0; 2 * self.ty.nodes], 0.0, 0.0); + for (_, a) in points { + volume += a.dv; + divergence(&mut g, a.g.get(..self.ty.nodes).unwrap_or_default(), a.dv); + if let Some(t) = thermal { + let e = Self::thermal_strain(t, a.shape.values()); + thermal_volume += a.dv * (e[0] + e[1] + e[2]); + } + } + if !(volume > 0.0) { + return Err(self.degenerate()); + } + Ok(Mixed2 { bulk, g, volume, thermal_volume }) + } + /// The consistent mass matrix per displacement component, n × n: m_ij = ∫ ρ Nᵢ Nⱼ dV. pub fn mass(&self, density: f64) -> Result { let n = self.ty.nodes; @@ -518,8 +579,18 @@ impl<'a> Plane<'a> { /// strains the hoop direction too); incompatible modes are condensed as in the stiffness. pub fn thermal_load(&self, d: &Elasticity, thermal: &ThermalStrain) -> Result> { thermal.check(self.ty.nodes)?; - let d2 = self.theory.elasticity(d)?; let n = 2 * self.ty.nodes; + if let Some(dev) = self.deviatoric(d)? { + // ∫ Bᵀ D_dev ε_th dV + κ g (∫ tr ε_th dV) / V. + let points = self.points()?; + let mixed = self.mixed(d, &points, Some(thermal))?; + let (mut fu, _) = Self::thermal_forces(&dev, thermal, &points, None, self.ty.nodes); + for (f, g) in fu.iter_mut().zip(&mixed.g) { + *f += mixed.bulk * g * mixed.thermal_volume / mixed.volume; + } + return Ok(fu.as_chunks::<2>().0.to_vec()); + } + let d2 = self.theory.elasticity(d)?; let mut fu = if self.formulation == Formulation::IncompatibleModes { let (full, points, modes) = self.full_stiffness(&d2)?; let (mut fu, mut fa) = Self::thermal_forces(&d2, thermal, &points, modes.as_ref(), self.ty.nodes); @@ -549,9 +620,21 @@ impl<'a> Plane<'a> { if let Some(t) = thermal { t.check(self.ty.nodes)?; } - let d2 = self.theory.elasticity(d)?; + let deviatoric = self.deviatoric(d)?; + let d2 = match deviatoric { + Some(dev) => dev, + None => self.theory.elasticity(d)?, + }; let (full, points, modes) = if self.formulation == Formulation::IncompatibleModes { self.full_stiffness(&d2)? } else { (Matrix::zeros(0), self.points()?, None) }; + // Hybrid: the element's pressure κ (gᵀ u − ∫ tr ε_th dV) / V. + let pressure = if deviatoric.is_some() { + let m = self.mixed(d, &points, thermal)?; + let gu: f64 = m.g.iter().zip(u.iter().flatten()).map(|(g, x)| g * x).sum(); + m.bulk * (gu - m.thermal_volume) / m.volume + } else { + 0.0 + }; // Incompatible modes: K_aa α = f_a − K_au u (f_a the thermal load on the modes). let alpha = if modes.is_some() { let n = 2 * self.ty.nodes; @@ -595,10 +678,11 @@ impl<'a> Plane<'a> { }; } // Plane stress: D̃ has no zz row or column, so the εzz just recovered adds nothing. - let s: [f64; 4] = match eth { + let mut s: [f64; 4] = match eth { Some(t) => std::array::from_fn(|r| (0..4).map(|c| d2[r][c] * (e[c] - t[c])).sum()), None => std::array::from_fn(|r| (0..4).map(|c| d2[r][c] * e[c]).sum()), }; + s.iter_mut().take(3).for_each(|v| *v += pressure); out.push(PointStress { xyz: a.x, strain: [e[0], e[1], e[2], e[3], 0.0, 0.0], stress: [s[0], s[1], s[2], s[3], 0.0, 0.0] }); } Ok(out) @@ -614,6 +698,32 @@ impl<'a> Plane<'a> { } } +/// The parts of a hybrid 2D element ([`Plane::mixed`]). +pub(crate) struct Mixed2 { + bulk: f64, + pub(crate) g: Vec, + volume: f64, + thermal_volume: f64, +} + +/// g += w Bᵀ m over strain columns `cols` (u takes ∂N/∂x and the hoop N/r, v takes ∂N/∂y): the +/// divergence of each nodal motion, which the hybrid elements' volume constraint acts on. +pub(crate) fn divergence(g: &mut [f64], cols: &[Column], w: f64) { + for (pair, col) in g.as_chunks_mut::<2>().0.iter_mut().zip(cols) { + pair[0] += w * (col[0] + col[2]); + pair[1] += w * col[1]; + } +} + +/// k += c g gᵀ. +pub(crate) fn add_outer(k: &mut Matrix, g: &[f64], c: f64) { + for (i, gi) in g.iter().enumerate() { + for (j, gj) in g.iter().enumerate() { + k.add(i, j, c * gi * gj); + } + } +} + /// The `size` × `size` block of `m` starting at (`row`, `col`). fn sub_block(m: &Matrix, row: usize, col: usize, size: usize) -> Matrix { let mut out = Matrix::zeros(size); @@ -661,6 +771,53 @@ mod tests { const THEORIES: [Theory; 3] = [Theory::PlaneStress { thickness: 0.7 }, Theory::PlaneStrain { thickness: 1.0 }, Theory::Axisymmetric]; + #[test] + fn hybrid_elements_reproduce_a_uniform_strain_and_do_not_lock() { + // u = a x, v = b y: a uniform strain (axisymmetric: the hoop strain a too), whose + // stress the mean dilatation reproduces exactly at every point. + let d = isotropic(1000.0, 0.4999).unwrap(); + let (a, b) = (0.002, -0.003); + for theory in [Theory::PlaneStrain { thickness: 1.5 }, Theory::Axisymmetric] { + for kind in [ElementKind::Quad4, ElementKind::Quad8, ElementKind::Tri6] { + let x = element(kind); + let plane = Plane::new(kind, &x, Formulation::Hybrid { pressures: 1 }, theory).unwrap(); + let u: Vec<[f64; 2]> = x.iter().map(|p| [a * p[0], b * p[1]]).collect(); + let hoop = if theory == Theory::Axisymmetric { a } else { 0.0 }; + let e = [a, b, hoop, 0.0, 0.0, 0.0]; + let exact = crate::material::stress(&d, &e); + for p in plane.stresses(&d, &u).unwrap() { + for c in 0..4 { + assert!((p.stress[c] - exact[c]).abs() < 1e-9 * exact[0].abs().max(1.0), "{theory:?} {kind:?} {c}: {:?} {exact:?}", p.stress); + } + } + // K u is the forces of that stress: the stiffness agrees with the stresses. + let k = plane.stiffness(&d).unwrap(); + let flat: Vec = u.iter().flatten().copied().collect(); + let full = Plane::new(kind, &x, Formulation::Full, theory).unwrap().stiffness(&d).unwrap(); + let (ku, kf): (Vec, Vec) = (0..flat.len()) + .map(|i| (0..flat.len()).fold((0.0, 0.0), |(s, t), j| (s + k.get(i, j) * flat[j], t + full.get(i, j) * flat[j]))) + .unzip(); + for (p, q) in ku.iter().zip(&kf) { + assert!((p - q).abs() < 1e-8 * kf.iter().fold(0.0f64, |m, v| m.max(v.abs())), "{p} {q}"); + } + } + } + // Nearly incompressible: the hybrid element's volume is one constraint, so a motion that + // keeps the element's volume (orthogonal to g) costs only shear, where the full element's + // four point constraints lock it. + let x = element(ElementKind::Quad4); + let hybrid = Plane::new(ElementKind::Quad4, &x, Formulation::Hybrid { pressures: 1 }, Theory::PlaneStrain { thickness: 1.0 }).unwrap(); + let g = hybrid.mixed(&d, &hybrid.points().unwrap(), None).unwrap().g; + let raw: Vec = (0..8).map(|i| ((i * 7 % 5) as f64 - 2.0) * 0.01).collect(); + let along = raw.iter().zip(&g).map(|(a, b)| a * b).sum::() / g.iter().map(|v| v * v).sum::(); + let shear: Vec = raw.iter().zip(&g).map(|(a, b)| a - along * b).collect(); + let energy = |f: Formulation| { + let k = Plane::new(ElementKind::Quad4, &x, f, Theory::PlaneStrain { thickness: 1.0 }).unwrap().stiffness(&d).unwrap(); + (0..8).map(|i| (0..8).map(|j| shear[i] * k.get(i, j) * shear[j]).sum::()).sum::() + }; + assert!(energy(Formulation::Hybrid { pressures: 1 }) < 0.01 * energy(Formulation::Full)); + } + #[test] fn rigid_motion_strains_nothing_and_matrices_are_symmetric() { let d = isotropic(1000.0, 0.3).unwrap(); @@ -783,6 +940,9 @@ mod tests { let on_axis = sq.map(|p| [p[0] - 1.0, p[1], 0.0]); let e = Plane::new(ElementKind::Quad4, &on_axis, Formulation::IncompatibleModes, Theory::Axisymmetric).unwrap(); assert!(e.stiffness(&d).is_ok()); - assert!(formulation(ElementKind::Quad4, &fieldcraft_model::Formulation { hybrid: true, ..Default::default() }).is_err()); + let hybrid = fieldcraft_model::Formulation { hybrid: true, ..Default::default() }; + assert_eq!(formulation(ElementKind::Quad4, &hybrid).unwrap(), Formulation::Hybrid { pressures: 1 }); + assert!(formulation(ElementKind::Tri3, &hybrid).is_err()); + assert!(formulation(ElementKind::Quad4, &fieldcraft_model::Formulation { incompatible_modes: true, ..hybrid }).is_err()); } } diff --git a/crates/model/src/sections.rs b/crates/model/src/sections.rs index 58a89d4..e779c60 100644 --- a/crates/model/src/sections.rs +++ b/crates/model/src/sections.rs @@ -278,8 +278,9 @@ impl Formulation { if !kind.is_continuum() && self.incompatible_modes { return Err(Error::Invalid(format!("incompatible modes are for solid and 2D sections, not {}", kind.name()))); } - if kind.dimension() != 3 && self.hybrid { - return Err(Error::Invalid(format!("hybrid formulations are for solid sections, not {}", kind.name()))); + let volume_constrained = kind.dimension() == 3 || matches!(kind, SectionKind::PlaneStrain { .. } | SectionKind::Axisymmetric {}); + if self.hybrid && !volume_constrained { + return Err(Error::Invalid(format!("hybrid formulations are for solid, plane strain and axisymmetric sections, not {}", kind.name()))); } Ok(()) } @@ -391,7 +392,10 @@ mod tests { let im = Formulation { incompatible_modes: true, ..Formulation::default() }; im.validate_for(&stress).unwrap(); im.validate_for(&axi).unwrap(); - assert!(Formulation { hybrid: true, ..Formulation::default() }.validate_for(&strain).is_err()); + let hybrid = Formulation { hybrid: true, ..Formulation::default() }; + hybrid.validate_for(&strain).unwrap(); + hybrid.validate_for(&axi).unwrap(); + assert!(hybrid.validate_for(&stress).is_err()); assert!(im.validate_for(&SectionKind::Shell { thickness: 1.0, offset: 0.0, points: 5 }).is_err()); } } diff --git a/docs/fem.md b/docs/fem.md index d026a15..c7ac898 100644 --- a/docs/fem.md +++ b/docs/fem.md @@ -146,7 +146,10 @@ pressure field. The choice of pressure space follows the inf-sup analysis (docs checkerboard pressures; linear tetrahedra and wedges lock even so and are refused. A thermal strain splits the same way: ∫ Bᵀ D_dev ε_th on the nodes and, through the pressure, K G M⁻¹ ∫ N_p tr ε_th, so a heated nearly incompressible solid takes its bulk stress −3KαΔT -when confined. A continuous +when confined. Plane strain and axisymmetric quad4, quad8 and tri6 sections take one constant +pressure (the mean dilatation, the hoop strain in the divergence): K_dev + κ g gᵀ / V with +g = ∫ Bᵀ m dV (`crates/fem/src/plane.rs`); the linear triangle locks and plane stress has no +volume constraint, so both are refused. A continuous nodal pressure (`DofKind::P`, Taylor–Hood) needs a symmetric indefinite factorization and is planned. @@ -251,6 +254,7 @@ differences of its forces. | Element | Nonlinear steps | |---|---| | tet4, tet10, wedge6, hex8, hex20 (full), hex20 (reduced) | total Lagrangian | +| hybrid quad4, quad8, tri6 (plane strain, axisymmetric) | mean dilatation: θ = ∫ J dV₀ / V₀ with the hoop stretch, p = U′(θ), stiffness U″ g gᵀ / V₀ | | hybrid hex8, hex20, tet10 | Simo–Taylor–Pister three-field form on the linear hybrid elements' pressure fields (θ the projection of J, p the projection of U′(θ)); needs the neo-Hookean option (St Venant–Kirchhoff has no volumetric split) | | hex8 and quad4 incompatible modes | refused: their finite-strain form (enhanced assumed strain, Simo–Armero 1992) needs the modes as unknowns carried between iterations and hourglasses in large compression; use full integration, hybrid or quadratic elements | | plane strain, axisymmetric (tri3, tri6, quad4, quad8) | total Lagrangian with F₃₃ = 1 or r/R (hoop column N/r, geometric term τ_θθ N_a N_b/r²) | diff --git a/tests/oracle/nl-mooney-rivlin-tube-axisymmetric/README.md b/tests/oracle/nl-mooney-rivlin-tube-axisymmetric/README.md new file mode 100644 index 0000000..e662b8f --- /dev/null +++ b/tests/oracle/nl-mooney-rivlin-tube-axisymmetric/README.md @@ -0,0 +1,8 @@ +# Inflation of a stretched Mooney–Rivlin tube, axisymmetric + +The tube of `nl-mooney-rivlin-tube` (A = 10, B = 20, one unit long; E = 3, ν = 0.4995, +C₀₁ = 0.2) as an axisymmetric model of 12 × 2 hybrid quad4 (one constant pressure each, the mean +dilatation with the hoop strain): the end y = 0 on a symmetry line, the end y = 1 stretched to +λ_z = 1.5 while the follower pressure p = 0.540850 inflates the bore. The closed form puts the +bore at a = 15 and the outer wall at b = 20.615528 (see `nl-mooney-rivlin-tube`); this case checks +the nearly incompressible 2D hybrid elements at large strain against it. diff --git a/tests/oracle/nl-mooney-rivlin-tube-axisymmetric/expected.json b/tests/oracle/nl-mooney-rivlin-tube-axisymmetric/expected.json new file mode 100644 index 0000000..7163415 --- /dev/null +++ b/tests/oracle/nl-mooney-rivlin-tube-axisymmetric/expected.json @@ -0,0 +1,53 @@ +{ + "title": "Inflation of a stretched Mooney–Rivlin tube (axisymmetric hybrid quad4)", + "kind": "theory", + "source": "Finite inflation and extension of an incompressible Mooney–Rivlin tube (R. S. Rivlin, Large elastic deformations of isotropic materials VI, Phil. Trans. R. Soc. A 242:173–195, 1949; R. W. Ogden, Non-Linear Elastic Deformations, Dover 1997, §5.3): at axial stretch λ_z, r² = a² + (R² − A²)/λ_z, σ_θ − σ_r = 2 (λ_θ² − λ_r²)(C₁₀ + C₀₁ λ_z²) with λ_θ = r/R, λ_r = 1/(λ_θ λ_z), integrated through the wall: p = ∫ (σ_θ − σ_r) dr/r (Simpson, 2·10⁵ intervals)", + "quantities": [ + { + "id": "ua", + "quantity": "radial displacement of the bore (10, 0, 0) at the pressure that inflates it to a = 15 at λ_z = 1.5", + "measure": { + "command": "solve.probe", + "params": { + "step": "Inflate", + "at": [ + [ + 10.0, + 0, + 0.0 + ] + ] + }, + "pointer": "/nodes/0/displacement/0" + }, + "value": 5.0, + "tolerance": { + "relative": 0.01 + }, + "reason": "16 × 6 hybrid hex20 around the quarter; ν = 0.4995 stands in for incompressibility (a volume change of order p/κ ≈ 10⁻³ shifts the radii by about 0.1 %), and the tube is stretched axially by its ends, which a nearly incompressible solid takes with larger pressures than the plane-strain neo-Hookean tube" + }, + { + "id": "ub", + "quantity": "radial displacement of the outer wall (20, 0, 0): b − B, a small difference of radii", + "measure": { + "command": "solve.probe", + "params": { + "step": "Inflate", + "at": [ + [ + 20.0, + 0, + 0.0 + ] + ] + }, + "pointer": "/nodes/0/displacement/0" + }, + "value": 0.6155281280883038, + "tolerance": { + "absolute": 0.04 + }, + "reason": "0.2 % of the outer radius b = 20.616: the same volume change of order 10⁻³ moves b by about 0.03, which is 4.5 % of b − B; the incompressible kinematics r² = a² + (R² − A²)/λ_z amplify it there" + } + ] +} \ No newline at end of file diff --git a/tests/oracle/nl-mooney-rivlin-tube-axisymmetric/fixture.json b/tests/oracle/nl-mooney-rivlin-tube-axisymmetric/fixture.json new file mode 100644 index 0000000..b578e4d --- /dev/null +++ b/tests/oracle/nl-mooney-rivlin-tube-axisymmetric/fixture.json @@ -0,0 +1,13 @@ +{ + "note": "FieldCraft's accepted results for this case, written by `cargo xtask oracle --bless`. A result outside its band fails CI.", + "values": { + "ua": { + "band": 5.03287906011305e-9, + "value": 5.03287906011305 + }, + "ub": { + "band": 0.000040000644728979273, + "value": 0.6447289792693816 + } + } +} diff --git a/tests/oracle/nl-mooney-rivlin-tube-axisymmetric/recipe.json b/tests/oracle/nl-mooney-rivlin-tube-axisymmetric/recipe.json new file mode 100644 index 0000000..387fde4 --- /dev/null +++ b/tests/oracle/nl-mooney-rivlin-tube-axisymmetric/recipe.json @@ -0,0 +1,14 @@ +[ + {"command": "mesh.block", "params": {"body": "Pipe", "kind": "quad4", "corners": [[10.0, 0, 0], [20.0, 0, 0], [20.0, 1, 0], [10.0, 1, 0]], "divisions": [12, 2], "edges": [{"corners": [1, 4], "name": "bore"}, {"corners": [1, 2], "name": "bottom"}, {"corners": [3, 4], "name": "top"}]}}, + {"command": "model.set.create", "params": {"name": "all", "geometry": [{"body": "Pipe", "face": "area"}]}}, + {"command": "model.set.create", "params": {"name": "bore", "geometry": [{"body": "Pipe", "edge": "bore"}]}}, + {"command": "model.set.create", "params": {"name": "bottom", "geometry": [{"body": "Pipe", "edge": "bottom"}]}}, + {"command": "model.set.create", "params": {"name": "top", "geometry": [{"body": "Pipe", "edge": "top"}]}}, + {"command": "model.material.create", "params": {"name": "Rubber", "elastic": {"type": "isotropic", "young": 3.0, "poisson": 0.4995}, "hyperelastic": {"type": "mooney_rivlin", "c01": 0.2}}}, + {"command": "model.section.create", "params": {"name": "pipe", "set": "all", "material": "Rubber", "kind": {"type": "axisymmetric"}, "formulation": {"hybrid": true}}}, + {"command": "model.boundary_condition.create", "params": {"name": "sb", "type": "symmetry", "target": {"set": "bottom"}, "plane": "y"}}, + {"command": "model.boundary_condition.create", "params": {"name": "stretch", "type": "displacement", "target": {"set": "top"}, "values": [null, 0.5, null, null, null, null]}}, + {"command": "model.load.create", "params": {"name": "p", "type": "pressure", "set": "bore", "magnitude": 0.5408503415684837}}, + {"command": "model.step.create", "params": {"name": "Inflate", "procedure": {"type": "nonlinear_static"}, "boundary_conditions": ["sb", "stretch"], "loads": ["p"]}}, + {"command": "solve.nonlinear", "params": {"step": "Inflate"}, "id": "solve"} +]