Isotropic Elastoplasticity in Logarithmic Strains Space

Note

[PDME92]’s formulation of rate-independent isotropic elastoplasticity at large strains and derived algorithms are reported below. Ratel implementation is in progress.

Linearized Virtual Work: Spatial Tangent Modulus

Being \(\chi(\Omega)\) the region occupied by the deformed body \(\Omega\) and \(\rho\bm{g}\) and \(\bar{\bm{t}}\) prescribed body forces and surface traction, respectively, the spatial (or current) version of the linearized virtual work equation at \(\bm{u}\) for the quasistatic case (i.e., neglecting inertia terms) reads:

(184)\[ \int_{\chi(\Omega)}{ \textbf{a} : \nabla_x \bm{u} : \nabla_x \bm{v}} \ + \int_{\chi(\Omega)}{(\bm{\sigma} : \nabla_x \bm{v} - \rho\bm{g} \cdot \bm{v})} \ - \int_{\chi(\partial \Omega)}{\bar{\bm{t}} \cdot \bm{v}} \ = 0, \quad \forall \bm{v} \in \mathcal{V}. \]

where \(\mathcal{V}\) is the space of admissible displacements \(\bm{v}\) and gradients \(\nabla_x\) are calculated in the spatial configuration. In elastoplasticity, the Cauchy stress \(\bm{\sigma}\) also depends on internal variable(s) associated with dissipative mechanisms. The resulting strain path dependency is generally dealt with by introducing sufficiently small time steps for the problem to be considered path independent within one time increment \([t, t_{n+1}]\) and updating \(\bm{\sigma}_{n+1}\) implicitly via plastic corrector (or return-mapping) algorithms. To ensure quadratic convergence, the tangent modulus \(\textbf{a}\) must be obtained by linearization of the algorithmic expression for the stress update procedure [JCRL85].

General Framework

Expressions for \(\bm{\sigma}\) and \(\textbf{a}\) in (184) depend on the specific constitutive model. Following the formalism of thermodynamics with internal variables, these laws are derived from postulated free energy potential \(\psi\), dissipation potential \(\Xi\), and yield function \(\bm{\phi}\) [GFA10].

Free Energy Potential

In elastoplasticity in logarithmic strains space, the free energy potential \(\psi=\psi(\textbf{e}^e, \bm{\alpha})\) is a function of the (Eulerian) logarithmic strain tensor \(\textbf{e}^e\) and a set of \(k\) internal variables \(\bm{\alpha}=\{\alpha_1, ..., \alpha_k\}\). As in (11), \(\textbf{e}^e\) can be equally defined via the (elastic) left stretch tensor \(\textbf{v}^e\) or the (elastic) left Cauchy-Green strain tensor \(\bm{b}^e=(\textbf{v}^e)^2=\bm{F}^e(\bm{F}^e)^T\)

(185)\[ \textbf{e}^e \equiv \log{\textbf{v}^e}=\frac{1}{2}\log\bm{b}^e \]

The (symmetric) Kirchhoff stress tensor \(\bm{\tau}\equiv2\rho_0\frac{\partial\psi}{\partial\bm{b}^e}\bm{b}^e= J\bm{\sigma}\) can thus be expressed as

(186)\[ \bm{\tau}=2\rho_0\frac{\partial\psi}{\partial\bm{b}^e}\bm{b}^e=\rho_0\frac{\partial\psi}{\partial\textbf{e}^e}:\frac{\partial\log\bm{b}^e}{\partial\bm{b}^e}\bm{b}^e \]

where \(\rho_0\) is the mass density in the reference configuration. Assuming isotropic elastic response and taking advantage of certain properties of derivatives of isotropic tensor functions, the following identity can be derived [dSNEARJ98]

(187)\[ \frac{\partial\psi}{\partial\textbf{e}^e}:\frac{\partial\log\bm{b}^e}{\partial\bm{b}^e}\bm{b}^e=\frac{\partial\psi}{\partial\textbf{e}^e} \]

Substituting (187)into (186) finally gives

\[ \bm{\tau}=\rho_0\frac{\partial\psi}{\partial\textbf{e}^e} \]

