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::ImplicitDaeRoot,
20        optimize::{EqualityConstraint, RootFinding},
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/// Root-finding methods for viscoelastic solid constitutive models.
144pub trait Root {
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 ImplicitDaeRoot<
154            FirstPiolaKirchhoffStress,
155            FirstPiolaKirchhoffRateTangentStiffness,
156            DeformationGradient,
157            DeformationGradients,
158            DeformationGradientRates,
159        >,
160        solver: impl RootFinding<
161            FirstPiolaKirchhoffStress,
162            FirstPiolaKirchhoffRateTangentStiffness,
163            DeformationGradientRate,
164        >,
165    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError>;
166}
167
168impl<T> Root for T
169where
170    T: Viscoelastic,
171{
172    fn root(
173        &self,
174        applied_load: AppliedLoad,
175        integrator: impl ImplicitDaeRoot<
176            FirstPiolaKirchhoffStress,
177            FirstPiolaKirchhoffRateTangentStiffness,
178            DeformationGradient,
179            DeformationGradients,
180            DeformationGradientRates,
181        >,
182        solver: impl RootFinding<
183            FirstPiolaKirchhoffStress,
184            FirstPiolaKirchhoffRateTangentStiffness,
185            DeformationGradientRate,
186        >,
187    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError> {
188        match applied_load {
189            AppliedLoad::UniaxialStress(deformation_gradient_rate_11, time) => {
190                let mut matrix = Matrix::zero(4, 9);
191                let mut vector = Vector::zero(4);
192                matrix[0][0] = 1.0;
193                matrix[1][1] = 1.0;
194                matrix[2][2] = 1.0;
195                matrix[3][5] = 1.0;
196                integrator.integrate(
197                    |_: Quantity<Time>,
198                     deformation_gradient: &DeformationGradient,
199                     deformation_gradient_rate: &DeformationGradientRate| {
200                        Ok(self.first_piola_kirchhoff_stress(
201                            deformation_gradient,
202                            deformation_gradient_rate,
203                        )?)
204                    },
205                    |_: Quantity<Time>,
206                     deformation_gradient: &DeformationGradient,
207                     deformation_gradient_rate: &DeformationGradientRate| {
208                        Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
209                            deformation_gradient,
210                            deformation_gradient_rate,
211                        )?)
212                    },
213                    solver,
214                    time,
215                    DeformationGradient::identity(),
216                    |t: Quantity<Time>| {
217                        vector[0] = deformation_gradient_rate_11(t);
218                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
219                    },
220                )
221            }
222            AppliedLoad::BiaxialStress(
223                deformation_gradient_rate_11,
224                deformation_gradient_rate_22,
225                time,
226            ) => {
227                let mut matrix = Matrix::zero(5, 9);
228                let mut vector = Vector::zero(5);
229                matrix[0][0] = 1.0;
230                matrix[1][1] = 1.0;
231                matrix[2][2] = 1.0;
232                matrix[3][5] = 1.0;
233                matrix[4][4] = 1.0;
234                integrator.integrate(
235                    |_: Quantity<Time>,
236                     deformation_gradient: &DeformationGradient,
237                     deformation_gradient_rate: &DeformationGradientRate| {
238                        Ok(self.first_piola_kirchhoff_stress(
239                            deformation_gradient,
240                            deformation_gradient_rate,
241                        )?)
242                    },
243                    |_: Quantity<Time>,
244                     deformation_gradient: &DeformationGradient,
245                     deformation_gradient_rate: &DeformationGradientRate| {
246                        Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
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        .map_err(|error| ConstitutiveError::upstream(error, self))
263    }
264}
conspire/constitutive/solid/viscoelastic/mod.rs

conspire/constitutive/solid/viscoelastic/
mod.rs
