conspire/constitutive/solid/viscoelastic/
mod.rs

1//! Viscoelastic solid constitutive models.
2//!
3//! ---
4//!
5//! Viscoelastic solid constitutive models cannot be defined by a Helmholtz free energy density and a 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, Vector,
19    integrate::{ImplicitDaeFirstOrderRoot, ImplicitDaeZerothOrderRoot},
20    optimize::{EqualityConstraint, FirstOrderRootFinding, ZerothOrderRootFinding},
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 solid 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}
139
140/// Zeroth-order root-finding methods for viscoelastic solid constitutive models.
141pub trait ZerothOrderRoot {
142    /// Solve for the unknown components of the deformation gradient and rate under an applied load.
143    ///
144    /// ```math
145    /// \mathbf{P}(\mathbf{F},\dot{\mathbf{F}}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
146    /// ```
147    fn root(
148        &self,
149        applied_load: AppliedLoad,
150        integrator: impl ImplicitDaeZerothOrderRoot<DeformationGradient, DeformationGradients>,
151        solver: impl ZerothOrderRootFinding<DeformationGradient>,
152    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError>;
153}
154
155/// Zeroth-order root-finding methods for viscoelastic solid constitutive models.
156pub trait FirstOrderRoot {
157    /// Solve for the unknown components of the deformation gradient and rate under an applied load.
158    ///
159    /// ```math
160    /// \mathbf{P}(\mathbf{F},\dot{\mathbf{F}}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
161    /// ```
162    fn root(
163        &self,
164        applied_load: AppliedLoad,
165        integrator: impl ImplicitDaeFirstOrderRoot<
166            FirstPiolaKirchhoffStress,
167            FirstPiolaKirchhoffRateTangentStiffness,
168            DeformationGradientRate,
169            DeformationGradientRates,
170        >,
171        solver: impl FirstOrderRootFinding<
172            FirstPiolaKirchhoffStress,
173            FirstPiolaKirchhoffRateTangentStiffness,
174            DeformationGradientRate,
175        >,
176    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError>;
177}
178
179impl<T> ZerothOrderRoot for T
180where
181    T: Viscoelastic,
182{
183    fn root(
184        &self,
185        applied_load: AppliedLoad,
186        integrator: impl ImplicitDaeZerothOrderRoot<DeformationGradient, DeformationGradients>,
187        solver: impl ZerothOrderRootFinding<DeformationGradientRate>,
188    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError> {
189        match match applied_load {
190            AppliedLoad::UniaxialStress(deformation_gradient_rate_11, time) => {
191                let mut matrix = Matrix::zero(4, 9);
192                let mut vector = Vector::zero(4);
193                matrix[0][0] = 1.0;
194                matrix[1][1] = 1.0;
195                matrix[2][2] = 1.0;
196                matrix[3][5] = 1.0;
197                integrator.integrate(
198                    |_: Scalar,
199                     deformation_gradient: &DeformationGradient,
200                     deformation_gradient_rate: &DeformationGradientRate| {
201                        Ok(self.first_piola_kirchhoff_stress(
202                            deformation_gradient,
203                            deformation_gradient_rate,
204                        )?)
205                    },
206                    solver,
207                    time,
208                    DeformationGradient::identity(),
209                    |t: Scalar| {
210                        vector[0] = deformation_gradient_rate_11(t);
211                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
212                    },
213                )
214            }
215            AppliedLoad::BiaxialStress(
216                deformation_gradient_rate_11,
217                deformation_gradient_rate_22,
218                time,
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                integrator.integrate(
228                    |_: Scalar,
229                     deformation_gradient: &DeformationGradient,
230                     deformation_gradient_rate: &DeformationGradientRate| {
231                        Ok(self.first_piola_kirchhoff_stress(
232                            deformation_gradient,
233                            deformation_gradient_rate,
234                        )?)
235                    },
236                    solver,
237                    time,
238                    DeformationGradient::identity(),
239                    |t: Scalar| {
240                        vector[0] = deformation_gradient_rate_11(t);
241                        vector[4] = deformation_gradient_rate_22(t);
242                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
243                    },
244                )
245            }
246        } {
247            Ok(results) => Ok(results),
248            Err(error) => Err(ConstitutiveError::Upstream(
249                format!("{error}"),
250                format!("{self:?}"),
251            )),
252        }
253    }
254}
255
256impl<T> FirstOrderRoot for T
257where
258    T: Viscoelastic,
259{
260    fn root(
261        &self,
262        applied_load: AppliedLoad,
263        integrator: impl ImplicitDaeFirstOrderRoot<
264            FirstPiolaKirchhoffStress,
265            FirstPiolaKirchhoffRateTangentStiffness,
266            DeformationGradientRate,
267            DeformationGradientRates,
268        >,
269        solver: impl FirstOrderRootFinding<
270            FirstPiolaKirchhoffStress,
271            FirstPiolaKirchhoffRateTangentStiffness,
272            DeformationGradientRate,
273        >,
274    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError> {
275        match match applied_load {
276            AppliedLoad::UniaxialStress(deformation_gradient_rate_11, time) => {
277                let mut matrix = Matrix::zero(4, 9);
278                let mut vector = Vector::zero(4);
279                matrix[0][0] = 1.0;
280                matrix[1][1] = 1.0;
281                matrix[2][2] = 1.0;
282                matrix[3][5] = 1.0;
283                integrator.integrate(
284                    |_: Scalar,
285                     deformation_gradient: &DeformationGradient,
286                     deformation_gradient_rate: &DeformationGradientRate| {
287                        Ok(self.first_piola_kirchhoff_stress(
288                            deformation_gradient,
289                            deformation_gradient_rate,
290                        )?)
291                    },
292                    |_: Scalar,
293                     deformation_gradient: &DeformationGradient,
294                     deformation_gradient_rate: &DeformationGradientRate| {
295                        Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
296                            deformation_gradient,
297                            deformation_gradient_rate,
298                        )?)
299                    },
300                    solver,
301                    time,
302                    DeformationGradient::identity(),
303                    |t: Scalar| {
304                        vector[0] = deformation_gradient_rate_11(t);
305                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
306                    },
307                )
308            }
309            AppliedLoad::BiaxialStress(
310                deformation_gradient_rate_11,
311                deformation_gradient_rate_22,
312                time,
313            ) => {
314                let mut matrix = Matrix::zero(5, 9);
315                let mut vector = Vector::zero(5);
316                matrix[0][0] = 1.0;
317                matrix[1][1] = 1.0;
318                matrix[2][2] = 1.0;
319                matrix[3][5] = 1.0;
320                matrix[4][4] = 1.0;
321                integrator.integrate(
322                    |_: Scalar,
323                     deformation_gradient: &DeformationGradient,
324                     deformation_gradient_rate: &DeformationGradientRate| {
325                        Ok(self.first_piola_kirchhoff_stress(
326                            deformation_gradient,
327                            deformation_gradient_rate,
328                        )?)
329                    },
330                    |_: Scalar,
331                     deformation_gradient: &DeformationGradient,
332                     deformation_gradient_rate: &DeformationGradientRate| {
333                        Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
334                            deformation_gradient,
335                            deformation_gradient_rate,
336                        )?)
337                    },
338                    solver,
339                    time,
340                    DeformationGradient::identity(),
341                    |t: Scalar| {
342                        vector[0] = deformation_gradient_rate_11(t);
343                        vector[4] = deformation_gradient_rate_22(t);
344                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
345                    },
346                )
347            }
348        } {
349            Ok(results) => Ok(results),
350            Err(error) => Err(ConstitutiveError::Upstream(
351                format!("{error}"),
352                format!("{self:?}"),
353            )),
354        }
355    }
356}
conspire/constitutive/solid/viscoelastic/mod.rs

conspire/constitutive/solid/viscoelastic/
mod.rs
