conspire/constitutive/solid/viscoelastic/
mod.rs

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

conspire/constitutive/solid/viscoelastic/
mod.rs
