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    super::fluid::viscous::Viscous,
30    viscoelastic::{AppliedLoad, Viscoelastic},
31    *,
32};
33use crate::math::{
34    Matrix, TensorVec, Vector,
35    integrate::{Explicit, IntegrationError},
36    optimize::{EqualityConstraint, NewtonRaphson, OptimizeError, SecondOrderOptimization},
37};
38use std::fmt::Debug;
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    /// Solve for the unknown components of the deformation gradient and rate under an applied load.
71    ///
72    /// ```math
73    /// \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}}
74    /// ```
75    fn minimize(
76        &self,
77        applied_load: AppliedLoad,
78        integrator: impl Explicit<DeformationGradientRate, DeformationGradientRates>,
79    ) -> Result<(Times, DeformationGradients, DeformationGradientRates), IntegrationError> {
80        match applied_load {
81            AppliedLoad::UniaxialStress(deformation_gradient_rate_11, time) => integrator
82                .integrate(
83                    |t: Scalar, deformation_gradient: &DeformationGradient| {
84                        Ok(self.minimize_uniaxial_inner(
85                            deformation_gradient,
86                            deformation_gradient_rate_11(t),
87                        )?)
88                    },
89                    time,
90                    DeformationGradient::identity(),
91                ),
92            AppliedLoad::BiaxialStress(
93                deformation_gradient_rate_11,
94                deformation_gradient_rate_22,
95                time,
96            ) => integrator.integrate(
97                |t: Scalar, deformation_gradient: &DeformationGradient| {
98                    Ok(self.minimize_biaxial_inner(
99                        deformation_gradient,
100                        deformation_gradient_rate_11(t),
101                        deformation_gradient_rate_22(t),
102                    )?)
103                },
104                time,
105                DeformationGradient::identity(),
106            ),
107        }
108    }
109    #[doc(hidden)]
110    fn minimize_uniaxial_inner(
111        &self,
112        deformation_gradient: &DeformationGradient,
113        deformation_gradient_rate_11: Scalar,
114    ) -> Result<DeformationGradientRate, OptimizeError> {
115        let solver = NewtonRaphson {
116            ..Default::default()
117        };
118        let mut matrix = Matrix::zero(4, 9);
119        let mut vector = Vector::zero(4);
120        matrix[0][0] = 1.0;
121        matrix[1][1] = 1.0;
122        matrix[2][2] = 1.0;
123        matrix[3][5] = 1.0;
124        vector[0] = deformation_gradient_rate_11;
125        solver.minimize(
126            |deformation_gradient_rate: &DeformationGradientRate| {
127                Ok(self.dissipation_potential(deformation_gradient, deformation_gradient_rate)?)
128            },
129            |deformation_gradient_rate: &DeformationGradientRate| {
130                Ok(self.first_piola_kirchhoff_stress(
131                    deformation_gradient,
132                    deformation_gradient_rate,
133                )?)
134            },
135            |deformation_gradient_rate: &DeformationGradientRate| {
136                Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
137                    deformation_gradient,
138                    deformation_gradient_rate,
139                )?)
140            },
141            DeformationGradientRate::zero(),
142            EqualityConstraint::Linear(matrix, vector),
143            None,
144        )
145    }
146    #[doc(hidden)]
147    fn minimize_biaxial_inner(
148        &self,
149        deformation_gradient: &DeformationGradient,
150        deformation_gradient_rate_11: Scalar,
151        deformation_gradient_rate_22: Scalar,
152    ) -> Result<DeformationGradientRate, OptimizeError> {
153        let solver = NewtonRaphson {
154            ..Default::default()
155        };
156        let mut matrix = Matrix::zero(5, 9);
157        let mut vector = Vector::zero(5);
158        matrix[0][0] = 1.0;
159        matrix[1][1] = 1.0;
160        matrix[2][2] = 1.0;
161        matrix[3][5] = 1.0;
162        matrix[4][4] = 1.0;
163        vector[0] = deformation_gradient_rate_11;
164        vector[4] = deformation_gradient_rate_22;
165        solver.minimize(
166            |deformation_gradient_rate: &DeformationGradientRate| {
167                Ok(self.dissipation_potential(deformation_gradient, deformation_gradient_rate)?)
168            },
169            |deformation_gradient_rate: &DeformationGradientRate| {
170                Ok(self.first_piola_kirchhoff_stress(
171                    deformation_gradient,
172                    deformation_gradient_rate,
173                )?)
174            },
175            |deformation_gradient_rate: &DeformationGradientRate| {
176                Ok(self.first_piola_kirchhoff_rate_tangent_stiffness(
177                    deformation_gradient,
178                    deformation_gradient_rate,
179                )?)
180            },
181            DeformationGradientRate::zero(),
182            EqualityConstraint::Linear(matrix, vector),
183            None,
184        )
185    }
186}
conspire/constitutive/solid/elastic_hyperviscous/mod.rs

conspire/constitutive/solid/elastic_hyperviscous/
mod.rs