The thermodynamic forces \(A_i\space\), work conjugates of the internal variables \(\alpha_i\space\) , are similarly defined as

\[ A_i=\rho_0\frac{\partial\psi}{\partial\alpha_i} \]

Multiplicative decomposition of deformation gradient

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 a plastic \(\bm{F}^p\) and elastic part \(\bm{F}^e\) such that [EH69]

(188)\[ \bm{F}=\bm{F}^e\bm{F}^p \]

Plastic deformation is assumed isochoric (i.e., \(\det\bm{F}^p=1\)) and mapping material points to an unstressed intermediate configuration, as shown in Multiplicative decomposition of the deformation gradient tensor.. 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. With (188), the spatial quantity describing the evolution of \(\bm{F}\), \(\bm{l}\equiv\bm{\dot{F}}\bm{F}^{-1}\), can be additively decomposed as

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

where \(\bm{l}^e\equiv\dot{\bm{F}^e}(\bm{F}^e)^{-1}\) and \(\bm{L}^p\equiv\dot{\bm{F}^p}(\bm{F}^p)^{-1}\), which lives in the intermediate configuration. To formulate the so-called (plastic) flow rule, it is convenient to further split \(\bm{L}^p\) into its symmetric \(\bm{D}^p\equiv\text{sym}(\bm{L}^p)\) and skew \(\bm{W}^p\equiv\text{skew}(\bm{L}^p)\) parts. So defined, the plastic stretching tensor \(\bm{D}^p\) measures the instantaneous rate of plastic stretching along the orthonormal eigenvectors in the intermediate configuration. On the other hand, \(\bm{W}^p\) represents the triad’s rigid spinning.

Fig. 5 Multiplicative decomposition of the deformation gradient tensor.

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]. For additive plasticity, the Lagrangian strain \(\bm{E}\) is broken into elastic and plastic parts:

(190)\[ \bm{E} = \bm{E}^e + \bm{E}^p \]

where \(\bm{E}\) is the strain (for the linear case we use linear strain); in logarithmic space, using a geometric preprocessor, \(\bm{E}^e\) becomes:

(191)\[ \bm{E} = \frac{1}{2} \log \left( \bm{C} \right). \]

From Miehe et al. [MAL02], \(\bm{E}^p\) is defined as follows:

(192)\[ \begin{aligned} \bm{E}^p &= \frac{1}{2} \log \left( \bm{C}^p \right), \\ \bm{C}^p &= (\bm{F}^p)^T\bm{F}^p \end{aligned} \]

where \(\bm{C}^p\) is the plastic part of the Cauchy-Green tensor (also refered to as the plastic metric).

The Second Piola-Kirchhoff stress, \(\bm{S}\), now becomes

(193)\[ \bm{S} \left( \bm{E} \right) = \frac{\partial \psi}{\partial \bm{E}^e} \left( \bm{E}^e, \bm{\alpha} \right). \]

\(\mathcal{I}\) contains the internal (history) variables of the plasticity model. For a minimal plasticity model \(\mathcal{I} = {\bar{\epsilon}^p}\), where \(\bar{\epsilon}^p\) is the accumulated plastic strain; for a plasticity model with hardening \(\mathcal{I} = {\bar{\epsilon}^p, \bm{\textbf{A}}, \bm{\alpha}}\). At the convergence of each solve iteration \(\mathcal{I}\) is updated with the most recent values.

Yield Function

The yield function \(\phi (\bm{\tau}, \bm{\textbf{A}}, \bm{\alpha})\) is used to define the set of admissible (i.e., physically possible given the current value of the internal variables) stress states \(\mathscr{S}\) as

