From b98178947758ef9fb0f839339c73502c3234f76f Mon Sep 17 00:00:00 2001 From: mrbuche Date: Mon, 26 Aug 2024 11:14:19 -0600 Subject: [PATCH] seems to work poorly on AH for some reason --- .../solid/elastic/almansi_hamel/test.rs | 2 - src/constitutive/solid/elastic/mod.rs | 44 --------- src/constitutive/solid/elastic/test.rs | 74 --------------- src/constitutive/solid/hyperelastic/mod.rs | 61 +++++++++++++ src/constitutive/solid/hyperelastic/test.rs | 91 ++++++++++++++++++- 5 files changed, 151 insertions(+), 121 deletions(-) diff --git a/src/constitutive/solid/elastic/almansi_hamel/test.rs b/src/constitutive/solid/elastic/almansi_hamel/test.rs index 8608341..3557751 100644 --- a/src/constitutive/solid/elastic/almansi_hamel/test.rs +++ b/src/constitutive/solid/elastic/almansi_hamel/test.rs @@ -6,5 +6,3 @@ test_solid_elastic_constitutive_model!( ALMANSIHAMELPARAMETERS, AlmansiHamel::new(ALMANSIHAMELPARAMETERS) ); - -test_solve_uniaxial!(AlmansiHamel::new(ALMANSIHAMELPARAMETERS)); diff --git a/src/constitutive/solid/elastic/mod.rs b/src/constitutive/solid/elastic/mod.rs index e52c0b7..0d60a5c 100644 --- a/src/constitutive/solid/elastic/mod.rs +++ b/src/constitutive/solid/elastic/mod.rs @@ -19,8 +19,6 @@ pub use almansi_hamel::AlmansiHamel; use super::*; -const MAXIMUM_STEPS: usize = 10_000; - /// Required methods for elastic constitutive models. pub trait Elastic<'a> where @@ -151,46 +149,4 @@ where &second_piola_kirchoff_stress, )) } - fn solve_biaxial( - &self, - _deformation_gradient_11: &Scalar, - _deformation_gradient_22: &Scalar, - ) -> Result<(Scalar, Scalar), ConstitutiveError> { - todo!() - } - fn solve_pure_shear(&self, _: &Scalar) -> Result { - todo!() - } - fn solve_simple_shear(&self, _: &Scalar) -> Result { - todo!() - } - fn solve_uniaxial( - &self, - deformation_gradient_11: &Scalar, - ) -> Result<(DeformationGradient, CauchyStress), ConstitutiveError> { - let mut cauchy_stress = ZERO; - let mut deformation_gradient = IDENTITY_10; - deformation_gradient[0][0] = *deformation_gradient_11; - deformation_gradient[1][1] = 1.0 / deformation_gradient_11.sqrt(); - deformation_gradient[2][2] = deformation_gradient[1][1]; - let mut residual = 0.0; - let mut residual_abs = 1.0; - let mut residual_rel = 1.0; - let mut steps: usize = 0; - while residual_abs >= ABS_TOL && residual_rel >= REL_TOL { - if steps > MAXIMUM_STEPS { - return Err(ConstitutiveError::SolveError); - } else { - deformation_gradient[1][1] -= residual - / self.calculate_cauchy_tangent_stiffness(&deformation_gradient)?[1][1][1][1]; - deformation_gradient[2][2] = deformation_gradient[1][1]; - cauchy_stress = self.calculate_cauchy_stress(&deformation_gradient)?; - residual = cauchy_stress[1][1]; - residual_abs = residual.abs(); - residual_rel = residual_abs / cauchy_stress[0][0].abs(); - steps += 1; - } - } - Ok((deformation_gradient, cauchy_stress)) - } } diff --git a/src/constitutive/solid/elastic/test.rs b/src/constitutive/solid/elastic/test.rs index 20c03ba..556add6 100644 --- a/src/constitutive/solid/elastic/test.rs +++ b/src/constitutive/solid/elastic/test.rs @@ -2,80 +2,6 @@ use crate::mechanics::Scalar; pub const ALMANSIHAMELPARAMETERS: &[Scalar; 2] = &[13.0, 3.0]; -macro_rules! test_solve_uniaxial { - ($constitutive_model_constructed: expr) => { - #[test] - fn solve_biaxial_compression() { - let (a, b) = $constitutive_model_constructed - .solve_biaxial(&0.55, &0.88) - .expect("the unexpected"); - assert!(a < 0.0); - assert!(b < 0.0); - assert!(a < b); - } - #[test] - fn solve_biaxial_mixed() { - let (a, b) = $constitutive_model_constructed - .solve_biaxial(&3.3, &0.44) - .expect("the unexpected"); - assert!(a > 0.0); - assert!(b < 0.0); - } - #[test] - fn solve_biaxial_tension() { - let (a, b) = $constitutive_model_constructed - .solve_biaxial(&3.3, &2.2) - .expect("the unexpected"); - assert!(a > 0.0); - assert!(b > 0.0); - assert!(a > b); - } - #[test] - fn solve_biaxial_undeformed() { - let (a, b) = $constitutive_model_constructed - .solve_biaxial(&1.0, &1.0) - .expect("the unexpected"); - assert_eq!(a, 0.0); - assert_eq!(b, 0.0); - } - #[test] - fn solve_uniaxial_compression() { - assert!( - $constitutive_model_constructed - .solve_uniaxial(&0.22) - .expect("the unexpected") - .1[0][0] - < 0.0 - ) - } - #[test] - fn solve_uniaxial_tension() { - let (deformation_gradient, cauchy_stress) = $constitutive_model_constructed - .solve_uniaxial(&4.4) - .expect("the unexpected"); - assert!(cauchy_stress[0][0] > 0.0); - crate::test::assert_eq_within_tols(&(cauchy_stress[1][1] / cauchy_stress[0][0]), &0.0); - crate::test::assert_eq_within_tols(&(cauchy_stress[2][2] / cauchy_stress[0][0]), &0.0); - assert!(cauchy_stress.is_diagonal()); - assert!(deformation_gradient[0][0] > deformation_gradient[1][1]); - crate::test::assert_eq_within_tols( - &deformation_gradient[1][1], - &deformation_gradient[2][2], - ); - assert!(deformation_gradient.is_diagonal()); - } - #[test] - fn solve_uniaxial_undeformed() { - let (deformation_gradient, cauchy_stress) = $constitutive_model_constructed - .solve_uniaxial(&1.0) - .expect("the unexpected"); - assert!(cauchy_stress.is_zero()); - assert!(deformation_gradient.is_identity()); - } - }; -} -pub(crate) use test_solve_uniaxial; - macro_rules! calculate_cauchy_stress_from_deformation_gradient { ($constitutive_model_constructed: expr, $deformation_gradient: expr) => { $constitutive_model_constructed diff --git a/src/constitutive/solid/hyperelastic/mod.rs b/src/constitutive/solid/hyperelastic/mod.rs index 3cd5bc4..071770b 100644 --- a/src/constitutive/solid/hyperelastic/mod.rs +++ b/src/constitutive/solid/hyperelastic/mod.rs @@ -35,6 +35,8 @@ pub use self::{ }; use super::{elastic::Elastic, *}; +const MAXIMUM_STEPS: usize = 10_000; + /// Required methods for hyperelastic constitutive models. pub trait Hyperelastic<'a> where @@ -49,4 +51,63 @@ where &self, deformation_gradient: &DeformationGradient, ) -> Result; + fn solve_biaxial( + &self, + deformation_gradient_11: &Scalar, + deformation_gradient_22: &Scalar, + ) -> Result<(DeformationGradient, CauchyStress), ConstitutiveError> { + let mut cauchy_stress = ZERO; + let mut deformation_gradient = IDENTITY_10; + deformation_gradient[0][0] = *deformation_gradient_11; + deformation_gradient[1][1] = *deformation_gradient_22; + deformation_gradient[2][2] = 1.0 / deformation_gradient_11 / deformation_gradient_22; + let mut residual = 0.0; + let mut residual_abs = 1.0; + let mut residual_rel = 1.0; + let mut steps: usize = 0; + while residual_abs >= ABS_TOL && residual_rel >= REL_TOL { + if steps > MAXIMUM_STEPS { + return Err(ConstitutiveError::SolveError); + } else { + deformation_gradient[2][2] -= residual + / self.calculate_cauchy_tangent_stiffness(&deformation_gradient)?[2][2][2][2]; + cauchy_stress = self.calculate_cauchy_stress(&deformation_gradient)?; + residual = cauchy_stress[2][2]; + residual_abs = residual.abs(); + residual_rel = residual_abs + / (cauchy_stress[0][0].powi(2) + cauchy_stress[1][1].powi(2)).sqrt(); + steps += 1; + } + } + Ok((deformation_gradient, cauchy_stress)) + } + fn solve_uniaxial( + &self, + deformation_gradient_11: &Scalar, + ) -> Result<(DeformationGradient, CauchyStress), ConstitutiveError> { + let mut cauchy_stress = ZERO; + let mut deformation_gradient = IDENTITY_10; + deformation_gradient[0][0] = *deformation_gradient_11; + deformation_gradient[1][1] = 1.0 / deformation_gradient_11.sqrt(); + deformation_gradient[2][2] = deformation_gradient[1][1]; + let mut residual = 0.0; + let mut residual_abs = 1.0; + let mut residual_rel = 1.0; + let mut steps: usize = 0; + while residual_abs >= ABS_TOL && residual_rel >= REL_TOL { + if steps > MAXIMUM_STEPS { + return Err(ConstitutiveError::SolveError); + } else { + deformation_gradient[1][1] -= residual + / self.calculate_cauchy_tangent_stiffness(&deformation_gradient)?[1][1][1][1]; + deformation_gradient[2][2] = deformation_gradient[1][1]; + cauchy_stress = self.calculate_cauchy_stress(&deformation_gradient)?; + residual = cauchy_stress[1][1]; + residual_abs = residual.abs(); + residual_rel = residual_abs / cauchy_stress[0][0].abs(); + steps += 1; + } + } + Ok((deformation_gradient, cauchy_stress)) + } } diff --git a/src/constitutive/solid/hyperelastic/test.rs b/src/constitutive/solid/hyperelastic/test.rs index 02cfb5f..488c8b0 100644 --- a/src/constitutive/solid/hyperelastic/test.rs +++ b/src/constitutive/solid/hyperelastic/test.rs @@ -25,6 +25,96 @@ pub const YEOHPARAMETERS: &[Scalar; 6] = &[ 1e-5, ]; +macro_rules! test_solve_uniaxial { + ($constitutive_model_constructed: expr) => { + #[test] + fn solve_biaxial_compression() { + let (deformation_gradient, cauchy_stress) = $constitutive_model_constructed + .solve_biaxial(&0.55, &0.88) + .expect("the unexpected"); + assert!(cauchy_stress[0][0] < 0.0); + assert!(cauchy_stress[1][1] < 0.0); + crate::test::assert_eq_within_tols( + &(cauchy_stress[2][2] + / (cauchy_stress[0][0].powi(2) + cauchy_stress[1][1].powi(2)).sqrt()), + &0.0, + ); + assert!(cauchy_stress.is_diagonal()); + assert!(deformation_gradient.is_diagonal()); + } + #[test] + fn solve_biaxial_mixed() { + let (deformation_gradient, cauchy_stress) = $constitutive_model_constructed + .solve_biaxial(&3.3, &0.44) + .expect("the unexpected"); + assert!(cauchy_stress[0][0] > cauchy_stress[1][1]); + crate::test::assert_eq_within_tols( + &(cauchy_stress[2][2] + / (cauchy_stress[0][0].powi(2) + cauchy_stress[1][1].powi(2)).sqrt()), + &0.0, + ); + assert!(cauchy_stress.is_diagonal()); + assert!(deformation_gradient.is_diagonal()); + } + #[test] + fn solve_biaxial_tension() { + let (deformation_gradient, cauchy_stress) = $constitutive_model_constructed + .solve_biaxial(&3.3, &2.2) + .expect("the unexpected"); + assert!(cauchy_stress[0][0] > cauchy_stress[1][1]); + assert!(cauchy_stress[1][1] > 0.0); + crate::test::assert_eq_within_tols( + &(cauchy_stress[2][2] + / (cauchy_stress[0][0].powi(2) + cauchy_stress[1][1].powi(2)).sqrt()), + &0.0, + ); + assert!(cauchy_stress.is_diagonal()); + assert!(deformation_gradient.is_diagonal()); + } + #[test] + fn solve_biaxial_undeformed() { + let (deformation_gradient, cauchy_stress) = $constitutive_model_constructed + .solve_biaxial(&1.0, &1.0) + .expect("the unexpected"); + assert!(cauchy_stress.is_zero()); + assert!(deformation_gradient.is_identity()); + } + #[test] + fn solve_uniaxial_compression() { + let (deformation_gradient, cauchy_stress) = $constitutive_model_constructed + .solve_uniaxial(&0.44) + .expect("the unexpected"); + assert!(cauchy_stress[0][0] < 0.0); + crate::test::assert_eq_within_tols(&(cauchy_stress[1][1] / cauchy_stress[0][0]), &0.0); + crate::test::assert_eq_within_tols(&(cauchy_stress[2][2] / cauchy_stress[0][0]), &0.0); + assert!(cauchy_stress.is_diagonal()); + assert_eq!(deformation_gradient[1][1], deformation_gradient[2][2],); + assert!(deformation_gradient.is_diagonal()); + } + #[test] + fn solve_uniaxial_tension() { + let (deformation_gradient, cauchy_stress) = $constitutive_model_constructed + .solve_uniaxial(&4.4) + .expect("the unexpected"); + assert!(cauchy_stress[0][0] > 0.0); + crate::test::assert_eq_within_tols(&(cauchy_stress[1][1] / cauchy_stress[0][0]), &0.0); + crate::test::assert_eq_within_tols(&(cauchy_stress[2][2] / cauchy_stress[0][0]), &0.0); + assert!(cauchy_stress.is_diagonal()); + assert_eq!(deformation_gradient[1][1], deformation_gradient[2][2],); + assert!(deformation_gradient.is_diagonal()); + } + #[test] + fn solve_uniaxial_undeformed() { + let (deformation_gradient, cauchy_stress) = $constitutive_model_constructed + .solve_uniaxial(&1.0) + .expect("the unexpected"); + assert!(cauchy_stress.is_zero()); + assert!(deformation_gradient.is_identity()); + } + }; +} +pub(crate) use test_solve_uniaxial; + macro_rules! calculate_helmholtz_free_energy_density_from_deformation_gradient_simple { ($constitutive_model_constructed: expr, $deformation_gradient: expr) => { $constitutive_model_constructed @@ -50,7 +140,6 @@ macro_rules! use_elastic_macros { calculate_second_piola_kirchoff_stress_from_deformation_gradient_rotated, calculate_second_piola_kirchoff_stress_from_deformation_gradient_simple, calculate_second_piola_kirchoff_tangent_stiffness_from_deformation_gradient, - test_solve_uniaxial, }; }; }