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_with_second(deformation_gradient, deformation_gradient)
71            / deformation_gradient.determinant()
72            - CauchyTangentStiffness::dyad_ij_kl(
73                &cauchy_stress,
74                &deformation_gradient_inverse_transpose,
75            )
76            + CauchyTangentStiffness::dyad_il_kj(&some_stress, &IDENTITY)
77            + CauchyTangentStiffness::dyad_ik_jl(&IDENTITY, &some_stress))
78    }
79    /// Calculates and returns the first Piola-Kirchhoff stress.
80    ///
81    /// ```math
82    /// \mathbf{P} = J\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}
83    /// ```
84    fn first_piola_kirchhoff_stress(
85        &self,
86        deformation_gradient: &DeformationGradient,
87    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
88        Ok(self.cauchy_stress(deformation_gradient)?
89            * deformation_gradient.inverse_transpose()
90            * deformation_gradient.determinant())
91    }
92    /// Calculates and returns the tangent stiffness associated with the first Piola-Kirchhoff stress.
93    ///
94    /// ```math
95    /// \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}
96    /// ```
97    fn first_piola_kirchhoff_tangent_stiffness(
98        &self,
99        deformation_gradient: &DeformationGradient,
100    ) -> Result<FirstPiolaKirchhoffTangentStiffness, ConstitutiveError> {
101        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
102        let first_piola_kirchhoff_stress =
103            self.first_piola_kirchhoff_stress(deformation_gradient)?;
104        Ok(self
105            .cauchy_tangent_stiffness(deformation_gradient)?
106            .contract_second_with_first(&deformation_gradient_inverse_transpose)
107            * deformation_gradient.determinant()
108            + FirstPiolaKirchhoffTangentStiffness::dyad_ij_kl(
109                &first_piola_kirchhoff_stress,
110                &deformation_gradient_inverse_transpose,
111            )
112            - FirstPiolaKirchhoffTangentStiffness::dyad_il_kj(
113                &first_piola_kirchhoff_stress,
114                &deformation_gradient_inverse_transpose,
115            ))
116    }
117    /// Calculates and returns the second Piola-Kirchhoff stress.
118    ///
119    /// ```math
120    /// \mathbf{S} = \mathbf{F}^{-1}\cdot\mathbf{P}
121    /// ```
122    fn second_piola_kirchhoff_stress(
123        &self,
124        deformation_gradient: &DeformationGradient,
125    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
126        Ok(deformation_gradient.inverse()
127            * self.first_piola_kirchhoff_stress(deformation_gradient)?)
128    }
129    /// Calculates and returns the tangent stiffness associated with the second Piola-Kirchhoff stress.
130    ///
131    /// ```math
132    /// \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}
133    /// ```
134    fn second_piola_kirchhoff_tangent_stiffness(
135        &self,
136        deformation_gradient: &DeformationGradient,
137    ) -> Result<SecondPiolaKirchhoffTangentStiffness, ConstitutiveError> {
138        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
139        let deformation_gradient_inverse = deformation_gradient_inverse_transpose.transpose();
140        let second_piola_kirchhoff_stress =
141            self.second_piola_kirchhoff_stress(deformation_gradient)?;
142        Ok(self
143            .cauchy_tangent_stiffness(deformation_gradient)?
144            .contract_first_second_with_second(
145                &deformation_gradient_inverse,
146                &deformation_gradient_inverse,
147            )
148            * deformation_gradient.determinant()
149            + SecondPiolaKirchhoffTangentStiffness::dyad_ij_kl(
150                &second_piola_kirchhoff_stress,
151                &deformation_gradient_inverse_transpose,
152            )
153            - SecondPiolaKirchhoffTangentStiffness::dyad_il_kj(
154                &second_piola_kirchhoff_stress,
155                &deformation_gradient_inverse_transpose,
156            )
157            - SecondPiolaKirchhoffTangentStiffness::dyad_ik_jl(
158                &deformation_gradient_inverse,
159                &second_piola_kirchhoff_stress,
160            ))
161    }
162}
163
164/// Zeroth-order root-finding methods for elastic solid constitutive models.
165pub trait ZerothOrderRoot {
166    /// Solve for the unknown components of the deformation gradient under an applied load.
167    ///
168    /// ```math
169    /// \mathbf{P}(\mathbf{F}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
170    /// ```
171    fn root(
172        &self,
173        applied_load: AppliedLoad,
174        solver: impl ZerothOrderRootFinding<DeformationGradient>,
175    ) -> Result<DeformationGradient, ConstitutiveError>;
176}
177
178/// First-order root-finding methods for elastic solid constitutive models.
179pub trait FirstOrderRoot {
180    /// Solve for the unknown components of the deformation gradient under an applied load.
181    ///
182    /// ```math
183    /// \mathbf{P}(\mathbf{F}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
184    /// ```
185    fn root(
186        &self,
187        applied_load: AppliedLoad,
188        solver: impl FirstOrderRootFinding<
189            FirstPiolaKirchhoffStress,
190            FirstPiolaKirchhoffTangentStiffness,
191            DeformationGradient,
192        >,
193    ) -> Result<DeformationGradient, ConstitutiveError>;
194}
195
196impl<T> ZerothOrderRoot for T
197where
198    T: Elastic,
199{
200    fn root(
201        &self,
202        applied_load: AppliedLoad,
203        solver: impl ZerothOrderRootFinding<DeformationGradient>,
204    ) -> Result<DeformationGradient, ConstitutiveError> {
205        let (matrix, vector) = bcs(applied_load);
206        match solver.root(
207            |deformation_gradient: &DeformationGradient| {
208                Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
209            },
210            DeformationGradient::identity(),
211            EqualityConstraint::Linear(matrix, vector),
212        ) {
213            Ok(deformation_gradient) => Ok(deformation_gradient),
214            Err(error) => Err(ConstitutiveError::Upstream(
215                format!("{error}"),
216                format!("{self:?}"),
217            )),
218        }
219    }
220}
221
222impl<T> FirstOrderRoot for T
223where
224    T: Elastic,
225{
226    fn root(
227        &self,
228        applied_load: AppliedLoad,
229        solver: impl FirstOrderRootFinding<
230            FirstPiolaKirchhoffStress,
231            FirstPiolaKirchhoffTangentStiffness,
232            DeformationGradient,
233        >,
234    ) -> Result<DeformationGradient, ConstitutiveError> {
235        let (matrix, vector) = bcs(applied_load);
236        match solver.root(
237            |deformation_gradient: &DeformationGradient| {
238                Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
239            },
240            |deformation_gradient: &DeformationGradient| {
241                Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)?)
242            },
243            DeformationGradient::identity(),
244            EqualityConstraint::Linear(matrix, vector),
245            None,
246        ) {
247            Ok(deformation_gradient) => Ok(deformation_gradient),
248            Err(error) => Err(ConstitutiveError::Upstream(
249                format!("{error}"),
250                format!("{self:?}"),
251            )),
252        }
253    }
254}
255
256#[doc(hidden)]
257pub fn bcs(applied_load: AppliedLoad) -> (Matrix, Vector) {
258    match applied_load {
259        AppliedLoad::UniaxialStress(deformation_gradient_11) => {
260            let mut matrix = Matrix::zero(4, 9);
261            let mut vector = Vector::zero(4);
262            matrix[0][0] = 1.0;
263            matrix[1][1] = 1.0;
264            matrix[2][2] = 1.0;
265            matrix[3][5] = 1.0;
266            vector[0] = deformation_gradient_11;
267            (matrix, vector)
268        }
269        AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => {
270            let mut matrix = Matrix::zero(5, 9);
271            let mut vector = Vector::zero(5);
272            matrix[0][0] = 1.0;
273            matrix[1][1] = 1.0;
274            matrix[2][2] = 1.0;
275            matrix[3][5] = 1.0;
276            matrix[4][4] = 1.0;
277            vector[0] = deformation_gradient_11;
278            vector[4] = deformation_gradient_22;
279            (matrix, vector)
280        }
281    }
282}
conspire/constitutive/solid/elastic/mod.rs

conspire/constitutive/solid/elastic/
mod.rs
