conspire/constitutive/solid/elastic/
mod.rs

1//! Elastic 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
13pub mod internal_variables;
14
15mod almansi_hamel;
16mod hencky;
17mod saint_venant_kirchhoff;
18
19pub use self::{
20    almansi_hamel::AlmansiHamel, hencky::Hencky, saint_venant_kirchhoff::SaintVenantKirchhoff,
21};
22
23use super::*;
24use crate::math::{
25    Matrix, Vector,
26    optimize::{EqualityConstraint, FirstOrderRootFinding, ZerothOrderRootFinding},
27};
28
29/// Possible applied loads.
30pub enum AppliedLoad {
31    /// Uniaxial stress given $`F_{11}`$.
32    UniaxialStress(Scalar),
33    /// Biaxial stress given $`F_{11}`$ and $`F_{22}`$.
34    BiaxialStress(Scalar, Scalar),
35}
36
37/// Required methods for elastic constitutive models.
38pub trait Elastic
39where
40    Self: Solid,
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    ) -> Result<CauchyStress, ConstitutiveError> {
51        Ok(deformation_gradient
52            * self.second_piola_kirchhoff_stress(deformation_gradient)?
53            * deformation_gradient.transpose()
54            / deformation_gradient.determinant())
55    }
56    /// Calculates and returns the tangent stiffness associated with the Cauchy stress.
57    ///
58    /// ```math
59    /// \mathcal{T}_{ijkL} = \frac{\partial\sigma_{ij}}{\partial F_{kL}} = J^{-1} \mathcal{G}_{MNkL} F_{iM} F_{jN} - \sigma_{ij} F_{kL}^{-T} + \left(\delta_{jk}\sigma_{is} + \delta_{ik}\sigma_{js}\right)F_{sL}^{-T}
60    /// ```
61    fn cauchy_tangent_stiffness(
62        &self,
63        deformation_gradient: &DeformationGradient,
64    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
65        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
66        let cauchy_stress = self.cauchy_stress(deformation_gradient)?;
67        let some_stress = &cauchy_stress * &deformation_gradient_inverse_transpose;
68        Ok(self
69            .second_piola_kirchhoff_tangent_stiffness(deformation_gradient)?
70            .contract_first_second_indices_with_second_indices_of(
71                deformation_gradient,
72                deformation_gradient,
73            )
74            / deformation_gradient.determinant()
75            - CauchyTangentStiffness::dyad_ij_kl(
76                &cauchy_stress,
77                &deformation_gradient_inverse_transpose,
78            )
79            + CauchyTangentStiffness::dyad_il_kj(&some_stress, &IDENTITY)
80            + CauchyTangentStiffness::dyad_ik_jl(&IDENTITY, &some_stress))
81    }
82    /// Calculates and returns the first Piola-Kirchhoff stress.
83    ///
84    /// ```math
85    /// \mathbf{P} = J\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}
86    /// ```
87    fn first_piola_kirchhoff_stress(
88        &self,
89        deformation_gradient: &DeformationGradient,
90    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
91        Ok(self.cauchy_stress(deformation_gradient)?
92            * deformation_gradient.inverse_transpose()
93            * deformation_gradient.determinant())
94    }
95    /// Calculates and returns the tangent stiffness associated with the first Piola-Kirchhoff stress.
96    ///
97    /// ```math
98    /// \mathcal{C}_{iJkL} = \frac{\partial P_{iJ}}{\partial F_{kL}} = J \mathcal{T}_{iskL} F_{sJ}^{-T} + P_{iJ} F_{kL}^{-T} - P_{iL} F_{kJ}^{-T}
99    /// ```
100    fn first_piola_kirchhoff_tangent_stiffness(
101        &self,
102        deformation_gradient: &DeformationGradient,
103    ) -> Result<FirstPiolaKirchhoffTangentStiffness, ConstitutiveError> {
104        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
105        let first_piola_kirchhoff_stress =
106            self.first_piola_kirchhoff_stress(deformation_gradient)?;
107        Ok(self
108            .cauchy_tangent_stiffness(deformation_gradient)?
109            .contract_second_index_with_first_index_of(&deformation_gradient_inverse_transpose)
110            * deformation_gradient.determinant()
111            + FirstPiolaKirchhoffTangentStiffness::dyad_ij_kl(
112                &first_piola_kirchhoff_stress,
113                &deformation_gradient_inverse_transpose,
114            )
115            - FirstPiolaKirchhoffTangentStiffness::dyad_il_kj(
116                &first_piola_kirchhoff_stress,
117                &deformation_gradient_inverse_transpose,
118            ))
119    }
120    /// Calculates and returns the second Piola-Kirchhoff stress.
121    ///
122    /// ```math
123    /// \mathbf{S} = \mathbf{F}^{-1}\cdot\mathbf{P}
124    /// ```
125    fn second_piola_kirchhoff_stress(
126        &self,
127        deformation_gradient: &DeformationGradient,
128    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
129        Ok(deformation_gradient.inverse()
130            * self.first_piola_kirchhoff_stress(deformation_gradient)?)
131    }
132    /// Calculates and returns the tangent stiffness associated with the second Piola-Kirchhoff stress.
133    ///
134    /// ```math
135    /// \mathcal{G}_{IJkL} = \frac{\partial S_{IJ}}{\partial F_{kL}} = \mathcal{C}_{mJkL}F_{mI}^{-T} - S_{LJ}F_{kI}^{-T} = J \mathcal{T}_{mnkL} F_{mI}^{-T} F_{nJ}^{-T} + S_{IJ} F_{kL}^{-T} - S_{IL} F_{kJ}^{-T} -S_{LJ} F_{kI}^{-T}
136    /// ```
137    fn second_piola_kirchhoff_tangent_stiffness(
138        &self,
139        deformation_gradient: &DeformationGradient,
140    ) -> Result<SecondPiolaKirchhoffTangentStiffness, ConstitutiveError> {
141        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
142        let deformation_gradient_inverse = deformation_gradient_inverse_transpose.transpose();
143        let second_piola_kirchhoff_stress =
144            self.second_piola_kirchhoff_stress(deformation_gradient)?;
145        Ok(self
146            .cauchy_tangent_stiffness(deformation_gradient)?
147            .contract_first_second_indices_with_second_indices_of(
148                &deformation_gradient_inverse,
149                &deformation_gradient_inverse,
150            )
151            * deformation_gradient.determinant()
152            + SecondPiolaKirchhoffTangentStiffness::dyad_ij_kl(
153                &second_piola_kirchhoff_stress,
154                &deformation_gradient_inverse_transpose,
155            )
156            - SecondPiolaKirchhoffTangentStiffness::dyad_il_kj(
157                &second_piola_kirchhoff_stress,
158                &deformation_gradient_inverse_transpose,
159            )
160            - SecondPiolaKirchhoffTangentStiffness::dyad_ik_jl(
161                &deformation_gradient_inverse,
162                &second_piola_kirchhoff_stress,
163            ))
164    }
165}
166
167/// Zeroth-order root-finding methods for elastic constitutive models.
168pub trait ZerothOrderRoot {
169    /// Solve for the unknown components of the deformation gradient under an applied load.
170    ///
171    /// ```math
172    /// \mathbf{P}(\mathbf{F}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
173    /// ```
174    fn root(
175        &self,
176        applied_load: AppliedLoad,
177        solver: impl ZerothOrderRootFinding<DeformationGradient>,
178    ) -> Result<DeformationGradient, ConstitutiveError>;
179}
180
181/// First-order root-finding methods for elastic constitutive models.
182pub trait FirstOrderRoot {
183    /// Solve for the unknown components of the deformation gradient under an applied load.
184    ///
185    /// ```math
186    /// \mathbf{P}(\mathbf{F}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
187    /// ```
188    fn root(
189        &self,
190        applied_load: AppliedLoad,
191        solver: impl FirstOrderRootFinding<
192            FirstPiolaKirchhoffStress,
193            FirstPiolaKirchhoffTangentStiffness,
194            DeformationGradient,
195        >,
196    ) -> Result<DeformationGradient, ConstitutiveError>;
197}
198
199impl<T> ZerothOrderRoot for T
200where
201    T: Elastic,
202{
203    fn root(
204        &self,
205        applied_load: AppliedLoad,
206        solver: impl ZerothOrderRootFinding<DeformationGradient>,
207    ) -> Result<DeformationGradient, ConstitutiveError> {
208        let (matrix, vector) = match applied_load {
209            AppliedLoad::UniaxialStress(deformation_gradient_11) => {
210                let mut matrix = Matrix::zero(4, 9);
211                let mut vector = Vector::zero(4);
212                matrix[0][0] = 1.0;
213                matrix[1][1] = 1.0;
214                matrix[2][2] = 1.0;
215                matrix[3][5] = 1.0;
216                vector[0] = deformation_gradient_11;
217                (matrix, vector)
218            }
219            AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => {
220                let mut matrix = Matrix::zero(5, 9);
221                let mut vector = Vector::zero(5);
222                matrix[0][0] = 1.0;
223                matrix[1][1] = 1.0;
224                matrix[2][2] = 1.0;
225                matrix[3][5] = 1.0;
226                matrix[4][4] = 1.0;
227                vector[0] = deformation_gradient_11;
228                vector[4] = deformation_gradient_22;
229                (matrix, vector)
230            }
231        };
232        match solver.root(
233            |deformation_gradient: &DeformationGradient| {
234                Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
235            },
236            DeformationGradient::identity(),
237            EqualityConstraint::Linear(matrix, vector),
238        ) {
239            Ok(deformation_gradient) => Ok(deformation_gradient),
240            Err(error) => Err(ConstitutiveError::Upstream(
241                format!("{error}"),
242                format!("{self:?}"),
243            )),
244        }
245    }
246}
247
248impl<T> FirstOrderRoot for T
249where
250    T: Elastic,
251{
252    fn root(
253        &self,
254        applied_load: AppliedLoad,
255        solver: impl FirstOrderRootFinding<
256            FirstPiolaKirchhoffStress,
257            FirstPiolaKirchhoffTangentStiffness,
258            DeformationGradient,
259        >,
260    ) -> Result<DeformationGradient, ConstitutiveError> {
261        let (matrix, vector) = match applied_load {
262            AppliedLoad::UniaxialStress(deformation_gradient_11) => {
263                let mut matrix = Matrix::zero(4, 9);
264                let mut vector = Vector::zero(4);
265                matrix[0][0] = 1.0;
266                matrix[1][1] = 1.0;
267                matrix[2][2] = 1.0;
268                matrix[3][5] = 1.0;
269                vector[0] = deformation_gradient_11;
270                (matrix, vector)
271            }
272            AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => {
273                let mut matrix = Matrix::zero(5, 9);
274                let mut vector = Vector::zero(5);
275                matrix[0][0] = 1.0;
276                matrix[1][1] = 1.0;
277                matrix[2][2] = 1.0;
278                matrix[3][5] = 1.0;
279                matrix[4][4] = 1.0;
280                vector[0] = deformation_gradient_11;
281                vector[4] = deformation_gradient_22;
282                (matrix, vector)
283            }
284        };
285        match solver.root(
286            |deformation_gradient: &DeformationGradient| {
287                Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
288            },
289            |deformation_gradient: &DeformationGradient| {
290                Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)?)
291            },
292            DeformationGradient::identity(),
293            EqualityConstraint::Linear(matrix, vector),
294        ) {
295            Ok(deformation_gradient) => Ok(deformation_gradient),
296            Err(error) => Err(ConstitutiveError::Upstream(
297                format!("{error}"),
298                format!("{self:?}"),
299            )),
300        }
301    }
302}
conspire/constitutive/solid/elastic/mod.rs

conspire/constitutive/solid/elastic/
mod.rs
