conspire/constitutive/solid/hyperviscoelastic/saint_venant_kirchhoff/
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            hyperviscoelastic::Hyperviscoelastic, viscoelastic::Viscoelastic,
11        },
12    },
13    math::{IDENTITY_00, Rank2},
14    mechanics::{
15        Deformation, DeformationGradient, DeformationGradientRate, Scalar,
16        SecondPiolaKirchhoffRateTangentStiffness, SecondPiolaKirchhoffStress,
17    },
18};
19
20/// The Saint Venant-Kirchhoff hyperviscoelastic 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 Green-Saint Venant strain measure is given by $`\mathbf{E}=\tfrac{1}{2}(\mathbf{C}-\mathbf{1})`$.
37#[derive(Debug)]
38pub struct SaintVenantKirchhoff {
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 SaintVenantKirchhoff {
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 SaintVenantKirchhoff {
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 SaintVenantKirchhoff {
68    /// Calculates and returns the second Piola-Kirchhoff stress.
69    ///
70    /// ```math
71    /// \mathbf{S}(\mathbf{F},\dot\mathbf{F}) = 2\mu\mathbf{E}' + \kappa\,\mathrm{tr}(\mathbf{E})\mathbf{1} + 2\eta\dot{\mathbf{E}}' + \zeta\,\mathrm{tr}(\dot{\mathbf{E}})\mathbf{1}
72    /// ```
73    fn second_piola_kirchhoff_stress(
74        &self,
75        deformation_gradient: &DeformationGradient,
76        deformation_gradient_rate: &DeformationGradientRate,
77    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
78        let _jacobian = self.jacobian(deformation_gradient)?;
79        let (deviatoric_strain, strain_trace) =
80            ((deformation_gradient.right_cauchy_green() - IDENTITY_00) * 0.5)
81                .deviatoric_and_trace();
82        let first_term = deformation_gradient_rate.transpose() * deformation_gradient;
83        let (deviatoric_strain_rate, strain_rate_trace) =
84            ((&first_term + first_term.transpose()) * 0.5).deviatoric_and_trace();
85        Ok(deviatoric_strain * (2.0 * self.shear_modulus())
86            + deviatoric_strain_rate * (2.0 * self.shear_viscosity())
87            + IDENTITY_00
88                * (self.bulk_modulus() * strain_trace + self.bulk_viscosity() * strain_rate_trace))
89    }
90    /// Calculates and returns the rate tangent stiffness associated with the second Piola-Kirchhoff stress.
91    ///
92    /// ```math
93    /// \mathcal{W}_{IJkL}(\mathbf{F}) = \eta\,\delta_{JL}F_{kI} + \eta\,\delta_{IL}F_{kJ} + \left(\zeta - \frac{2}{3}\,\eta\right)\delta_{IJ}F_{kL}
94    /// ```
95    fn second_piola_kirchhoff_rate_tangent_stiffness(
96        &self,
97        deformation_gradient: &DeformationGradient,
98        _: &DeformationGradientRate,
99    ) -> Result<SecondPiolaKirchhoffRateTangentStiffness, ConstitutiveError> {
100        let _jacobian = self.jacobian(deformation_gradient)?;
101        let scaled_deformation_gradient_transpose =
102            deformation_gradient.transpose() * self.shear_viscosity();
103        Ok(SecondPiolaKirchhoffRateTangentStiffness::dyad_ik_jl(
104            &scaled_deformation_gradient_transpose,
105            &IDENTITY_00,
106        ) + SecondPiolaKirchhoffRateTangentStiffness::dyad_il_jk(
107            &IDENTITY_00,
108            &scaled_deformation_gradient_transpose,
109        ) + SecondPiolaKirchhoffRateTangentStiffness::dyad_ij_kl(
110            &(IDENTITY_00 * (self.bulk_viscosity() - TWO_THIRDS * self.shear_viscosity())),
111            deformation_gradient,
112        ))
113    }
114}
115
116impl ElasticHyperviscous for SaintVenantKirchhoff {
117    /// Calculates and returns the viscous dissipation.
118    ///
119    /// ```math
120    /// \phi(\mathbf{F},\dot{\mathbf{F}}) = \eta\,\mathrm{tr}(\dot{\mathbf{E}}^2) + \frac{1}{2}\left(\zeta - \frac{2}{3}\,\eta\right)\mathrm{tr}(\dot{\mathbf{E}})^2
121    /// ```
122    fn viscous_dissipation(
123        &self,
124        deformation_gradient: &DeformationGradient,
125        deformation_gradient_rate: &DeformationGradientRate,
126    ) -> Result<Scalar, ConstitutiveError> {
127        let _jacobian = self.jacobian(deformation_gradient)?;
128        let first_term = deformation_gradient_rate.transpose() * deformation_gradient;
129        let strain_rate = (&first_term + first_term.transpose()) * 0.5;
130        Ok(self.shear_viscosity() * strain_rate.squared_trace()
131            + 0.5
132                * (self.bulk_viscosity() - TWO_THIRDS * self.shear_viscosity())
133                * strain_rate.trace().powi(2))
134    }
135}
136
137impl Hyperviscoelastic for SaintVenantKirchhoff {
138    /// Calculates and returns the Helmholtz free energy density.
139    ///
140    /// ```math
141    /// a(\mathbf{F}) = \mu\,\mathrm{tr}(\mathbf{E}^2) + \frac{1}{2}\left(\kappa - \frac{2}{3}\,\mu\right)\mathrm{tr}(\mathbf{E})^2
142    /// ```
143    fn helmholtz_free_energy_density(
144        &self,
145        deformation_gradient: &DeformationGradient,
146    ) -> Result<Scalar, ConstitutiveError> {
147        let _jacobian = self.jacobian(deformation_gradient)?;
148        let strain = (deformation_gradient.right_cauchy_green() - IDENTITY_00) * 0.5;
149        Ok(self.shear_modulus() * strain.squared_trace()
150            + 0.5
151                * (self.bulk_modulus() - TWO_THIRDS * self.shear_modulus())
152                * strain.trace().powi(2))
153    }
154}