Return Mapping, Consistent Tangent, and Rate-Dependent Extensions
📽 Slides: Open presentation
Regime. This chapter assumes small-strain plasticity: \(\boldsymbol{\varepsilon} = \tfrac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T)\) with additive decomposition \(\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}^e + \boldsymbol{\varepsilon}^p\).
Over a time step \([t_n, t_{n+1}]\) with strain increment \(\Delta\boldsymbol{\varepsilon}\):
\[ \boldsymbol{\varepsilon}_{n+1} = \boldsymbol{\varepsilon}_n + \Delta\boldsymbol{\varepsilon}, \qquad \boldsymbol{\varepsilon}^e_{n+1} = \boldsymbol{\varepsilon}^e_n + \Delta\boldsymbol{\varepsilon} - \Delta\boldsymbol{\varepsilon}^p. \]
Elastic predictor — plastic corrector paradigm:
where \(f^\text{tr} := f(\boldsymbol{\sigma}^\text{tr}, q_n)\) is the trial yield function value.
For J2 plasticity the return mapping has a closed-form solution (radial return in deviatoric space):
Given \(\mathbf{s}^\text{tr}\) (trial deviatoric stress) and \(\bar{\varepsilon}^p_n\):
\[ \Delta\gamma = \frac{\sigma_\text{eq}^\text{tr} - \sigma_{y0} - H\bar{\varepsilon}^p_n}{2\mu + \tfrac{2}{3}H} \cdot \frac{1}{1} > 0 \]
Updated variables: \[ \mathbf{s}_{n+1} = \left(1 - \frac{3\mu\Delta\gamma}{\sigma_\text{eq}^\text{tr}}\right)\mathbf{s}^\text{tr}, \qquad \bar{\varepsilon}^p_{n+1} = \bar{\varepsilon}^p_n + \Delta\gamma. \]
def j2_return_mapping(s_tr, ep_bar_n, mu, sigma_y0, H):
"""Radial return for J2 plasticity with linear isotropic hardening."""
import numpy as np
seq_tr = np.sqrt(1.5 * np.tensordot(s_tr, s_tr))
f_tr = seq_tr - (sigma_y0 + H * ep_bar_n)
if f_tr <= 0.0: # elastic step
return s_tr, ep_bar_n
# plastic step
dgamma = f_tr / (2 * mu + 2/3 * H)
s_new = (1 - 3 * mu * dgamma / seq_tr) * s_tr
ep_bar_new = ep_bar_n + dgamma
return s_new, ep_bar_newThe continuum tangent \(\mathbb{C}^{ep}\) is not the right modulus for FEM — it loses quadratic convergence. The consistent (algorithmic) tangent \(\mathbb{C}^\text{alg} = \frac{d\boldsymbol{\sigma}_{n+1}}{d\boldsymbol{\varepsilon}_{n+1}}\) is required.
For J2 with radial return: \[ \mathbb{C}^\text{alg} = \kappa\,\mathbf{I}\otimes\mathbf{I} + 2\mu\theta_1\,\mathbb{P}_\text{dev} - 2\mu\theta_2\,\mathbf{n}\otimes\mathbf{n} \]
where \(\mathbf{n} = \mathbf{s}^\text{tr}/\|\mathbf{s}^\text{tr}\|\), \(\theta_1 = 1 - 3\mu\Delta\gamma/\sigma_\text{eq}^\text{tr}\), \(\theta_2 = 1/(1 + H/3\mu) - (1-\theta_1)\).
General return mapping for any yield surface (not only J2):
Minimise \(\tfrac{1}{2}\|\boldsymbol{\sigma} - \boldsymbol{\sigma}^\text{tr}\|_{\mathbb{C}^{e-1}}^2\) subject to \(f(\boldsymbol{\sigma}) = 0\).
Iterative Newton loop: \[ \begin{bmatrix}\boldsymbol{\sigma}\\\Delta\gamma\end{bmatrix}^{k+1} = \begin{bmatrix}\boldsymbol{\sigma}\\\Delta\gamma\end{bmatrix}^k - \mathbf{J}^{-1}\begin{bmatrix}\mathbf{r}_\sigma \\ f\end{bmatrix}^k \]
where \(\mathbf{J}\) is the Jacobian of the residual system.
Perzyna-type: allows stress to lie outside the yield surface — the yield function serves as an overstress measure:
\[ \dot{\bar{\varepsilon}}^p = \frac{1}{\eta}\langle f / \sigma_y \rangle^N, \qquad \langle x \rangle = \max(0, x). \]
Consistency viscoplasticity (Wang et al.): extends standard KKT conditions to rate-dependent case with a viscoplastic yield surface \(f^{vp}(\boldsymbol{\sigma}, \kappa, \dot{\kappa}) = 0\).
Update algorithm: implicit scheme, Newton solve for \(\Delta\gamma\).
Drucker-Prager yield function: \[ f = \alpha I_1 + \sqrt{J_2} - k \leq 0 \]
Return mapping: no closed form → cut-plane or linearized iterations.
Special cases: apex return (if trial stress maps to the apex cone).
Mohr-Coulomb in principal-stress space (Koiter’s multi-surface plasticity): \[ f_i = \sigma_j - \sigma_k - (\sigma_j + \sigma_k)\sin\phi - 2c\cos\phi \leq 0, \quad i,j,k = 1,2,3 \]
Six yield planes form the hexagonal pyramid.
Return mapping: depends on which face(s) are active — requires case distinction (edge/apex/face returns).
General approach for any yield surface using backward Euler (fully implicit):
Newton-Raphson system: \[ \boldsymbol{\sigma}_{n+1} = \mathbf{C}^e:(\boldsymbol{\varepsilon}_{n+1} - \boldsymbol{\varepsilon}^p_n - \Delta\gamma\mathbf{n}_{n+1}), \quad f(\boldsymbol{\sigma}_{n+1}, q_{n+1}) = 0, \]
where \(\Delta\gamma\) and state variables are unknowns. Iterations via: \[ \mathbf{J}\begin{bmatrix}\Delta\boldsymbol{\sigma}\\\Delta\Delta\gamma\end{bmatrix}^{(k)} = -\begin{bmatrix}\text{stress residual}\\f\end{bmatrix}^{(k)} \]
Jacobian combines elasticity, yield surface curvature, and hardening effects.
The consistent (algorithmic) tangent for general plasticity: \[ \mathbb{C}^\text{alg} = \mathbb{C}^e - \frac{\mathbb{C}^e:\mathbf{n} \otimes \mathbf{n}:\mathbb{C}^e}{h + \mathbf{n}:\mathbb{C}^e:\mathbf{n}} \]
where: - \(\mathbf{n} = \partial f/\partial\boldsymbol{\sigma}\) (flow direction)1 - \(h = -\partial f/\partial q \cdot dq/d\bar{\varepsilon}^p\) (hardening modulus, negative for softening)
Critical: Use \(\mathbb{C}^\text{alg}\) for FEM assembly. Continuum \(\mathbb{C}^{ep}\) loses quadratic convergence.
Drucker-Prager apex: if trial stress maps “inside” the cone apex, return directly to apex (constrain deviatoric flow).
Mohr-Coulomb edges: hexagonal pyramid has 12 edges. Returning to an edge means active constraints on two yield planes. Requires 2D projection in (\(I_1, \sqrt{J_2}\)) space.
Robust strategy: Check all candidate faces/edges; select the closest admissible surface.
Perzyna model: allows stress overstress beyond yield surface: \[ \dot{\bar{\varepsilon}}^p = \frac{1}{\eta}\left\langle \frac{f}{\sigma_y} \right\rangle^N, \quad \langle x \rangle = \max(0, x) \]
Overstress relaxes via viscous flow; material “creeps” above yield.
Penalty-type viscoplasticity: embedded directly in time-stepping (implicit method can use smaller steps for rate effects).
Key parameters: \(\eta\) (viscosity) and \(N\) (power law exponent).
For combined isotropic-kinematic hardening in incremental form: \[ \Delta\boldsymbol{\sigma} = \mathbb{C}^e:\Delta\boldsymbol{\varepsilon}^e = \mathbb{C}^e:(\Delta\boldsymbol{\varepsilon} - \Delta\gamma\mathbf{n}) \]
Yield condition at step \(n+1\): \[ f = \|\mathbf{s}_{n+1} - \boldsymbol{\alpha}_{n+1}\| - \sqrt{\frac{2}{3}}[\sigma_{y0} + (1-\beta)H\bar{\varepsilon}^p_{n+1}] = 0 \]
Back stress evolves: \[ \boldsymbol{\alpha}_{n+1} = \boldsymbol{\alpha}_n + \frac{2}{3}\beta H\Delta\gamma\mathbf{n} \]
Effective plastic strain: \[ \bar{\varepsilon}^p_{n+1} = \bar{\varepsilon}^p_n + \sqrt{\frac{2}{3}}\Delta\gamma \]
For a fully worked example combining linear kinematic (Prager) and power-law isotropic hardening — including the non-linear scalar residual \(R(\Delta\gamma) = 0\), its local Newton solve, and the consistent-tangent outline — see L08 Appendix — Worked Examples.
For J2 with linear hardening, the consistent tangent has explicit form: \[ \mathbb{C}^\text{alg} = \kappa(\mathbf{I}\otimes\mathbf{I}) + 2\mu\left[ \theta_1\mathbb{P}_\text{dev} - \theta_2\mathbf{n}\otimes\mathbf{n} \right] \]
where: - \(\mathbb{P}_\text{dev} = \mathbb{I} - \frac{1}{3}(\mathbf{I}\otimes\mathbf{I})\) is deviatoric projector - \(\theta_1 = 1 - 3\mu\Delta\gamma/\sigma_\text{eq}^\text{tr}\) captures deviatoric scaling - \(\theta_2 = 1 - \theta_1 + H/(3\mu)\) balances elastic vs. plastic stiffness
Ensures quadratic N-R convergence in FEM.
Voigt convention. Throughout this chapter, the 4th-order elasticity tensor \(\mathbb{C}^e\) is stored as a 6×6 matrix using Voigt ordering: index pairs (i,j) map to a single index as [1,1]→1, [2,2]→2, [3,3]→3, [2,3]=[3,2]→4, [1,3]=[3,1]→5, [1,2]=[2,1]→6. Stress and strain vectors are stored with the same ordering. Note that strain components 4–6 are engineering shears (\(\gamma_{ij} = 2\varepsilon_{ij}\) for \(i \neq j\)), so Voigt is not an isometry — use Mandel ordering (scaled by \(\sqrt{2}\)) if isometric inner products are required.
State variables to store at each Gauss point: - Cauchy stress \(\boldsymbol{\sigma}\) - Effective plastic strain \(\bar{\varepsilon}^p\) (or cumulative plastic multiplier \(\sum\Delta\gamma\)) - Back stress \(\boldsymbol{\alpha}\) (if kinematic hardening) - Accumulated damage \(D\) (if damage is coupled)
Algorithm summary: 1. Compute trial stress from elastic predictor 2. Check yield; if elastic, update stresses only 3. If plastic: implicit return mapping (Newton loop) → updated \(\boldsymbol{\sigma}, \bar{\varepsilon}^p, \boldsymbol{\alpha}, q\) 4. Compute algorithmic tangent for FEM assembly 5. Assemble element stiffness and internal force; proceed to next global iteration
Common pitfalls: - Using continuum tangent instead of algorithmic → poor convergence - Stress drift from yield surface → use consistent/accurate return mapping - Integration error accumulation → monitor time step size
📖 Worked example: L08 Appendix — Worked Examples carries out the full RRM derivation for the combined linear kinematic (Prager) + power-law isotropic hardening case, including the non-linear scalar residual \(R(\Delta\gamma) = 0\) and its local Newton solve.
Constitutive Models — Technion | Website