Inelasticity

Overview

Ratel provides a general framework for constructing inelastic constitutive models of isotropic materials through the parallel assembly of Hooke (elastic), Maxwell (viscoelastic), and Prandtl (elasto-plastic) rheological elements. These elements share common helper functions for the elastic predictor, the update of stored kinematic quantities, and the consistent tangent, while differing in their specific stress update algorithms [PD03]. By adopting Hencky strain energy functionals and restricting updates to deviatoric components in principal directions, Ratel ensures an efficient implementation of models that can capture a wide range of isotropic mechanical responses with a minimal set of material parameters.

Thermodynamic basis

Supported by experimental observations of plastic deformation in crystalline materials (e.g., metals), the deformation gradient tensor \(\bm{F}= \bm{I}+\nabla_X\bm{u}\) can be conveniently decomposed into inelastic \(\bm{F}^i\) and elastic \(\bm{F}^e\) parts such that [EH69]

(252)\[ \bm{F}=\bm{F}^e\bm{F}^i \]

Viscoelastic (\(\bm{F}^v=\bm{F}^i\)) or elastoplastic (\(\bm{F}^p=\bm{F}^i\)) deformations are assumed isochoric (i.e., \(\det\bm{F}^i=1\)) and mapping material points to an unstressed intermediate configuration. Hence, \(\bm{F}^e = \textbf{v}^e \bm{R}^e\) maps the intermediate configuration to the spatial configuration via rigid rotations \(\bm{R}^e\) and elastic distortions, \(\textbf{v}^e\). With (252), the spatial velocity gradient \(\bm{l}\equiv \frac{\partial\bm{v}}{\partial\bm{x}} =\bm{\dot{F}}\bm{F}^{-1}\), can be additively decomposed as

(253)\[ \bm{l}=\bm{l}^e+\bm{F}^e\bm{L}^i(\bm{F}^e)^{-1} \]

where \(\bm{l}^e\equiv\dot{\bm{F}^e}(\bm{F}^e)^{-1}\) and \(\bm{L}^i\equiv\dot{\bm{F}^i}(\bm{F}^i)^{-1}\), which lives in the intermediate configuration. The velocity gradient tensor \(\bm{l}\) can be decomposed into its symmetric \(\bm{d}\equiv\text{sym}(\bm{l})\) and skew \(\bm{w}\equiv\text{skew}(\bm{l})\) parts. So defined, the elastic stretch tensor \(\bm{d}\) measures the instantaneous rate of stretching along the orthonormal eigenvectors, whereas \(\bm{w}\) represents the triad’s rigid spinning, which is usually neglected (e.g. by setting \(\bm{w} = \bm{0}\)) in inelastic isotropy.

With \(\bm{l}\) in (253), the Lie derivative of the left Cauchy-Green tensor \(\bm{b}^e=\bm{F}^e\bm{F}^{eT}\) with respect to the velocity field \(\bm{v}\) can be expressed as,

(254)\[ \mathcal{L}_{\bm{v}} \bm{b}^e \coloneqq \bm{F}\frac{\partial}{\partial t} \left( \bm{F}^{-1} \bm{b}^e \bm{F}^{-T} \right) \bm{F}^T = \dot{\bm{b}}^e - \bm{l}\bm{b}^e - \bm{b}^e\bm{l}^T \]

Having described the relevant kinematic, we can now consider the (isothermal) Clausius-Duhem inequality

(255)\[ \bm{\tau}: \bm{d} - \dot{\psi} \geq 0 \]

where \(\bm{\tau}\) is the Kirchhoff stress work conjugate of \(\bm{d}\) and \(\psi(\bm{b^e}, \bm{\alpha})\) a free energy function of elastic strains and a set of internal variables \(\bm{\alpha}\) associated with dissipative mechanisms [GFA10]. As derived in [JCRL85], the rate of change can be reformulated using (254) as

(256)\[ \dot{\psi}(\bm{b}^e, \bm{\xi}) = \frac{\partial \psi}{\partial \bm{b}^e} : \dot{\bm{b}}^e + \frac{\partial \psi}{\partial \bm{\xi}^e} \cdot \dot{\bm{\xi}} = 2\frac{\partial \psi}{\partial \bm{b}^e} \bm{b}^e : (\bm{d} + \frac{1}{2}(\mathcal{L}_{\bm{v}} \bm{b}^e) \bm{b}^{e-1}) + \frac{\partial \psi}{\partial \bm{\alpha}} \cdot \dot{\bm{\alpha}} \]

Substituting (256) into (255) and rearranging terms gives

