conspire/constitutive/solid/elastic_viscoplastic/
mod.rs

1//! Elastic-viscoplastic solid constitutive models.
2//!
3//! ---
4//!
5#![doc = include_str!("doc.md")]
6
7#[cfg(feature = "doc")]
8pub mod doc;
9
10mod canonical;
11
12use crate::{
13    constitutive::{
14        ConstitutiveError,
15        fluid::viscoplastic::{
16            Viscoplastic, ViscoplasticEvolution, ViscoplasticEvolutionHistory,
17            ViscoplasticStateVariables, ViscoplasticStateVariablesHistory,
18        },
19    },
20    math::{
21        ContractWith, Derivative, Differentiable, Intermediate, Quantity, Rank2, Reference, Scalar,
22        Tensor, TensorArray, TensorRank2, TensorTuple, TensorVec, Vector,
23        integrate::{
24            ButcherTableau, EmbeddedTableau, EvolvedIncrement, ExplicitDaeFirstOrderRoot,
25            ExplicitDaeZerothOrderRoot, Flat, Integrable, Product, StateEvolution, Unimodular,
26            integrate_rkmk_dae_adaptive_first_order_root, rkmk_dae_step_first_order_root,
27        },
28        optimize::{EqualityConstraint, FirstOrderRootFinding, ZerothOrderRootFinding},
29    },
30    mechanics::{
31        DeformationGradient, DeformationGradients, FirstPiolaKirchhoffStress,
32        FirstPiolaKirchhoffTangentStiffness, Times,
33    },
34    units::{Dissipation, Rate, Time},
35};
36use std::ops::{Add, Mul};
37
38use crate::constitutive::solid::elastic_plastic::bcs;
39pub use crate::constitutive::solid::elastic_plastic::{
40    AppliedLoad, ElasticPlasticOrViscoplastic, PlasticTangents,
41};
42
43/// Required methods for elastic-viscoplastic solid constitutive models.
44pub trait ElasticViscoplastic<Y>
45where
46    Self: ElasticPlasticOrViscoplastic + Viscoplastic<Y>,
47    Y: Differentiable + Tensor,
48{
49    /// Calculates and returns the internal dissipation.
50    ///
51    /// ```math
52    /// T\dot{s} = \mathbf{M}_\mathrm{e}':\mathbf{D}_\mathrm{p}
53    /// ```
54    fn internal_dissipation(
55        &self,
56        deformation_gradient: &DeformationGradient,
57        state_variables: &ViscoplasticStateVariables<Y>,
58    ) -> Result<Quantity<Dissipation>, ConstitutiveError> {
59        let deformation_gradient_p = &state_variables.0;
60        let plastic_stretching_rate = self
61            .state_variables_evolution(deformation_gradient, state_variables)?
62            .0
63            * deformation_gradient_p.inverse();
64        Ok(self
65            .mandel_stress(deformation_gradient, deformation_gradient_p)?
66            .deviatoric()
67            .contract_with(&plastic_stretching_rate))
68    }
69    /// Calculates and returns the evolution of the state variables.
70    fn state_variables_evolution(
71        &self,
72        deformation_gradient: &DeformationGradient,
73        state_variables: &ViscoplasticStateVariables<Y>,
74    ) -> Result<ViscoplasticEvolution<Y>, ConstitutiveError> {
75        self.plastic_evolution(
76            self.mandel_stress(deformation_gradient, &state_variables.0)?,
77            state_variables,
78        )
79    }
80}
81
82/// Zeroth-order root-finding methods for elastic-viscoplastic solid constitutive models.
83pub trait ZerothOrderRoot<Y>
84where
85    Y: Differentiable + Tensor,
86{
87    /// Solve for the unknown components of the deformation gradients under an applied load.
88    ///
89    /// ```math
90    /// \mathbf{P}(\mathbf{F},\mathbf{F}_\mathrm{p}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
91    /// ```
92    fn root(
93        &self,
94        applied_load: AppliedLoad,
95        integrator: impl ExplicitDaeZerothOrderRoot<
96            FirstPiolaKirchhoffStress,
97            ViscoplasticStateVariables<Y>,
98            DeformationGradient,
99            ViscoplasticStateVariablesHistory<Y>,
100            DeformationGradients,
101            ViscoplasticEvolutionHistory<Y>,
102        >,
103        solver: impl ZerothOrderRootFinding<FirstPiolaKirchhoffStress, DeformationGradient>,
104    ) -> Result<
105        (
106            Times,
107            DeformationGradients,
108            ViscoplasticStateVariablesHistory<Y>,
109        ),
110        ConstitutiveError,
111    >;
112}
113
114/// First-order root-finding methods for elastic-viscoplastic solid constitutive models.
115pub trait FirstOrderRoot<Y>
116where
117    Y: Differentiable + Tensor,
118{
119    /// Solve for the unknown components of the deformation gradients under an applied load.
120    ///
121    /// ```math
122    /// \mathbf{P}(\mathbf{F},\mathbf{F}_\mathrm{p}) - \boldsymbol{\lambda} - \mathbf{P}_0 = \mathbf{0}
123    /// ```
124    fn root(
125        &self,
126        applied_load: AppliedLoad,
127        integrator: impl ExplicitDaeFirstOrderRoot<
128            FirstPiolaKirchhoffStress,
129            FirstPiolaKirchhoffTangentStiffness,
130            ViscoplasticStateVariables<Y>,
131            DeformationGradient,
132            ViscoplasticStateVariablesHistory<Y>,
133            DeformationGradients,
134            ViscoplasticEvolutionHistory<Y>,
135        >,
136        solver: impl FirstOrderRootFinding<
137            FirstPiolaKirchhoffStress,
138            FirstPiolaKirchhoffTangentStiffness,
139            DeformationGradient,
140        >,
141    ) -> Result<
142        (
143            Times,
144            DeformationGradients,
145            ViscoplasticStateVariablesHistory<Y>,
146        ),
147        ConstitutiveError,
148    >;
149}
150
151impl<C, Y> ZerothOrderRoot<Y> for C
152where
153    C: ElasticViscoplastic<Y>,
154    Y: Differentiable + Tensor,
155{
156    fn root(
157        &self,
158        applied_load: AppliedLoad,
159        integrator: impl ExplicitDaeZerothOrderRoot<
160            FirstPiolaKirchhoffStress,
161            ViscoplasticStateVariables<Y>,
162            DeformationGradient,
163            ViscoplasticStateVariablesHistory<Y>,
164            DeformationGradients,
165            ViscoplasticEvolutionHistory<Y>,
166        >,
167        solver: impl ZerothOrderRootFinding<FirstPiolaKirchhoffStress, DeformationGradient>,
168    ) -> Result<
169        (
170            Times,
171            DeformationGradients,
172            ViscoplasticStateVariablesHistory<Y>,
173        ),
174        ConstitutiveError,
175    > {
176        let (matrix, prescribed, time) = bcs(applied_load);
177        let mut vector = Vector::zero(matrix.len());
178        let (times, state_variables, _, deformation_gradients) = integrator
179            .integrate(
180                |_: Quantity<Time>,
181                 state_variables: &ViscoplasticStateVariables<Y>,
182                 deformation_gradient: &DeformationGradient| {
183                    Ok(self.state_variables_evolution(deformation_gradient, state_variables)?)
184                },
185                |_: Quantity<Time>,
186                 state_variables: &ViscoplasticStateVariables<Y>,
187                 deformation_gradient: &DeformationGradient| {
188                    let deformation_gradient_p = &state_variables.0;
189                    Ok(self.first_piola_kirchhoff_stress(
190                        deformation_gradient,
191                        deformation_gradient_p,
192                    )?)
193                },
194                solver,
195                time,
196                (self.initial_state(), DeformationGradient::identity()),
197                |t: Quantity<Time>| {
198                    prescribed
199                        .iter()
200                        .for_each(|(index, function)| vector[*index] = function(t));
201                    EqualityConstraint::Linear(matrix.clone(), vector.clone())
202                },
203            )
204            .map_err(|error| ConstitutiveError::upstream(error, self))?;
205        Ok((times, deformation_gradients, state_variables))
206    }
207}
208
209impl<C, Y> FirstOrderRoot<Y> for C
210where
211    C: ElasticViscoplastic<Y>,
212    Y: Differentiable + Tensor,
213{
214    fn root(
215        &self,
216        applied_load: AppliedLoad,
217        integrator: impl ExplicitDaeFirstOrderRoot<
218            FirstPiolaKirchhoffStress,
219            FirstPiolaKirchhoffTangentStiffness,
220            ViscoplasticStateVariables<Y>,
221            DeformationGradient,
222            ViscoplasticStateVariablesHistory<Y>,
223            DeformationGradients,
224            ViscoplasticEvolutionHistory<Y>,
225        >,
226        solver: impl FirstOrderRootFinding<
227            FirstPiolaKirchhoffStress,
228            FirstPiolaKirchhoffTangentStiffness,
229            DeformationGradient,
230        >,
231    ) -> Result<
232        (
233            Times,
234            DeformationGradients,
235            ViscoplasticStateVariablesHistory<Y>,
236        ),
237        ConstitutiveError,
238    > {
239        let (matrix, prescribed, time) = bcs(applied_load);
240        let mut vector = Vector::zero(matrix.len());
241        let (times, state_variables, _, deformation_gradients) = integrator
242            .integrate(
243                |_: Quantity<Time>,
244                 state_variables: &ViscoplasticStateVariables<Y>,
245                 deformation_gradient: &DeformationGradient| {
246                    Ok(self.state_variables_evolution(deformation_gradient, state_variables)?)
247                },
248                |_: Quantity<Time>,
249                 state_variables: &ViscoplasticStateVariables<Y>,
250                 deformation_gradient: &DeformationGradient| {
251                    let deformation_gradient_p = &state_variables.0;
252                    Ok(self.first_piola_kirchhoff_stress(
253                        deformation_gradient,
254                        deformation_gradient_p,
255                    )?)
256                },
257                |_: Quantity<Time>,
258                 state_variables: &ViscoplasticStateVariables<Y>,
259                 deformation_gradient: &DeformationGradient| {
260                    let deformation_gradient_p = &state_variables.0;
261                    Ok(self.first_piola_kirchhoff_tangent_stiffness(
262                        deformation_gradient,
263                        deformation_gradient_p,
264                    )?)
265                },
266                solver,
267                time,
268                (self.initial_state(), DeformationGradient::identity()),
269                |t: Quantity<Time>| {
270                    prescribed
271                        .iter()
272                        .for_each(|(index, function)| vector[*index] = function(t));
273                    EqualityConstraint::Linear(matrix.clone(), vector.clone())
274                },
275            )
276            .map_err(|error| ConstitutiveError::upstream(error, self))?;
277        Ok((times, deformation_gradients, state_variables))
278    }
279}
280
281/// The internal state `(F_p, Y)` evolves as `F_p` on the unimodular group
282/// (`Reference → Intermediate`, so its algebra element `D_p Δt` is
283/// `Intermediate → Intermediate`) and the hardening variable `Y` additively. The
284/// rate is `(D_p, Ẏ)`, from the model's [`plastic_evolution`] (`D_p` recovered as
285/// `Ḟ_p F_p⁻¹`), driven by the total deformation gradient through the Mandel
286/// stress. Blanket over any [`ElasticViscoplastic`] model — not
287/// `Canonical`-specific — so a hybrid composition gets it automatically as
288/// soon as it implements [`ElasticViscoplastic<Y>`].
289///
290/// [`plastic_evolution`]: Viscoplastic::plastic_evolution
291impl<C, Y> StateEvolution<Time, Y> for C
292where
293    C: ElasticViscoplastic<Y>,
294    Y: Clone + Differentiable<Time> + Tensor,
295    for<'a> Y: Add<&'a Y, Output = Y>,
296    TensorTuple<TensorRank2<3, Intermediate, Intermediate>, Y>: Differentiable<
297            Time,
298            Derivative = TensorTuple<
299                TensorRank2<3, Intermediate, Intermediate, Rate>,
300                Derivative<Y>,
301            >,
302        >,
303{
304    type Field = Product<Unimodular<Intermediate, Reference>, Flat<Y>>;
305    type Drive = DeformationGradient;
306    fn initial_state(&self) -> ViscoplasticStateVariables<Y> {
307        <Self as Viscoplastic<Y>>::initial_state(self)
308    }
309    fn state_rate(
310        &self,
311        _time: Quantity<Time>,
312        deformation_gradient: &DeformationGradient,
313        state: &ViscoplasticStateVariables<Y>,
314    ) -> Result<Derivative<<Self::Field as Integrable>::Increment, Time>, String> {
315        let mandel_stress = self.mandel_stress(deformation_gradient, &state.0)?;
316        let evolution = self.plastic_evolution(mandel_stress, state)?;
317        let plastic_stretching_rate = evolution.0 * state.0.inverse();
318        Ok(TensorTuple(plastic_stretching_rate, evolution.1))
319    }
320}
321
322/// RKMK-DAE stage-equilibrium methods for elastic-viscoplastic solid constitutive
323/// models. The sibling of [`FirstOrderRoot`] that keeps `F_p` on its manifold
324/// instead of marching it additively, by resolving `F` from equilibrium at
325/// every stage abscissa rather than freezing it across the window — so the
326/// coupling is the tableau's own order, not first order. Blanket over any
327/// [`ElasticViscoplastic`] model, same as [`FirstOrderRoot`] itself.
328pub trait RootRkmkDae<Y>
329where
330    Y: Differentiable + Tensor,
331{
332    /// `F` is re-solved from equilibrium at every stage abscissa of the
333    /// window while `F_p` advances on its group, generic over the hardening
334    /// variable `Y`. Stage `i` reconstructs `F_p` at `σᵢ`, solves
335    /// `P(F, F_p^i) - λ(t + cᵢ Δt) - P_0 = 0` for `F` there, and evaluates the
336    /// plastic rate at that consistent pair.
337    ///
338    /// This is the half-explicit RK treatment of the index-1 DAE that
339    /// [`FirstOrderRoot::root`] already performs, with the state leg moved off
340    /// the additive march onto `expm`/`dexpinv` — so `det F_p = 1` is kept
341    /// rather than drifting.
342    fn root_rkmk_dae<Tab: ButcherTableau>(
343        &self,
344        applied_load: AppliedLoad,
345        solver: impl FirstOrderRootFinding<
346            FirstPiolaKirchhoffStress,
347            FirstPiolaKirchhoffTangentStiffness,
348            DeformationGradient,
349        >,
350    ) -> Result<
351        (
352            Times,
353            DeformationGradients,
354            ViscoplasticStateVariablesHistory<Y>,
355        ),
356        ConstitutiveError,
357    >;
358    /// As [`Self::root_rkmk_dae`], but the whole span is stepped under
359    /// embedded (`Tab::D`) error control rather than on the supplied load
360    /// grid.
361    ///
362    /// Two times in `applied_load` give only the span, and the controller's
363    /// own accepted steps are reported. More than two are requested report
364    /// times — the convention of the flat DAE loop — and `F_p` is served at
365    /// each from the geodesic `HermiteSegment` of the accepted step
366    /// containing it, so it is on the unimodular group at every reported time
367    /// and not just at the accepted ones; `F` is then re-solved from
368    /// equilibrium there.
369    fn root_rkmk_dae_adaptive<Tab: EmbeddedTableau>(
370        &self,
371        applied_load: AppliedLoad,
372        solver: impl FirstOrderRootFinding<
373            FirstPiolaKirchhoffStress,
374            FirstPiolaKirchhoffTangentStiffness,
375            DeformationGradient,
376        >,
377        abs_tol: Scalar,
378        rel_tol: Scalar,
379    ) -> Result<
380        (
381            Times,
382            DeformationGradients,
383            ViscoplasticStateVariablesHistory<Y>,
384        ),
385        ConstitutiveError,
386    >;
387}
388
389impl<C, Y> RootRkmkDae<Y> for C
390where
391    C: ElasticViscoplastic<Y>
392        + StateEvolution<
393            Time,
394            Y,
395            Drive = DeformationGradient,
396            Field: Integrable<Point = ViscoplasticStateVariables<Y>>,
397        >,
398    Y: Differentiable + Tensor,
399    EvolvedIncrement<C, Time, Y>: Clone + Differentiable<Time>,
400    for<'a> &'a Derivative<EvolvedIncrement<C, Time, Y>, Time>:
401        Mul<Quantity<Time>, Output = EvolvedIncrement<C, Time, Y>>,
402{
403    #[allow(clippy::type_complexity)]
404    fn root_rkmk_dae<Tab: ButcherTableau>(
405        &self,
406        applied_load: AppliedLoad,
407        solver: impl FirstOrderRootFinding<
408            FirstPiolaKirchhoffStress,
409            FirstPiolaKirchhoffTangentStiffness,
410            DeformationGradient,
411        >,
412    ) -> Result<
413        (
414            Times,
415            DeformationGradients,
416            ViscoplasticStateVariablesHistory<Y>,
417        ),
418        ConstitutiveError,
419    > {
420        let (matrix, prescribed, time) = bcs(applied_load);
421        let mut state = <Self as StateEvolution<Time, Y>>::initial_state(self);
422        let mut scratch = Vec::new();
423        let equality_constraint = |t: Quantity<Time>| {
424            let mut vector = Vector::zero(matrix.len());
425            prescribed
426                .iter()
427                .for_each(|(index, function)| vector[*index] = function(t));
428            EqualityConstraint::Linear(matrix.clone(), vector)
429        };
430        let function = |_: Quantity<Time>,
431                        state: &ViscoplasticStateVariables<Y>,
432                        deformation_gradient: &DeformationGradient|
433         -> Result<FirstPiolaKirchhoffStress, String> {
434            Ok(self.first_piola_kirchhoff_stress(deformation_gradient, &state.0)?)
435        };
436        let jacobian = |_: Quantity<Time>,
437                        state: &ViscoplasticStateVariables<Y>,
438                        deformation_gradient: &DeformationGradient|
439         -> Result<FirstPiolaKirchhoffTangentStiffness, String> {
440            Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient, &state.0)?)
441        };
442        let mut deformation_gradient = solver
443            .root(
444                |deformation_gradient: &DeformationGradient| {
445                    function(time[0], &state, deformation_gradient)
446                },
447                |deformation_gradient: &DeformationGradient| {
448                    jacobian(time[0], &state, deformation_gradient)
449                },
450                DeformationGradient::identity(),
451                equality_constraint(time[0]),
452                None,
453            )
454            .map_err(|error| ConstitutiveError::upstream(String::from(error), self))?;
455        let mut times = Times::new();
456        let mut deformation_gradients = DeformationGradients::new();
457        let mut state_variables = ViscoplasticStateVariablesHistory::new();
458        let mut carry = None;
459        times.push(time[0]);
460        deformation_gradients.push(deformation_gradient.clone());
461        state_variables.push(state.clone());
462        for step in time.windows(2) {
463            let advanced = rkmk_dae_step_first_order_root::<
464                <Self as StateEvolution<Time, Y>>::Field,
465                Tab,
466                FirstPiolaKirchhoffStress,
467                FirstPiolaKirchhoffTangentStiffness,
468                DeformationGradient,
469                Time,
470            >(
471                &mut |t, state, deformation_gradient| {
472                    self.state_rate(t, deformation_gradient, state)
473                },
474                function,
475                jacobian,
476                &solver,
477                &state,
478                &deformation_gradient,
479                step[0],
480                step[1] - step[0],
481                &mut scratch,
482                carry.as_ref(),
483                equality_constraint,
484            )
485            .map_err(|error| ConstitutiveError::upstream(error, self))?;
486            state = advanced.0;
487            deformation_gradient = advanced.1;
488            carry = advanced.2;
489            times.push(step[1]);
490            deformation_gradients.push(deformation_gradient.clone());
491            state_variables.push(state.clone());
492        }
493        Ok((times, deformation_gradients, state_variables))
494    }
495    #[allow(clippy::type_complexity)]
496    fn root_rkmk_dae_adaptive<Tab: EmbeddedTableau>(
497        &self,
498        applied_load: AppliedLoad,
499        solver: impl FirstOrderRootFinding<
500            FirstPiolaKirchhoffStress,
501            FirstPiolaKirchhoffTangentStiffness,
502            DeformationGradient,
503        >,
504        abs_tol: Scalar,
505        rel_tol: Scalar,
506    ) -> Result<
507        (
508            Times,
509            DeformationGradients,
510            ViscoplasticStateVariablesHistory<Y>,
511        ),
512        ConstitutiveError,
513    > {
514        let (matrix, prescribed, time) = bcs(applied_load);
515        let state = <Self as StateEvolution<Time, Y>>::initial_state(self);
516        let equality_constraint = |t: Quantity<Time>| {
517            let mut vector = Vector::zero(matrix.len());
518            prescribed
519                .iter()
520                .for_each(|(index, function)| vector[*index] = function(t));
521            EqualityConstraint::Linear(matrix.clone(), vector)
522        };
523        let function = |_: Quantity<Time>,
524                        state: &ViscoplasticStateVariables<Y>,
525                        deformation_gradient: &DeformationGradient|
526         -> Result<FirstPiolaKirchhoffStress, String> {
527            Ok(self.first_piola_kirchhoff_stress(deformation_gradient, &state.0)?)
528        };
529        let jacobian = |_: Quantity<Time>,
530                        state: &ViscoplasticStateVariables<Y>,
531                        deformation_gradient: &DeformationGradient|
532         -> Result<FirstPiolaKirchhoffTangentStiffness, String> {
533            Ok(self.first_piola_kirchhoff_tangent_stiffness(deformation_gradient, &state.0)?)
534        };
535        let deformation_gradient = solver
536            .root(
537                |deformation_gradient: &DeformationGradient| {
538                    function(time[0], &state, deformation_gradient)
539                },
540                |deformation_gradient: &DeformationGradient| {
541                    jacobian(time[0], &state, deformation_gradient)
542                },
543                DeformationGradient::identity(),
544                equality_constraint(time[0]),
545                None,
546            )
547            .map_err(|error| ConstitutiveError::upstream(String::from(error), self))?;
548        let (times, state_variables, deformation_gradients) =
549            integrate_rkmk_dae_adaptive_first_order_root::<
550                <Self as StateEvolution<Time, Y>>::Field,
551                Tab,
552                FirstPiolaKirchhoffStress,
553                FirstPiolaKirchhoffTangentStiffness,
554                DeformationGradient,
555                ViscoplasticStateVariablesHistory<Y>,
556                DeformationGradients,
557                Time,
558            >(
559                |t, state, deformation_gradient| self.state_rate(t, deformation_gradient, state),
560                function,
561                jacobian,
562                &solver,
563                time,
564                (state, deformation_gradient),
565                abs_tol,
566                rel_tol,
567                equality_constraint,
568            )
569            .map_err(|error| ConstitutiveError::upstream(error, self))?;
570        Ok((times, deformation_gradients, state_variables))
571    }
572}
conspire/constitutive/solid/elastic_viscoplastic/mod.rs

conspire/constitutive/solid/elastic_viscoplastic/
mod.rs
