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
13mod almansi_hamel;
14mod hencky;
15mod saint_venant_kirchhoff;
16
17pub use self::{
18    almansi_hamel::AlmansiHamel, hencky::Hencky, saint_venant_kirchhoff::SaintVenantKirchhoff,
19};
20
21use super::*;
22use crate::math::{
23    Matrix, TensorVec, Vector,
24    optimize::{self, EqualityConstraint, OptimizeError},
25};
26
27/// Possible applied loads.
28pub enum AppliedLoad {
29    /// Uniaxial stress given $`F_{11}`$.
30    UniaxialStress(Scalar),
31    /// Biaxial stress given $`F_{11}`$ and $`F_{22}`$.
32    BiaxialStress(Scalar, Scalar),
33}
34
35/// Required methods for elastic constitutive models.
36pub trait Elastic
37where
38    Self: Solid,
39{
40    /// Calculates and returns the Cauchy stress.
41    ///
42    /// ```math
43    /// \boldsymbol{\sigma} = J^{-1}\mathbf{P}\cdot\mathbf{F}^T
44    /// ```
45    fn cauchy_stress(
46        &self,
47        deformation_gradient: &DeformationGradient,
48    ) -> Result<CauchyStress, ConstitutiveError> {
49        Ok(deformation_gradient
50            * self.second_piola_kirchhoff_stress(deformation_gradient)?
51            * deformation_gradient.transpose()
52            / deformation_gradient.determinant())
53    }
54    /// Calculates and returns the tangent stiffness associated with the Cauchy stress.
55    ///
56    /// ```math
57    /// \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}
58    /// ```
59    fn cauchy_tangent_stiffness(
60        &self,
61        deformation_gradient: &DeformationGradient,
62    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
63        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
64        let cauchy_stress = self.cauchy_stress(deformation_gradient)?;
65        let some_stress = &cauchy_stress * &deformation_gradient_inverse_transpose;
66        Ok(self
67            .second_piola_kirchhoff_tangent_stiffness(deformation_gradient)?
68            .contract_first_second_indices_with_second_indices_of(
69                deformation_gradient,
70                deformation_gradient,
71            )
72            / deformation_gradient.determinant()
73            - CauchyTangentStiffness::dyad_ij_kl(
74                &cauchy_stress,
75                &deformation_gradient_inverse_transpose,
76            )
77            + CauchyTangentStiffness::dyad_il_kj(&some_stress, &IDENTITY)
78            + CauchyTangentStiffness::dyad_ik_jl(&IDENTITY, &some_stress))
79    }
80    /// Calculates and returns the first Piola-Kirchhoff stress.
81    ///
82    /// ```math
83    /// \mathbf{P} = J\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}
84    /// ```
85    fn first_piola_kirchhoff_stress(
86        &self,
87        deformation_gradient: &DeformationGradient,
88    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
89        Ok(self.cauchy_stress(deformation_gradient)?
90            * deformation_gradient.inverse_transpose()
91            * deformation_gradient.determinant())
92    }
93    /// Calculates and returns the tangent stiffness associated with the first Piola-Kirchhoff stress.
94    ///
95    /// ```math
96    /// \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}
97    /// ```
98    fn first_piola_kirchhoff_tangent_stiffness(
99        &self,
100        deformation_gradient: &DeformationGradient,
101    ) -> Result<FirstPiolaKirchhoffTangentStiffness, ConstitutiveError> {
102        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
103        let first_piola_kirchhoff_stress =
104            self.first_piola_kirchhoff_stress(deformation_gradient)?;
105        Ok(self
106            .cauchy_tangent_stiffness(deformation_gradient)?
107            .contract_second_index_with_first_index_of(&deformation_gradient_inverse_transpose)
108            * deformation_gradient.determinant()
109            + FirstPiolaKirchhoffTangentStiffness::dyad_ij_kl(
110                &first_piola_kirchhoff_stress,
111                &deformation_gradient_inverse_transpose,
112            )
113            - FirstPiolaKirchhoffTangentStiffness::dyad_il_kj(
114                &first_piola_kirchhoff_stress,
115                &deformation_gradient_inverse_transpose,
116            ))
117    }
118    /// Calculates and returns the second Piola-Kirchhoff stress.
119    ///
120    /// ```math
121    /// \mathbf{S} = \mathbf{F}^{-1}\cdot\mathbf{P}
122    /// ```
123    fn second_piola_kirchhoff_stress(
124        &self,
125        deformation_gradient: &DeformationGradient,
126    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
127        Ok(deformation_gradient.inverse()
128            * self.cauchy_stress(deformation_gradient)?
129            * deformation_gradient.inverse_transpose()
130            * deformation_gradient.determinant())
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 optimize::ZerothOrderRootFinding<DeformationGradient>,
178    ) -> Result<DeformationGradient, OptimizeError>;
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 optimize::FirstOrderRootFinding<
192            FirstPiolaKirchhoffStress,
193            FirstPiolaKirchhoffTangentStiffness,
194            DeformationGradient,
195        >,
196    ) -> Result<DeformationGradient, OptimizeError>;
197}
198
199impl<T> ZerothOrderRoot for T
200where
201    T: Elastic,
202{
203    fn root(
204        &self,
205        applied_load: AppliedLoad,
206        solver: impl optimize::ZerothOrderRootFinding<DeformationGradient>,
207    ) -> Result<DeformationGradient, OptimizeError> {
208        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                solver.root(
218                    |deformation_gradient: &DeformationGradient| {
219                        Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
220                    },
221                    DeformationGradient::identity(),
222                    EqualityConstraint::Linear(matrix, vector),
223                )
224            }
225            AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => {
226                let mut matrix = Matrix::zero(5, 9);
227                let mut vector = Vector::zero(5);
228                matrix[0][0] = 1.0;
229                matrix[1][1] = 1.0;
230                matrix[2][2] = 1.0;
231                matrix[3][5] = 1.0;
232                matrix[4][4] = 1.0;
233                vector[0] = deformation_gradient_11;
234                vector[4] = deformation_gradient_22;
235                solver.root(
236                    |deformation_gradient: &DeformationGradient| {
237                        Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
238                    },
239                    DeformationGradient::identity(),
240                    EqualityConstraint::Linear(matrix, vector),
241                )
242            }
243        }
244    }
245}
246
247impl<T> FirstOrderRoot for T
248where
249    T: Elastic,
250{
251    fn root(
252        &self,
253        applied_load: AppliedLoad,
254        solver: impl optimize::FirstOrderRootFinding<
255            FirstPiolaKirchhoffStress,
256            FirstPiolaKirchhoffTangentStiffness,
257            DeformationGradient,
258        >,
259    ) -> Result<DeformationGradient, OptimizeError> {
260        match applied_load {
261            AppliedLoad::UniaxialStress(deformation_gradient_11) => {
262                let mut matrix = Matrix::zero(4, 9);
263                let mut vector = Vector::zero(4);
264                matrix[0][0] = 1.0;
265                matrix[1][1] = 1.0;
266                matrix[2][2] = 1.0;
267                matrix[3][5] = 1.0;
268                vector[0] = deformation_gradient_11;
269                solver.root(
270                    |deformation_gradient: &DeformationGradient| {
271                        Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
272                    },
273                    |deformation_gradient: &DeformationGradient| {
274                        Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)?)
275                    },
276                    DeformationGradient::identity(),
277                    EqualityConstraint::Linear(matrix, vector),
278                )
279            }
280            AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => {
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                vector[0] = deformation_gradient_11;
289                vector[4] = deformation_gradient_22;
290                solver.root(
291                    |deformation_gradient: &DeformationGradient| {
292                        Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
293                    },
294                    |deformation_gradient: &DeformationGradient| {
295                        Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)?)
296                    },
297                    DeformationGradient::identity(),
298                    EqualityConstraint::Linear(matrix, vector),
299                )
300            }
301        }
302    }
303}
conspire/constitutive/solid/elastic/mod.rs

conspire/constitutive/solid/elastic/
mod.rs
