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

1#[cfg(test)]
2mod test;
3
4use crate::{
5    constitutive::{
6        Constitutive, ConstitutiveError,
7        solid::{FIVE_THIRDS, Solid, TWO_THIRDS, elastic::Elastic, hyperelastic::Hyperelastic},
8    },
9    math::{IDENTITY, Rank2},
10    mechanics::{CauchyStress, CauchyTangentStiffness, Deformation, DeformationGradient, Scalar},
11};
12use std::iter::once;
13
14const SEVEN_THIRDS: Scalar = 7.0 / 3.0;
15
16#[doc = include_str!("doc.md")]
17#[derive(Debug)]
18pub struct Yeoh<const N: usize> {
19    /// The bulk modulus $`\kappa`$.
20    pub bulk_modulus: Scalar,
21    /// The shear modulus $`\mu`$.
22    pub shear_modulus: Scalar,
23    /// The extra moduli $`\mu_n`$ for $`n=2\ldots N`$.
24    pub extra_moduli: [Scalar; N],
25}
26
27impl<const N: usize> Yeoh<N> {
28    /// Returns an array of the extra moduli.
29    pub fn extra_moduli(&self) -> &[Scalar] {
30        &self.extra_moduli
31    }
32}
33
34impl<const N: usize> Solid for Yeoh<N> {
35    fn bulk_modulus(&self) -> &Scalar {
36        &self.bulk_modulus
37    }
38    fn shear_modulus(&self) -> &Scalar {
39        &self.shear_modulus
40    }
41}
42
43impl<const N: usize> Elastic for Yeoh<N> {
44    #[doc = include_str!("cauchy_stress.md")]
45    fn cauchy_stress(
46        &self,
47        deformation_gradient: &DeformationGradient,
48    ) -> Result<CauchyStress, ConstitutiveError> {
49        let jacobian = self.jacobian(deformation_gradient)?;
50        let (deviatoric_left_cauchy_green_deformation, left_cauchy_green_deformation_trace) =
51            deformation_gradient
52                .left_cauchy_green()
53                .deviatoric_and_trace();
54        let scalar_term = left_cauchy_green_deformation_trace / jacobian.powf(TWO_THIRDS) - 3.0;
55        Ok(deviatoric_left_cauchy_green_deformation
56            * once(self.shear_modulus())
57                .chain(self.extra_moduli().iter())
58                .enumerate()
59                .map(|(n, modulus)| ((n as Scalar) + 1.0) * modulus * scalar_term.powi(n as i32))
60                .sum::<Scalar>()
61            / jacobian.powf(FIVE_THIRDS)
62            + IDENTITY * self.bulk_modulus() * 0.5 * (jacobian - 1.0 / jacobian))
63    }
64    #[doc = include_str!("cauchy_tangent_stiffness.md")]
65    fn cauchy_tangent_stiffness(
66        &self,
67        deformation_gradient: &DeformationGradient,
68    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
69        let jacobian = self.jacobian(deformation_gradient)?;
70        let inverse_transpose_deformation_gradient = deformation_gradient.inverse_transpose();
71        let left_cauchy_green_deformation = deformation_gradient.left_cauchy_green();
72        let scalar_term = left_cauchy_green_deformation.trace() / jacobian.powf(TWO_THIRDS) - 3.0;
73        let scaled_modulus = once(self.shear_modulus())
74            .chain(self.extra_moduli().iter())
75            .enumerate()
76            .map(|(n, modulus)| ((n as Scalar) + 1.0) * modulus * scalar_term.powi(n as i32))
77            .sum::<Scalar>()
78            / jacobian.powf(FIVE_THIRDS);
79        let deviatoric_left_cauchy_green_deformation = left_cauchy_green_deformation.deviatoric();
80        let last_term = CauchyTangentStiffness::dyad_ij_kl(
81            &deviatoric_left_cauchy_green_deformation,
82            &((left_cauchy_green_deformation.deviatoric()
83                * &inverse_transpose_deformation_gradient)
84                * (2.0
85                    * self
86                        .extra_moduli()
87                        .iter()
88                        .enumerate()
89                        .map(|(n, modulus)| {
90                            ((n as Scalar) + 2.0)
91                                * ((n as Scalar) + 1.0)
92                                * modulus
93                                * scalar_term.powi(n as i32)
94                        })
95                        .sum::<Scalar>()
96                    / jacobian.powf(SEVEN_THIRDS))),
97        );
98        Ok(
99            (CauchyTangentStiffness::dyad_ik_jl(&IDENTITY, deformation_gradient)
100                + CauchyTangentStiffness::dyad_il_jk(deformation_gradient, &IDENTITY)
101                - CauchyTangentStiffness::dyad_ij_kl(&IDENTITY, deformation_gradient)
102                    * (TWO_THIRDS))
103                * scaled_modulus
104                + CauchyTangentStiffness::dyad_ij_kl(
105                    &(IDENTITY * (0.5 * self.bulk_modulus() * (jacobian + 1.0 / jacobian))
106                        - deviatoric_left_cauchy_green_deformation
107                            * (scaled_modulus * FIVE_THIRDS)),
108                    &inverse_transpose_deformation_gradient,
109                )
110                + last_term,
111        )
112    }
113}
114
115impl<const N: usize> Hyperelastic for Yeoh<N> {
116    #[doc = include_str!("helmholtz_free_energy_density.md")]
117    fn helmholtz_free_energy_density(
118        &self,
119        deformation_gradient: &DeformationGradient,
120    ) -> Result<Scalar, ConstitutiveError> {
121        let jacobian = self.jacobian(deformation_gradient)?;
122        let scalar_term =
123            deformation_gradient.left_cauchy_green().trace() / jacobian.powf(TWO_THIRDS) - 3.0;
124        Ok(0.5
125            * (once(self.shear_modulus())
126                .chain(self.extra_moduli().iter())
127                .enumerate()
128                .map(|(n, modulus)| modulus * scalar_term.powi((n + 1) as i32))
129                .sum::<Scalar>()
130                + self.bulk_modulus() * (0.5 * (jacobian.powi(2) - 1.0) - jacobian.ln())))
131    }
132}