Isotropic Elastoplasticity in Logarithmic Strains Space¶
Note
The formulation of rate-independent isotropic elastoplasticity at large strains and derived algorithms are based on the pioneer work by [EB90] and [WA90].
General Framework¶
In elastoplasticity, the stress state depends on internal variable(s) associated with dissipative deformation 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 stresses implicitly via plastic corrector (or return-mapping) algorithms. To ensure quadratic convergence of the Newton’s method, the tangent modulus must be obtained by linearization of the algorithmic expression for the stress update procedure [JCRL85], which depends 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 hyperelasticity in logarithmic strain space, the free energy potential \(\psi\) can be expressed as a function of the (Eulerian) logarithmic strain tensor \(\textbf{e}^e\) using the same functional format of the small-strains counterpart
where \(\textbf{D}\) the fourth order elasticity tensor of small-strains (linear) elasticity. 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\)
The (symmetric) Kirchhoff stress tensor \(\bm{\tau}\equiv2\frac{\partial\psi}{\partial\bm{b}^e}\bm{b}^e\) can thus be expressed as
Assuming isotropic elastic response and taking advantage of certain properties of derivatives of isotropic tensor functions, the following identity can be derived [dSNEARJ98]
Substituting (187) into (186) finally gives
which similarly retains the functional format of linear elasticity.
In elastoplasticity, the free energy potential also depends on a set of \(k\) internal variables \(\bm{\alpha}=\{\alpha_1, ..., \alpha_k\}\) such that \(\psi=\psi(\textbf{e}^e, \bm{\alpha})\). Differentiation with respect to these variables defines so-called thermodynamic forces \(A_i\space\) that drive plastic deformation and are work conjugates of the internal variables \(\alpha_i\space\),
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]
Plastic deformation is assumed isochoric (i.e., \(\det\bm{F}^p=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. With (190), the spatial quantity describing the evolution of \(\bm{F}\), \(\bm{l}\equiv\bm{\dot{F}}\bm{F}^{-1}\), can be additively decomposed as
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.
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] and 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]. For additive plasticity, the Lagrangian strain \(\bm{E}\) is broken into elastic and plastic parts:
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:
From Miehe et al. [MAL02], \(\bm{E}^p\) is defined as follows:
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
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
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.
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]:
Where \(\dot\gamma\) is the plastic multiplier calculated through optimality conditions
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
Note that assumption of zero plastic spin in (199) 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 (199) 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 (190), (191) and (199), \(\bm{F}^p_{n+1}\) at time \((t+\Delta{t})\) can be calculated as
Analogously, \(\bm{F}^e_{n+1}\) can be expressed from (200) 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 the isotropic property of exponential function, it gives
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 tr}_{n+1}=\bm{F}_\Delta\bm{F}_n^e\) and (201) rewritten as:
The update of kinematic quantities can be further simplified with manipulation via the stretch tensor \(\bm{v}^e_{n+1}\) and (185) to obtain
(203) is completed with the internal variables update
and constrains
This system of equations for \(\Delta\gamma\) has the same functional format as its small-strains counterpart.
Rate-independent Isotropic von Mises Model¶
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
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:
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:
with (scalar) \(q=\sqrt{3J_2({\bm{s}})}= \sqrt{\frac{3}{2}}||\bm{s}||\) termed the von Mises (effective or equivalent) stress.
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 (\(\bm{N}^p\)) is thus normal to the yield surface. For von Mises plasticity, we can then write
It also follows from (200) that
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 kinematic state variables \(\bar\varepsilon^{p}_n\) and \(\bm{C}^{p-1}_n\) is given below.
Note
Choosing \(\bm{C}^{p-1}_n\) as additional state variable alongside \(\bar\varepsilon^{p}_n\) minimizes the stored arrays, as \(\bm{C}^{p-1}_n\) has length 6 in Voigt notation. However, other kinematic state variables could be considered. For instance, one could use \(\bm{F}^{p-1}_n\) and update via (210). This approach would nevertheless increase storage and introduce a significant amount of additional computation, e.g. calculating \(\textbf{v}^{e}_{n+1}\) needed in \(\bm{R}^{e \space tr}_{n+1} = \bm{R}^e_{n+1} = (\textbf{v}^{e}_{n+1})^{-1} \bm{F}^{e}_{n+1}\).
Elastic Predictor¶
The trial effective stress \(q^{tr}\) can be calculated at \(t+\Delta t\) as:
Admissibility can thus be verified after the initial guess \(\bar\varepsilon^{p}_{n+1} =\bar{\varepsilon}^p_n\) :
Note
The isotropic tensor function for \(\textbf{e}^{e \space tr}_{\space n+1}\) in (211)can be conveniently computed via the spectral decomposition of \(\textbf{b}^{e \space tr}\) as
with eigenvalues \(\lambda^b_i=\lambda^2_i\) and eigenvectors \(\hat{\bm{N}}_i\).
Return Mapping¶
In the von Mises model, the system of equations composed of (203), (204) and (205) can be reduced to a single non-linear equation
which can be solved efficiently using the Newton method with tolerance \(\tilde{\bm{\phi}}\leq tol\),
Note
In linear hardening, \(H\) is constant so a closed form expression for \(\Delta\gamma\) is obtained
Update¶
With \(\Delta\gamma\), deviatoric stresses can be updated as
and used to compute \(\bm{\tau}_{n+1}\). The update of state variables is also straightforward. The current accumulated plastic strain is computed as
To update \(\bm{C}^{p-1}_n\), on the other hand, we first compute \(\bm{b}^e_{n+1} = \exp(2\textbf{e}^{e}_{n+1})\) and then
Note
Ratel also implements a principal stress space-based formulation whereby the stress update is directly carried out on the eigenvalues of \(\textbf{s}\)
and \(\textbf{s}\) retrieved as
Consistent Tangent¶
The (consistent) linearization of \(\bm \tau\) with respect to \(\bm F\) in (121) is derived from its algorithmic expression \(\tilde\tau=\tilde\tau(\textbf{e}^{e \space tr}_{n+1}(\bm{B}^{e \space tr}_{n+1}(F_n^e,F_{n+1})), \bm{\alpha}_n)\) and is conveniently expressed as the double contraction of three fourth-order tensors [dSNEARJ98]
with
\(\textbf{D}^{ep}\) is symmetric for associative plasticity and the only quantity in (224) that depends on the specific plasticity model. Also note that \(\textbf{L}\) is the derivative of tensor-valued function of symmetric tensor, for details see [C98] [BJ16],
Expression for \(\textbf{B}\)¶
The stored \(\mathbf{C_n^{p-1}}\) can be retrieved to compute the change of \(\bm{b}^{e \space tr}_{n+1}\) with \(\bm{F}\) as,
Expression for \(\textbf{L}\)¶
To calculate the change \(\diff\textbf{e}^{e \space tr}_{n+1}\) due to a perturbation \(\diff\bm{b}^{e \space tr}_{n+1}\), it is convenient to express the former in term of eigenvalues \(\lambda^2_i\) and eigenbasis \(\hat{\bm{N}}_i\) of \(\bm{b}^{e \space tr}_{n+1}\) as
Note
Helper Qfunctions are implemented in Ratel to perform the spectral decomposition of a \(3\times3\) symmetric matrix and compute the increments \(\diff\lambda^b_i\) and \(\diff\hat{\bm{N}}_i\) {eq}``, as described for mixed elasticity (see derivation of terms in (164)).
Expression for \(\textbf D^{ep}\) : von Mises case¶
The elastoplastic consistent tangent tensor \(\textbf{D}^{ep}\) in (223) finally links \(\diff\tilde{\bm \tau}\) to the change in logarithmic elastic strains \(\diff\textbf{e}^{e}_{\space n+1} = \diff\textbf{e}^{e}_{v \space n+1} + \diff\textbf{e}^{e}_{d \space n+1}\). For the von Mises model, only the deviatoric part is updated. Differentiating the expression for \(\diff\textbf{e}^{e}_{d \space n+1}\) in (218)and making use of the relationships in (216), we get:
where (215)gives \(\diff\tilde{\phi}=\diff q^{tr}_{n+1}\). Recalling \(q^{tr}_{n+1}=\sqrt{\frac{3}{2}s^{tr}_{n+1}: s^{tr}_{n+1}}\) and \(s^{tr}_{n+1}=2\mu\textbf{e}^{e \space tr}_{d \space n+1}\), differentiating \(q^{tr}_{n+1}\) with respect to \(\textbf{e}^{e \space tr}_{d \space n+1}\) gives,
Note
In the principal axes implementation, the linearization of \(\bm{\tau}\) in (227) is expressed in terms of its eigenvalues as
where