conspire/constitutive/solid/elastic_plastic/
mod.rs

1//! Elastic-plastic solid constitutive models.
2
3use crate::{
4    constitutive::{ConstitutiveError, fluid::plastic::Plastic, solid::Solid},
5    math::{
6        ContractFirstSecondWithSecond, ContractSecondWithFirst, IDENTITY, Matrix, Quantity, Rank2,
7        TensorArray, TensorRank2, TensorRank4,
8    },
9    mechanics::{
10        CauchyStress, CauchyTangentStiffness, CauchyTangentStiffnessPlastic, DeformationGradient,
11        DeformationGradientPlastic, FirstPiolaKirchhoffStress, FirstPiolaKirchhoffTangentStiffness,
12        FirstPiolaKirchhoffTangentStiffnessPlastic, MandelStressElastic,
13        MandelStressTangentElastic, MandelStressTangentElasticPlastic, Scalar,
14        SecondPiolaKirchhoffStress, SecondPiolaKirchhoffTangentStiffness,
15    },
16    units::Time,
17};
18use std::array::from_fn;
19
20/// Possible applied loads.
21pub enum AppliedLoad<'a> {
22    /// Uniaxial stress given $`F_{11}`$.
23    UniaxialStress(fn(Quantity<Time>) -> Scalar, &'a [Quantity<Time>]),
24    /// Biaxial stress given $`F_{11}`$ and $`F_{22}`$.
25    BiaxialStress(
26        fn(Quantity<Time>) -> Scalar,
27        fn(Quantity<Time>) -> Scalar,
28        &'a [Quantity<Time>],
29    ),
30}
31
32type Prescribed = Vec<(usize, fn(Quantity<Time>) -> Scalar)>;
33
34#[doc(hidden)]
35pub fn bcs(applied_load: AppliedLoad<'_>) -> (Matrix, Prescribed, &'_ [Quantity<Time>]) {
36    let (mut matrix, prescribed, time) = match applied_load {
37        AppliedLoad::UniaxialStress(deformation_gradient_11, time) => {
38            (Matrix::zero(4, 9), vec![(0, deformation_gradient_11)], time)
39        }
40        AppliedLoad::BiaxialStress(deformation_gradient_11, deformation_gradient_22, time) => (
41            Matrix::zero(5, 9),
42            vec![(0, deformation_gradient_11), (4, deformation_gradient_22)],
43            time,
44        ),
45    };
46    matrix[0][0] = 1.0;
47    matrix[1][1] = 1.0;
48    matrix[2][2] = 1.0;
49    matrix[3][5] = 1.0;
50    if matrix.len() == 5 {
51        matrix[4][4] = 1.0
52    }
53    (matrix, prescribed, time)
54}
55
56/// Required methods for elastic-plastic or elastic-viscoplastic solid constitutive models.
57pub trait ElasticPlasticOrViscoplastic
58where
59    Self: Solid + Plastic,
60{
61    /// Calculates and returns the Cauchy stress.
62    ///
63    /// ```math
64    /// \boldsymbol{\sigma} = \boldsymbol{\sigma}_\mathrm{e}
65    /// ```
66    fn cauchy_stress(
67        &self,
68        deformation_gradient: &DeformationGradient,
69        deformation_gradient_p: &DeformationGradientPlastic,
70    ) -> Result<CauchyStress, ConstitutiveError> {
71        Ok(deformation_gradient
72            * self.second_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_p)?
73            * deformation_gradient.transpose()
74            / deformation_gradient.determinant())
75    }
76    /// Calculates and returns the tangent stiffness associated with the Cauchy stress.
77    ///
78    /// ```math
79    /// \boldsymbol{\mathcal{T}} = \boldsymbol{\mathcal{T}}_\mathrm{e}\cdot\mathbf{F}_\mathrm{p}^{-T}
80    /// ```
81    fn cauchy_tangent_stiffness(
82        &self,
83        deformation_gradient: &DeformationGradient,
84        deformation_gradient_p: &DeformationGradientPlastic,
85    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
86        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
87        let cauchy_stress = self.cauchy_stress(deformation_gradient, deformation_gradient_p)?;
88        let some_stress = &cauchy_stress * &deformation_gradient_inverse_transpose;
89        Ok(self
90            .second_piola_kirchhoff_tangent_stiffness(deformation_gradient, deformation_gradient_p)?
91            .contract_first_second_with_second(deformation_gradient, deformation_gradient)
92            / deformation_gradient.determinant()
93            - CauchyTangentStiffness::dyad_ij_kl(
94                &cauchy_stress,
95                &deformation_gradient_inverse_transpose,
96            )
97            + CauchyTangentStiffness::dyad_il_kj(&some_stress, &IDENTITY)
98            + CauchyTangentStiffness::dyad_ik_jl(&IDENTITY, &some_stress))
99    }
100    /// Calculates and returns the first Piola-Kirchhoff stress.
101    ///
102    /// ```math
103    /// \mathbf{P} = \mathbf{P}_\mathrm{e}\cdot\mathbf{F}_\mathrm{p}^{-T}
104    /// ```
105    fn first_piola_kirchhoff_stress(
106        &self,
107        deformation_gradient: &DeformationGradient,
108        deformation_gradient_p: &DeformationGradientPlastic,
109    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
110        Ok(
111            self.cauchy_stress(deformation_gradient, deformation_gradient_p)?
112                * deformation_gradient.inverse_transpose()
113                * deformation_gradient.determinant(),
114        )
115    }
116    /// Calculates and returns the tangent stiffness associated with the first Piola-Kirchhoff stress.
117    ///
118    /// ```math
119    /// \mathcal{C}_{iJkL} = \mathcal{C}^\mathrm{e}_{iMkN} F_{MJ}^{\mathrm{p}-T} F_{NL}^{\mathrm{p}-T}
120    /// ```
121    fn first_piola_kirchhoff_tangent_stiffness(
122        &self,
123        deformation_gradient: &DeformationGradient,
124        deformation_gradient_p: &DeformationGradientPlastic,
125    ) -> Result<FirstPiolaKirchhoffTangentStiffness, ConstitutiveError> {
126        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
127        let first_piola_kirchhoff_stress =
128            self.first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_p)?;
129        Ok(self
130            .cauchy_tangent_stiffness(deformation_gradient, deformation_gradient_p)?
131            .contract_second_with_first(&deformation_gradient_inverse_transpose)
132            * deformation_gradient.determinant()
133            + FirstPiolaKirchhoffTangentStiffness::dyad_ij_kl(
134                &first_piola_kirchhoff_stress,
135                &deformation_gradient_inverse_transpose,
136            )
137            - FirstPiolaKirchhoffTangentStiffness::dyad_il_kj(
138                &first_piola_kirchhoff_stress,
139                &deformation_gradient_inverse_transpose,
140            ))
141    }
142    /// Calculates and returns the Mandel stress.
143    ///
144    /// ```math
145    /// \mathbf{M}_\mathrm{e} = J\mathbf{F}_\mathrm{e}^T\cdot\boldsymbol{\sigma}\cdot\mathbf{F}_\mathrm{e}^{-T}
146    /// ```
147    fn mandel_stress(
148        &self,
149        deformation_gradient: &DeformationGradient,
150        deformation_gradient_p: &DeformationGradientPlastic,
151    ) -> Result<MandelStressElastic, ConstitutiveError> {
152        let jacobian = self.jacobian(deformation_gradient)?;
153        let deformation_gradient_e = deformation_gradient * deformation_gradient_p.inverse();
154        let cauchy_stress = self.cauchy_stress(deformation_gradient, deformation_gradient_p)?;
155        Ok((deformation_gradient_e.transpose()
156            * cauchy_stress
157            * deformation_gradient_e.inverse_transpose())
158            * jacobian)
159    }
160    /// Calculates and returns the second Piola-Kirchhoff stress.
161    ///
162    /// ```math
163    /// \mathbf{S} = \mathbf{F}_\mathrm{p}^{-1}\cdot\mathbf{S}_\mathrm{e}\cdot\mathbf{F}_\mathrm{p}^{-T}
164    /// ```
165    fn second_piola_kirchhoff_stress(
166        &self,
167        deformation_gradient: &DeformationGradient,
168        deformation_gradient_p: &DeformationGradientPlastic,
169    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
170        Ok(deformation_gradient.inverse()
171            * self.first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_p)?)
172    }
173    /// Calculates and returns the tangent stiffness associated with the second Piola-Kirchhoff stress.
174    ///
175    /// ```math
176    /// \mathcal{G}_{IJkL} = \mathcal{G}^\mathrm{e}_{MNkO} F_{MI}^{\mathrm{p}-T} F_{NJ}^{\mathrm{p}-T} F_{OL}^{\mathrm{p}-T}
177    /// ```
178    fn second_piola_kirchhoff_tangent_stiffness(
179        &self,
180        deformation_gradient: &DeformationGradient,
181        deformation_gradient_p: &DeformationGradientPlastic,
182    ) -> Result<SecondPiolaKirchhoffTangentStiffness, ConstitutiveError> {
183        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
184        let deformation_gradient_inverse = deformation_gradient_inverse_transpose.transpose();
185        let second_piola_kirchhoff_stress =
186            self.second_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_p)?;
187        Ok(self
188            .cauchy_tangent_stiffness(deformation_gradient, deformation_gradient_p)?
189            .contract_first_second_with_second(
190                &deformation_gradient_inverse,
191                &deformation_gradient_inverse,
192            )
193            * deformation_gradient.determinant()
194            + SecondPiolaKirchhoffTangentStiffness::dyad_ij_kl(
195                &second_piola_kirchhoff_stress,
196                &deformation_gradient_inverse_transpose,
197            )
198            - SecondPiolaKirchhoffTangentStiffness::dyad_il_kj(
199                &second_piola_kirchhoff_stress,
200                &deformation_gradient_inverse_transpose,
201            )
202            - SecondPiolaKirchhoffTangentStiffness::dyad_ik_jl(
203                &deformation_gradient_inverse,
204                &second_piola_kirchhoff_stress,
205            ))
206    }
207}
208
209pub(crate) type Matrix3 = [[Scalar; 3]; 3];
210pub(crate) type Entries4 = [[[[Scalar; 3]; 3]; 3]; 3];
211
212pub(crate) fn matrix_3<I, J, U>(tensor: &TensorRank2<3, I, J, U>) -> Matrix3 {
213    from_fn(|i| from_fn(|j| tensor[i][j].value()))
214}
215
216pub(crate) fn entries_4<I, J, K, L, U>(tensor: &TensorRank4<3, I, J, K, L, U>) -> Entries4 {
217    from_fn(|i| from_fn(|j| from_fn(|k| from_fn(|l| tensor[i][j][k][l].value()))))
218}
219
220pub(crate) fn rank_4<I, J, K, L, U>(entries: &Entries4) -> TensorRank4<3, I, J, K, L, U> {
221    let mut tensor = TensorRank4::zero();
222    (0..3).for_each(|i| {
223        (0..3).for_each(|j| {
224            (0..3).for_each(|k| {
225                (0..3).for_each(|l| tensor[i][j][k][l] = Quantity::new(entries[i][j][k][l]))
226            })
227        })
228    });
229    tensor
230}
231
232//
233// M_ij = J F^e_ki sigma_kl F^{e-1}_jl, differentiated once for a direction (a, b) in
234// which dF^e_ki = g_ka F^{p-1}_bi: g = 1 for F itself and g = -F^e for F^p, with the
235// Jacobian term present only for F.
236//
237#[allow(clippy::too_many_arguments)]
238fn mandel_stress_tangent_entries(
239    jacobian: Scalar,
240    mandel_stress: &Matrix3,
241    cauchy_stress: &Matrix3,
242    deformation_gradient_e: &Matrix3,
243    deformation_gradient_e_inverse: &Matrix3,
244    deformation_gradient_p_inverse: &Matrix3,
245    g: &Matrix3,
246    d_jacobian: Option<&Matrix3>,
247    d_cauchy_stress: &Entries4,
248) -> Entries4 {
249    let w: Matrix3 = from_fn(|a| {
250        from_fn(|j| {
251            (0..3)
252                .map(|k| {
253                    (0..3)
254                        .map(|l| {
255                            g[k][a] * cauchy_stress[k][l] * deformation_gradient_e_inverse[j][l]
256                        })
257                        .sum::<Scalar>()
258                })
259                .sum()
260        })
261    });
262    let v: Matrix3 = from_fn(|j| {
263        from_fn(|a| {
264            (0..3)
265                .map(|m| deformation_gradient_e_inverse[j][m] * g[m][a])
266                .sum()
267        })
268    });
269    let u: Matrix3 = from_fn(|i| {
270        from_fn(|b| {
271            (0..3)
272                .map(|n| mandel_stress[i][n] * deformation_gradient_p_inverse[b][n])
273                .sum()
274        })
275    });
276    from_fn(|i| {
277        from_fn(|j| {
278            from_fn(|a| {
279                from_fn(|b| {
280                    let elastic = (0..3)
281                        .map(|k| {
282                            (0..3)
283                                .map(|l| {
284                                    deformation_gradient_e[k][i]
285                                        * d_cauchy_stress[k][l][a][b]
286                                        * deformation_gradient_e_inverse[j][l]
287                                })
288                                .sum::<Scalar>()
289                        })
290                        .sum::<Scalar>();
291                    d_jacobian.map_or(0.0, |d_jacobian| d_jacobian[a][b] * mandel_stress[i][j])
292                        + jacobian * (deformation_gradient_p_inverse[b][i] * w[a][j] + elastic)
293                        - u[i][b] * v[j][a]
294                })
295            })
296        })
297    })
298}
299
300/// Tangents with respect to the plastic deformation gradient, and of the Mandel
301/// stress with respect to both deformation gradients — the pieces a coupled
302/// (condensed) return map needs beyond the elastic tangents.
303pub trait PlasticTangents
304where
305    Self: ElasticPlasticOrViscoplastic,
306{
307    /// Calculates and returns the tangent stiffness of the Cauchy stress with
308    /// respect to the plastic deformation gradient.
309    ///
310    /// ```math
311    /// \frac{\partial\sigma_{ij}}{\partial F^\mathrm{p}_{NO}} = -\mathcal{T}^\mathrm{e}_{ijmA} F^\mathrm{e}_{mN} F^{\mathrm{p}-1}_{OA}
312    /// ```
313    fn cauchy_tangent_stiffness_p(
314        &self,
315        deformation_gradient: &DeformationGradient,
316        deformation_gradient_p: &DeformationGradientPlastic,
317    ) -> Result<CauchyTangentStiffnessPlastic, ConstitutiveError>;
318    /// Calculates and returns the tangent stiffness of the first Piola-Kirchhoff
319    /// stress with respect to the plastic deformation gradient.
320    ///
321    /// ```math
322    /// \frac{\partial P_{iJ}}{\partial F^\mathrm{p}_{NO}} = -\mathcal{C}^\mathrm{e}_{iAmB} F^\mathrm{e}_{mN} F^{\mathrm{p}-1}_{OB} F^{\mathrm{p}-1}_{JA} - P_{iO} F^{\mathrm{p}-1}_{JN}
323    /// ```
324    fn first_piola_kirchhoff_tangent_stiffness_p(
325        &self,
326        deformation_gradient: &DeformationGradient,
327        deformation_gradient_p: &DeformationGradientPlastic,
328    ) -> Result<FirstPiolaKirchhoffTangentStiffnessPlastic, ConstitutiveError>;
329    /// Calculates and returns the tangent stiffness of the Mandel stress with
330    /// respect to the deformation gradient.
331    fn mandel_stress_tangent(
332        &self,
333        deformation_gradient: &DeformationGradient,
334        deformation_gradient_p: &DeformationGradientPlastic,
335    ) -> Result<MandelStressTangentElastic, ConstitutiveError> {
336        let deformation_gradient_p_inverse = deformation_gradient_p.inverse();
337        let deformation_gradient_e = deformation_gradient * &deformation_gradient_p_inverse;
338        Ok(rank_4(&mandel_stress_tangent_entries(
339            self.jacobian(deformation_gradient)?,
340            &matrix_3(&self.mandel_stress(deformation_gradient, deformation_gradient_p)?),
341            &matrix_3(&self.cauchy_stress(deformation_gradient, deformation_gradient_p)?),
342            &matrix_3(&deformation_gradient_e),
343            &matrix_3(&deformation_gradient_e.inverse()),
344            &matrix_3(&deformation_gradient_p_inverse),
345            &matrix_3(&IDENTITY),
346            Some(&matrix_3(&deformation_gradient.inverse_transpose())),
347            &entries_4(
348                &self.cauchy_tangent_stiffness(deformation_gradient, deformation_gradient_p)?,
349            ),
350        )))
351    }
352    /// Calculates and returns the tangent stiffness of the Mandel stress with
353    /// respect to the plastic deformation gradient.
354    fn mandel_stress_tangent_p(
355        &self,
356        deformation_gradient: &DeformationGradient,
357        deformation_gradient_p: &DeformationGradientPlastic,
358    ) -> Result<MandelStressTangentElasticPlastic, ConstitutiveError> {
359        let deformation_gradient_p_inverse = deformation_gradient_p.inverse();
360        let deformation_gradient_e = deformation_gradient * &deformation_gradient_p_inverse;
361        Ok(rank_4(&mandel_stress_tangent_entries(
362            self.jacobian(deformation_gradient)?,
363            &matrix_3(&self.mandel_stress(deformation_gradient, deformation_gradient_p)?),
364            &matrix_3(&self.cauchy_stress(deformation_gradient, deformation_gradient_p)?),
365            &matrix_3(&deformation_gradient_e),
366            &matrix_3(&deformation_gradient_e.inverse()),
367            &matrix_3(&deformation_gradient_p_inverse),
368            &matrix_3(&(&deformation_gradient_e * -1.0)),
369            None,
370            &entries_4(
371                &self.cauchy_tangent_stiffness_p(deformation_gradient, deformation_gradient_p)?,
372            ),
373        )))
374    }
375}
376
377/// Required methods for elastic-plastic solid constitutive models.
378pub trait ElasticPlastic
379where
380    Self: ElasticPlasticOrViscoplastic,
381{
382}
conspire/constitutive/solid/elastic_plastic/mod.rs

conspire/constitutive/solid/elastic_plastic/
mod.rs
