conspire/constitutive/solid/elastic_hyperviscous/almansi_hamel/
mod.rs

1#[cfg(test)]
2mod test;
3
4use crate::{
5    constitutive::{
6        Constitutive, ConstitutiveError,
7        fluid::viscous::Viscous,
8        solid::{
9            Solid, TWO_THIRDS, elastic_hyperviscous::ElasticHyperviscous,
10            viscoelastic::Viscoelastic,
11        },
12    },
13    math::{IDENTITY, Rank2},
14    mechanics::{
15        CauchyRateTangentStiffness, CauchyStress, DeformationGradient, DeformationGradientRate,
16        Scalar,
17    },
18};
19
20/// The Almansi-Hamel viscoelastic constitutive model.
21///
22/// **Parameters**
23/// - The bulk modulus $`\kappa`$.
24/// - The shear modulus $`\mu`$.
25/// - The bulk viscosity $`\zeta`$.
26/// - The shear viscosity $`\eta`$.
27///
28/// **External variables**
29/// - The deformation gradient $`\mathbf{F}`$.
30/// - The deformation gradient rate $`\dot{\mathbf{F}}`$.
31///
32/// **Internal variables**
33/// - None.
34///
35/// **Notes**
36/// - The Almansi-Hamel strain measure is given by $`\mathbf{e}=\tfrac{1}{2}(\mathbf{1}-\mathbf{B}^{-1})`$.
37#[derive(Debug)]
38pub struct AlmansiHamel {
39    /// The bulk modulus $`\kappa`$.
40    pub bulk_modulus: Scalar,
41    /// The shear modulus $`\mu`$.
42    pub shear_modulus: Scalar,
43    /// The bulk viscosity $`\zeta`$.
44    pub bulk_viscosity: Scalar,
45    /// The shear viscosity $`\eta`$.
46    pub shear_viscosity: Scalar,
47}
48
49impl Solid for AlmansiHamel {
50    fn bulk_modulus(&self) -> &Scalar {
51        &self.bulk_modulus
52    }
53    fn shear_modulus(&self) -> &Scalar {
54        &self.shear_modulus
55    }
56}
57
58impl Viscous for AlmansiHamel {
59    fn bulk_viscosity(&self) -> &Scalar {
60        &self.bulk_viscosity
61    }
62    fn shear_viscosity(&self) -> &Scalar {
63        &self.shear_viscosity
64    }
65}
66
67impl Viscoelastic for AlmansiHamel {
68    /// Calculates and returns the Cauchy stress.
69    ///
70    /// ```math
71    /// \boldsymbol{\sigma}(\mathbf{F},\dot\mathbf{F}) = 2\mu\mathbf{e}' + \kappa\,\mathrm{tr}(\mathbf{e})\mathbf{1} + 2\eta\mathbf{D}' + \zeta\,\mathrm{tr}(\mathbf{D})\mathbf{1}
72    /// ```
73    fn cauchy_stress(
74        &self,
75        deformation_gradient: &DeformationGradient,
76        deformation_gradient_rate: &DeformationGradientRate,
77    ) -> Result<CauchyStress, ConstitutiveError> {
78        let jacobian = self.jacobian(deformation_gradient)?;
79        let inverse_deformation_gradient = deformation_gradient.inverse();
80        let strain = (IDENTITY
81            - inverse_deformation_gradient.transpose() * &inverse_deformation_gradient)
82            * 0.5;
83        let (deviatoric_strain, strain_trace) = strain.deviatoric_and_trace();
84        let velocity_gradient = deformation_gradient_rate * inverse_deformation_gradient;
85        let strain_rate = (&velocity_gradient + velocity_gradient.transpose()) * 0.5;
86        let (deviatoric_strain_rate, strain_rate_trace) = strain_rate.deviatoric_and_trace();
87        Ok(deviatoric_strain * (2.0 * self.shear_modulus() / jacobian)
88            + deviatoric_strain_rate * (2.0 * self.shear_viscosity() / jacobian)
89            + IDENTITY
90                * ((self.bulk_modulus() * strain_trace
91                    + self.bulk_viscosity() * strain_rate_trace)
92                    / jacobian))
93    }
94    /// Calculates and returns the rate tangent stiffness associated with the Cauchy stress.
95    ///
96    /// ```math
97    /// \mathcal{V}_{IJkL}(\mathbf{F}) = \eta\,\delta_{ik}F_{jL}^{-T} + \eta\,\delta_{jk}F_{iL}^{-T} + \left(\zeta - \frac{2}{3}\,\eta\right)\delta_{ij}F_{kL}^{-T}
98    /// ```
99    fn cauchy_rate_tangent_stiffness(
100        &self,
101        deformation_gradient: &DeformationGradient,
102        _: &DeformationGradientRate,
103    ) -> Result<CauchyRateTangentStiffness, ConstitutiveError> {
104        let jacobian = self.jacobian(deformation_gradient)?;
105        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
106        let scaled_deformation_gradient_inverse_transpose =
107            &deformation_gradient_inverse_transpose * self.shear_viscosity() / jacobian;
108        Ok(CauchyRateTangentStiffness::dyad_ik_jl(
109            &IDENTITY,
110            &scaled_deformation_gradient_inverse_transpose,
111        ) + CauchyRateTangentStiffness::dyad_il_jk(
112            &scaled_deformation_gradient_inverse_transpose,
113            &IDENTITY,
114        ) + CauchyRateTangentStiffness::dyad_ij_kl(
115            &(IDENTITY
116                * ((self.bulk_viscosity() - TWO_THIRDS * self.shear_viscosity()) / jacobian)),
117            &deformation_gradient_inverse_transpose,
118        ))
119    }
120}
121
122impl ElasticHyperviscous for AlmansiHamel {
123    /// Calculates and returns the viscous dissipation.
124    ///
125    /// ```math
126    /// \phi(\mathbf{F},\dot{\mathbf{F}}) = \eta\,\mathrm{tr}(\mathbf{D}^2) + \frac{1}{2}\left(\zeta - \frac{2}{3}\,\eta\right)\mathrm{tr}(\mathbf{D})^2
127    /// ```
128    fn viscous_dissipation(
129        &self,
130        deformation_gradient: &DeformationGradient,
131        deformation_gradient_rate: &DeformationGradientRate,
132    ) -> Result<Scalar, ConstitutiveError> {
133        let _jacobian = self.jacobian(deformation_gradient)?;
134        let velocity_gradient = deformation_gradient_rate * deformation_gradient.inverse();
135        let strain_rate = (&velocity_gradient + velocity_gradient.transpose()) * 0.5;
136        Ok(self.shear_viscosity() * strain_rate.squared_trace()
137            + 0.5
138                * (self.bulk_viscosity() - TWO_THIRDS * self.shear_viscosity())
139                * strain_rate.trace().powi(2))
140    }
141}