\[ \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 as illustrated for the von Mises case in principal deviatoric stress space (the \(\pi\)-plane is a yield surface projection normal to the hydrostatic axis) Yield surface and representation of the return-mapping scheme for von Mises plasticity.. 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 as described below.

Fig. 6 Yield surface and representation of the return-mapping scheme for von Mises plasticity.

Dissipation Potential

The non-negative and convex dissipation potential \(\Xi(A, \bm\alpha)\) is used to state constitutive laws for the evolution of internal variables and the flow rule. To ensure thermodynamic principles are satisfied, a normality rule is assumed [BN75]:

\[ \dot{\bm{\alpha}}=-\dot\gamma\frac{\partial\Xi}{\partial\bm{\textbf{A}}} \]

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

\[ \dot\gamma\geq0 \qquad \bm{\phi}\leq0 \qquad \dot\gamma\bm{\phi}=0 \]

The plastic multiplier also appears in the flow rule, which is conveniently stated as constitutive laws for \(\bm{D}^p\) (rotated to the spatial configuration via \(\bm{R}^e\)) and \(\bm{W}^p\) as

(194)\[ \bm{R}^e\bm{D}^p(\bm{R}^e)^T=\dot\gamma\frac{\partial\Xi}{\partial\bm{\tau}} \medspace, \qquad \bm{W}^p=\bm{0} \]

Note that assumption of zero plastic spin in (194) is consistent with that of plastic isotropy. Plastic isotropy dictates coaxiality, i.e. sharing of principal directions, between the flow vector \({\partial\Xi}/{\partial\bm{\tau}}\) in (194) and \(\bm{\tau}\). Assumption of isochoric plastic flow further implies that \({\partial\Xi}/{\partial\bm{\tau}}\) is a traceless tensor [GFA10]. These properties help to greatly simplify the update of kinematic quantities as described next.

Exponential Map Integration

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]. Using backward integration and via substitution of (188), (189) and (194), \(\bm{F}^p_{n+1}\) at time \((t+\Delta{t})\) can be calculated as

(195)\[ \bm{F}^p_{n+1}=\exp(\Delta\gamma\bm{R}_{n+1}^{e \space T}\frac{\partial\bm{\Xi}}{\partial\bm{\tau}}\vert_{n+1}\bm{R}_{n+1}^e)\bm{F}_n^p \]

Analogously, \(\bm{F}^e_{n+1}\) can be expressed from (195) by introducing the deformation gradient increment \(\bm{F}_\Delta=\bm{F}_{n+1}\bm{F}_n^{-1}\), which expresses the “difference” between deformations at different times. Taking also advantage of isotropic property of exponential function, it gives

(196)\[ \bm{F}^e_{n+1}=\bm{F}_\Delta\bm{F}_n^e \bm{R}_{n+1}^{e \space T}\exp(-\Delta\gamma\frac{\partial\bm{\Xi}}{\partial\bm{\tau}}\vert_{n+1})\bm{R}_{n+1}^e \]

Tentatively considering \(\bm{F}_\Delta\) as purely elastic would generally lead to inadmissible stress states. A trial elastic deformation increment is thus introduced as \(\bm{F}^{e \space \text{trial}}_{n+1}=\bm{F}_\Delta\bm{F}_n^e\) and (196) rewritten as:

\[ \bm{F}^e_{n+1}=\bm{F}^{e \space \text{trial}}_{n+1} \bm{R}_{n+1}^{e \space T}\exp(-\Delta\gamma\frac{\partial\bm{\Xi}}{\partial\bm{\tau}}\vert_{n+1})\bm{R}_{n+1}^e \]

The update of kinematic quantities can be further simplified with manipulation via the stretch tensor \(\bm{v}^e_{n+1}\) and (185) to obtain

(197)\[ \textbf{e}^e_{n+1} =\textbf{e}^{e \space \text{trial}}_{n+1} - \Delta\gamma\frac{\partial\bm{\Xi}}{\partial\bm{\tau}}\vert_{n+1} \]

(197)is completed with the internal variables update

(198)\[ \bm{\alpha}_{n+1} =\bm{\alpha}_n - \Delta\gamma\frac{\partial\bm{\Xi}}{\partial\bm{\textbf{A}}}\vert_{n+1} \]

and constrains

(199)\[ \Delta\gamma\geq0 \medspace, \qquad \bm{\phi}(\bm{\tau}_{n+1} ,\bm{\textbf{A}}_{n+1}) \leq 0 \medspace, \qquad \Delta\gamma \bm{\phi}(\bm{\tau}_{n+1}, \bm{\textbf{A}}_{n+1})=0 \]

