conspire/constitutive/solid/elastic_hyperviscous/
mod.rs

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

conspire/constitutive/solid/elastic_hyperviscous/
mod.rs
