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, RootFinding},
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/// Root-finding methods for elastic solid constitutive models.
176pub trait Root {
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 RootFinding<
186            FirstPiolaKirchhoffStress,
187            FirstPiolaKirchhoffTangentStiffness,
188            DeformationGradient,
189        >,
190    ) -> Result<DeformationGradient, ConstitutiveError>;
191}
192
193impl<T> Root for T
194where
195    T: Elastic,
196{
197    fn root(
198        &self,
199        applied_load: AppliedLoad,
200        solver: impl RootFinding<
201            FirstPiolaKirchhoffStress,
202            FirstPiolaKirchhoffTangentStiffness,
203            DeformationGradient,
204        >,
205    ) -> Result<DeformationGradient, ConstitutiveError> {
206        let (matrix, vector) = bcs(applied_load);
207        solver
208            .root(
209                |deformation_gradient: &DeformationGradient| {
210                    Ok(self.first_piola_kirchhoff_stress(deformation_gradient)?)
211                },
212                |deformation_gradient: &DeformationGradient| {
213                    Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient)?)
214                },
215                DeformationGradient::identity(),
216                vec![EqualityConstraint::Linear(matrix, vector)],
217                None,
218            )
219            .map_err(|error| ConstitutiveError::upstream(error, self))
220    }
221}
222
223#[doc(hidden)]
224pub fn bcs(applied_load: AppliedLoad) -> (Matrix, Vector) {
225    match applied_load {
226        AppliedLoad::UniaxialStress(deformation_gradient_11) => {
227            let mut matrix = Matrix::zero(4, 9);
228            let mut vector = Vector::zero(4);
229            matrix[0][0] = 1.0;
230            matrix[1][1] = 1.0;
231            matrix[2][2] = 1.0;
232            matrix[3][5] = 1.0;
233            vector[0] = deformation_gradient_11;
234            (matrix, vector)
235        }
236        AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22) => {
237            let mut matrix = Matrix::zero(5, 9);
238            let mut vector = Vector::zero(5);
239            matrix[0][0] = 1.0;
240            matrix[1][1] = 1.0;
241            matrix[2][2] = 1.0;
242            matrix[3][5] = 1.0;
243            matrix[4][4] = 1.0;
244            vector[0] = deformation_gradient_11;
245            vector[4] = deformation_gradient_22;
246            (matrix, vector)
247        }
248    }
249}
conspire/constitutive/solid/elastic/mod.rs

conspire/constitutive/solid/elastic/
mod.rs
