conspire/constitutive/solid/hyperelastic/ogden/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{
5    constitutive::{
6        ConstitutiveError,
7        solid::{Solid, elastic::Elastic, hyperelastic::Hyperelastic},
8    },
9    math::{
10        ContractThirdFourthWithFirstSecond, Current, IDENTITY, Quantity, Rank2, Tensor,
11        TensorRank0List, TensorRank2, TensorRank4,
12    },
13    mechanics::{
14        CauchyStress, CauchyTangentStiffness, Deformation, DeformationGradient,
15        LeftCauchyGreenDeformation, Scalar,
16    },
17    units::{EnergyDensity, Stress},
18};
19
20#[doc = include_str!("doc.md")]
21#[derive(Clone, Debug)]
22pub struct Ogden {
23    /// The bulk modulus $`\kappa`$.
24    pub bulk_modulus: Quantity<Stress>,
25    /// The shear moduli $`\mu_n`$.
26    pub shear_moduli: Vec<Quantity<Stress>>,
27    /// The exponents $`\alpha_n`$.
28    pub exponents: Vec<Scalar>,
29}
30
31enum Spectrum {
32    Eigen(TensorRank0List<3>, LeftCauchyGreenDeformation),
33    Fallback(LeftCauchyGreenDeformation),
34}
35
36impl Spectrum {
37    fn new(tensor: &LeftCauchyGreenDeformation, model: &Ogden) -> Result<Self, ConstitutiveError> {
38        if tensor.is_diagonal() || (tensor - &IDENTITY).norm() < 1e-2 {
39            Ok(Self::Fallback(tensor.clone()))
40        } else {
41            let (eigenvalues, eigenvectors) = tensor
42                .eigen()
43                .map_err(|error| ConstitutiveError::upstream(error, model))?;
44            Ok(Self::Eigen(eigenvalues, eigenvectors))
45        }
46    }
47    fn powm(
48        &self,
49        exponent: Scalar,
50        model: &Ogden,
51    ) -> Result<LeftCauchyGreenDeformation, ConstitutiveError> {
52        match self {
53            Self::Eigen(eigenvalues, eigenvectors) => {
54                TensorRank2::powm_from_eigen(eigenvalues, eigenvectors, exponent)
55            }
56            Self::Fallback(tensor) => tensor.powm(exponent),
57        }
58        .map_err(|error| ConstitutiveError::upstream(error, model))
59    }
60    fn dpowm(
61        &self,
62        exponent: Scalar,
63        model: &Ogden,
64    ) -> Result<TensorRank4<3, Current, Current, Current, Current>, ConstitutiveError> {
65        match self {
66            Self::Eigen(eigenvalues, eigenvectors) => {
67                TensorRank2::dpowm_from_eigen(eigenvalues, eigenvectors, exponent)
68            }
69            Self::Fallback(tensor) => tensor.dpowm(exponent),
70        }
71        .map_err(|error| ConstitutiveError::upstream(error, model))
72    }
73}
74
75impl Solid for Ogden {
76    fn bulk_modulus(&self) -> Quantity<Stress> {
77        self.bulk_modulus
78    }
79    fn shear_modulus(&self) -> Quantity<Stress> {
80        self.shear_moduli
81            .iter()
82            .zip(self.exponents.iter())
83            .map(|(modulus, exponent)| *modulus * exponent)
84            .sum::<Quantity<Stress>>()
85            * 0.5
86    }
87}
88
89impl Elastic for Ogden {
90    #[doc = include_str!("cauchy_stress.md")]
91    fn cauchy_stress(
92        &self,
93        deformation_gradient: &DeformationGradient,
94    ) -> Result<CauchyStress, ConstitutiveError> {
95        let jacobian = self.jacobian(deformation_gradient)?;
96        let left_cauchy_green = deformation_gradient.left_cauchy_green();
97        let spectrum = Spectrum::new(&left_cauchy_green, self)?;
98        let mut cauchy_stress = IDENTITY * self.bulk_modulus() * 0.5 * (jacobian - 1.0 / jacobian);
99        for (modulus, exponent) in self.shear_moduli.iter().zip(self.exponents.iter()) {
100            let scaling = *modulus / jacobian.powf(exponent / 3.0 + 1.0);
101            cauchy_stress += spectrum.powm(exponent / 2.0, self)?.deviatoric() * scaling;
102        }
103        Ok(cauchy_stress)
104    }
105    #[doc = include_str!("cauchy_tangent_stiffness.md")]
106    fn cauchy_tangent_stiffness(
107        &self,
108        deformation_gradient: &DeformationGradient,
109    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
110        let jacobian = self.jacobian(deformation_gradient)?;
111        let left_cauchy_green = deformation_gradient.left_cauchy_green();
112        let spectrum = Spectrum::new(&left_cauchy_green, self)?;
113        let inverse_transpose_deformation_gradient = deformation_gradient.inverse_transpose();
114        let mut cauchy_tangent_stiffness = TensorRank4::dyad_ij_kl(
115            &(IDENTITY * (self.bulk_modulus() * 0.5 * (jacobian + 1.0 / jacobian))),
116            &inverse_transpose_deformation_gradient,
117        );
118        for (modulus, exponent) in self.shear_moduli.iter().zip(self.exponents.iter()) {
119            let half_exponent = exponent / 2.0;
120            let scaling = *modulus / jacobian.powf(exponent / 3.0 + 1.0);
121            let scaled_deformation_gradient = deformation_gradient * scaling;
122            let raw = spectrum
123                .dpowm(half_exponent, self)?
124                .contract_third_fourth_with_first_second(
125                    &(TensorRank4::dyad_il_jk(&scaled_deformation_gradient, &IDENTITY)
126                        + TensorRank4::dyad_ik_jl(&IDENTITY, &scaled_deformation_gradient)),
127                );
128            let trace_term = spectrum.powm(half_exponent - 1.0, self)?
129                * deformation_gradient
130                * (scaling * half_exponent * 2.0);
131            let cauchy_stress_n = spectrum.powm(half_exponent, self)?.deviatoric() * scaling;
132            cauchy_tangent_stiffness += raw
133                - TensorRank4::dyad_ij_kl(&(IDENTITY * (1.0 / 3.0)), &trace_term)
134                - TensorRank4::dyad_ij_kl(
135                    &(cauchy_stress_n * (exponent / 3.0 + 1.0)),
136                    &inverse_transpose_deformation_gradient,
137                );
138        }
139        Ok(cauchy_tangent_stiffness)
140    }
141}
142
143impl Hyperelastic for Ogden {
144    #[doc = include_str!("helmholtz_free_energy_density.md")]
145    fn helmholtz_free_energy_density(
146        &self,
147        deformation_gradient: &DeformationGradient,
148    ) -> Result<Quantity<EnergyDensity>, ConstitutiveError> {
149        let jacobian = self.jacobian(deformation_gradient)?;
150        let left_cauchy_green = deformation_gradient.left_cauchy_green();
151        let spectrum = Spectrum::new(&left_cauchy_green, self)?;
152        let mut helmholtz_free_energy_density =
153            self.bulk_modulus() * 0.5 * (0.5 * (jacobian.powi(2) - 1.0) - jacobian.ln());
154        for (modulus, exponent) in self.shear_moduli.iter().zip(self.exponents.iter()) {
155            helmholtz_free_energy_density += (*modulus / *exponent)
156                * (spectrum.powm(exponent / 2.0, self)?.trace() / jacobian.powf(exponent / 3.0)
157                    - 3.0);
158        }
159        Ok(helmholtz_free_energy_density)
160    }
161}
conspire/constitutive/solid/hyperelastic/ogden/mod.rs

conspire/constitutive/solid/hyperelastic/ogden/
mod.rs
