conspire/constitutive/solid/elastic/
mod.rs

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

conspire/constitutive/solid/elastic/
mod.rs
