10 L08 — Plasticity: Algorithms

Return Mapping, Consistent Tangent, and Rate-Dependent Extensions

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\).

10.1 Incremental Elasto-Plasticity

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:

Freeze plastic strain: assume \(\Delta\boldsymbol{\varepsilon}^p = \mathbf{0}\).
Compute trial stress: \(\boldsymbol{\sigma}^\text{tr} = \boldsymbol{\sigma}_n + \mathbb{C}^e:\Delta\boldsymbol{\varepsilon}\).

where \(f^\text{tr} := f(\boldsymbol{\sigma}^\text{tr}, q_n)\) is the trial yield function value.

Check yield: if \(f^\text{tr} \leq 0\) → elastic step, accept trial.
If \(f^\text{tr} > 0\) → plastic step, correct by return mapping.

10.2 Radial Return Algorithm (J2)

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_new

10.3 Consistent (Algorithmic) Tangent

The 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)\).

10.4 Cut-Plane (Closest-Point Projection) Algorithm

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.

10.5 Rate-Dependent Plasticity (Viscoplasticity)

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\).

10.6 Pressure-Dependent Plasticity: Drucker-Prager

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).

10.7 Pressure-Dependent Plasticity: Mohr-Coulomb

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).

10.8 Multi-Dimensional Return Mapping: Implicit Integration

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.

10.9 Consistent Tangent: Multi-Dimensional

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)¹ - \(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.

10.10 Apex and Edge Returns

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.

10.11 Viscoplasticity: Rate-Dependent Plasticity

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).

10.12 Incremental Form: General Hardening

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.

10.13 Consistent Tangent: J2 Plasticity Explicit Form

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.

10.14 Implementation Notes

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.

For associated flow rules, the flow direction equals the yield-normal: \(\mathbf{n} = \partial f/\partial\boldsymbol{\sigma}\). For non-associated flow rules, a separate plastic potential \(g(\boldsymbol{\sigma})\) defines the direction \(\mathbf{n} = \partial g/\partial\boldsymbol{\sigma}\) (see L07 for the general setup). Unless stated otherwise, this chapter assumes associated flow (\(g = f\)).↩︎