conspire/constitutive/solid/elastic/
mod.rs

1//! Elastic solid 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 solid 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 solid 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 solid 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) = bcs(applied_load);
209        match solver.root(
210            |deformation_gradient: &DeformationGradient| {
211                Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
212            },
213            DeformationGradient::identity(),
214            EqualityConstraint::Linear(matrix, vector),
215        ) {
216            Ok(deformation_gradient) => Ok(deformation_gradient),
217            Err(error) => Err(ConstitutiveError::Upstream(
218                format!("{error}"),
219                format!("{self:?}"),
220            )),
221        }
222    }
223}
224
225impl<T> FirstOrderRoot for T
226where
227    T: Elastic,
228{
229    fn root(
230        &self,
231        applied_load: AppliedLoad,
232        solver: impl FirstOrderRootFinding<
233            FirstPiolaKirchhoffStress,
234            FirstPiolaKirchhoffTangentStiffness,
235            DeformationGradient,
236        >,
237    ) -> Result<DeformationGradient, ConstitutiveError> {
238        let (matrix, vector) = bcs(applied_load);
239        match solver.root(
240            |deformation_gradient: &DeformationGradient| {
241                Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
242            },
243            |deformation_gradient: &DeformationGradient| {
244                Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)?)
245            },
246            DeformationGradient::identity(),
247            EqualityConstraint::Linear(matrix, vector),
248        ) {
249            Ok(deformation_gradient) => Ok(deformation_gradient),
250            Err(error) => Err(ConstitutiveError::Upstream(
251                format!("{error}"),
252                format!("{self:?}"),
253            )),
254        }
255    }
256}
257
258#[doc(hidden)]
259pub fn bcs(applied_load: AppliedLoad) -> (Matrix, Vector) {
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            (matrix, vector)
270        }
271        AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => {
272            let mut matrix = Matrix::zero(5, 9);
273            let mut vector = Vector::zero(5);
274            matrix[0][0] = 1.0;
275            matrix[1][1] = 1.0;
276            matrix[2][2] = 1.0;
277            matrix[3][5] = 1.0;
278            matrix[4][4] = 1.0;
279            vector[0] = deformation_gradient_11;
280            vector[4] = deformation_gradient_22;
281            (matrix, vector)
282        }
283    }
284}
conspire/constitutive/solid/elastic/mod.rs

conspire/constitutive/solid/elastic/
mod.rs