This system of equations for \(\Delta\gamma\) has the same functional format as its small-strains counterpart.

Consistent Tangent

The forth order tensor \(\textbf{a}\) in (184) of major symmetry \(\text{a}_{ijkl}=\text{a}_{klij}\) is obtained by differentiating the algorithmic expression for the Kirchhoff stress \(\tilde\tau=\tilde\tau(\textbf{e}^{e \space \text{trial}}_{n+1}, \bm{\alpha}_n)\) in

(200)\[ \text{a}_{ijkl}=\frac{1}{J}\frac{\partial\tilde{\tau}_{ij}}{\partial{F}_{kq}}F_{lq}-\sigma_{il}\delta_{jk} \]

where \(\delta_{ij}\) is the Kronecker delta (\(\delta_{ij}=1 \text{ if } i=j, \text{ else } \delta_{ij}=0)\) [dSNEARJ98]. In logarithmic strains space, (200) can be expressed as a function of three forth order tensors (See [dSNEARJ98] for derivation)

\[ \text{a}_{ijkl}=\frac{1}{2J}(\textbf{D}:\textbf{L}:\textbf{B})_{ijkl}-\sigma_{il}\delta_{jk} \]

with

(201)\[ \textbf{D}=\frac{\partial\tilde{\tau}}{\partial\textbf{e}^{e \space \text{trial}}_{n+1}} \quad, \quad \textbf{L}=\frac{\partial\log(\bm{B}^{e \space \text{trial}}_{n+1})}{\partial\bm{B}^{e \space \text{trial}}_{n+1}} \quad, \space \text{and} \quad \textbf{B}=\delta_{ik}(\bm{b}^{e \space \text{trial}}_{n+1})_{jl}+\delta_{jk}(\bm{b}^{e \space \text{trial}}_{n+1})_{il} \]

\(\textbf{D}\) is symmetric for associative plasticity and the only quantity in (201) that is affected by the specific material model. Its expression depends on whether plastic loading or elastic loading/unloading occur:

(202)\[ \textbf{D}= \begin{cases} \textbf{D}^e \qquad \text{for elastic loading/unloading} \\ \textbf{D}^{ep} \qquad \text{for plastic loading} \end{cases} \]

Whereas the expression for \(\textbf{D}^{ep}\) depends on the particular choice of plasticity model (e.g. von Mises plasticity), \(\textbf{D}^e\) is the more familiar elasticity tensor.

\[ \textbf{D}^{e} = 2\mu\textbf{I}_d +K \bm{I} \otimes \bm{I} \]

where \(K\) is the bulk modulus and \(\textbf{I}_d=\textbf{I}-\frac{1}{3}\bm{I}\otimes\bm{I}\) is the forth order projection tensor with symmetric identity forth order tensor \(\textbf{I}_{ijkl}= \frac{1}{2}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\). This result follows from the definition of Hencky strain energy \(\psi^e \equiv \frac{1}{2} \textbf{e}^e:\textbf{D}^e:\textbf{e}^e\) describing the hyperelastic response.

Note

Several methods have been proposed to calculate the derivative of tensor-valued function of symmetric tensor \(\textbf{L}\) in a computationally efficient way [C98] [BJ16].

Rate-independent Isotropic von Mises Model

The principle of maximum plastic dissipation leads to associative plasticity [BN75], in which the dissipation potential coincides with the yield function and the flow vector is thus normal to the yield surface, as shown in Yield surface and representation of the return-mapping scheme for von Mises plasticity.. 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 conveniently expressed in terms of the second invariant

\[ J_2=\frac{1}{2}\bm{\hat\sigma}:\bm{\hat\sigma} \]