(257)\[ (\bm{\tau} - 2\frac{\partial \psi}{\partial \bm{b}^e} \bm{b}^e) : \bm{d} + 2\frac{\partial \psi}{\partial \bm{b}^e} \bm{b}^e : (-\frac{1}{2}(\mathcal{L}_{\bm{v}} \bm{b}^e) \bm{b}^{e-1}) - \frac{\partial \psi}{\partial \bm{\alpha}} \cdot \dot{\bm{\alpha}}\geq 0 \]

Standard arguments then provide the constitutive equation for \(\bm\tau\) and simplify the inequality constraint,

(258)\[ \bm{\tau}=2\frac{\partial\psi}{\partial\bm{b}^e}\bm{b}^e, \qquad \bm{\tau} : (-\frac{1}{2}(\mathcal{L}_{\bm{v}} \bm{b}^e) \bm{b}^{e-1}) - \frac{\partial \psi}{\partial \bm{\alpha}} \cdot \dot{\bm{\alpha}}\geq 0 \]

Finally, following the pioneer work in [RG98], a isotropic fourth order and positive definite tensor \(\mathsf{V}\) can be introduced such that

(259)\[ -\frac{1}{2}(\mathcal{L}_{\bm{v}} \bm{b}^e) \bm{b}^{e-1} \coloneqq \mathsf{V}:\bm{\tau} \]

Additive decomposition of strain

As an alternative to the multiplicative decomposition, we can take an additive decomposition of the strain tensor which has comparable accuracy at small strains [MAL02]. However, this is true only for co-axial loading conditions (e.g. pure shear). For studies on the softening effect of additive plasticity in non-coaxial loading conditions see [PDME92] and the more recent review in [FMS22].

Stress Update

As \(\bm{F}\) is assumed constant at each Newton iteration of the non-linear solve over the time increment \(\Delta t_{n+1} = t_{n+1} - t_n\), \(\bm{l} = \bm{0}\) and (254) is simply \(\mathcal{L}_{\bm{v}} \bm{b}^e = \dot{\bm{b}}^e\). Substituting in (259), and rearranging terms, gives

(260)\[ \dot{\bm{b}^e} = -2 (\mathsf{V}:\bm{\tau})\bm{b}^e \]

which can be solved using a backward Euler approximation based on the exponential mapping integrator, [EB90] and [WA90],

(261)\[ \bm{b}^e= \exp \left( -2 \Delta t\mathsf{V}:\bm{\tau} \right) \cdot \bm{b}^{e \space tr} \]

Where \(tr\) indicate trial quantities associated with the inadmissible state that generally results from taking a purely elastic increment as initial guess. The trail left-Cauchy strain tensor \(\bm{b}^{e \space tr}\) can be computed after retrieving the inverse of the (symmetric) right Cauchy strain tensor from the previously converged time step, \(\bm{C}^{i-1}_n\),

(262)\[ \bm{b}^{e \space tr}=\bm{F}\bm{C}^{i-1}_n \bm{F}^{T} \]

The correction in (261) further simplifies when applied to the logarithmic strain. In particular, this can be limited to the deviatoric component \(\textbf{e}^{e}_{d}\), as we only consider deviatoric inelastic deformations,

(263)\[ \textbf{e}_{d}^e= \textbf{e}^{e \space tr}_{d} - \Delta t \mathsf{V}:\bm{\tau} \]

Using (132), the Kirchhoff stress can be finally computed as

(264)\[ \bm{\tau} = 2\mu \textbf{e}_{d}^e + \bulk\trace(\textbf{e}) \]

Note

The exponential tensor function has desirable properties in the numerical integration of constitutive laws of (both isotropic and anisotropic) elastoplasticity [GFA10]. In particular, it maps traceless tensor,\(\space\bm{X}\), to unimodular tensors for which \(\det{\bm{X}}=1\). The function therefore preserves the isochoric properties of plastic deformation, the loss of which would otherwise require (nonphysical) ad-hoc corrections [BHL+16].

Assembly of Rheological Elements

The above results can be generalized for rheological components arranged in parallel as the total free energy can be expressed as the sum of the free energies of the single components \(\psi = \sum_k \psi_k (\bm{b}^e_k, \bm{\alpha}_k)\) and \(\bm{\tau} = \sum_k \bm{\tau}_k\) ,

