conspire/constitutive/solid/viscoelastic/
mod.rs

1//! Viscoelastic constitutive models.
2//!
3//! ---
4//!
5//! Viscoelastic constitutive models cannot be defined by a Helmholtz free energy density function and viscous dissipation function.
6//! These constitutive models are therefore defined by a relation for the stress as a function of the deformation gradient and rate.
7//! Consequently, the rate tangent stiffness associated with the first Piola-Kirchhoff stress is not symmetric for these models.
8//!
9//! ```math
10//! \mathcal{U}_{iJkL} \neq \mathcal{U}_{kLiJ}
11//! ```
12
13#[cfg(test)]
14pub mod test;
15
16use super::{super::fluid::viscous::Viscous, *};
17use crate::math::{
18    Matrix, TensorVec, Vector,
19    integrate::{Explicit, IntegrationError},
20    optimize::{EqualityConstraint, FirstOrderRootFinding, NewtonRaphson, OptimizeError},
21};
22
23/// Possible applied loads.
24pub enum AppliedLoad<'a> {
25    /// Uniaxial stress given $`\dot{F}_{11}`$.
26    UniaxialStress(fn(Scalar) -> Scalar, &'a [Scalar]),
27    /// Biaxial stress given $`\dot{F}_{11}`$ and $`\dot{F}_{22}`$.
28    BiaxialStress(fn(Scalar) -> Scalar, fn(Scalar) -> Scalar, &'a [Scalar]),
29}
30
31/// Required methods for viscoelastic constitutive models.
32pub trait Viscoelastic
33where
34    Self: Solid + Viscous,
35{
36    /// Calculates and returns the Cauchy stress.
37    ///
38    /// ```math
39    /// \boldsymbol{\sigma} = J^{-1}\mathbf{P}\cdot\mathbf{F}^T
40    /// ```
41    fn cauchy_stress(
42        &self,
43        deformation_gradient: &DeformationGradient,
44        deformation_gradient_rate: &DeformationGradientRate,
45    ) -> Result<CauchyStress, ConstitutiveError> {
46        Ok(deformation_gradient
47            * self
48                .second_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_rate)?
49            * deformation_gradient.transpose()
50            / deformation_gradient.determinant())
51    }
52    /// Calculates and returns the rate tangent stiffness associated with the Cauchy stress.
53    ///
54    /// ```math
55    /// \mathcal{V}_{ijkL} = \frac{\partial\sigma_{ij}}{\partial\dot{F}_{kL}} = J^{-1} \mathcal{W}_{MNkL} F_{iM} F_{jN}
56    /// ```
57    fn cauchy_rate_tangent_stiffness(
58        &self,
59        deformation_gradient: &DeformationGradient,
60        deformation_gradient_rate: &DeformationGradientRate,
61    ) -> Result<CauchyRateTangentStiffness, ConstitutiveError> {
62        Ok(self
63            .second_piola_kirchhoff_rate_tangent_stiffness(
64                deformation_gradient,
65                deformation_gradient_rate,
66            )?
67            .contract_first_second_indices_with_second_indices_of(
68                deformation_gradient,
69                deformation_gradient,
70            )
71            / deformation_gradient.determinant())
72    }
73    /// Calculates and returns the first Piola-Kirchhoff stress.
74    ///
75    /// ```math
76    /// \mathbf{P} = J\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}
77    /// ```
78    fn first_piola_kirchhoff_stress(
79        &self,
80        deformation_gradient: &DeformationGradient,
81        deformation_gradient_rate: &DeformationGradientRate,
82    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
83        Ok(
84            self.cauchy_stress(deformation_gradient, deformation_gradient_rate)?
85                * deformation_gradient.inverse_transpose()
86                * deformation_gradient.determinant(),
87        )
88    }
89    /// Calculates and returns the rate tangent stiffness associated with the first Piola-Kirchhoff stress.
90    ///
91    /// ```math
92    /// \mathcal{U}_{iJkL} = \frac{\partial P_{iJ}}{\partial\dot{F}_{kL}} = J \mathcal{V}_{iskL} F_{sJ}^{-T}
93    /// ```
94    fn first_piola_kirchhoff_rate_tangent_stiffness(
95        &self,
96        deformation_gradient: &DeformationGradient,
97        deformation_gradient_rate: &DeformationGradientRate,
98    ) -> Result<FirstPiolaKirchhoffRateTangentStiffness, ConstitutiveError> {
99        Ok(self
100            .cauchy_rate_tangent_stiffness(deformation_gradient, deformation_gradient_rate)?
101            .contract_second_index_with_first_index_of(&deformation_gradient.inverse_transpose())
102            * deformation_gradient.determinant())
103    }
104    /// Calculates and returns the second Piola-Kirchhoff stress.
105    ///
106    /// ```math
107    /// \mathbf{S} = \mathbf{F}^{-1}\cdot\mathbf{P}
108    /// ```
109    fn second_piola_kirchhoff_stress(
110        &self,
111        deformation_gradient: &DeformationGradient,
112        deformation_gradient_rate: &DeformationGradientRate,
113    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
114        Ok(deformation_gradient.inverse()
115            * self.cauchy_stress(deformation_gradient, deformation_gradient_rate)?
116            * deformation_gradient.inverse_transpose()
117            * deformation_gradient.determinant())
118    }
119    /// Calculates and returns the rate tangent stiffness associated with the second Piola-Kirchhoff stress.
120    ///
121    /// ```math
122    /// \mathcal{W}_{IJkL} = \frac{\partial S_{IJ}}{\partial\dot{F}_{kL}} = \mathcal{U}_{mJkL}F_{mI}^{-T} = J \mathcal{V}_{mnkL} F_{mI}^{-T} F_{nJ}^{-T}
123    /// ```
124    fn second_piola_kirchhoff_rate_tangent_stiffness(
125        &self,
126        deformation_gradient: &DeformationGradient,
127        deformation_gradient_rate: &DeformationGradientRate,
128    ) -> Result<SecondPiolaKirchhoffRateTangentStiffness, ConstitutiveError> {
129        let deformation_gradient_inverse = deformation_gradient.inverse();
130        Ok(self
131            .cauchy_rate_tangent_stiffness(deformation_gradient, deformation_gradient_rate)?
132            .contract_first_second_indices_with_second_indices_of(
133                &deformation_gradient_inverse,
134                &deformation_gradient_inverse,
135            )
136            * deformation_gradient.determinant())
137    }
138    /// Solve for the unknown components of the deformation gradient and rate under an applied load.
139    ///
140    /// ```math
141    /// \mathbf{P}(\mathbf{F},\dot{\mathbf{F}}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
142    /// ```
143    fn root(
144        &self,
145        applied_load: AppliedLoad,
146        integrator: impl Explicit<DeformationGradientRate, DeformationGradientRates>,
147    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), IntegrationError> {
148        match applied_load {
149            AppliedLoad::UniaxialStress(deformation_gradient_rate_11, time) => integrator
150                .integrate(
151                    |t: Scalar, deformation_gradient: &DeformationGradient| {
152                        Ok(self.root_uniaxial_inner(
153                            deformation_gradient,
154                            deformation_gradient_rate_11(t),
155                        )?)
156                    },
157                    time,
158                    DeformationGradient::identity(),
159                ),
160            AppliedLoad::BiaxialStress(
161                deformation_gradient_rate_11,
162                deformation_gradient_rate_22,
163                time,
164            ) => integrator.integrate(
165                |t: Scalar, deformation_gradient: &DeformationGradient| {
166                    Ok(self.root_biaxial_inner(
167                        deformation_gradient,
168                        deformation_gradient_rate_11(t),
169                        deformation_gradient_rate_22(t),
170                    )?)
171                },
172                time,
173                DeformationGradient::identity(),
174            ),
175        }
176    }
177    #[doc(hidden)]
178    fn root_uniaxial_inner(
179        &self,
180        deformation_gradient: &DeformationGradient,
181        deformation_gradient_rate_11: Scalar,
182    ) -> Result<DeformationGradientRate, OptimizeError> {
183        let solver = NewtonRaphson {
184            ..Default::default()
185        };
186        let mut matrix = Matrix::zero(4, 9);
187        let mut vector = Vector::zero(4);
188        matrix[0][0] = 1.0;
189        matrix[1][1] = 1.0;
190        matrix[2][2] = 1.0;
191        matrix[3][5] = 1.0;
192        vector[0] = deformation_gradient_rate_11;
193        solver.root(
194            |deformation_gradient_rate: &DeformationGradientRate| {
195                Ok(self.first_piola_kirchhoff_stress(
196                    deformation_gradient,
197                    deformation_gradient_rate,
198                )?)
199            },
200            |deformation_gradient_rate: &DeformationGradientRate| {
201                Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
202                    deformation_gradient,
203                    deformation_gradient_rate,
204                )?)
205            },
206            DeformationGradientRate::zero(),
207            EqualityConstraint::Linear(matrix, vector),
208        )
209    }
210    #[doc(hidden)]
211    fn root_biaxial_inner(
212        &self,
213        deformation_gradient: &DeformationGradient,
214        deformation_gradient_rate_11: Scalar,
215        deformation_gradient_rate_22: Scalar,
216    ) -> Result<DeformationGradientRate, OptimizeError> {
217        let solver = NewtonRaphson {
218            ..Default::default()
219        };
220        let mut matrix = Matrix::zero(5, 9);
221        let mut vector = Vector::zero(5);
222        matrix[0][0] = 1.0;
223        matrix[1][1] = 1.0;
224        matrix[2][2] = 1.0;
225        matrix[3][5] = 1.0;
226        matrix[4][4] = 1.0;
227        vector[0] = deformation_gradient_rate_11;
228        vector[4] = deformation_gradient_rate_22;
229        solver.root(
230            |deformation_gradient_rate: &DeformationGradientRate| {
231                Ok(self.first_piola_kirchhoff_stress(
232                    deformation_gradient,
233                    deformation_gradient_rate,
234                )?)
235            },
236            |deformation_gradient_rate: &DeformationGradientRate| {
237                Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
238                    deformation_gradient,
239                    deformation_gradient_rate,
240                )?)
241            },
242            DeformationGradientRate::zero(),
243            EqualityConstraint::Linear(matrix, vector),
244        )
245    }
246}
conspire/constitutive/solid/viscoelastic/mod.rs

conspire/constitutive/solid/viscoelastic/
mod.rs
