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

conspire/constitutive/solid/viscoelastic/
mod.rs