(265)\[ \bm{\tau}_k = 2\frac{\partial\psi_k}{\partial\bm{b}^e_k}\bm{b}^e_k, \qquad \bm{\tau}_k : (-\frac{1}{2}(\mathcal{L}_{\bm{v}} \bm{b}^e_k) \bm{b}^{e-1}_k) - \frac{\partial \psi_k}{\partial \bm{\alpha}_k} \cdot \dot{\bm{\alpha}}_k\geq 0 \]

Evolutions laws for strain and internal variables specific to each rheological element must be prescribed to ensure thermodynamic consistency. Linearization of \(\bm\tau_k\) with respect to \(\bm{F}\) is also straightforward as it can be obtained from the contribution of single components. These share the format of the Hooke (spring) element in (133)

(266)\[ \qquad \frac{\partial\bm{\tau}_k}{\partial\bm{F}}=\mathsf{C}^i_k:\mathsf{N}_k:\mathsf{B}_k \]

\(\mathsf C^i\) is the only quantity in (266) that depends on the specific material model. In the cases considered below (e.g. associative plasticity), it retains its symmetry. To ensure quadratic convergence of the Newton’s method, it is computed by applying the linearization to the algorithmic expression for the stress update [JCRL85].

Maxwell element (viscoelasticity)

To trivially satisfy thermodynamic consistency in (258), [RG98] prescribe no internal variables (i.e. \(\bm{\alpha} = \bm{0}\)) and

(267)\[ \mathsf{V}:\bm{\tau} = \frac{1}{2\eta_d}\bm{\tau}_d \]

where \(\eta_d\) is the deviatoric viscosity. Equation (263) becomes

(268)\[ \textbf{e}_{d}^e= \textbf{e}^{e \space tr}_{d} - \frac{\Delta t}{\eta_d}\bm{\tau}_d \]

As \(\bm{\tau}_d\) is generally a non-linear function of \(\textbf{e}_{d}^e\), (268) is commonly solved via local Newton iterations. Nevertheless, (268) reduces to a closed-form update thanks to the adoption of the Hencky model [FPS06]

(269)\[ \textbf{e}_{d}^e= \left( 1 + \frac{\Delta t}{\eta_d} \right)^{-1} \textbf{e}^{e \space tr}_{d} \]

and, consequently, to a straightforward computation of \(\diff\textbf{e}_{d}^e\) in (137)

(270)\[ \diff\textbf{e}_{d}^e= \left( 1 + \frac{\Delta t}{\eta_d} \right)^{-1} \diff\textbf{e}^{e \space tr}_{d} \]

Command-line interface

To enable the viscoelastic Hencky model, use the model option -model viscoelasticity-hencky-current-principal or viscoelasticity-hencky-initial-principal and set the material parameters listed in Viscoelastic Hencky model options. Any parameter without a default option is required.

Table 31 Viscoelastic Hencky model options

Option	Description
-model viscoelasticity-hencky-current-principal | viscoelasticity-hencky-initial-principal	Required to enable the neo-Hookean model, | denotes an option for the user.
-E [real]	Young’s modulus, \(E > 0\)
-nu [real]	Poisson’s ratio, \(\nu \leq 0.5\).
deviatoric_viscosity [real]	Deviatoric viscosity, \(\eta_d > 0\)

An example using the viscoelastic Hencky model can be run via

$ ./bin/ratel-quasistatic -options_file ex02-quasistatic-viscoelasticity-hencky-current-principal.yml

Prandtl element

Plasticity relies on the definition of a non-negative, rank-one convex dissipation potential \(\Xi(\bm\tau, A, \bm\alpha)\) to formulate constitutive laws governing the evolution of internal variables and plastic flow. In associative plasticity, the dissipation potential coincides with the yield function, and the principle of maximum plastic dissipation implies that the plastic flow vector (\(\bm{N}^p\)) is normal to the yield surface [BN75]. We can then write

(271)\[ \bm{N}^p \coloneqq \frac{\partial\Xi}{\partial\bm{\tau}} = \frac{\partial\phi}{\partial\bm{\tau}} \]

To ensure thermodynamic principles in (258) are satisfied, we can set [PD03]

(272)\[ -\frac{1}{2}(\mathcal{L}_{\bm{v}} \bm{b}^e) \bm{b}^{e-1} =\dot\gamma\frac{\partial\phi}{\partial\bm{\tau}}=\dot\gamma\bm{N}^p, \qquad \dot{\bm{\alpha}}=\dot\gamma\frac{\partial\phi}{\partial\bm{A}}, \]

where \(\dot\gamma\) is the plastic multiplier calculated through optimality conditions,

(273)\[ \dot\gamma\geq0 \qquad \phi\leq0 \qquad \dot\gamma\phi=0 \]

