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, 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        solver: impl FirstOrderRootFinding<
148            FirstPiolaKirchhoffStress,
149            FirstPiolaKirchhoffTangentStiffness,
150            DeformationGradient,
151        >,
152    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), IntegrationError> {
153        let mut solution = DeformationGradientRate::zero();
154        match applied_load {
155            AppliedLoad::UniaxialStress(deformation_gradient_rate_11, time) => integrator
156                .integrate(
157                    |t: Scalar, deformation_gradient: &DeformationGradient| {
158                        solution = self.root_uniaxial_inner(
159                            deformation_gradient,
160                            deformation_gradient_rate_11(t),
161                            &solver,
162                            &solution,
163                        )?;
164                        Ok(solution.clone())
165                    },
166                    time,
167                    DeformationGradient::identity(),
168                ),
169            AppliedLoad::BiaxialStress(
170                deformation_gradient_rate_11,
171                deformation_gradient_rate_22,
172                time,
173            ) => integrator.integrate(
174                |t: Scalar, deformation_gradient: &DeformationGradient| {
175                    solution = self.root_biaxial_inner(
176                        deformation_gradient,
177                        deformation_gradient_rate_11(t),
178                        deformation_gradient_rate_22(t),
179                        &solver,
180                        &solution,
181                    )?;
182                    Ok(solution.clone())
183                },
184                time,
185                DeformationGradient::identity(),
186            ),
187        }
188    }
189    #[doc(hidden)]
190    fn root_uniaxial_inner(
191        &self,
192        deformation_gradient: &DeformationGradient,
193        deformation_gradient_rate_11: Scalar,
194        solver: &impl FirstOrderRootFinding<
195            FirstPiolaKirchhoffStress,
196            FirstPiolaKirchhoffTangentStiffness,
197            DeformationGradient,
198        >,
199        initial_guess: &DeformationGradientRate,
200    ) -> Result<DeformationGradientRate, OptimizeError> {
201        let mut matrix = Matrix::zero(4, 9);
202        let mut vector = Vector::zero(4);
203        matrix[0][0] = 1.0;
204        matrix[1][1] = 1.0;
205        matrix[2][2] = 1.0;
206        matrix[3][5] = 1.0;
207        vector[0] = deformation_gradient_rate_11;
208        solver.root(
209            |deformation_gradient_rate: &DeformationGradientRate| {
210                Ok(self.first_piola_kirchhoff_stress(
211                    deformation_gradient,
212                    deformation_gradient_rate,
213                )?)
214            },
215            |deformation_gradient_rate: &DeformationGradientRate| {
216                Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
217                    deformation_gradient,
218                    deformation_gradient_rate,
219                )?)
220            },
221            initial_guess.clone(),
222            EqualityConstraint::Linear(matrix, vector),
223        )
224    }
225    #[doc(hidden)]
226    fn root_biaxial_inner(
227        &self,
228        deformation_gradient: &DeformationGradient,
229        deformation_gradient_rate_11: Scalar,
230        deformation_gradient_rate_22: Scalar,
231        solver: &impl FirstOrderRootFinding<
232            FirstPiolaKirchhoffStress,
233            FirstPiolaKirchhoffTangentStiffness,
234            DeformationGradient,
235        >,
236        initial_guess: &DeformationGradientRate,
237    ) -> Result<DeformationGradientRate, OptimizeError> {
238        let mut matrix = Matrix::zero(5, 9);
239        let mut vector = Vector::zero(5);
240        matrix[0][0] = 1.0;
241        matrix[1][1] = 1.0;
242        matrix[2][2] = 1.0;
243        matrix[3][5] = 1.0;
244        matrix[4][4] = 1.0;
245        vector[0] = deformation_gradient_rate_11;
246        vector[4] = deformation_gradient_rate_22;
247        solver.root(
248            |deformation_gradient_rate: &DeformationGradientRate| {
249                Ok(self.first_piola_kirchhoff_stress(
250                    deformation_gradient,
251                    deformation_gradient_rate,
252                )?)
253            },
254            |deformation_gradient_rate: &DeformationGradientRate| {
255                Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
256                    deformation_gradient,
257                    deformation_gradient_rate,
258                )?)
259            },
260            initial_guess.clone(),
261            EqualityConstraint::Linear(matrix, vector),
262        )
263    }
264}
conspire/constitutive/solid/viscoelastic/mod.rs

conspire/constitutive/solid/viscoelastic/
mod.rs
