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::{
18    math::{
19        Matrix, Quantity, Vector,
20        integrate::{ImplicitDaeFirstOrderRoot, ImplicitDaeZerothOrderRoot},
21        optimize::{EqualityConstraint, FirstOrderRootFinding, ZerothOrderRootFinding},
22    },
23    units::Time,
24};
25
26/// Possible applied loads.
27pub enum AppliedLoad<'a> {
28    /// Uniaxial stress given $`\dot{F}_{11}`$.
29    UniaxialStress(fn(Quantity<Time>) -> Scalar, &'a [Quantity<Time>]),
30    /// Biaxial stress given $`\dot{F}_{11}`$ and $`\dot{F}_{22}`$.
31    BiaxialStress(
32        fn(Quantity<Time>) -> Scalar,
33        fn(Quantity<Time>) -> Scalar,
34        &'a [Quantity<Time>],
35    ),
36}
37
38/// Required methods for viscoelastic solid constitutive models.
39pub trait Viscoelastic
40where
41    Self: Solid + Viscous,
42{
43    /// Calculates and returns the Cauchy stress.
44    ///
45    /// ```math
46    /// \boldsymbol{\sigma} = J^{-1}\mathbf{P}\cdot\mathbf{F}^T
47    /// ```
48    fn cauchy_stress(
49        &self,
50        deformation_gradient: &DeformationGradient,
51        deformation_gradient_rate: &DeformationGradientRate,
52    ) -> Result<CauchyStress, ConstitutiveError> {
53        Ok(deformation_gradient
54            * self
55                .second_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_rate)?
56            * deformation_gradient.transpose()
57            / deformation_gradient.determinant())
58    }
59    /// Calculates and returns the rate tangent stiffness associated with the Cauchy stress.
60    ///
61    /// ```math
62    /// \mathcal{V}_{ijkL} = \frac{\partial\sigma_{ij}}{\partial\dot{F}_{kL}} = J^{-1} \mathcal{W}_{MNkL} F_{iM} F_{jN}
63    /// ```
64    fn cauchy_rate_tangent_stiffness(
65        &self,
66        deformation_gradient: &DeformationGradient,
67        deformation_gradient_rate: &DeformationGradientRate,
68    ) -> Result<CauchyRateTangentStiffness, ConstitutiveError> {
69        Ok(self
70            .second_piola_kirchhoff_rate_tangent_stiffness(
71                deformation_gradient,
72                deformation_gradient_rate,
73            )?
74            .contract_first_second_with_second(deformation_gradient, deformation_gradient)
75            / deformation_gradient.determinant())
76    }
77    /// Calculates and returns the first Piola-Kirchhoff stress.
78    ///
79    /// ```math
80    /// \mathbf{P} = J\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}
81    /// ```
82    fn first_piola_kirchhoff_stress(
83        &self,
84        deformation_gradient: &DeformationGradient,
85        deformation_gradient_rate: &DeformationGradientRate,
86    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
87        Ok(
88            self.cauchy_stress(deformation_gradient, deformation_gradient_rate)?
89                * deformation_gradient.inverse_transpose()
90                * deformation_gradient.determinant(),
91        )
92    }
93    /// Calculates and returns the rate tangent stiffness associated with the first Piola-Kirchhoff stress.
94    ///
95    /// ```math
96    /// \mathcal{U}_{iJkL} = \frac{\partial P_{iJ}}{\partial\dot{F}_{kL}} = J \mathcal{V}_{iskL} F_{sJ}^{-T}
97    /// ```
98    fn first_piola_kirchhoff_rate_tangent_stiffness(
99        &self,
100        deformation_gradient: &DeformationGradient,
101        deformation_gradient_rate: &DeformationGradientRate,
102    ) -> Result<FirstPiolaKirchhoffRateTangentStiffness, ConstitutiveError> {
103        Ok(self
104            .cauchy_rate_tangent_stiffness(deformation_gradient, deformation_gradient_rate)?
105            .contract_second_with_first(&deformation_gradient.inverse_transpose())
106            * deformation_gradient.determinant())
107    }
108    /// Calculates and returns the second Piola-Kirchhoff stress.
109    ///
110    /// ```math
111    /// \mathbf{S} = \mathbf{F}^{-1}\cdot\mathbf{P}
112    /// ```
113    fn second_piola_kirchhoff_stress(
114        &self,
115        deformation_gradient: &DeformationGradient,
116        deformation_gradient_rate: &DeformationGradientRate,
117    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
118        Ok(deformation_gradient.inverse()
119            * self.cauchy_stress(deformation_gradient, deformation_gradient_rate)?
120            * deformation_gradient.inverse_transpose()
121            * deformation_gradient.determinant())
122    }
123    /// Calculates and returns the rate tangent stiffness associated with the second Piola-Kirchhoff stress.
124    ///
125    /// ```math
126    /// \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}
127    /// ```
128    fn second_piola_kirchhoff_rate_tangent_stiffness(
129        &self,
130        deformation_gradient: &DeformationGradient,
131        deformation_gradient_rate: &DeformationGradientRate,
132    ) -> Result<SecondPiolaKirchhoffRateTangentStiffness, ConstitutiveError> {
133        let deformation_gradient_inverse = deformation_gradient.inverse();
134        Ok(self
135            .cauchy_rate_tangent_stiffness(deformation_gradient, deformation_gradient_rate)?
136            .contract_first_second_with_second(
137                &deformation_gradient_inverse,
138                &deformation_gradient_inverse,
139            )
140            * deformation_gradient.determinant())
141    }
142}
143
144/// Zeroth-order root-finding methods for viscoelastic solid constitutive models.
145pub trait ZerothOrderRoot {
146    /// Solve for the unknown components of the deformation gradient and rate under an applied load.
147    ///
148    /// ```math
149    /// \mathbf{P}(\mathbf{F},\dot{\mathbf{F}}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
150    /// ```
151    fn root(
152        &self,
153        applied_load: AppliedLoad,
154        integrator: impl ImplicitDaeZerothOrderRoot<
155            FirstPiolaKirchhoffStress,
156            DeformationGradient,
157            DeformationGradients,
158            DeformationGradientRates,
159        >,
160        solver: impl ZerothOrderRootFinding<FirstPiolaKirchhoffStress, DeformationGradientRate>,
161    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError>;
162}
163
164/// Zeroth-order root-finding methods for viscoelastic solid constitutive models.
165pub trait FirstOrderRoot {
166    /// Solve for the unknown components of the deformation gradient and rate under an applied load.
167    ///
168    /// ```math
169    /// \mathbf{P}(\mathbf{F},\dot{\mathbf{F}}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
170    /// ```
171    fn root(
172        &self,
173        applied_load: AppliedLoad,
174        integrator: impl ImplicitDaeFirstOrderRoot<
175            FirstPiolaKirchhoffStress,
176            FirstPiolaKirchhoffRateTangentStiffness,
177            DeformationGradient,
178            DeformationGradients,
179            DeformationGradientRates,
180        >,
181        solver: impl FirstOrderRootFinding<
182            FirstPiolaKirchhoffStress,
183            FirstPiolaKirchhoffRateTangentStiffness,
184            DeformationGradientRate,
185        >,
186    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError>;
187}
188
189impl<T> ZerothOrderRoot for T
190where
191    T: Viscoelastic,
192{
193    fn root(
194        &self,
195        applied_load: AppliedLoad,
196        integrator: impl ImplicitDaeZerothOrderRoot<
197            FirstPiolaKirchhoffStress,
198            DeformationGradient,
199            DeformationGradients,
200            DeformationGradientRates,
201        >,
202        solver: impl ZerothOrderRootFinding<FirstPiolaKirchhoffStress, DeformationGradientRate>,
203    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError> {
204        match match applied_load {
205            AppliedLoad::UniaxialStress(deformation_gradient_rate_11, time) => {
206                let mut matrix = Matrix::zero(4, 9);
207                let mut vector = Vector::zero(4);
208                matrix[0][0] = 1.0;
209                matrix[1][1] = 1.0;
210                matrix[2][2] = 1.0;
211                matrix[3][5] = 1.0;
212                integrator.integrate(
213                    |_: Quantity<Time>,
214                     deformation_gradient: &DeformationGradient,
215                     deformation_gradient_rate: &DeformationGradientRate| {
216                        Ok(self.first_piola_kirchhoff_stress(
217                            deformation_gradient,
218                            deformation_gradient_rate,
219                        )?)
220                    },
221                    solver,
222                    time,
223                    DeformationGradient::identity(),
224                    |t: Quantity<Time>| {
225                        vector[0] = deformation_gradient_rate_11(t);
226                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
227                    },
228                )
229            }
230            AppliedLoad::BiaxialStress(
231                deformation_gradient_rate_11,
232                deformation_gradient_rate_22,
233                time,
234            ) => {
235                let mut matrix = Matrix::zero(5, 9);
236                let mut vector = Vector::zero(5);
237                matrix[0][0] = 1.0;
238                matrix[1][1] = 1.0;
239                matrix[2][2] = 1.0;
240                matrix[3][5] = 1.0;
241                matrix[4][4] = 1.0;
242                integrator.integrate(
243                    |_: Quantity<Time>,
244                     deformation_gradient: &DeformationGradient,
245                     deformation_gradient_rate: &DeformationGradientRate| {
246                        Ok(self.first_piola_kirchhoff_stress(
247                            deformation_gradient,
248                            deformation_gradient_rate,
249                        )?)
250                    },
251                    solver,
252                    time,
253                    DeformationGradient::identity(),
254                    |t: Quantity<Time>| {
255                        vector[0] = deformation_gradient_rate_11(t);
256                        vector[4] = deformation_gradient_rate_22(t);
257                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
258                    },
259                )
260            }
261        } {
262            Ok(results) => Ok(results),
263            Err(error) => Err(ConstitutiveError::Upstream(
264                format!("{error}"),
265                format!("{self:?}"),
266            )),
267        }
268    }
269}
270
271impl<T> FirstOrderRoot for T
272where
273    T: Viscoelastic,
274{
275    fn root(
276        &self,
277        applied_load: AppliedLoad,
278        integrator: impl ImplicitDaeFirstOrderRoot<
279            FirstPiolaKirchhoffStress,
280            FirstPiolaKirchhoffRateTangentStiffness,
281            DeformationGradient,
282            DeformationGradients,
283            DeformationGradientRates,
284        >,
285        solver: impl FirstOrderRootFinding<
286            FirstPiolaKirchhoffStress,
287            FirstPiolaKirchhoffRateTangentStiffness,
288            DeformationGradientRate,
289        >,
290    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError> {
291        match match applied_load {
292            AppliedLoad::UniaxialStress(deformation_gradient_rate_11, time) => {
293                let mut matrix = Matrix::zero(4, 9);
294                let mut vector = Vector::zero(4);
295                matrix[0][0] = 1.0;
296                matrix[1][1] = 1.0;
297                matrix[2][2] = 1.0;
298                matrix[3][5] = 1.0;
299                integrator.integrate(
300                    |_: Quantity<Time>,
301                     deformation_gradient: &DeformationGradient,
302                     deformation_gradient_rate: &DeformationGradientRate| {
303                        Ok(self.first_piola_kirchhoff_stress(
304                            deformation_gradient,
305                            deformation_gradient_rate,
306                        )?)
307                    },
308                    |_: Quantity<Time>,
309                     deformation_gradient: &DeformationGradient,
310                     deformation_gradient_rate: &DeformationGradientRate| {
311                        Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
312                            deformation_gradient,
313                            deformation_gradient_rate,
314                        )?)
315                    },
316                    solver,
317                    time,
318                    DeformationGradient::identity(),
319                    |t: Quantity<Time>| {
320                        vector[0] = deformation_gradient_rate_11(t);
321                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
322                    },
323                )
324            }
325            AppliedLoad::BiaxialStress(
326                deformation_gradient_rate_11,
327                deformation_gradient_rate_22,
328                time,
329            ) => {
330                let mut matrix = Matrix::zero(5, 9);
331                let mut vector = Vector::zero(5);
332                matrix[0][0] = 1.0;
333                matrix[1][1] = 1.0;
334                matrix[2][2] = 1.0;
335                matrix[3][5] = 1.0;
336                matrix[4][4] = 1.0;
337                integrator.integrate(
338                    |_: Quantity<Time>,
339                     deformation_gradient: &DeformationGradient,
340                     deformation_gradient_rate: &DeformationGradientRate| {
341                        Ok(self.first_piola_kirchhoff_stress(
342                            deformation_gradient,
343                            deformation_gradient_rate,
344                        )?)
345                    },
346                    |_: Quantity<Time>,
347                     deformation_gradient: &DeformationGradient,
348                     deformation_gradient_rate: &DeformationGradientRate| {
349                        Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
350                            deformation_gradient,
351                            deformation_gradient_rate,
352                        )?)
353                    },
354                    solver,
355                    time,
356                    DeformationGradient::identity(),
357                    |t: Quantity<Time>| {
358                        vector[0] = deformation_gradient_rate_11(t);
359                        vector[4] = deformation_gradient_rate_22(t);
360                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
361                    },
362                )
363            }
364        } {
365            Ok(results) => Ok(results),
366            Err(error) => Err(ConstitutiveError::Upstream(
367                format!("{error}"),
368                format!("{self:?}"),
369            )),
370        }
371    }
372}
conspire/constitutive/solid/viscoelastic/mod.rs

conspire/constitutive/solid/viscoelastic/
mod.rs
