flavio/constitutive/solid/hyperelastic/gent/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::*;
5
6/// The Gent hyperelastic constitutive model.[^cite]
7///
8/// [^cite]: A.N. Gent, [Rubber Chem. Technol. **69**, 59 (1996)](https://doi.org/10.5254/1.3538357).
9///
10/// **Parameters**
11/// - The bulk modulus $`\kappa`$.
12/// - The shear modulus $`\mu`$.
13/// - The extensibility $`J_m`$.
14///
15/// **External variables**
16/// - The deformation gradient $`\mathbf{F}`$.
17///
18/// **Internal variables**
19/// - None.
20///
21/// **Notes**
22/// - The Gent model reduces to the [Neo-Hookean model](NeoHookean) when $`J_m\to\infty`$.
23#[derive(Debug)]
24pub struct Gent<'a> {
25    parameters: Parameters<'a>,
26}
27
28impl Gent<'_> {
29    /// Returns the extensibility.
30    fn get_extensibility(&self) -> &Scalar {
31        &self.parameters[2]
32    }
33}
34
35impl<'a> Constitutive<'a> for Gent<'a> {
36    fn new(parameters: Parameters<'a>) -> Self {
37        Self { parameters }
38    }
39}
40
41impl<'a> Solid<'a> for Gent<'a> {
42    fn get_bulk_modulus(&self) -> &Scalar {
43        &self.parameters[0]
44    }
45    fn get_shear_modulus(&self) -> &Scalar {
46        &self.parameters[1]
47    }
48}
49
50impl<'a> Elastic<'a> for Gent<'a> {
51    /// Calculates and returns the Cauchy stress.
52    ///
53    /// ```math
54    /// \boldsymbol{\sigma}(\mathbf{F}) = \frac{J^{-1}\mu J_m {\mathbf{B}^* }'}{J_m - \mathrm{tr}(\mathbf{B}^* ) + 3} + \frac{\kappa}{2}\left(J - \frac{1}{J}\right)\mathbf{1}
55    /// ```
56    fn calculate_cauchy_stress(
57        &self,
58        deformation_gradient: &DeformationGradient,
59    ) -> Result<CauchyStress, ConstitutiveError> {
60        let jacobian = deformation_gradient.determinant();
61        if jacobian > 0.0 {
62            let isochoric_left_cauchy_green_deformation = self
63                .calculate_left_cauchy_green_deformation(deformation_gradient)
64                / jacobian.powf(TWO_THIRDS);
65            let (
66                deviatoric_isochoric_left_cauchy_green_deformation,
67                isochoric_left_cauchy_green_deformation_trace,
68            ) = isochoric_left_cauchy_green_deformation.deviatoric_and_trace();
69            let denominator =
70                self.get_extensibility() - isochoric_left_cauchy_green_deformation_trace + 3.0;
71            if denominator <= 0.0 {
72                Err(ConstitutiveError::Custom(
73                    "Maximum extensibility reached.".to_string(),
74                    deformation_gradient.copy(),
75                    format!("{:?}", &self),
76                ))
77            } else {
78                Ok((deviatoric_isochoric_left_cauchy_green_deformation
79                    * self.get_shear_modulus()
80                    * self.get_extensibility()
81                    / jacobian)
82                    / denominator
83                    + IDENTITY * self.get_bulk_modulus() * 0.5 * (jacobian - 1.0 / jacobian))
84            }
85        } else {
86            Err(ConstitutiveError::InvalidJacobian(
87                jacobian,
88                deformation_gradient.copy(),
89                format!("{:?}", &self),
90            ))
91        }
92    }
93    /// Calculates and returns the tangent stiffness associated with the Cauchy stress.
94    ///
95    /// ```math
96    /// \mathcal{T}_{ijkL}(\mathbf{F}) = \frac{J^{-5/3}\mu J_m}{J_m - \mathrm{tr}(\mathbf{B}^* ) + 3}\Bigg[ \delta_{ik}F_{jL} + \delta_{jk}F_{iL} - \frac{2}{3}\,\delta_{ij}F_{kL} + \frac{2{B_{ij}^* }' F_{kL}}{J_m - \mathrm{tr}(\mathbf{B}^* ) + 3} - \left(\frac{5}{3} + \frac{2}{3}\frac{\mathrm{tr}(\mathbf{B}^* )}{J_m - \mathrm{tr}(\mathbf{B}^* ) + 3}\right) J^{2/3} {B_{ij}^* }' F_{kL}^{-T} \Bigg] + \frac{\kappa}{2} \left(J + \frac{1}{J}\right)\delta_{ij}F_{kL}^{-T}
97    /// ```
98    fn calculate_cauchy_tangent_stiffness(
99        &self,
100        deformation_gradient: &DeformationGradient,
101    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
102        let jacobian = deformation_gradient.determinant();
103        if jacobian > 0.0 {
104            let inverse_transpose_deformation_gradient = deformation_gradient.inverse_transpose();
105            let isochoric_left_cauchy_green_deformation = self
106                .calculate_left_cauchy_green_deformation(deformation_gradient)
107                / jacobian.powf(TWO_THIRDS);
108            let (
109                deviatoric_isochoric_left_cauchy_green_deformation,
110                isochoric_left_cauchy_green_deformation_trace,
111            ) = isochoric_left_cauchy_green_deformation.deviatoric_and_trace();
112            let denominator =
113                self.get_extensibility() - isochoric_left_cauchy_green_deformation_trace + 3.0;
114            if denominator <= 0.0 {
115                Err(ConstitutiveError::Custom(
116                    "Maximum extensibility reached.".to_string(),
117                    deformation_gradient.copy(),
118                    format!("{:?}", &self),
119                ))
120            } else {
121                let prefactor =
122                    self.get_shear_modulus() * self.get_extensibility() / jacobian / denominator;
123                Ok(
124                    (CauchyTangentStiffness::dyad_ik_jl(&IDENTITY, deformation_gradient)
125                        + CauchyTangentStiffness::dyad_il_jk(deformation_gradient, &IDENTITY)
126                        - CauchyTangentStiffness::dyad_ij_kl(&IDENTITY, deformation_gradient)
127                            * (TWO_THIRDS)
128                        + CauchyTangentStiffness::dyad_ij_kl(
129                            &deviatoric_isochoric_left_cauchy_green_deformation,
130                            deformation_gradient,
131                        ) * (2.0 / denominator))
132                        * (prefactor / jacobian.powf(TWO_THIRDS))
133                        + CauchyTangentStiffness::dyad_ij_kl(
134                            &(IDENTITY
135                                * (0.5 * self.get_bulk_modulus() * (jacobian + 1.0 / jacobian))
136                                - deviatoric_isochoric_left_cauchy_green_deformation
137                                    * prefactor
138                                    * ((5.0
139                                        + 2.0 * isochoric_left_cauchy_green_deformation_trace
140                                            / denominator)
141                                        / 3.0)),
142                            &inverse_transpose_deformation_gradient,
143                        ),
144                )
145            }
146        } else {
147            Err(ConstitutiveError::InvalidJacobian(
148                jacobian,
149                deformation_gradient.copy(),
150                format!("{:?}", &self),
151            ))
152        }
153    }
154}
155
156impl<'a> Hyperelastic<'a> for Gent<'a> {
157    /// Calculates and returns the Helmholtz free energy density.
158    ///
159    /// ```math
160    /// a(\mathbf{F}) = -\frac{\mu J_m}{2}\,\ln\left[1 - \frac{\mathrm{tr}(\mathbf{B}^* ) - 3}{J_m}\right] + \frac{\kappa}{2}\left[\frac{1}{2}\left(J^2 - 1\right) - \ln J\right]
161    /// ```
162    fn calculate_helmholtz_free_energy_density(
163        &self,
164        deformation_gradient: &DeformationGradient,
165    ) -> Result<Scalar, ConstitutiveError> {
166        let jacobian = deformation_gradient.determinant();
167        if jacobian > 0.0 {
168            let factor = (self
169                .calculate_left_cauchy_green_deformation(deformation_gradient)
170                .trace()
171                / jacobian.powf(TWO_THIRDS)
172                - 3.0)
173                / self.get_extensibility();
174            if factor >= 1.0 {
175                Err(ConstitutiveError::Custom(
176                    "Maximum extensibility reached.".to_string(),
177                    deformation_gradient.copy(),
178                    format!("{:?}", &self),
179                ))
180            } else {
181                Ok(0.5
182                    * (-self.get_shear_modulus() * self.get_extensibility() * (1.0 - factor).ln()
183                        + self.get_bulk_modulus()
184                            * (0.5 * (jacobian.powi(2) - 1.0) - jacobian.ln())))
185            }
186        } else {
187            Err(ConstitutiveError::InvalidJacobian(
188                jacobian,
189                deformation_gradient.copy(),
190                format!("{:?}", &self),
191            ))
192        }
193    }
194}