and \(\bm{A}\) the set of conjugate forces to the internal variables,

(274)\[ \bm{A}=-\frac{\partial\psi}{\partial\bm{\alpha}} \]

Note

The second condition in (273) defines the set of admissible (i.e., physically possible given the current value of the internal variables) stress states \(\mathscr{S}\) as

(275)\[ \mathscr{S} = \{ \phi (\bm{\tau}, \bm{\textbf{A}}, \bm{\alpha}) \leq 0 \} \]

In stress space, \(\phi = 0\) defines a yield surface that delimits the elastic domain. Plastic deformation occurs when stress states lie on such surface and the flow vector points outwards. However, the incremental nature of applied load/displacements may cause stress states to fall outside \(\mathscr{S}\). To ensure admissibility, these are projected back to the yield surface.

Motivated by experimental observations of pressure insensitivity of yielding in metals, plastic deformation in von Mises plasticity occurs when the deviatoric strain energy exceeds a critical value. This energy can be expressed in terms of the second invariant

(276)\[ J_2=\frac{1}{2}\bm{s}:\bm{s} \]

of the (symmetric) deviatoric stress tensor \(\bm{s}=I_d(\bm{\sigma})\). The critical value is given by the strain hardening law \(\sigma_y=\sigma_y(\bar{\varepsilon}^p_n)\) that is function of a single (scalar) internal variable termed the accumulated (effective or equivalent) plastic strain \(\bar{\varepsilon}^p\) and defined as:

(277)\[ \bar{\varepsilon}^p=\int_0^t |\dot{\gamma}|dt \]

As uniaxial tensile testing is the preferred method to calibrate plasticity models, the von Mises yield function is generally stated in terms of the yield stress \(\sigma_y\) as:

(278)\[ \phi=q-\sigma_y \]

with (scalar) \(q=\sqrt{3J_2({\bm{s}})}= \sqrt{\frac{3}{2}}||\bm{s}||\) termed the von Mises (effective or equivalent) stress. If plastic deformation occurs, the flow vector takes the form,

(279)\[ \bm{N}^p =\sqrt{\frac{3}{2}}\frac{\bm{s}}{||\bm{s}||} = \frac{3}{2}\frac{\bm{s}}{q} \]

Note

The implemented strain hardening law combines a linear hardening with Voce-type nonlinear hardening as

(280)\[ \sigma_y= \sigma_0 + H\bar{\varepsilon}^p + (\sigma_\infty - \sigma_0) (1-\exp(-\beta\bar{\varepsilon}^p)) \]

where \(\sigma_0\) is the yield stress and {\(H\), \(\sigma_\infty\), \(\beta\)} are the three strain hardening parameters respectively representing the linear hardening factor, the saturation (flow) stress, and the hardening decay parameter.

Stress Update

The stress updating procedure reuses the returning mapping algorithm for small-strains with additional preprocessor and postprocessor steps to map strains and stresses from Euclidean to logarithmic strain spaces and vice versa [C94]. An overview of the integration algorithm with state variable \(\bar\varepsilon^{p}_n\) and stored kinematic quantity \(\bm{C}^{p-1}_n\) is given below.

The trial effective stress \(q^{tr}\) can be calculated at \(t+\Delta t\) as:

(281)\[ \bm{b}^{e \space tr}=\bm{F}\bm{C}^{p-1}_n \bm{F}^{T} \space \rightarrow \space \textbf{e}^{e \space tr}=0.5\log(\bm{b}^{e \space tr}) \space \rightarrow \bm{s}^{tr} = 2\mu \textbf{e}^{e \space tr}_{d }\space \rightarrow \space q^{tr}= \sqrt{3J_2({\bm{s}^{tr}})} \]

Admissibility can thus be verified after the initial guess \(\bar\varepsilon^{p} =\bar{\varepsilon}^p_n\) :

(282)\[ \begin{aligned} & \text{if:} \qquad q^{tr}-\sigma_y(\bar\varepsilon^{p}_{n})\leq0 \qquad \text{elastic step} \\ &\text {else: \quad elastoplastic step} \rightarrow \text{return mapping} \end{aligned} \]

Note

The logarithmic strain in (281) is computed in a numerically stable manner, as shown in (129).

In the von Mises model, the equations of the return mapping can be reduced to a single non-linear equation [Sim92]

(283)\[ \tilde{\phi}=q^{tr}-3\mu\Delta\gamma-\sigma_y(\bar\varepsilon^{p}_{n} + \Delta\gamma)=0 \]

