conspire/constitutive/solid/viscoelastic/
mod.rs

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

conspire/constitutive/solid/viscoelastic/
mod.rs
