conspire/constitutive/solid/elastic_viscoplastic/
mod.rs

1//! Elastic-viscoplastic solid constitutive models.
2
3use crate::{
4    constitutive::{ConstitutiveError, fluid::viscoplastic::Viscoplastic, solid::Solid},
5    math::{
6        ContractFirstSecondIndicesWithSecondIndicesOf, ContractSecondIndexWithFirstIndexOf,
7        IDENTITY, Matrix, Rank2, TensorArray, TensorTuple, TensorTupleVec, Vector,
8        integrate::ExplicitInternalVariables,
9        optimize::{
10            EqualityConstraint, FirstOrderRootFinding, OptimizationError, ZerothOrderRootFinding,
11        },
12    },
13    mechanics::{
14        CauchyStress, CauchyTangentStiffness, Deformation, DeformationGradient,
15        DeformationGradientPlastic, DeformationGradients, FirstPiolaKirchhoffStress,
16        FirstPiolaKirchhoffTangentStiffness, Scalar, SecondPiolaKirchhoffStress,
17        SecondPiolaKirchhoffTangentStiffness, Times,
18    },
19};
20
21/// Elastic-viscoplastic state variables.
22pub type StateVariables = TensorTuple<DeformationGradientPlastic, Scalar>;
23
24/// Elastic-viscoplastic state variables history.
25pub type StateVariablesHistory = TensorTupleVec<DeformationGradientPlastic, Scalar>;
26
27/// Possible applied loads.
28pub enum AppliedLoad<'a> {
29    /// Uniaxial stress given $`F_{11}`$.
30    UniaxialStress(fn(Scalar) -> Scalar, &'a [Scalar]),
31    // /// Biaxial stress given $`F_{11}`$ and $`F_{22}`$.
32    // BiaxialStress(fn(Scalar) -> Scalar, fn(Scalar) -> Scalar, &'a [Scalar]),
33}
34
35/// Required methods for elastic-viscoplastic solid constitutive models.
36pub trait ElasticViscoplastic
37where
38    Self: Solid + Viscoplastic,
39{
40    /// Calculates and returns the Cauchy stress.
41    ///
42    /// ```math
43    /// \boldsymbol{\sigma} = \boldsymbol{\sigma}_\mathrm{e}
44    /// ```
45    fn cauchy_stress(
46        &self,
47        deformation_gradient: &DeformationGradient,
48        deformation_gradient_p: &DeformationGradientPlastic,
49    ) -> Result<CauchyStress, ConstitutiveError> {
50        Ok(deformation_gradient
51            * self.second_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_p)?
52            * deformation_gradient.transpose()
53            / deformation_gradient.determinant())
54    }
55    /// Calculates and returns the tangent stiffness associated with the Cauchy stress.
56    ///
57    /// ```math
58    /// \boldsymbol{\mathcal{T}} = \boldsymbol{\mathcal{T}}_\mathrm{e}\cdot\mathbf{F}_\mathrm{p}^{-T}
59    /// ```
60    fn cauchy_tangent_stiffness(
61        &self,
62        deformation_gradient: &DeformationGradient,
63        deformation_gradient_p: &DeformationGradientPlastic,
64    ) -> Result<CauchyTangentStiffness, ConstitutiveError> {
65        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
66        let cauchy_stress = self.cauchy_stress(deformation_gradient, deformation_gradient_p)?;
67        let some_stress = &cauchy_stress * &deformation_gradient_inverse_transpose;
68        Ok(self
69            .second_piola_kirchhoff_tangent_stiffness(deformation_gradient, deformation_gradient_p)?
70            .contract_first_second_indices_with_second_indices_of(
71                deformation_gradient,
72                deformation_gradient,
73            )
74            / deformation_gradient.determinant()
75            - CauchyTangentStiffness::dyad_ij_kl(
76                &cauchy_stress,
77                &deformation_gradient_inverse_transpose,
78            )
79            + CauchyTangentStiffness::dyad_il_kj(&some_stress, &IDENTITY)
80            + CauchyTangentStiffness::dyad_ik_jl(&IDENTITY, &some_stress))
81    }
82    /// Calculates and returns the first Piola-Kirchhoff stress.
83    ///
84    /// ```math
85    /// \mathbf{P} = \mathbf{P}_\mathrm{e}\cdot\mathbf{F}_\mathrm{p}^{-T}
86    /// ```
87    fn first_piola_kirchhoff_stress(
88        &self,
89        deformation_gradient: &DeformationGradient,
90        deformation_gradient_p: &DeformationGradientPlastic,
91    ) -> Result<FirstPiolaKirchhoffStress, ConstitutiveError> {
92        Ok(
93            self.cauchy_stress(deformation_gradient, deformation_gradient_p)?
94                * deformation_gradient.inverse_transpose()
95                * deformation_gradient.determinant(),
96        )
97    }
98    /// Calculates and returns the tangent stiffness associated with the first Piola-Kirchhoff stress.
99    ///
100    /// ```math
101    /// \mathcal{C}_{iJkL} = \mathcal{C}^\mathrm{e}_{iMkN} F_{MJ}^{\mathrm{p}-T} F_{NL}^{\mathrm{p}-T}
102    /// ```
103    fn first_piola_kirchhoff_tangent_stiffness(
104        &self,
105        deformation_gradient: &DeformationGradient,
106        deformation_gradient_p: &DeformationGradientPlastic,
107    ) -> Result<FirstPiolaKirchhoffTangentStiffness, ConstitutiveError> {
108        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
109        let first_piola_kirchhoff_stress =
110            self.first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_p)?;
111        Ok(self
112            .cauchy_tangent_stiffness(deformation_gradient, deformation_gradient_p)?
113            .contract_second_index_with_first_index_of(&deformation_gradient_inverse_transpose)
114            * deformation_gradient.determinant()
115            + FirstPiolaKirchhoffTangentStiffness::dyad_ij_kl(
116                &first_piola_kirchhoff_stress,
117                &deformation_gradient_inverse_transpose,
118            )
119            - FirstPiolaKirchhoffTangentStiffness::dyad_il_kj(
120                &first_piola_kirchhoff_stress,
121                &deformation_gradient_inverse_transpose,
122            ))
123    }
124    /// Calculates and returns the second Piola-Kirchhoff stress.
125    ///
126    /// ```math
127    /// \mathbf{S} = \mathbf{F}_\mathrm{p}^{-1}\cdot\mathbf{S}_\mathrm{e}\cdot\mathbf{F}_\mathrm{p}^{-T}
128    /// ```
129    fn second_piola_kirchhoff_stress(
130        &self,
131        deformation_gradient: &DeformationGradient,
132        deformation_gradient_p: &DeformationGradientPlastic,
133    ) -> Result<SecondPiolaKirchhoffStress, ConstitutiveError> {
134        Ok(deformation_gradient.inverse()
135            * self.first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_p)?)
136    }
137    /// Calculates and returns the tangent stiffness associated with the second Piola-Kirchhoff stress.
138    ///
139    /// ```math
140    /// \mathcal{G}_{IJkL} = \mathcal{G}^\mathrm{e}_{MNkO} F_{MI}^{\mathrm{p}-T} F_{NJ}^{\mathrm{p}-T} F_{OL}^{\mathrm{p}-T}
141    /// ```
142    fn second_piola_kirchhoff_tangent_stiffness(
143        &self,
144        deformation_gradient: &DeformationGradient,
145        deformation_gradient_p: &DeformationGradientPlastic,
146    ) -> Result<SecondPiolaKirchhoffTangentStiffness, ConstitutiveError> {
147        let deformation_gradient_inverse_transpose = deformation_gradient.inverse_transpose();
148        let deformation_gradient_inverse = deformation_gradient_inverse_transpose.transpose();
149        let second_piola_kirchhoff_stress =
150            self.second_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_p)?;
151        Ok(self
152            .cauchy_tangent_stiffness(deformation_gradient, deformation_gradient_p)?
153            .contract_first_second_indices_with_second_indices_of(
154                &deformation_gradient_inverse,
155                &deformation_gradient_inverse,
156            )
157            * deformation_gradient.determinant()
158            + SecondPiolaKirchhoffTangentStiffness::dyad_ij_kl(
159                &second_piola_kirchhoff_stress,
160                &deformation_gradient_inverse_transpose,
161            )
162            - SecondPiolaKirchhoffTangentStiffness::dyad_il_kj(
163                &second_piola_kirchhoff_stress,
164                &deformation_gradient_inverse_transpose,
165            )
166            - SecondPiolaKirchhoffTangentStiffness::dyad_ik_jl(
167                &deformation_gradient_inverse,
168                &second_piola_kirchhoff_stress,
169            ))
170    }
171    /// Calculates and returns the evolution of the state variables.
172    ///
173    /// ```math
174    /// \dot{\mathbf{F}}_\mathrm{p} = \mathbf{D}_\mathrm{p}\cdot\mathbf{F}_\mathrm{p}
175    /// ```
176    fn state_variables_evolution(
177        &self,
178        deformation_gradient: &DeformationGradient,
179        state_variables: &StateVariables,
180    ) -> Result<StateVariables, ConstitutiveError> {
181        let (deformation_gradient_p, yield_stress) = state_variables.into();
182        let jacobian = deformation_gradient.jacobian().unwrap();
183        let deformation_gradient_e = deformation_gradient * deformation_gradient_p.inverse();
184        let cauchy_stress = self.cauchy_stress(deformation_gradient, deformation_gradient_p)?;
185        let mandel_stress_e = (deformation_gradient_e.transpose()
186            * cauchy_stress
187            * deformation_gradient_e.inverse_transpose())
188            * jacobian;
189        let plastic_stretching_rate =
190            self.plastic_stretching_rate(mandel_stress_e.deviatoric(), *yield_stress)?;
191        Ok(StateVariables::from((
192            &plastic_stretching_rate * deformation_gradient_p,
193            self.yield_stress_evolution(&plastic_stretching_rate)?,
194        )))
195    }
196}
197
198/// Zeroth-order root-finding methods for elastic-viscoplastic solid constitutive models.
199pub trait ZerothOrderRoot {
200    /// Solve for the unknown components of the deformation gradients under an applied load.
201    ///
202    /// ```math
203    /// \mathbf{P}(\mathbf{F},\mathbf{F}_\mathrm{p}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
204    /// ```
205    fn root(
206        &self,
207        applied_load: AppliedLoad,
208        integrator: impl ExplicitInternalVariables<
209            StateVariables,
210            DeformationGradient,
211            StateVariablesHistory,
212            DeformationGradients,
213        >,
214        solver: impl ZerothOrderRootFinding<DeformationGradient>,
215    ) -> Result<(Times, DeformationGradients, StateVariablesHistory), ConstitutiveError>;
216    #[doc(hidden)]
217    fn root_inner_0(
218        &self,
219        deformation_gradient_p: &DeformationGradientPlastic,
220        equality_constraint: EqualityConstraint,
221        solver: &impl ZerothOrderRootFinding<DeformationGradient>,
222        initial_guess: &DeformationGradient,
223    ) -> Result<DeformationGradient, OptimizationError>;
224}
225
226/// First-order root-finding methods for elastic-viscoplastic solid constitutive models.
227pub trait FirstOrderRoot {
228    /// Solve for the unknown components of the deformation gradients under an applied load.
229    ///
230    /// ```math
231    /// \mathbf{P}(\mathbf{F},\mathbf{F}_\mathrm{p}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
232    /// ```
233    fn root(
234        &self,
235        applied_load: AppliedLoad,
236        integrator: impl ExplicitInternalVariables<
237            StateVariables,
238            DeformationGradient,
239            StateVariablesHistory,
240            DeformationGradients,
241        >,
242        solver: impl FirstOrderRootFinding<
243            FirstPiolaKirchhoffStress,
244            FirstPiolaKirchhoffTangentStiffness,
245            DeformationGradient,
246        >,
247    ) -> Result<(Times, DeformationGradients, StateVariablesHistory), ConstitutiveError>;
248    #[doc(hidden)]
249    fn root_inner_1(
250        &self,
251        deformation_gradient_p: &DeformationGradientPlastic,
252        equality_constraint: EqualityConstraint,
253        solver: &impl FirstOrderRootFinding<
254            FirstPiolaKirchhoffStress,
255            FirstPiolaKirchhoffTangentStiffness,
256            DeformationGradient,
257        >,
258        initial_guess: &DeformationGradient,
259    ) -> Result<DeformationGradient, OptimizationError>;
260}
261
262impl<T> ZerothOrderRoot for T
263where
264    T: ElasticViscoplastic,
265{
266    fn root(
267        &self,
268        applied_load: AppliedLoad,
269        integrator: impl ExplicitInternalVariables<
270            StateVariables,
271            DeformationGradient,
272            StateVariablesHistory,
273            DeformationGradients,
274        >,
275        solver: impl ZerothOrderRootFinding<DeformationGradient>,
276    ) -> Result<(Times, DeformationGradients, StateVariablesHistory), ConstitutiveError> {
277        match match applied_load {
278            AppliedLoad::UniaxialStress(deformation_gradient_11, time) => {
279                let mut matrix = Matrix::zero(4, 9);
280                let mut vector = Vector::zero(4);
281                matrix[0][0] = 1.0;
282                matrix[1][1] = 1.0;
283                matrix[2][2] = 1.0;
284                matrix[3][5] = 1.0;
285                integrator.integrate_and_evaluate(
286                    |_: Scalar,
287                     state_variables: &StateVariables,
288                     deformation_gradient: &DeformationGradient| {
289                        Ok(self.state_variables_evolution(deformation_gradient, state_variables)?)
290                    },
291                    |t: Scalar,
292                     state_variables: &StateVariables,
293                     deformation_gradient: &DeformationGradient| {
294                        let (deformation_gradient_p, _) = state_variables.into();
295                        vector[0] = deformation_gradient_11(t);
296                        Ok(self.root_inner_0(
297                            deformation_gradient_p,
298                            EqualityConstraint::Linear(matrix.clone(), vector.clone()),
299                            &solver,
300                            deformation_gradient,
301                        )?)
302                    },
303                    time,
304                    StateVariables::from((
305                        DeformationGradientPlastic::identity(),
306                        self.initial_yield_stress(),
307                    )),
308                    DeformationGradient::identity(),
309                )
310            }
311        } {
312            Ok((times, state_variables, _, deformation_gradients)) => {
313                Ok((times, deformation_gradients, state_variables))
314            }
315            Err(error) => Err(ConstitutiveError::Upstream(
316                format!("{error}"),
317                format!("{self:?}"),
318            )),
319        }
320    }
321    #[doc(hidden)]
322    fn root_inner_0(
323        &self,
324        deformation_gradient_p: &DeformationGradientPlastic,
325        equality_constraint: EqualityConstraint,
326        solver: &impl ZerothOrderRootFinding<DeformationGradient>,
327        initial_guess: &DeformationGradient,
328    ) -> Result<DeformationGradient, OptimizationError> {
329        solver.root(
330            |deformation_gradient: &DeformationGradient| {
331                Ok(self
332                    .first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_p)?)
333            },
334            initial_guess.clone(),
335            equality_constraint,
336        )
337    }
338}
339
340impl<T> FirstOrderRoot for T
341where
342    T: ElasticViscoplastic,
343{
344    fn root(
345        &self,
346        applied_load: AppliedLoad,
347        integrator: impl ExplicitInternalVariables<
348            StateVariables,
349            DeformationGradient,
350            StateVariablesHistory,
351            DeformationGradients,
352        >,
353        solver: impl FirstOrderRootFinding<
354            FirstPiolaKirchhoffStress,
355            FirstPiolaKirchhoffTangentStiffness,
356            DeformationGradient,
357        >,
358    ) -> Result<(Times, DeformationGradients, StateVariablesHistory), ConstitutiveError> {
359        match match applied_load {
360            AppliedLoad::UniaxialStress(deformation_gradient_11, time) => {
361                let mut matrix = Matrix::zero(4, 9);
362                let mut vector = Vector::zero(4);
363                matrix[0][0] = 1.0;
364                matrix[1][1] = 1.0;
365                matrix[2][2] = 1.0;
366                matrix[3][5] = 1.0;
367                integrator.integrate_and_evaluate(
368                    |_: Scalar,
369                     state_variables: &StateVariables,
370                     deformation_gradient: &DeformationGradient| {
371                        Ok(self.state_variables_evolution(deformation_gradient, state_variables)?)
372                    },
373                    |t: Scalar,
374                     state_variables: &StateVariables,
375                     deformation_gradient: &DeformationGradient| {
376                        let (deformation_gradient_p, _) = state_variables.into();
377                        vector[0] = deformation_gradient_11(t);
378                        Ok(self.root_inner_1(
379                            deformation_gradient_p,
380                            EqualityConstraint::Linear(matrix.clone(), vector.clone()),
381                            &solver,
382                            deformation_gradient,
383                        )?)
384                    },
385                    time,
386                    StateVariables::from((
387                        DeformationGradientPlastic::identity(),
388                        self.initial_yield_stress(),
389                    )),
390                    DeformationGradient::identity(),
391                )
392            }
393        } {
394            Ok((times, state_variables, _, deformation_gradients)) => {
395                Ok((times, deformation_gradients, state_variables))
396            }
397            Err(error) => Err(ConstitutiveError::Upstream(
398                format!("{error}"),
399                format!("{self:?}"),
400            )),
401        }
402    }
403    #[doc(hidden)]
404    fn root_inner_1(
405        &self,
406        deformation_gradient_p: &DeformationGradientPlastic,
407        equality_constraint: EqualityConstraint,
408        solver: &impl FirstOrderRootFinding<
409            FirstPiolaKirchhoffStress,
410            FirstPiolaKirchhoffTangentStiffness,
411            DeformationGradient,
412        >,
413        initial_guess: &DeformationGradient,
414    ) -> Result<DeformationGradient, OptimizationError> {
415        solver.root(
416            |deformation_gradient: &DeformationGradient| {
417                Ok(self
418                    .first_piola_kirchhoff_stress(deformation_gradient, deformation_gradient_p)?)
419            },
420            |deformation_gradient: &DeformationGradient| {
421                Ok(self.first_piola_kirchhoff_tangent_stiffness(
422                    deformation_gradient,
423                    deformation_gradient_p,
424                )?)
425            },
426            initial_guess.clone(),
427            equality_constraint,
428        )
429    }
430}
conspire/constitutive/solid/elastic_viscoplastic/mod.rs

conspire/constitutive/solid/elastic_viscoplastic/
mod.rs