which can be solved efficiently for \(\Delta\gamma\) via Newton-Raphson iterations to tolerance \(\tilde{\phi}\leq tol\),

(284)\[ H:=\frac{\diff\sigma_y}{\diff\bar{\varepsilon}^p}\vert_{\bar{\varepsilon}^p_n+\Delta\gamma} \space \rightarrow \frac{\diff\tilde{\phi}}{\diff\Delta\gamma}:=-3\mu-H \space \rightarrow \Delta\gamma:=\Delta\gamma+\frac{\tilde{\phi}}{3\mu +H} \]

With \(\Delta\gamma\), we can update the deviatoric strain component in (263)

(285)\[ \textbf{e}_{d}^e= \textbf{e}^{e \space tr}_{d} - \Delta t \left( \frac{\Delta\gamma}{\Delta t}\bm{N}^{p \space tr} \right) = \textbf{e}^{e \space tr}_{d} - \Delta\gamma \left(\frac{3}{{2}}\frac{\bm{s^{tr}}}{q^{tr}} \right) = \left( 1 -3 \mu\frac{\Delta\gamma}{q^{tr}} \right) \textbf{e}^{e \space tr}_{d } \]

and the accumulated plastic strain

(286)\[ \bar\varepsilon^{p} =\bar{\varepsilon}^p_n+\Delta\gamma \]

To update \(\bm{C}^{p-1}_n\), on the other hand, we first update \(\bm{b}^e = \exp(2\textbf{e}^{e})\) and then compute

(287)\[ \bm{C}^{p-1} = \bm{F}^{-1} \bm{b}^e \bm{F}^{-T} \]

Expression for \(\mathsf{C}^i\)

The linearization of the Kirchhoff stress in (137) for von Mises plasticity is obtained by differentiating the expression for \(\diff\textbf{e}^{e}_{d }\) in (285) and making use of the relationships in (284)

(288)\[ \diff\bm{\tau} = \diff(2\mu (1-3\mu\frac{\Delta\gamma}{q^{tr}})\textbf{e}^{e}_{d } + \bulk\trace(\textbf{e}))= 2\mu(1-3\mu\frac{\Delta\gamma}{q^{tr}})\diff\textbf{e}^{e \space tr}_{d } - 6 \mu^2 \frac{\space \diff\tilde{\phi}}{q^{tr}}( \frac{1}{3\mu+H}- \frac{\Delta\gamma}{q^{tr}})\textbf{e}^{e \space tr}_{d } + \bulk\trace(\diff\textbf{e}) \]

where (283) gives \(\diff\tilde{\phi}=\diff q^{tr}\). Recalling \(q^{tr}=\sqrt{\frac{3}{2}s^{tr}: s^{tr}}\) and \(s^{tr}=2\mu\textbf{e}^{e \space tr}_{d }\), differentiating \(q^{tr}\) with respect to \(\textbf{e}^{e \space tr}_{d }\) gives,

(289)\[ \diff q^{tr} = \sqrt\frac{3}{2} \frac{\diff(s^{tr}: s^{tr})}{2\sqrt{s^{tr}: s^{tr}}} = \sqrt\frac{3}{2} \frac{2s^{tr} : \diff s^{tr}} {2\sqrt\frac{2}{3} {q^{tr}}}=\frac{3}{2} \frac{s^{tr} : \diff s^{tr}} {{q^{tr}}} \]

Command-line interface

To enable the elastoplastic Hencky model, use the model option -model plasticity-hencky-current or model plasticity-hencky-initial and set the material parameters listed in Hencky elastoplasticity model options. Any parameter without a default option is required.

Table 32 Hencky elastoplasticity model options

Option	Description	Default Value
-model plasticity-hencky-current | plasticity-hencky-initial	Required to enable the Hencky elastoplasticity model, | denotes an option for the user.
-E [real]	Young’s modulus, \(E > 0\)
-nu [real]	Poisson’s ratio, \(\nu \leq 0.5\).
-yield_stress [real]	Initial yield stress threshold, \(\sigma_0 > 0\)	Plasticity
-linear_hardening [real]	Isotropic linear hardening parameter, \(H >= 0\)	Plasticity
-saturation_stress [real]	Saturation stress for Voce hardening component, \(\sigma_\infty >= \sigma_0\)	Plasticity
-hardening_decay [real]	Hardening decay parameter for Voce hardening component, \(\beta >= 0\)	Plasticity

An example using the Hencky elastoplasticity model can be run via

$ ./bin/ratel-quasistatic -options_file examples/ymls/ex02-quasistatic-plasticity-hencky-initial-principal.yml