of the (symmetric) deviatoric stress tensor \(\bm{\hat\sigma}=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 [needs double check]:

\[ \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 can be conveniently expressed in terms of uniaxial stress as:

\[ \bm{\phi}=q-\sigma_y \]

with (scalar) \(q=\sqrt{3J_2({\bm{s}})}\) termed the von Mises (effective or equivalent) stress.

Stress Update

In comparison with the returning mapping algorithm for the small-strain formulation, the one for logarithmic strains only requires 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 is given below.

Elastic Predictor

Having calculated \(\bm{F}_\Delta\) and recovered the elastic strain tensor obtained in the previous time step, \(\textbf{e}^e_n\), the trial effective stress \(q^{\text{\text{trial}}}\) can be calculated at \(t+\Delta t\) as:

\[ \bm{b}^e_n=\exp(2\textbf{e}^e_n) \space \rightarrow \space \bm{b}^{e \space \text{trial}}_{n+1}=\bm{F}_\Delta \bm{b}^e_n \bm{F}_\Delta^T \space \rightarrow \space \bm{\varepsilon}^{e \space \text{trial}}_{n+1}=0.5\log(\bm{b}^{e \space \text{trial}}_{n+1}) \space \rightarrow \]

\[ \space \rightarrow \space \textbf{e}^{e \space \text{trial}}_{d \space n+1}= \space \rightarrow \space \bm{\hat\sigma}^{\text{\text{trial}}}_{n+1} = 2\mu \textbf{e}^{e \space \text{trial}}_{d \space n+1}\space \rightarrow \space q^{\text{\text{trial}}}_{n+1}= \sqrt{3J_2({\bm{\hat\sigma}^{\text{\text{trial}}}_{n+1}})} \]

Admissibility can thus be verified after the initial guess \(\bar\varepsilon^{p \space \text{trial}}_{n+1} =\bar{\varepsilon}^p_n\) :

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

Return Mapping

In the von Mises model , the system of equations composed of (197), (198) and (199) can be reduced to a single non-linear equation

\[ \tilde{\bm{\phi}}=q^{\text{\text{trial}}}_{n+1}-3\mu\Delta\gamma-\sigma_y=0 \]

which can be solved efficiently using the Raphson-Newton method with tolerance \(\tilde{\bm{\phi}}\leq tol\) ,

\[ H:=\frac{d\sigma_y}{d\bar{\varepsilon}^p_n}\vert_{\bar{\varepsilon}^p_n+\Delta\gamma} \space \rightarrow d:=-3\mu-H \space \rightarrow \Delta\gamma:=\Delta\gamma-\frac{\tilde{\bm{\phi}}}{d} \]

Note

Linear Hardening slope is constant and this can be closed form that non-linear as piecewise linear hardening at time step values

\[ \Delta\gamma=\frac{q^{\text{\text{trial}}}}{3\mu+H} \]

Update

Once \(\Delta\gamma\) is calculated, the deviatoric stress tensor is corrected

\[ \bm{\hat\sigma}_{n+1} =(1-3\mu\frac{\Delta\gamma}{q^{\text{\text{trial}}}_{n+1}})\bm{\hat\sigma}^{\text{\text{trial}}}_{n+1} \]

and finally \(\bm{\sigma}\) in (184) updated with \(p_{n+1}=p_{n}\) as

\[ \bm{\tau}_{n+1} =\bm{s}_{n+1}+p_{n+1}\bm{I} \space \rightarrow \space \bm{\sigma}_{n+1} =J^{-1}\bm{\tau}_{n+1} \]

The update of the accumulated plastic strain is also straightforward:

\[ \bar\varepsilon^{p}_{n+1} =\bar{\varepsilon}^p_n+\Delta\gamma \]

Consistent Tangent for von Mises Model

The elastoplastic consistent tangent tensor \(\textbf{D}^{ep}\) in (202) for the von Mises model has the compact form (See [dSNEARJ98] for derivation)

\[ \textbf{D}^{ep}= \textbf{D}^e - 6\mu^2 \frac{\Delta\gamma}{q^{\text{\text{trial}}}_{n+1}} \textbf{I}_d+ 6\mu^2 \left(\frac{\Delta\gamma}{q^{\text{\text{trial}}}_{n+1}}-\frac{1}{3\mu+H}\right)\bar{N}_{n+1}\otimes\bar{N}_{n+1} \]

with unit flow vector \(\bar{N}_{n+1} = \textbf{e}^{e \space \text{trial}}_{d \space n+1}/{\Vert\textbf{e}^{e \space \text{trial}}_{d \space n+1}\Vert }\).