From 3f9e9e30d868711da2e042593b7d4a442faf5be3 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Wed, 14 Oct 2026 07:04:28 -0400 Subject: [PATCH] fem, solver: thermal strain in nonlinear static steps Temperatures now strain the material in nonlinear steps instead of being refused. At finite strain the deformation gradient splits multiplicatively, F = F_m F_theta with F_theta = diag(1 + alpha_i dT) at each point: the material sees F_m and the loops run on the thermally expanded volume J_theta dV0 (mixed elements project det F_m), so the tangent stays consistent. At small strain the split is additive, and elements on their linear stiffness take the equivalent nodal forces as a dead load. The step's temperatures are its temperature loads, or the steady conduction of its temperature conditions and heat loads (one-way coupling, as the linear coupled step); both at once are refused. They ramp with the load factor from each material's reference temperature. Johnson-Cook plasticity softens at the point's temperature. Verified: free expansion and a fully constrained neo-Hookean block at 20 % thermal stretch, full and hybrid (exact); a bar held at its ends heated past yield (exact); the linear limit and the conducted temperatures against the linear solvers; tangents by central differences with a temperature field. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm --- crates/fem/src/finite.rs | 274 ++++++++++++++++++----- crates/fem/src/plasticity.rs | 5 + crates/solver/src/analysis.rs | 6 +- crates/solver/src/heat.rs | 7 + crates/solver/src/nonlinear.rs | 96 ++++++-- crates/solver/src/system.rs | 42 ++-- crates/solver/tests/nonlinear.rs | 9 +- crates/solver/tests/nonlinear_thermal.rs | 147 ++++++++++++ 8 files changed, 492 insertions(+), 94 deletions(-) create mode 100644 crates/solver/tests/nonlinear_thermal.rs diff --git a/crates/fem/src/finite.rs b/crates/fem/src/finite.rs index d65fe10..68de9fe 100644 --- a/crates/fem/src/finite.rs +++ b/crates/fem/src/finite.rs @@ -29,6 +29,17 @@ //! elements of an elastic material are refused in nonlinear steps; J2 plasticity has it exactly //! (the return is deviatoric, the volume Hencky's K/2 (ln J)², [`crate::elastoplastic`]). //! +//! **Thermal strain** ([`Temperature`]): at finite strain the multiplicative split F = F_m F_θ +//! with the thermal stretches F_θ = diag(1 + αᵢ (T − T_ref)) at each point (T interpolated from +//! the nodes; V. A. Lubarda, *Constitutive theories based on the multiplicative decomposition of +//! deformation gradient: thermoelasticity, elastoplasticity, and biomechanics*, Appl. Mech. Rev. +//! 57(2):95–108, 2004): the material sees +//! F_m, the stress per reference volume is J_θ τ_m(F_m), and since F_θ is fixed over an +//! iteration the spatial tangent is J_θ c_m, so the loops run unchanged on the thermally +//! expanded volume J_θ dV₀ (the mixed elements project det F_m on it). The linear stretch +//! 1 + α ΔT is the engineering definition of α, so a free body grows by exactly α ΔT. At small +//! strain the split is additive, ε_m = ε − α ΔT. Reported strains are the total ones. +//! //! **Small strain** ([`Kinematics::Small`], a step without nonlinear geometry with an inelastic //! material): the same loops on the linear kinematics, ε = sym ∇u with the reference gradients, //! f = ∫ Bᵀ σ dV₀ and K = ∫ Bᵀ 𝔻 B dV₀ with the material's consistent tangent 𝔻, and no @@ -73,14 +84,76 @@ pub enum Kinematics { Small, } +/// The temperatures of an element in a nonlinear step: the linear expansion coefficients along +/// X, Y, Z (zero for a material that does not expand), the reference temperature at which the +/// thermal strain is zero, and the temperature at each node. +#[derive(Debug, Clone, PartialEq)] +pub struct Temperature { + pub alpha: [f64; 3], + pub reference: f64, + pub nodal: Vec, +} + +impl Temperature { + fn check(&self, nodes: usize) -> Result<()> { + if self.nodal.len() != nodes { + return Err(Error::Invalid(format!("{} temperatures for a {nodes}-node element", self.nodal.len()))); + } + if self.alpha.iter().chain(&self.nodal).chain([&self.reference]).any(|v| !v.is_finite()) { + return Err(Error::Invalid("temperatures and expansion coefficients must be finite".into())); + } + Ok(()) + } +} + +/// The temperature at a point and its thermal stretches ϑᵢ = 1 + αᵢ (T − T_ref). +#[derive(Debug, Clone, Copy)] +struct Heat { + temperature: Option, + stretch: [f64; 3], +} + +impl Heat { + /// At a point of shape function values `n`. + fn at(t: Option<&Temperature>, n: &[f64]) -> Result { + let Some(t) = t else { return Ok(Heat { temperature: None, stretch: [1.0; 3] }) }; + let temperature: f64 = t.nodal.iter().zip(n).map(|(a, b)| a * b).sum(); + let stretch = t.alpha.map(|a| 1.0 + a * (temperature - t.reference)); + if !stretch.iter().all(|v| *v > 0.0 && v.is_finite()) { + return Err(Error::Degenerate(format!("the thermal stretches {stretch:?} at T = {temperature} are not positive"))); + } + Ok(Heat { temperature: Some(temperature), stretch }) + } + + /// J_θ = ϑ₁ ϑ₂ ϑ₃. + fn volume(&self) -> f64 { + self.stretch.iter().product() + } + + /// The mechanical part F_m = F F_θ⁻¹ of a deformation gradient. + fn mechanical(&self, f: &Mat3) -> Mat3 { + std::array::from_fn(|i| std::array::from_fn(|j| f[i][j] / self.stretch[j])) + } + + /// The mechanical part ε − α (T − T_ref) of a small strain. + fn mechanical_strain(&self, e: &[f64; 6]) -> [f64; 6] { + std::array::from_fn(|k| if k < 3 { e[k] - (self.stretch[k] - 1.0) } else { e[k] }) + } +} + /// What one integration point gives the element. struct PointKinematics { + /// The point's share of the reference volume (of the thermally expanded volume J_θ dV₀ at + /// finite strain). weight: f64, + /// The deformation gradient, and its mechanical part F_m = F F_θ⁻¹ the material sees. + total: Mat3, f: Mat3, - /// det F, or the dilatation 1 + tr ε of small strain (what a mixed element projects). + /// det F_m, or the dilatation 1 + tr ε_m of small strain (what a mixed element projects). j: f64, - /// The small strain ε = sym ∇u (engineering shears); zero at finite strain. + /// The mechanical small strain ε − ε_th (engineering shears); zero at finite strain. strain: [f64; 6], + temperature: Option, /// ∇ₓN of every node (∇₀N at small strain). gradients: Vec<[f64; 3]>, /// Reference position. @@ -100,8 +173,8 @@ fn small_strain(f: &Mat3) -> [f64; 6] { /// The response of one point's material with its history slices. fn point_response(m: &MaterialModel, kin: Kinematics, p: &PointKinematics, committed: &[f64], updated: &mut [f64]) -> Result { match kin { - Kinematics::Finite => m.response(&p.f, None, committed, updated), - Kinematics::Small => m.small_strain(&p.strain, None, committed, updated, true), + Kinematics::Finite => m.response(&p.f, p.temperature, committed, updated), + Kinematics::Small => m.small_strain(&p.strain, p.temperature, committed, updated, true), } } @@ -182,43 +255,52 @@ impl Solid<'_> { /// F, J and the gradients at every stiffness point: the current ones at finite strain, the /// reference ones at small strain. - fn kinematics(&self, u: &[[f64; 3]], kin: Kinematics) -> Result> { + fn kinematics(&self, u: &[[f64; 3]], kin: Kinematics, thermal: Option<&Temperature>) -> Result> { let n = self.ty.nodes; + if let Some(t) = thermal { + t.check(n)?; + } let mut out = Vec::new(); for p in self.rule().points() { let g = self.gradients(p.xi)?; let reference = g.dx.get(..n).unwrap_or_default(); - let f = crate::constitutive::deformation_gradient(u, reference); + let total = crate::constitutive::deformation_gradient(u, reference); let (weight, x0) = (g.det * p.weight, self.position_of(g.shape.values())); + let heat = Heat::at(thermal, g.shape.values())?; if kin == Kinematics::Small { - let strain = small_strain(&f); + let strain = heat.mechanical_strain(&small_strain(&total)); let j = 1.0 + strain[0] + strain[1] + strain[2]; - out.push(PointKinematics { weight, f, j, strain, gradients: reference.to_vec(), x0 }); + let gradients = reference.to_vec(); + out.push(PointKinematics { weight, total, f: total, j, strain, temperature: heat.temperature, gradients, x0 }); continue; } - let (inv, j) = inverse3(&f).ok_or_else(|| inverted(self.ty.kind.name(), 0.0))?; - if !(j > 0.0 && j.is_finite()) { - return Err(inverted(self.ty.kind.name(), j)); + let (inv, det) = inverse3(&total).ok_or_else(|| inverted(self.ty.kind.name(), 0.0))?; + let j = det / heat.volume(); + if !(det > 0.0 && j.is_finite()) { + return Err(inverted(self.ty.kind.name(), det)); } let gradients = reference.iter().map(|d| std::array::from_fn(|i| (0..3).map(|jj| d[jj] * inv[jj][i]).sum())).collect(); - out.push(PointKinematics { weight, f, j, strain: [0.0; 6], gradients, x0 }); + let f = heat.mechanical(&total); + out.push(PointKinematics { weight: weight * heat.volume(), total, f, j, strain: [0.0; 6], temperature: heat.temperature, gradients, x0 }); } Ok(out) } /// The effective response (Kirchhoff stress and tangent, or the stress of small strain) at /// every point, the forces, and the tangent when `tangent`. + #[allow(clippy::too_many_arguments)] fn nonlinear_points( &self, m: &MaterialModel, u: &[[f64; 3]], kin: Kinematics, + thermal: Option<&Temperature>, committed: &[f64], updated: &mut [f64], tangent: bool, ) -> Result<(Finite, Vec<(PointKinematics, Response)>)> { self.check_nonlinear(m, u, kin)?; - let points = self.kinematics(u, kin)?; + let points = self.kinematics(u, kin, thermal)?; let len = check_history(m, points.len(), committed, updated)?; let ndof = 3 * self.ty.nodes; let mut force = vec![0.0; ndof]; @@ -287,11 +369,11 @@ impl Solid<'_> { let pr: f64 = (0..np).map(|a| nb[a] * pressure[a]).sum(); 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, None, c, u, false)?; + let mut r = m.small_strain(&p.strain, p.temperature, c, u, false)?; r.stress.iter_mut().take(3).for_each(|s| *s += pr); r } else { - add_pressure(&m.isochoric(&p.f, None, c, u)?, p.j, pr, 0.0) + add_pressure(&m.isochoric(&p.f, p.temperature, c, u)?, p.j, pr, 0.0) }); } if let Some(k) = k { @@ -331,48 +413,57 @@ impl Solid<'_> { /// with the material history `committed` at the start of the increment (points × the /// material's history length) written, updated, to `updated`. pub fn finite_strain(&self, m: &MaterialModel, u: &[[f64; 3]], committed: &[f64], updated: &mut [f64], tangent: bool) -> Result { - Ok(self.nonlinear_points(m, u, Kinematics::Finite, committed, updated, tangent)?.0) + Ok(self.nonlinear_points(m, u, Kinematics::Finite, None, committed, updated, tangent)?.0) } - /// As [`Self::finite_strain`] with either kinematics; with small strain an inelastic - /// material on the linear kinematics, f = ∫ Bᵀ σ dV and K = ∫ Bᵀ 𝔻 B dV with the - /// material's consistent tangent 𝔻. + /// As [`Self::finite_strain`] with either kinematics and the element's temperatures when + /// it has some; with small strain an inelastic material on the linear kinematics, + /// f = ∫ Bᵀ σ dV and K = ∫ Bᵀ 𝔻 B dV with the material's consistent tangent 𝔻. + #[allow(clippy::too_many_arguments)] pub fn nonlinear( &self, m: &MaterialModel, u: &[[f64; 3]], kin: Kinematics, + thermal: Option<&Temperature>, committed: &[f64], updated: &mut [f64], tangent: bool, ) -> Result { - Ok(self.nonlinear_points(m, u, kin, committed, updated, tangent)?.0) + Ok(self.nonlinear_points(m, u, kin, thermal, committed, updated, tangent)?.0) } - /// The Cauchy stress and the strain (logarithmic at finite strain, ε at small strain) at - /// every integration point (at its reference position). - pub fn nonlinear_stresses(&self, m: &MaterialModel, u: &[[f64; 3]], kin: Kinematics, committed: &[f64]) -> Result> { + /// The Cauchy stress and the total strain (logarithmic at finite strain, ε at small strain) + /// at every integration point (at its reference position). + pub fn nonlinear_stresses( + &self, + m: &MaterialModel, + u: &[[f64; 3]], + kin: Kinematics, + thermal: Option<&Temperature>, + committed: &[f64], + ) -> Result> { let mut scratch = committed.to_vec(); - let (_, points) = self.nonlinear_points(m, u, kin, committed, &mut scratch, false)?; + let (_, points) = self.nonlinear_points(m, u, kin, thermal, committed, &mut scratch, false)?; points .iter() .map(|(p, r)| { if kin == Kinematics::Small { - return Ok(PointStress { xyz: p.x0, strain: p.strain, stress: r.stress }); + return Ok(PointStress { xyz: p.x0, strain: small_strain(&p.total), stress: r.stress }); } - let strain = logarithmic_strain(&p.f).ok_or_else(|| inverted(self.ty.kind.name(), p.j))?; + let strain = logarithmic_strain(&p.total).ok_or_else(|| inverted(self.ty.kind.name(), p.j))?; Ok(PointStress { xyz: p.x0, strain, stress: r.stress.map(|v| v / p.j) }) }) .collect() } } -/// One integration point of a 2D element in a nonlinear step. +/// 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 { f: Mat3, j: f64, x: [f64; 3], - strain: [f64; 6], response: Response, } @@ -399,17 +490,22 @@ impl Plane<'_> { /// Per point: F, J, the reference position, the small strain and the response; the forces /// and the tangent. + #[allow(clippy::too_many_arguments)] fn nonlinear_points( &self, m: &MaterialModel, u: &[[f64; 2]], kin: Kinematics, + thermal: Option<&Temperature>, committed: &[f64], updated: &mut [f64], tangent: bool, ) -> Result<(Finite, Vec)> { self.check_nonlinear(u, kin)?; let n = self.ty.nodes; + if let Some(t) = thermal { + t.check(n)?; + } let rule = self.rule().points(); let len = check_history(m, rule.len(), committed, updated)?; let ndof = 2 * n; @@ -433,10 +529,11 @@ impl Plane<'_> { } } let updated_point = updated.get_mut(i * len..(i + 1) * len).unwrap_or_default(); - let (j, strain, current, r) = if small { - let strain = small_strain(&f); - let r = m.small_strain(&strain, None, point_history(committed, len, i), updated_point, true)?; - (1.0 + strain[0] + strain[1] + strain[2], strain, cols.to_vec(), r) + let heat = Heat::at(thermal, at.shape.values())?; + let (j, dv, current, r) = 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) } else { let det2 = f[0][0] * f[1][1] - f[0][1] * f[1][0]; let j = det2 * f[2][2]; @@ -446,32 +543,32 @@ 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(&f, None, point_history(committed, len, i), updated_point)?; - (j, [0.0; 6], current, r) + 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 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) { - *x += at.dv * (g[0] * tau[0] + g[2] * tau[2] + g[1] * tau[3]); + *x += dv * (g[0] * tau[0] + g[2] * tau[2] + g[1] * tau[3]); } if let Some(x) = force.get_mut(2 * a + 1) { - *x += at.dv * (g[1] * tau[1] + g[0] * tau[3]); + *x += dv * (g[1] * tau[1] + g[0] * tau[3]); } } 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, at.dv); + crate::plane::accumulate(&mut k.data, ¤t, &c4, dv); if !small { for (a, ga) in current.iter().enumerate() { for (b, gb) in current.iter().enumerate() { let s = ga[0] * gb[0] * tau[0] + ga[1] * gb[1] * tau[1] + (ga[0] * gb[1] + ga[1] * gb[0]) * tau[3]; - k.add(2 * a, 2 * b, at.dv * (s + tau[2] * ga[2] * gb[2])); - k.add(2 * a + 1, 2 * b + 1, at.dv * s); + k.add(2 * a, 2 * b, dv * (s + tau[2] * ga[2] * gb[2])); + k.add(2 * a + 1, 2 * b + 1, dv * s); } } } } - out.push(PlanePoint { f, j, x: at.x, strain, response: r }); + out.push(PlanePoint { f, j, x: at.x, response: r }); } Ok((Finite { force, tangent: k }, out)) } @@ -479,32 +576,41 @@ impl Plane<'_> { /// 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 { - Ok(self.nonlinear_points(m, u, Kinematics::Finite, committed, updated, tangent)?.0) + Ok(self.nonlinear_points(m, u, Kinematics::Finite, None, committed, updated, tangent)?.0) } /// As [`Solid::nonlinear`], on the 2D element. + #[allow(clippy::too_many_arguments)] pub fn nonlinear( &self, m: &MaterialModel, u: &[[f64; 2]], kin: Kinematics, + thermal: Option<&Temperature>, committed: &[f64], updated: &mut [f64], tangent: bool, ) -> Result { - Ok(self.nonlinear_points(m, u, kin, committed, updated, tangent)?.0) + Ok(self.nonlinear_points(m, u, kin, thermal, committed, updated, tangent)?.0) } - /// The Cauchy stress and the strain at every integration point, as the solid elements' - /// 6-vectors (axisymmetric: rr, axial, hoop, r-axial). - pub fn nonlinear_stresses(&self, m: &MaterialModel, u: &[[f64; 2]], kin: Kinematics, committed: &[f64]) -> Result> { + /// The Cauchy stress and the total strain at every integration point, as the solid + /// elements' 6-vectors (axisymmetric: rr, axial, hoop, r-axial). + pub fn nonlinear_stresses( + &self, + m: &MaterialModel, + u: &[[f64; 2]], + kin: Kinematics, + thermal: Option<&Temperature>, + committed: &[f64], + ) -> Result> { let mut scratch = committed.to_vec(); - let (_, points) = self.nonlinear_points(m, u, kin, committed, &mut scratch, false)?; + let (_, points) = self.nonlinear_points(m, u, kin, thermal, committed, &mut scratch, false)?; points .iter() .map(|p| { if kin == Kinematics::Small { - return Ok(PointStress { xyz: p.x, strain: p.strain, stress: p.response.stress }); + return Ok(PointStress { xyz: p.x, strain: small_strain(&p.f), stress: p.response.stress }); } let strain = logarithmic_strain(&p.f).ok_or_else(|| inverted(&format!("2D {}", self.ty.kind.name()), p.j))?; Ok(PointStress { xyz: p.x, strain, stress: p.response.stress.map(|v| v / p.j) }) @@ -552,7 +658,7 @@ mod tests { /// The same from a committed history, with either kinematics. fn check_tangent_from(solid: &Solid, m: &MaterialModel, u: &[[f64; 3]], kin: Kinematics, committed: &[f64], tol: f64) { let mut scratch = committed.to_vec(); - let k = solid.nonlinear_points(m, u, kin, committed, &mut scratch, true).unwrap().0.tangent.unwrap(); + let k = solid.nonlinear_points(m, u, kin, None, committed, &mut scratch, true).unwrap().0.tangent.unwrap(); let scale = k.data.iter().fold(0.0f64, |a, v| a.max(v.abs())); let n = 3 * u.len(); let h = if kin == Kinematics::Small { 1e-9 } else { 1e-6 }; @@ -560,7 +666,7 @@ mod tests { let shifted = |s: f64| { let mut v = u.to_vec(); v[j / 3][j % 3] += s; - solid.nonlinear_points(m, &v, kin, committed, &mut committed.to_vec(), false).unwrap().0.force + solid.nonlinear_points(m, &v, kin, None, committed, &mut committed.to_vec(), false).unwrap().0.force }; let (p, q) = (shifted(h), shifted(-h)); for i in 0..n { @@ -658,7 +764,7 @@ mod tests { // A first increment from the virgin state, committed, then a second from it. let mut committed = vec![0.0; points * HISTORY]; let first = displacement(x.len(), 0.5 * scale); - solid.nonlinear_points(&m, &first, kin, &vec![0.0; points * HISTORY], &mut committed, false).unwrap(); + solid.nonlinear_points(&m, &first, kin, None, &vec![0.0; points * HISTORY], &mut committed, false).unwrap(); assert!((0..points).all(|p| committed[p * HISTORY + PLASTIC_STRAIN] > 0.0), "{kind:?} {f:?} {kin:?} yields everywhere"); check_tangent_from(&solid, &m, &displacement(x.len(), scale), kin, &committed, 1e-6); } @@ -666,7 +772,7 @@ mod tests { // Incompatible modes have no inelastic form. let im = Solid::new(ElementKind::Hex8, &x8, Formulation::IncompatibleModes).unwrap(); assert!(matches!( - im.nonlinear(&m, &displacement(8, 0.0), Kinematics::Small, &[0.0; 8 * HISTORY], &mut [0.0; 8 * HISTORY], true), + im.nonlinear(&m, &displacement(8, 0.0), Kinematics::Small, None, &[0.0; 8 * HISTORY], &mut [0.0; 8 * HISTORY], true), Err(Error::NotSupported(_)) )); } @@ -699,17 +805,17 @@ mod tests { for kin in [Kinematics::Finite, Kinematics::Small] { let u = |scale: f64| -> Vec<[f64; 2]> { displacement(x.len(), scale).iter().map(|v| [v[0], v[1]]).collect() }; let mut committed = vec![0.0; points * HISTORY]; - plane.nonlinear_points(&m, &u(0.03), kin, &vec![0.0; points * HISTORY], &mut committed, false).unwrap(); + plane.nonlinear_points(&m, &u(0.03), kin, None, &vec![0.0; points * HISTORY], &mut committed, false).unwrap(); assert!((0..points).all(|p| committed[p * HISTORY + PLASTIC_STRAIN] > 0.0), "{theory:?} {kind:?} {kin:?}"); let at = u(0.06); - let k = plane.nonlinear_points(&m, &at, kin, &committed, &mut committed.clone(), true).unwrap().0.tangent.unwrap(); + let k = plane.nonlinear_points(&m, &at, kin, None, &committed, &mut committed.clone(), true).unwrap().0.tangent.unwrap(); let scale = k.data.iter().fold(0.0f64, |a, v| a.max(v.abs())); let h = 1e-8; for j in 0..2 * x.len() { let shifted = |s: f64| { let mut v = at.clone(); v[j / 2][j % 2] += s; - plane.nonlinear_points(&m, &v, kin, &committed, &mut committed.clone(), false).unwrap().0.force + plane.nonlinear_points(&m, &v, kin, None, &committed, &mut committed.clone(), false).unwrap().0.force }; let (p, q) = (shifted(h), shifted(-h)); for i in 0..2 * x.len() { @@ -726,11 +832,67 @@ mod tests { } let stress = Plane::new(ElementKind::Quad4, &quad8[..4], Formulation::Full, Theory::PlaneStress { thickness: 1.0 }).unwrap(); assert!(matches!( - stress.nonlinear(&m, &[[0.0; 2]; 4], Kinematics::Small, &[0.0; 4 * HISTORY], &mut [0.0; 4 * HISTORY], true), + stress.nonlinear(&m, &[[0.0; 2]; 4], Kinematics::Small, None, &[0.0; 4 * HISTORY], &mut [0.0; 4 * HISTORY], true), Err(Error::NotSupported(_)) )); } + /// A temperature field: the tangent is the derivative of the forces (finite and small strain, + /// full and hybrid, elastic and plastic), and a free expansion u = (F_θ − I) X of a uniform + /// temperature stresses nothing. + #[test] + fn thermal_tangents_and_free_expansion() { + use crate::elastoplastic::HISTORY; + use fieldcraft_model::{Elastic, Hardening, Plastic}; + let j2 = crate::plasticity::J2::from_model( + &Elastic::Isotropic { young: 1000.0, poisson: 0.3 }, + Some(&Plastic::Bilinear { yield_stress: 5.0, hardening_modulus: 50.0, hardening: Hardening::Isotropic }), + None, + ) + .unwrap(); + let plastic = MaterialModel::Plastic(std::sync::Arc::new(j2)); + let neo = MaterialModel::NeoHookean { shear: 300.0, bulk: 3000.0 }; + let x20 = hex20(); + let field = |n: usize| Temperature { + alpha: [1e-3, 2e-3, 1.5e-3], + reference: 20.0, + nodal: (0..n).map(|i| 20.0 + 40.0 * (i as f64 * 0.37).sin()).collect(), + }; + for f in [Formulation::Full, Formulation::Hybrid { pressures: 4 }] { + let solid = Solid::new(ElementKind::Hex20, &x20, f).unwrap(); + let points = solid.rule().points().len(); + let t = field(20); + for (m, kin, len) in [(&neo, Kinematics::Finite, 0), (&plastic, Kinematics::Finite, HISTORY), (&plastic, Kinematics::Small, HISTORY)] { + let u = displacement(20, 0.03); + let committed = vec![0.0; points * len]; + let k = solid.nonlinear(m, &u, kin, Some(&t), &committed, &mut committed.clone(), true).unwrap().tangent.unwrap(); + let scale = k.data.iter().fold(0.0f64, |a, v| a.max(v.abs())); + let h = if kin == Kinematics::Small { 1e-9 } else { 1e-7 }; + for j in 0..60 { + let shifted = |s: f64| { + let mut v = u.clone(); + v[j / 3][j % 3] += s; + solid.nonlinear(m, &v, kin, Some(&t), &committed, &mut committed.clone(), false).unwrap().force + }; + let (p, q) = (shifted(h), shifted(-h)); + for i in 0..60 { + let numeric = (p[i] - q[i]) / (2.0 * h); + assert!((numeric - k.get(i, j)).abs() < 1e-6 * scale, "{f:?} {} {kin:?} {i} {j}: {numeric} vs {}", m.name(), k.get(i, j)); + } + } + } + // Free expansion by a uniform 150 degrees: no force, no stress. + let uniform = Temperature { alpha: [1e-3, 2e-3, 1.5e-3], reference: 20.0, nodal: vec![170.0; 20] }; + let stretch = uniform.alpha.map(|a| 1.0 + 150.0 * a); + let u: Vec<[f64; 3]> = x20.iter().map(|x| std::array::from_fn(|k| (stretch[k] - 1.0) * x[k])).collect(); + let r = solid.nonlinear(&neo, &u, Kinematics::Finite, Some(&uniform), &[], &mut [], false).unwrap(); + assert!(r.force.iter().all(|v| v.abs() < 1e-10 * 3000.0), "{f:?}: {:?}", r.force); + let s = solid.nonlinear_stresses(&neo, &u, Kinematics::Finite, Some(&uniform), &[]).unwrap(); + assert!(s.iter().all(|p| p.stress.iter().all(|v| v.abs() < 1e-10 * 3000.0)), "{f:?}"); + assert!((s[0].strain[1] - stretch[1].ln()).abs() < 1e-12, "total logarithmic strain: {:?}", s[0].strain); + } + } + #[test] fn rigid_rotation_strains_and_stresses_nothing() { let x20 = hex20(); @@ -739,7 +901,7 @@ mod tests { let neo = MaterialModel::NeoHookean { shear: 300.0, bulk: 3000.0 }; for f in [Formulation::Full, Formulation::Hybrid { pressures: 4 }] { let solid = Solid::new(ElementKind::Hex20, &x20, f).unwrap(); - let s = solid.nonlinear_stresses(&neo, &u, Kinematics::Finite, &[]).unwrap(); + let s = solid.nonlinear_stresses(&neo, &u, Kinematics::Finite, None, &[]).unwrap(); assert!(s.iter().all(|p| p.stress.iter().chain(&p.strain).all(|v| v.abs() < 1e-12 * 3000.0)), "{f:?}: {:?}", s[0]); let r = solid.finite_strain(&neo, &u, &[], &mut [], false).unwrap(); assert!(r.force.iter().all(|v| v.abs() < 1e-11 * 3000.0)); diff --git a/crates/fem/src/plasticity.rs b/crates/fem/src/plasticity.rs index 7f13484..d11dcc2 100644 --- a/crates/fem/src/plasticity.rs +++ b/crates/fem/src/plasticity.rs @@ -229,6 +229,11 @@ impl J2 { Ok(J2 { shear, bulk, flow, damage, heating }) } + /// Whether the flow stress depends on the temperature (Johnson–Cook softening). + pub fn is_temperature_dependent(&self) -> bool { + matches!(self.flow, Some(Flow::JohnsonCook { .. })) + } + /// The temperature a point starts at. pub fn initial_temperature(&self) -> f64 { match self.flow { diff --git a/crates/solver/src/analysis.rs b/crates/solver/src/analysis.rs index 9aa1bb3..dbbcb23 100644 --- a/crates/solver/src/analysis.rs +++ b/crates/solver/src/analysis.rs @@ -353,9 +353,9 @@ pub fn solve_static(model: &Model, step: Option<&str>, limits: &Limits) -> Resul // loads. let (field, temperature) = match step.procedure { Procedure::CoupledTemperatureDisplacement {} => { - let (labels, t) = crate::heat::temperatures(model, &resolved, step, limits)?; - let at = |l: &u32| labels.binary_search(l).ok().and_then(|i| t.get(i)).copied(); - (Some(mesh.labels.iter().map(at).collect::>()), Some(mesh.labels.iter().map(|l| at(l).unwrap_or(0.0)).collect())) + let field = crate::heat::temperatures_at(model, &resolved, step, limits, &mesh.labels)?; + let temperature = field.iter().map(|t| t.unwrap_or(0.0)).collect(); + (Some(field), Some(temperature)) } _ => (temperature_loads(model, &resolved, &mesh, step)?, None), }; diff --git a/crates/solver/src/heat.rs b/crates/solver/src/heat.rs index 11f489e..5df9298 100644 --- a/crates/solver/src/heat.rs +++ b/crates/solver/src/heat.rs @@ -414,6 +414,13 @@ pub(crate) fn temperatures(model: &Model, resolved: &Resolved, step: &Step, limi Ok((s.mesh.labels, t)) } +/// The steady temperatures of a step's thermal conditions at the nodes `labels` (sorted), `None` +/// at the nodes conduction does not reach. +pub(crate) fn temperatures_at(model: &Model, resolved: &Resolved, step: &Step, limits: &Limits, labels: &[u32]) -> Result>> { + let (solved, t) = temperatures(model, resolved, step, limits)?; + Ok(labels.iter().map(|l| solved.binary_search(l).ok().and_then(|i| t.get(i)).copied()).collect()) +} + /// Solve a heat transfer step: the step named `step`, or the model's only one. pub fn solve_heat(model: &Model, step: Option<&str>, limits: &Limits) -> Result { let resolved = model.resolve()?; diff --git a/crates/solver/src/nonlinear.rs b/crates/solver/src/nonlinear.rs index d2ac241..8b2dc36 100644 --- a/crates/solver/src/nonlinear.rs +++ b/crates/solver/src/nonlinear.rs @@ -55,8 +55,18 @@ //! ([`fieldcraft_fem::finite`]); without nonlinear geometry an element of a plastic material //! takes the same loops with small strain, and an elastic one its linear stiffness (K u). //! Springs act along their fixed directions. Trusses, shells and beams (large rotations), plane -//! stress, incompatible modes (with nonlinear geometry, or with an inelastic material), -//! temperature and centrifugal loads are refused by name. +//! stress, incompatible modes (with nonlinear geometry, or with an inelastic material) and +//! centrifugal loads are refused by name. +//! +//! **Thermal strain.** The step's temperatures are its temperature loads, or, when it has +//! temperature boundary conditions or heat loads, their steady conduction solved first (one-way +//! coupling, as the linear coupled temperature–displacement step). They ramp with the load +//! factor from each material's reference temperature, ΔT(λ) = λ (T − T_ref). With nonlinear +//! geometry the deformation gradient splits multiplicatively, F = F_m F_θ with +//! F_θ = diag(1 + αᵢ ΔT), the material sees F_m and the stress is J_θ τ_m per reference volume +//! (fem's `Temperature`); with small strain the split is additive, ε_m = ε − α ΔT; elements on +//! their linear stiffness take the equivalent nodal forces, a dead load. A Johnson–Cook material +//! softens at the point's temperature. //! //! **Plasticity** (`Material::plastic`, J2 with the radial return explicit dynamics uses, //! [`fieldcraft_fem::elastoplastic`]): Simo's exponential return map at finite strain, the @@ -75,11 +85,11 @@ use fieldcraft_fem::constitutive::MaterialModel; use fieldcraft_fem::dense::Matrix; use fieldcraft_fem::element::{Extrapolation, extrapolation}; -use fieldcraft_fem::finite::Kinematics; +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::{Extent, Hardening, Load, Model, NonlinearStatic, Plastic, Procedure, SectionKind, Step}; +use fieldcraft_model::{BoundaryCondition, Extent, Hardening, Load, Model, NonlinearStatic, Plastic, Procedure, SectionKind, Step}; use fieldcraft_parallel as parallel; use serde::Serialize; @@ -90,7 +100,7 @@ use crate::analysis::{ use crate::cholesky::FactorError; use crate::constraints::Reduction; use crate::sparse::{Csr, add_block}; -use crate::system::{Mesh, loads, translations}; +use crate::system::{Mesh, add_element_thermal_loads, loads, temperature_loads, translations}; use crate::{Error, Limits, Result}; /// One Newton iteration of an increment. @@ -194,6 +204,9 @@ struct Problem<'a> { element_kinds: Vec, /// Each element's material model. materials: Vec, + /// The step's temperature field at λ = 1, per mesh node (`None` where the field leaves a + /// node at its material's reference temperature); `None` without one. + field: Option>>, /// The tangent's pattern, zero. pattern: Csr, } @@ -257,9 +270,27 @@ impl Problem<'_> { if self.finite { Some(Kinematics::Finite) } else { self.materials.get(e).filter(|m| m.is_path_dependent()).map(|_| Kinematics::Small) } } + /// The temperatures of element `e` at load factor λ, which ramps them from the reference + /// temperature to the field's (J2 Johnson–Cook softening and thermal strain): the + /// material's expansion and reference temperature, or no expansion from the Johnson–Cook + /// room temperature; `None` without a field or for a material neither affects. + fn temperature(&self, e: usize, lambda: f64) -> Option { + let field = self.field.as_ref()?; + let el = self.mesh.elements.get(e)?; + let softens = el.plastic.as_ref().is_some_and(|j2| j2.is_temperature_dependent()); + let (alpha, reference) = match el.expansion { + Some(x) => x, + None if softens => ([0.0; 3], el.plastic.as_ref().map_or(0.0, |j2| j2.initial_temperature())), + None => return None, + }; + let nodal: Vec = + el.nodes.iter().map(|n| field.get(*n).copied().flatten().map_or(reference, |t| reference + lambda * (t - reference))).collect(); + (softens || nodal.iter().any(|t| *t != reference)).then_some(Temperature { alpha, reference, nodal }) + } + /// Internal forces and the tangent of element `e` at displacements `u` (all degrees of - /// freedom). - fn element(&self, e: usize, u: &[f64], tangent: bool, committed: &[f64]) -> Result<(Vec, Option, Vec)> { + /// freedom) and load factor λ. + fn element(&self, e: usize, u: &[f64], lambda: f64, tangent: bool, committed: &[f64]) -> Result<(Vec, Option, Vec)> { let el = self.mesh.elements.get(e).ok_or_else(|| Error::Internal("element index".into()))?; let at = |err| Error::element(el.label, err); let ue = self.element_vector(e, u); @@ -271,7 +302,7 @@ impl Problem<'_> { }; let material = self.materials.get(e).ok_or_else(|| Error::Internal("element material".into()))?; let mut updated = committed.to_vec(); - let r = kernel.nonlinear(material, &ue, kin, committed, &mut updated, tangent).map_err(at)?; + let r = kernel.nonlinear(material, &ue, kin, self.temperature(e, lambda).as_ref(), committed, &mut updated, tangent).map_err(at)?; Ok((r.force, r.tangent, updated)) } @@ -305,7 +336,7 @@ impl Problem<'_> { for round in parallel::blocks(self.mesh.elements.len(), ROUND) { let start = round.start; let computed = parallel::map_grained(round.len(), GRAIN, |i| { - self.element(start + i, u, tangent, committed.get(start + i).map(Vec::as_slice).unwrap_or_default()) + self.element(start + i, u, lambda, tangent, committed.get(start + i).map(Vec::as_slice).unwrap_or_default()) }); for (i, c) in computed.into_iter().enumerate() { let (f, ke, updated) = c?; @@ -647,12 +678,8 @@ fn check_supported(model: &Model, mesh: &Mesh, step: &Step, settings: &Nonlinear return Err(at(format!("a hyperelastic material without nonlinear geometry (element {})", e.label))); } } - for name in &step.loads { - match model.load(name) { - Some(Load::Temperature { .. }) => return Err(at(format!("temperature loads (load `{name}`): thermal strain"))), - Some(Load::Centrifugal { .. }) => return Err(at(format!("centrifugal loads (load `{name}`), which depend on the deformed position"))), - _ => {} - } + if let Some(name) = step.loads.iter().find(|n| matches!(model.load(n), Some(Load::Centrifugal { .. }))) { + return Err(at(format!("centrifugal loads (load `{name}`), which depend on the deformed position"))); } plasticity(model, step, settings, warnings)?; if settings.nonlinear_geometry @@ -704,6 +731,23 @@ fn plasticity(model: &Model, step: &Step, settings: &NonlinearStatic, warnings: Ok(()) } +/// The step's temperature field per mesh node: the steady heat solution of its temperature +/// boundary conditions and heat loads when it has some (one-way coupling, as the coupled +/// temperature–displacement step of the linear solver), or its temperature loads. +fn temperatures(model: &Model, resolved: &fieldcraft_model::Resolved, mesh: &Mesh, step: &Step, limits: &Limits) -> Result>>> { + let conducts = step.boundary_conditions.iter().any(|b| model.boundary_conditions().get(b).is_some_and(BoundaryCondition::is_thermal)) + || step.loads.iter().any(|l| model.load(l).is_some_and(Load::is_thermal)); + let loads = temperature_loads(model, resolved, mesh, step)?; + match (conducts, loads) { + (true, Some(_)) => Err(Error::Invalid(format!( + "step `{}` gives temperatures both by temperature loads and by heat conditions (temperature boundary conditions, heat loads); use one", + step.name + ))), + (true, None) => Ok(Some(crate::heat::temperatures_at(model, resolved, step, limits, &mesh.labels)?)), + (false, field) => Ok(field), + } +} + /// The follower pressures of a step (with nonlinear geometry), and the step without them for /// the dead loads. fn followers(model: &Model, resolved: &fieldcraft_model::Resolved, mesh: &Mesh, step: &Step) -> Result<(Vec, Step)> { @@ -769,6 +813,7 @@ pub fn solve_nonlinear(model: &Model, step: Option<&str>, limits: &Limits) -> Re let (follower_loads, dead_step) = if settings.nonlinear_geometry { followers(model, &resolved, &mesh, step)? } else { (Vec::new(), step.clone()) }; let dead = loads(model, &resolved, &mesh, &dead_step, &red.active)?; + let field = temperatures(model, &resolved, &mesh, step, limits)?; let springs = if mesh.discrete.springs.iter().any(|s| !s.dashpot) { let mut s = graph.zero_matrix(&mesh.layout); mesh.discrete.add_stiffness(&mut s)?; @@ -785,7 +830,7 @@ pub fn solve_nonlinear(model: &Model, step: Option<&str>, limits: &Limits) -> Re mesh.layout.element_dofs(&e.nodes, *kinds).ok_or_else(|| Error::Internal(format!("element {} degrees of freedom", e.label))) }) .collect::>>()?; - let problem = Problem { + let mut problem = Problem { mesh: &mesh, red: &red, dead, @@ -802,8 +847,16 @@ pub fn solve_nonlinear(model: &Model, step: Option<&str>, limits: &Limits) -> Re None => MaterialModel::from_elasticity(&e.d, e.hyperelastic.as_ref()), }) .collect(), + field, pattern: graph.zero_matrix(&mesh.layout), }; + // Elements on their linear stiffness take the thermal strain as the equivalent nodal + // forces, a dead load proportional to the load factor like the temperatures. + if let Some(field) = problem.field.take() { + let linear: Vec = (0..mesh.elements.len()).map(|i| problem.kinematics(i).is_none()).collect(); + add_element_thermal_loads(&mut problem.dead, &mesh, &field, |i| linear.get(i).copied().unwrap_or(false))?; + problem.field = Some(field); + } let (lo, hi) = mesh.xyz.iter().fold(([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]), |(lo, hi), p| { (std::array::from_fn(|k| lo[k].min(p[k])), std::array::from_fn(|k| hi[k].max(p[k]))) }); @@ -939,7 +992,7 @@ pub fn solve_nonlinear(model: &Model, step: Option<&str>, limits: &Limits) -> Re "follower pressure: the skew part of its load stiffness ({a:.1e} of it, over the free degrees of freedom) is left out of the tangent; the solution is exact, Newton may converge more slowly" )); } - finish(&mesh, step, &red, &problem, u, evaluation, &run.committed, history, run.factor_entries, warnings) + finish(&mesh, step, &red, &problem, u, lambda, evaluation, &run.committed, history, run.factor_entries, warnings) } /// The final state as a static result: displacements, reactions, Cauchy stresses and @@ -951,6 +1004,7 @@ fn finish( red: &Reduction, problem: &Problem, u: Vec, + lambda: f64, e: Evaluation, committed: &[Vec], history: History, @@ -984,10 +1038,14 @@ fn finish( let history = committed.get(i).map(Vec::as_slice).unwrap_or_default(); let material = problem.materials.get(i).ok_or_else(|| Error::Internal("element material".into()))?; let points = match problem.kinematics(i) { - Some(kin) => kernel.nonlinear_stresses(material, &ue, kin, history), + Some(kin) => kernel.nonlinear_stresses(material, &ue, kin, problem.temperature(i, lambda).as_ref(), history), None => { let per_node: Vec<[f64; 3]> = el.nodes.iter().map(|n| displacement.get(*n).copied().unwrap_or_default()).collect(); - kernel.stresses_with(&el.d, &per_node, None) + let thermal = problem.field.as_deref().and_then(|f| el.thermal_strain(f)).map(|mut t| { + t.delta.iter_mut().for_each(|d| *d *= lambda); + t + }); + kernel.stresses_with(&el.d, &per_node, thermal.as_ref()) } } .map_err(at)?; diff --git a/crates/solver/src/system.rs b/crates/solver/src/system.rs index 0b34469..7657768 100644 --- a/crates/solver/src/system.rs +++ b/crates/solver/src/system.rs @@ -11,7 +11,7 @@ use fieldcraft_fem::beam::{Beam, BeamSection}; use fieldcraft_fem::constitutive::MaterialModel; use fieldcraft_fem::dense::Matrix; use fieldcraft_fem::element::{Formulation, Solid, ThermalStrain}; -use fieldcraft_fem::finite::Kinematics; +use fieldcraft_fem::finite::{Kinematics, Temperature}; use fieldcraft_fem::plane::{Plane, Theory}; use fieldcraft_fem::plasticity::J2; use fieldcraft_fem::shell::{Shell, ShellSection}; @@ -122,27 +122,36 @@ impl Kernel<'_> { /// Nonlinear internal forces and tangent ([`fieldcraft_fem::finite`]) for the element's /// displacements `u` (its degrees of freedom in order: 3 per node, or 2 for 2D elements), /// total Lagrangian or with small strain. + #[allow(clippy::too_many_arguments)] pub fn nonlinear( &self, m: &MaterialModel, u: &[f64], kin: Kinematics, + thermal: Option<&Temperature>, committed: &[f64], updated: &mut [f64], tangent: bool, ) -> fieldcraft_fem::Result { match self { - Kernel::Solid(s) => s.nonlinear(m, u.as_chunks::<3>().0, kin, committed, updated, tangent), - Kernel::Plane(p) => p.nonlinear(m, u.as_chunks::<2>().0, kin, committed, updated, tangent), + Kernel::Solid(s) => s.nonlinear(m, u.as_chunks::<3>().0, kin, thermal, committed, updated, tangent), + Kernel::Plane(p) => p.nonlinear(m, u.as_chunks::<2>().0, kin, thermal, committed, updated, tangent), } } /// Cauchy stress and strain (logarithmic at finite strain) at the integration points, as /// [`Self::nonlinear`]. - pub fn nonlinear_stresses(&self, m: &MaterialModel, u: &[f64], kin: Kinematics, committed: &[f64]) -> fieldcraft_fem::Result> { + pub fn nonlinear_stresses( + &self, + m: &MaterialModel, + u: &[f64], + kin: Kinematics, + thermal: Option<&Temperature>, + committed: &[f64], + ) -> fieldcraft_fem::Result> { match self { - Kernel::Solid(s) => s.nonlinear_stresses(m, u.as_chunks::<3>().0, kin, committed), - Kernel::Plane(p) => p.nonlinear_stresses(m, u.as_chunks::<2>().0, kin, committed), + Kernel::Solid(s) => s.nonlinear_stresses(m, u.as_chunks::<3>().0, kin, thermal, committed), + Kernel::Plane(p) => p.nonlinear_stresses(m, u.as_chunks::<2>().0, kin, thermal, committed), } } @@ -936,13 +945,7 @@ pub(crate) fn hold(slot: Option<&mut Option>, value: f64, conflict: impl Fn /// Add the equivalent nodal forces of the thermal strains of a temperature field. pub fn add_thermal_loads(f: &mut [f64], mesh: &Mesh, field: &[Option]) -> Result<()> { - for e in &mesh.elements { - let Some(th) = e.thermal_strain(field) else { continue }; - let forces = e.kernel()?.thermal_load(&e.d, &th).map_err(|err| Error::element(e.label, err))?; - for (node, v) in e.nodes.iter().zip(forces) { - add(f, &translations(mesh, *node)?, v); - } - } + add_element_thermal_loads(f, mesh, field, |_| true)?; for t in &mesh.trusses { let Some(th) = t.thermal_strain(field) else { continue }; let forces = t.truss()?.thermal_load(t.young, t.area, &th).map_err(|err| Error::element(t.label, err))?; @@ -953,6 +956,19 @@ pub fn add_thermal_loads(f: &mut [f64], mesh: &Mesh, field: &[Option]) -> R Ok(()) } +/// Add the equivalent nodal forces of the thermal strains of the solid and 2D elements whose +/// index `keep` accepts. +pub fn add_element_thermal_loads(f: &mut [f64], mesh: &Mesh, field: &[Option], keep: impl Fn(usize) -> bool) -> Result<()> { + for (_, e) in mesh.elements.iter().enumerate().filter(|(i, _)| keep(*i)) { + let Some(th) = e.thermal_strain(field) else { continue }; + let forces = e.kernel()?.thermal_load(&e.d, &th).map_err(|err| Error::element(e.label, err))?; + for (node, v) in e.nodes.iter().zip(forces) { + add(f, &translations(mesh, *node)?, v); + } + } + Ok(()) +} + /// The uniform force per volume a step's gravity and body forces put on each shell and beam /// element ([`Mesh::structural`] order): beam end forces are recovered with it. pub fn structural_body_forces(model: &Model, resolved: &Resolved, mesh: &Mesh, step: &Step) -> Result> { diff --git a/crates/solver/tests/nonlinear.rs b/crates/solver/tests/nonlinear.rs index 0946da9..bf548c7 100644 --- a/crates/solver/tests/nonlinear.rs +++ b/crates/solver/tests/nonlinear.rs @@ -387,7 +387,7 @@ fn nonlinear_steps_fail_cleanly() { step.set_procedure(Procedure::NonlinearStatic(bad.clone())); assert!(capped.model.replace_step(step).is_err(), "{bad:?}"); } - // Refused by name: incompatible modes, a hybrid St Venant–Kirchhoff element, temperature. + // Refused by name: incompatible modes, a hybrid St Venant–Kirchhoff element, centrifugal loads. let mut im = grid([2, 1, 1], ElementKind::Hex8, |q| q); im.define(iso(100.0, 0.3), None, incompatible(), nonlinear()); im.fixed("root", |q| near(q[0], 0.0)); @@ -402,10 +402,13 @@ fn nonlinear_steps_fail_cleanly() { hy.load("p", Load::Pressure { set: "end".into(), magnitude: 1.0 }); let err = solve_nonlinear(&hy.model, None, &Limits::default()).unwrap_err().to_string(); assert!(err.contains("hybrid elements of a St Venant–Kirchhoff material"), "{err}"); - hy.load("hot", Load::Temperature { target: Target::Set("root".into()), value: 50.0 }); + hy.load("spin", Load::Centrifugal { angular_velocity: 10.0, point: [0.0; 3], axis: [0.0, 0.0, 1.0], set: None }); neo_hookean(&mut hy, 100.0, 0.3); + let mut m = hy.model.material("Mat").unwrap().clone(); + m.density = Some(1.0); + hy.model.replace_material("Mat", m).unwrap(); let err = solve_nonlinear(&hy.model, None, &Limits::default()).unwrap_err().to_string(); - assert!(err.contains("temperature loads"), "{err}"); + assert!(err.contains("centrifugal loads"), "{err}"); // A model that can move freely: the first tangent is singular, named like the linear solver. let mut free = grid([2, 1, 1], ElementKind::Hex8, |q| q); free.define(iso(100.0, 0.3), None, full(), nonlinear()); diff --git a/crates/solver/tests/nonlinear_thermal.rs b/crates/solver/tests/nonlinear_thermal.rs new file mode 100644 index 0000000..8d51ce1 --- /dev/null +++ b/crates/solver/tests/nonlinear_thermal.rs @@ -0,0 +1,147 @@ +//! Thermal strain in nonlinear static steps against closed forms: free expansion and a fully +//! constrained heated block at finite strain (the multiplicative split F = F_m F_θ), a heated +//! bar held at its ends past yield at small strain (the additive split), the linear limit +//! against the linear solver, and a step's own heat conduction feeding its temperatures. + +mod common; + +use common::*; +use fieldcraft_model::*; +use fieldcraft_solver::{Limits, StaticResult, solve_nonlinear, solve_static}; + +fn expanding(g: &mut Grid, e: f64, nu: f64, alpha: f64, hyperelastic: Option, plastic: Option) { + let thermal = Thermal { + expansion: Some(Directional::Isotropic(alpha)), + reference_temperature: 20.0, + conductivity: Some(Directional::Isotropic(45.0)), + ..Thermal::default() + }; + let m = Material { elastic: Some(iso(e, nu)), hyperelastic, plastic, thermal: Some(thermal), ..Material::default() }; + g.model.replace_material("Mat", m).unwrap(); +} + +fn heat_all(g: &mut Grid, value: f64) { + g.node_set("everything", |_| true); + g.load("heat", Load::Temperature { target: Target::Set("everything".into()), value }); +} + +fn solve(g: &Grid) -> StaticResult { + solve_nonlinear(&g.model, Some("Step"), &Limits::default()).unwrap() +} + +fn symmetric(g: &mut Grid) { + 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)); +} + +/// A free neo-Hookean cube heated by 200 with α = 10⁻³ grows by the thermal stretch 1.2 along +/// every edge, stress-free; held on all six faces it takes the pressure of the mechanical +/// volume ratio ϑ⁻³: σ = κ (ϑ⁻³ − 1), full and hybrid elements. +#[test] +fn finite_strain_free_and_constrained_expansion() { + let (e, nu, alpha) = (10.0, 0.3, 1e-3); + let kappa = e / (3.0 * (1.0 - 2.0 * nu)); + let theta: f64 = 1.0 + alpha * 200.0; + for formulation in [Formulation::default(), Formulation { hybrid: true, ..Formulation::default() }] { + let mut g = grid([2, 2, 2], ElementKind::Hex20, |x| x); + g.define(iso(e, nu), None, formulation, Procedure::NonlinearStatic(NonlinearStatic::default())); + expanding(&mut g, e, nu, alpha, Some(Hyperelastic::NeoHookean {}), None); + symmetric(&mut g); + heat_all(&mut g, 220.0); + let r = solve(&g); + for (label, u) in r.nodes.iter().zip(&r.displacement) { + let x = g.model.node(*label).unwrap().xyz; + assert!((0..3).all(|k| (u[k] - (theta - 1.0) * x[k]).abs() < 1e-9), "{formulation:?} node {label}: {u:?}"); + } + assert!(r.stress.iter().flatten().all(|s| s.abs() < 1e-9 * kappa), "{formulation:?}: {:?}", r.stress[0]); + assert!((r.elements[0].points[0].strain[0] - theta.ln()).abs() < 1e-9, "the total logarithmic strain"); + + g.symmetry("x1", Symmetry::X, |x| near(x[0], 1.0)); + g.symmetry("y1", Symmetry::Y, |x| near(x[1], 1.0)); + g.symmetry("z1", Symmetry::Z, |x| near(x[2], 1.0)); + let r = solve(&g); + let want = kappa * (theta.powi(-3) - 1.0); + for s in &r.stress { + assert!((0..3).all(|k| (s[k] - want).abs() < 1e-9 * kappa) && s[3].abs() < 1e-9 * kappa, "{formulation:?}: {s:?} vs {want}"); + } + eprintln!("{formulation:?}: constrained σ = {:.9} (closed form {want:.9})", r.stress[0][0]); + } +} + +/// A steel bar held axially at both ends (lateral faces free) and heated at small strain: the +/// mechanical strain is −αΔT, so the stress is −EαΔT until it yields, then +/// −(σ_y + E_t (αΔT − σ_y/E)) with linear hardening. +#[test] +fn small_strain_heated_bar_held_at_its_ends_yields() { + let (e, nu, alpha, y0, h) = (200_000.0, 0.3, 1.2e-5, 250.0, 2000.0); + let small = NonlinearStatic { nonlinear_geometry: false, ..NonlinearStatic::default() }; + let plastic = Plastic::Bilinear { yield_stress: y0, hardening_modulus: h, hardening: Hardening::Isotropic }; + let et = e * h / (e + h); + for (dt, want) in [(80.0, -e * alpha * 80.0), (300.0, -(y0 + et * (alpha * 300.0 - y0 / e)))] { + let mut g = grid([4, 1, 1], ElementKind::Hex8, |x| [4.0 * x[0], x[1], x[2]]); + g.define(iso(e, nu), None, Formulation::default(), Procedure::NonlinearStatic(small.clone())); + expanding(&mut g, e, nu, alpha, None, Some(plastic.clone())); + symmetric(&mut g); + g.symmetry("x1", Symmetry::X, |x| near(x[0], 1.0)); + heat_all(&mut g, 20.0 + dt); + let r = solve(&g); + for s in &r.stress { + assert!((s[0] - want).abs() < 1e-8 * want.abs() && s[1].abs() < 1e-8 * want.abs(), "ΔT = {dt}: {s:?} vs {want}"); + } + eprintln!("ΔT = {dt}: σ = {:.6} (closed form {want:.6})", r.stress[0][0]); + } +} + +/// Without nonlinear geometry an elastic model is the linear solution, thermal strain +/// included; a step's own temperature boundary conditions give its temperatures by conduction, +/// as in the linear coupled step. +#[test] +fn thermal_linear_limit_and_conducted_temperatures() { + let (e, nu, alpha) = (200_000.0, 0.3, 1.2e-5); + let build = |procedure: Procedure, plastic: Option, conducted: bool| { + let mut g = grid([6, 1, 2], ElementKind::Hex20, |q| [6.0 * q[0], q[1], q[2]]); + g.define(iso(e, nu), None, Formulation::default(), procedure); + expanding(&mut g, e, nu, alpha, None, plastic); + g.fixed("root", |q| near(q[0], 0.0)); + if conducted { + g.node_set("cold", |q| near(q[0], 0.0)); + g.bc("cold_t", BoundaryCondition::Temperature { target: Target::Set("cold".into()), value: 20.0 }); + g.node_set("hot", |q| near(q[0], 1.0)); + g.bc("hot_t", BoundaryCondition::Temperature { target: Target::Set("hot".into()), value: 320.0 }); + } else { + // Hot on top, cold below: the bar bends. + g.node_set("top", |q| near(q[2], 1.0)); + g.load("hot", Load::Temperature { target: Target::Set("top".into()), value: 120.0 }); + } + g + }; + let small = Procedure::NonlinearStatic(NonlinearStatic { nonlinear_geometry: false, ..NonlinearStatic::default() }); + let linear = solve_static(&build(Procedure::Static {}, None, false).model, None, &Limits::default()).unwrap(); + let r = solve(&build(small.clone(), None, false)); + let scale = linear.displacement.iter().flatten().fold(0.0f64, |m, v| m.max(v.abs())); + for (a, b) in r.displacement.iter().zip(&linear.displacement) { + assert!((0..3).all(|k| (a[k] - b[k]).abs() < 1e-9 * scale), "{a:?} vs {b:?}"); + } + for (a, b) in r.stress.iter().zip(&linear.stress) { + assert!((0..6).all(|k| (a[k] - b[k]).abs() < 1e-6 * e * alpha * 100.0), "{a:?} vs {b:?}"); + } + // Conducted: T linear in x from 20 to 320 over a bar fixed at its cold end; the step's + // temperature conditions are solved first, as the coupled linear step does, and a plastic + // material that stays below yield gives that step's solution. + let plastic = Plastic::Bilinear { yield_stress: 250.0, hardening_modulus: 0.0, hardening: Hardening::Isotropic }; + let r = solve(&build(small, Some(plastic), true)); + let coupled = solve_static(&build(Procedure::CoupledTemperatureDisplacement {}, None, true).model, None, &Limits::default()).unwrap(); + let largest = coupled.von_mises.iter().fold(0.0f64, |m, v| m.max(*v)); + assert!(largest < 250.0 && largest > 1.0, "{largest}"); + let scale = coupled.displacement.iter().flatten().fold(0.0f64, |m, v| m.max(v.abs())); + for (a, b) in r.displacement.iter().zip(&coupled.displacement) { + assert!((0..3).all(|k| (a[k] - b[k]).abs() < 1e-8 * scale), "{a:?} vs {b:?}"); + } + assert!(r.elements.iter().flat_map(|e| &e.plastic_strain).all(|p| *p == 0.0)); + // Both ways of giving temperatures at once is refused. + let mut both = build(Procedure::NonlinearStatic(NonlinearStatic::default()), None, true); + heat_all(&mut both, 50.0); + let err = solve_nonlinear(&both.model, None, &Limits::default()).unwrap_err().to_string(); + assert!(err.contains("both by temperature loads and by heat conditions"), "{err}"); +}