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

conspire/constitutive/solid/viscoelastic/
mod.rs
