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::{EqualityConstraint, FirstOrderRootFinding, NewtonRaphson, 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    /// Solve for the unknown components of the deformation gradient under an applied load.
159    ///
160    /// ```math
161    /// \mathbf{P}(\mathbf{F}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
162    /// ```
163    fn root(&self, applied_load: AppliedLoad) -> Result<DeformationGradient, OptimizeError> {
164        let solver = NewtonRaphson {
165            ..Default::default()
166        };
167        match applied_load {
168            AppliedLoad::UniaxialStress(deformation_gradient_11) => {
169                let mut matrix = Matrix::zero(4, 9);
170                let mut vector = Vector::zero(4);
171                matrix[0][0] = 1.0;
172                matrix[1][1] = 1.0;
173                matrix[2][2] = 1.0;
174                matrix[3][5] = 1.0;
175                vector[0] = deformation_gradient_11;
176                solver.root(
177                    |deformation_gradient: &DeformationGradient| {
178                        Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
179                    },
180                    |deformation_gradient: &DeformationGradient| {
181                        Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)?)
182                    },
183                    DeformationGradient::identity(),
184                    EqualityConstraint::Linear(matrix, vector),
185                )
186            }
187            AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => {
188                let mut matrix = Matrix::zero(5, 9);
189                let mut vector = Vector::zero(5);
190                matrix[0][0] = 1.0;
191                matrix[1][1] = 1.0;
192                matrix[2][2] = 1.0;
193                matrix[3][5] = 1.0;
194                matrix[4][4] = 1.0;
195                vector[0] = deformation_gradient_11;
196                vector[4] = deformation_gradient_22;
197                solver.root(
198                    |deformation_gradient: &DeformationGradient| {
199                        Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
200                    },
201                    |deformation_gradient: &DeformationGradient| {
202                        Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)?)
203                    },
204                    DeformationGradient::identity(),
205                    EqualityConstraint::Linear(matrix, vector),
206                )
207            }
208        }
209    }
210}
conspire/constitutive/solid/elastic/mod.rs

conspire/constitutive/solid/elastic/
mod.rs
