conspire/constitutive/solid/elastic_hyperviscous/
mod.rs

1//! Elastic-hyperviscous solid constitutive models.
2//!
3//! ---
4//!
5//! Elastic-hyperviscous solid constitutive models are defined by an elastic stress tensor function and a viscous dissipation function.
6//!
7//! ```math
8//! \mathbf{P}:\dot{\mathbf{F}} - \mathbf{P}^e(\mathbf{F}):\dot{\mathbf{F}} - \phi(\mathbf{F},\dot{\mathbf{F}}) \geq 0
9//! ```
10//! Satisfying the second law of thermodynamics though a minimum viscous dissipation principal yields a relation for the stress.
11//!
12//! ```math
13//! \mathbf{P} = \mathbf{P}^e + \frac{\partial\phi}{\partial\dot{\mathbf{F}}}
14//! ```
15//! Consequently, the rate tangent stiffness associated with the first Piola-Kirchhoff stress is symmetric for these constitutive models.
16//!
17//! ```math
18//! \mathcal{U}_{iJkL} = \mathcal{U}_{kLiJ}
19//! ```
20
21#[cfg(test)]
22pub mod test;
23
24mod almansi_hamel;
25
26pub use almansi_hamel::AlmansiHamel;
27
28use super::{
29    viscoelastic::{AppliedLoad, Viscoelastic},
30    *,
31};
32use crate::math::{
33    Matrix, Vector,
34    integrate::{ImplicitDaeFirstOrderMinimize, ImplicitDaeSecondOrderMinimize},
35    optimize::{EqualityConstraint, FirstOrderOptimization, SecondOrderOptimization},
36};
37
38/// Required methods for elastic-hyperviscous solid constitutive models.
39pub trait ElasticHyperviscous
40where
41    Self: Viscoelastic,
42{
43    /// Calculates and returns the dissipation potential.
44    ///
45    /// ```math
46    /// \mathbf{P}^e(\mathbf{F}):\dot{\mathbf{F}} + \phi(\mathbf{F},\dot{\mathbf{F}})
47    /// ```
48    fn dissipation_potential(
49        &self,
50        deformation_gradient: &DeformationGradient,
51        deformation_gradient_rate: &DeformationGradientRate,
52    ) -> Result<Scalar, ConstitutiveError> {
53        Ok(self
54            .first_piola_kirchhoff_stress(deformation_gradient, &ZERO_10)?
55            .full_contraction(deformation_gradient_rate)
56            + self.viscous_dissipation(deformation_gradient, deformation_gradient_rate)?)
57    }
58    /// Calculates and returns the viscous dissipation.
59    ///
60    /// ```math
61    /// \phi = \phi(\mathbf{F},\dot{\mathbf{F}})
62    /// ```
63    fn viscous_dissipation(
64        &self,
65        deformation_gradient: &DeformationGradient,
66        deformation_gradient_rate: &DeformationGradientRate,
67    ) -> Result<Scalar, ConstitutiveError>;
68}
69
70/// First-order optimization methods for elastic-hyperviscous solid constitutive models.
71pub trait FirstOrderMinimize {
72    /// Solve for the unknown components of the deformation gradient and rate under an applied load.
73    ///
74    /// ```math
75    /// \Pi(\mathbf{F},\dot{\mathbf{F}},\boldsymbol{\lambda}) = \mathbf{P}^e(\mathbf{F}):\dot{\mathbf{F}} + \phi(\mathbf{F},\dot{\mathbf{F}}) - \boldsymbol{\lambda}:(\dot{\mathbf{F}} - \dot{\mathbf{F}}_0) - \mathbf{P}_0:\dot{\mathbf{F}}
76    /// ```
77    fn minimize(
78        &self,
79        applied_load: AppliedLoad,
80        integrator: impl ImplicitDaeFirstOrderMinimize<
81            Scalar,
82            DeformationGradientRate,
83            DeformationGradientRates,
84        >,
85        solver: impl FirstOrderOptimization<Scalar, DeformationGradient>,
86    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError>;
87}
88
89/// Second-order optimization methods for elastic-hyperviscous solid constitutive models.
90pub trait SecondOrderMinimize {
91    /// Solve for the unknown components of the deformation gradient and rate under an applied load.
92    ///
93    /// ```math
94    /// \Pi(\mathbf{F},\dot{\mathbf{F}},\boldsymbol{\lambda}) = \mathbf{P}^e(\mathbf{F}):\dot{\mathbf{F}} + \phi(\mathbf{F},\dot{\mathbf{F}}) - \boldsymbol{\lambda}:(\dot{\mathbf{F}} - \dot{\mathbf{F}}_0) - \mathbf{P}_0:\dot{\mathbf{F}}
95    /// ```
96    fn minimize(
97        &self,
98        applied_load: AppliedLoad,
99        integrator: impl ImplicitDaeSecondOrderMinimize<
100            Scalar,
101            FirstPiolaKirchhoffStress,
102            FirstPiolaKirchhoffRateTangentStiffness,
103            DeformationGradientRate,
104            DeformationGradientRates,
105        >,
106        solver: impl SecondOrderOptimization<
107            Scalar,
108            FirstPiolaKirchhoffStress,
109            FirstPiolaKirchhoffRateTangentStiffness,
110            DeformationGradient,
111        >,
112    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError>;
113}
114
115impl<T> FirstOrderMinimize for T
116where
117    T: ElasticHyperviscous,
118{
119    fn minimize(
120        &self,
121        applied_load: AppliedLoad,
122        integrator: impl ImplicitDaeFirstOrderMinimize<
123            Scalar,
124            DeformationGradientRate,
125            DeformationGradientRates,
126        >,
127        solver: impl FirstOrderOptimization<Scalar, DeformationGradientRate>,
128    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError> {
129        match match applied_load {
130            AppliedLoad::UniaxialStress(deformation_gradient_rate_11, time) => {
131                let mut matrix = Matrix::zero(4, 9);
132                let mut vector = Vector::zero(4);
133                matrix[0][0] = 1.0;
134                matrix[1][1] = 1.0;
135                matrix[2][2] = 1.0;
136                matrix[3][5] = 1.0;
137                integrator.integrate(
138                    |_: Scalar,
139                     deformation_gradient: &DeformationGradient,
140                     deformation_gradient_rate: &DeformationGradientRate| {
141                        Ok(self.dissipation_potential(
142                            deformation_gradient,
143                            deformation_gradient_rate,
144                        )?)
145                    },
146                    |_: Scalar,
147                     deformation_gradient: &DeformationGradient,
148                     deformation_gradient_rate: &DeformationGradientRate| {
149                        Ok(self.first_piola_kirchhoff_stress(
150                            deformation_gradient,
151                            deformation_gradient_rate,
152                        )?)
153                    },
154                    solver,
155                    time,
156                    DeformationGradient::identity(),
157                    |t: Scalar| {
158                        vector[0] = deformation_gradient_rate_11(t);
159                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
160                    },
161                )
162            }
163            AppliedLoad::BiaxialStress(
164                deformation_gradient_rate_11,
165                deformation_gradient_rate_22,
166                time,
167            ) => {
168                let mut matrix = Matrix::zero(5, 9);
169                let mut vector = Vector::zero(5);
170                matrix[0][0] = 1.0;
171                matrix[1][1] = 1.0;
172                matrix[2][2] = 1.0;
173                matrix[3][5] = 1.0;
174                matrix[4][4] = 1.0;
175                integrator.integrate(
176                    |_: Scalar,
177                     deformation_gradient: &DeformationGradient,
178                     deformation_gradient_rate: &DeformationGradientRate| {
179                        Ok(self.dissipation_potential(
180                            deformation_gradient,
181                            deformation_gradient_rate,
182                        )?)
183                    },
184                    |_: Scalar,
185                     deformation_gradient: &DeformationGradient,
186                     deformation_gradient_rate: &DeformationGradientRate| {
187                        Ok(self.first_piola_kirchhoff_stress(
188                            deformation_gradient,
189                            deformation_gradient_rate,
190                        )?)
191                    },
192                    solver,
193                    time,
194                    DeformationGradient::identity(),
195                    |t: Scalar| {
196                        vector[0] = deformation_gradient_rate_11(t);
197                        vector[4] = deformation_gradient_rate_22(t);
198                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
199                    },
200                )
201            }
202        } {
203            Ok(results) => Ok(results),
204            Err(error) => Err(ConstitutiveError::Upstream(
205                format!("{error}"),
206                format!("{self:?}"),
207            )),
208        }
209    }
210}
211
212impl<T> SecondOrderMinimize for T
213where
214    T: ElasticHyperviscous,
215{
216    fn minimize(
217        &self,
218        applied_load: AppliedLoad,
219        integrator: impl ImplicitDaeSecondOrderMinimize<
220            Scalar,
221            FirstPiolaKirchhoffStress,
222            FirstPiolaKirchhoffRateTangentStiffness,
223            DeformationGradientRate,
224            DeformationGradientRates,
225        >,
226        solver: impl SecondOrderOptimization<
227            Scalar,
228            FirstPiolaKirchhoffStress,
229            FirstPiolaKirchhoffRateTangentStiffness,
230            DeformationGradientRate,
231        >,
232    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), ConstitutiveError> {
233        match match applied_load {
234            AppliedLoad::UniaxialStress(deformation_gradient_rate_11, time) => {
235                let mut matrix = Matrix::zero(4, 9);
236                let mut vector = Vector::zero(4);
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                integrator.integrate(
242                    |_: Scalar,
243                     deformation_gradient: &DeformationGradient,
244                     deformation_gradient_rate: &DeformationGradientRate| {
245                        Ok(self.dissipation_potential(
246                            deformation_gradient,
247                            deformation_gradient_rate,
248                        )?)
249                    },
250                    |_: Scalar,
251                     deformation_gradient: &DeformationGradient,
252                     deformation_gradient_rate: &DeformationGradientRate| {
253                        Ok(self.first_piola_kirchhoff_stress(
254                            deformation_gradient,
255                            deformation_gradient_rate,
256                        )?)
257                    },
258                    |_: Scalar,
259                     deformation_gradient: &DeformationGradient,
260                     deformation_gradient_rate: &DeformationGradientRate| {
261                        Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
262                            deformation_gradient,
263                            deformation_gradient_rate,
264                        )?)
265                    },
266                    solver,
267                    time,
268                    DeformationGradient::identity(),
269                    |t: Scalar| {
270                        vector[0] = deformation_gradient_rate_11(t);
271                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
272                    },
273                    None,
274                )
275            }
276            AppliedLoad::BiaxialStress(
277                deformation_gradient_rate_11,
278                deformation_gradient_rate_22,
279                time,
280            ) => {
281                let mut matrix = Matrix::zero(5, 9);
282                let mut vector = Vector::zero(5);
283                matrix[0][0] = 1.0;
284                matrix[1][1] = 1.0;
285                matrix[2][2] = 1.0;
286                matrix[3][5] = 1.0;
287                matrix[4][4] = 1.0;
288                integrator.integrate(
289                    |_: Scalar,
290                     deformation_gradient: &DeformationGradient,
291                     deformation_gradient_rate: &DeformationGradientRate| {
292                        Ok(self.dissipation_potential(
293                            deformation_gradient,
294                            deformation_gradient_rate,
295                        )?)
296                    },
297                    |_: Scalar,
298                     deformation_gradient: &DeformationGradient,
299                     deformation_gradient_rate: &DeformationGradientRate| {
300                        Ok(self.first_piola_kirchhoff_stress(
301                            deformation_gradient,
302                            deformation_gradient_rate,
303                        )?)
304                    },
305                    |_: Scalar,
306                     deformation_gradient: &DeformationGradient,
307                     deformation_gradient_rate: &DeformationGradientRate| {
308                        Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
309                            deformation_gradient,
310                            deformation_gradient_rate,
311                        )?)
312                    },
313                    solver,
314                    time,
315                    DeformationGradient::identity(),
316                    |t: Scalar| {
317                        vector[0] = deformation_gradient_rate_11(t);
318                        vector[4] = deformation_gradient_rate_22(t);
319                        EqualityConstraint::Linear(matrix.clone(), vector.clone())
320                    },
321                    None,
322                )
323            }
324        } {
325            Ok(results) => Ok(results),
326            Err(error) => Err(ConstitutiveError::Upstream(
327                format!("{error}"),
328                format!("{self:?}"),
329            )),
330        }
331    }
332}
conspire/constitutive/solid/elastic_hyperviscous/mod.rs

conspire/constitutive/solid/elastic_hyperviscous/
mod.rs
