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

conspire/constitutive/solid/elastic/
mod.rs
