22 L08 Appendix — Worked Examples
Combined Hardening: RRM with Linear Kinematic + Power-Law Isotropic
22.1 Purpose
L08 presents the radial return method (RRM) for J2 plasticity with linear isotropic hardening, where the plastic multiplier update has a closed form. This appendix carries out the full derivation for a richer combined-hardening model — linear kinematic (Prager) plus power-law isotropic — where the update becomes a non-linear scalar equation and a local Newton solve is required. Every step is written out so the derivation can be followed without prior knowledge beyond L08. Brief recap of standard material is included where it aids the flow; see L08 for the underlying theory.
The running conventions are the same as L08: \(\boldsymbol{\varepsilon}\) for total strain, \(\boldsymbol{\varepsilon}^p\) for plastic strain, \(\bar{\varepsilon}^p\) for equivalent plastic strain, \(\boldsymbol{\sigma}\) for Cauchy stress, \(\mathbf{s} = \text{dev}(\boldsymbol{\sigma})\) for the deviatoric stress, \(\boldsymbol{\alpha}\) for the backstress, \(\Delta\gamma\) for the discrete plastic multiplier, \(\mathbb{C}^e\) for the fourth-order elasticity tensor, and \(\mu\), \(\kappa\) for the shear and bulk moduli.
22.2 Model Definition
Yield function (combined hardening, J2): \[ f(\boldsymbol{\sigma}, \boldsymbol{\alpha}, \bar{\varepsilon}^p) = \|\mathbf{s} - \boldsymbol{\alpha}\| - \sqrt{\tfrac{2}{3}}\,\sigma_y(\bar{\varepsilon}^p) \;\le\; 0, \]
where
- \(\mathbf{s} = \text{dev}(\boldsymbol{\sigma})\) is the deviatoric Cauchy stress;
- \(\boldsymbol{\alpha}\) is the backstress tensor, deviatoric in J2 plasticity, marking the centre of the yield surface in deviatoric stress space;
- \(\bar{\varepsilon}^p\) is the scalar equivalent plastic strain (isotropic hardening variable);
- \(\|\cdot\|\) denotes the Frobenius norm, \(\|\mathbf{X}\| = \sqrt{\mathbf{X}:\mathbf{X}}\);
- \(\sigma_y(\bar{\varepsilon}^p)\) is the current yield stress.
Isotropic hardening — power law: \[ \sigma_y(\bar{\varepsilon}^p) \;=\; \sigma_{y0} \;+\; H_\text{iso}\,(\bar{\varepsilon}^p)^{m}, \]
with initial yield \(\sigma_{y0}\), isotropic hardening coefficient \(H_\text{iso}\), and exponent \(m \in [0, 1]\).
Kinematic hardening — linear Prager rule: \[ \dot{\boldsymbol{\alpha}} \;=\; H_\text{kin}\,\dot{\boldsymbol{\varepsilon}}^p \quad\Longrightarrow\quad \Delta\boldsymbol{\alpha} \;=\; H_\text{kin}\,\Delta\boldsymbol{\varepsilon}^p, \]
with constant scalar kinematic hardening modulus \(H_\text{kin}\). In J2 plasticity \(\dot{\boldsymbol{\varepsilon}}^p\) is deviatoric, so \(\dot{\boldsymbol{\alpha}}\) is deviatoric; if \(\boldsymbol{\alpha}_0\) is deviatoric, \(\boldsymbol{\alpha}\) remains deviatoric for all time.
Effective (shifted) stress. Introduce \[ \boldsymbol{\xi} \;=\; \boldsymbol{\sigma} - \boldsymbol{\alpha}, \qquad \mathbf{s}_\xi \;=\; \mathbf{s} - \boldsymbol{\alpha} \]
(the second identity uses that \(\boldsymbol{\alpha}\) is deviatoric). The yield function rewrites as \[ f(\mathbf{s}_\xi, \bar{\varepsilon}^p) \;=\; \|\mathbf{s}_\xi\| - \sqrt{\tfrac{2}{3}}\,\sigma_y(\bar{\varepsilon}^p) \;\le\; 0, \]
which makes the algorithm essentially a radial return in \(\mathbf{s}_\xi\)-space.
22.3 Step 1 — Elastic Predictor
Given the state at \(t_n\), namely \(\boldsymbol{\sigma}_n,\ \boldsymbol{\varepsilon}^p_n,\ \boldsymbol{\alpha}_n,\ \bar{\varepsilon}^p_n\), and the prescribed total strain increment \(\Delta\boldsymbol{\varepsilon}\), freeze plastic flow and form the trial state (denoted by a “tr” superscript):
- Trial stress: \[ \boldsymbol{\sigma}^\text{tr}_{n+1} \;=\; \boldsymbol{\sigma}_n + \mathbb{C}^e : \Delta\boldsymbol{\varepsilon}. \]
- Trial backstress (unchanged over an elastic step): \[ \boldsymbol{\alpha}^\text{tr}_{n+1} \;=\; \boldsymbol{\alpha}_n. \]
- Trial equivalent plastic strain (unchanged): \[ \bar{\varepsilon}^{p,\text{tr}}_{n+1} \;=\; \bar{\varepsilon}^p_n. \]
- Trial deviatoric effective stress: \[ \mathbf{s}^\text{tr}_\xi \;=\; \text{dev}(\boldsymbol{\sigma}^\text{tr}_{n+1}) - \boldsymbol{\alpha}^\text{tr}_{n+1}. \]
22.4 Step 2 — Yield Check
Evaluate the yield function at the trial state: \[ f^\text{tr} \;=\; \|\mathbf{s}^\text{tr}_\xi\| - \sqrt{\tfrac{2}{3}}\,\sigma_y(\bar{\varepsilon}^{p,\text{tr}}_{n+1}) \;=\; \|\mathbf{s}^\text{tr}_\xi\| - \sqrt{\tfrac{2}{3}}\,\bigl(\sigma_{y0} + H_\text{iso}\,(\bar{\varepsilon}^p_n)^m\bigr). \]
- If \(f^\text{tr} \le \text{TOL}\): the step is elastic. Accept the trial state, \[ \boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}^\text{tr}_{n+1},\quad \boldsymbol{\varepsilon}^p_{n+1} = \boldsymbol{\varepsilon}^p_n,\quad \boldsymbol{\alpha}_{n+1} = \boldsymbol{\alpha}_n,\quad \bar{\varepsilon}^p_{n+1} = \bar{\varepsilon}^p_n, \] and exit.
- If \(f^\text{tr} > \text{TOL}\): plastic yielding occurs — proceed to Step 3.
22.5 Step 3 — Plastic Corrector: Core Equations
Flow rule (associative, J2): \[ \Delta\boldsymbol{\varepsilon}^p \;=\; \Delta\gamma\,\frac{\partial f}{\partial \boldsymbol{\sigma}} \;=\; \Delta\gamma\,\frac{\mathbf{s}_{\xi, n+1}}{\|\mathbf{s}_{\xi, n+1}\|}, \qquad \Delta\gamma \ge 0. \]
Isotropic hardening update: \[ \bar{\varepsilon}^p_{n+1} \;=\; \bar{\varepsilon}^p_n + \sqrt{\tfrac{2}{3}}\,\Delta\gamma. \]
Kinematic hardening update: \[ \boldsymbol{\alpha}_{n+1} \;=\; \boldsymbol{\alpha}_n + H_\text{kin}\,\Delta\boldsymbol{\varepsilon}^p. \]
Stress update: \[ \boldsymbol{\sigma}_{n+1} \;=\; \boldsymbol{\sigma}^\text{tr}_{n+1} - \mathbb{C}^e : \Delta\boldsymbol{\varepsilon}^p. \]
Deviatoric effective stress update. Expand \(\mathbf{s}_{\xi, n+1} = \text{dev}(\boldsymbol{\sigma}_{n+1}) - \boldsymbol{\alpha}_{n+1}\) and substitute the two update equations above: \[ \mathbf{s}_{\xi, n+1} \;=\; \bigl(\text{dev}(\boldsymbol{\sigma}^\text{tr}_{n+1}) - 2\mu\,\Delta\boldsymbol{\varepsilon}^p\bigr) - \bigl(\boldsymbol{\alpha}_n + H_\text{kin}\,\Delta\boldsymbol{\varepsilon}^p\bigr). \]
Here \(\mathbb{C}^e : \Delta\boldsymbol{\varepsilon}^p = 2\mu\,\Delta\boldsymbol{\varepsilon}^p\) because \(\Delta\boldsymbol{\varepsilon}^p\) is deviatoric (so the volumetric \(\kappa\,\text{tr}(\Delta\boldsymbol{\varepsilon}^p)\,\mathbf{I}\) term vanishes). Rearranging gives \[ \mathbf{s}_{\xi, n+1} \;=\; \mathbf{s}^\text{tr}_\xi - (2\mu + H_\text{kin})\,\Delta\boldsymbol{\varepsilon}^p, \]
and substituting the flow rule for \(\Delta\boldsymbol{\varepsilon}^p\): \[ \mathbf{s}_{\xi, n+1} \;=\; \mathbf{s}^\text{tr}_\xi - (2\mu + H_\text{kin})\,\Delta\gamma\,\frac{\mathbf{s}_{\xi, n+1}}{\|\mathbf{s}_{\xi, n+1}\|}. \]
Radial return in effective-stress space. The previous equation shows that \(\mathbf{s}_{\xi, n+1}\) is a positive scalar multiple of \(\mathbf{s}^\text{tr}_\xi\), so they share a common direction: \[ \frac{\mathbf{s}_{\xi, n+1}}{\|\mathbf{s}_{\xi, n+1}\|} \;=\; \frac{\mathbf{s}^\text{tr}_\xi}{\|\mathbf{s}^\text{tr}_\xi\|} \;\equiv\; \mathbf{n}^\text{tr}_\xi. \]
Taking norms, \[ \|\mathbf{s}_{\xi, n+1}\| \;=\; \|\mathbf{s}^\text{tr}_\xi\| - (2\mu + H_\text{kin})\,\Delta\gamma. \tag{$\star$} \]
This is the “radial return” relation for the effective deviatoric stress.
22.6 Step 4 — Solving for the Plastic Multiplier \(\Delta\gamma\)
Consistency at \(t_{n+1}\) requires \(f(\mathbf{s}_{\xi, n+1}, \bar{\varepsilon}^p_{n+1}) = 0\): \[ \|\mathbf{s}_{\xi, n+1}\| \;=\; \sqrt{\tfrac{2}{3}}\,\sigma_y(\bar{\varepsilon}^p_{n+1}) \;=\; \sqrt{\tfrac{2}{3}}\,\Bigl(\sigma_{y0} + H_\text{iso}\,\bigl(\bar{\varepsilon}^p_n + \sqrt{\tfrac{2}{3}}\,\Delta\gamma\bigr)^m\Bigr). \]
Equating this with \((\star)\) and bringing everything onto one side gives a non-linear scalar equation in \(\Delta\gamma\) (non-linear because of the exponent \(m\)): \[ R(\Delta\gamma) \;=\; (2\mu + H_\text{kin})\,\Delta\gamma \;+\; \sqrt{\tfrac{2}{3}}\,\Bigl(\sigma_{y0} + H_\text{iso}\,\bigl(\bar{\varepsilon}^p_n + \sqrt{\tfrac{2}{3}}\,\Delta\gamma\bigr)^m\Bigr) \;-\; \|\mathbf{s}^\text{tr}_\xi\| \;=\; 0. \]
The trial overshoot \(f^\text{tr} = \|\mathbf{s}^\text{tr}_\xi\| - \sqrt{2/3}\,\sigma_y(\bar{\varepsilon}^p_n)\) is the “gap” the corrector must close; in particular \(R(0) = -f^\text{tr} < 0\) and \(R\) is monotonically increasing in \(\Delta\gamma\) for \(m \in [0, 1]\), so there is a unique root \(\Delta\gamma > 0\).
Local Newton–Raphson. At iterate \(j\), \[ \Delta\gamma^{(j+1)} \;=\; \Delta\gamma^{(j)} - \frac{R(\Delta\gamma^{(j)})}{R'(\Delta\gamma^{(j)})}, \]
with derivative \[ R'(\Delta\gamma) \;=\; (2\mu + H_\text{kin}) \;+\; \sqrt{\tfrac{2}{3}}\, H_\text{iso}\, m\,\bigl(\bar{\varepsilon}^p_n + \sqrt{\tfrac{2}{3}}\,\Delta\gamma\bigr)^{m-1}\cdot\sqrt{\tfrac{2}{3}} \] \[ \;=\; (2\mu + H_\text{kin}) \;+\; \tfrac{2}{3}\, H_\text{iso}\, m\,\bigl(\bar{\varepsilon}^p_n + \sqrt{\tfrac{2}{3}}\,\Delta\gamma\bigr)^{m-1}. \]
Practical notes:
- A good initial guess is \(\Delta\gamma^{(0)} = f^\text{tr} / (2\mu + H_\text{kin} + \tfrac{2}{3}H_\text{iso}\, m\,(\bar{\varepsilon}^p_n)^{m-1})\) — this is the first Newton step from zero with the argument frozen. A simpler fallback is \(\Delta\gamma^{(0)} = 0\).
- Stopping criterion: iterate until \(|R(\Delta\gamma^{(j+1)})| < \text{TOL}_R\) or \(|\Delta\gamma^{(j+1)} - \Delta\gamma^{(j)}| < \text{TOL}_\gamma\).
- Edge cases. If \(m = 0\) the isotropic term is constant (perfectly-plastic-isotropic with kinematic hardening still active), \(R\) is linear in \(\Delta\gamma\), and the closed-form solution is \(\Delta\gamma = f^\text{tr} / (2\mu + H_\text{kin})\). If \(m < 1\) and \(\bar{\varepsilon}^p_n = 0\), the derivative \(R'\) diverges at \(\Delta\gamma \to 0\) — regularise the first iterate (e.g. replace \(\bar{\varepsilon}^p_n\) by \(\max(\bar{\varepsilon}^p_n, \varepsilon_\text{reg})\) with \(\varepsilon_\text{reg} \sim 10^{-12}\)) or switch to a small-\(\bar{\varepsilon}^p\) closed-form starter.
- Because \(R\) is monotone and convex-enough on \([0, \Delta\gamma^\star]\) for \(m \in [0, 1]\), Newton converges quadratically from \(\Delta\gamma^{(0)} = 0\) in practice.
22.7 Step 5 — Final Updates
Once \(\Delta\gamma\) has converged to within tolerance:
- Equivalent plastic strain: \[ \bar{\varepsilon}^p_{n+1} \;=\; \bar{\varepsilon}^p_n + \sqrt{\tfrac{2}{3}}\,\Delta\gamma. \]
- Flow direction and plastic-strain increment: \[ \mathbf{n}^\text{tr}_\xi \;=\; \frac{\mathbf{s}^\text{tr}_\xi}{\|\mathbf{s}^\text{tr}_\xi\|}, \qquad \Delta\boldsymbol{\varepsilon}^p \;=\; \Delta\gamma\,\mathbf{n}^\text{tr}_\xi. \] (Well-defined because \(f^\text{tr} > 0\) implies \(\|\mathbf{s}^\text{tr}_\xi\| \neq 0\).)
- Plastic strain tensor: \[ \boldsymbol{\varepsilon}^p_{n+1} \;=\; \boldsymbol{\varepsilon}^p_n + \Delta\boldsymbol{\varepsilon}^p. \]
- Backstress: \[ \boldsymbol{\alpha}_{n+1} \;=\; \boldsymbol{\alpha}_n + H_\text{kin}\,\Delta\boldsymbol{\varepsilon}^p. \]
- Cauchy stress: \[ \boldsymbol{\sigma}_{n+1} \;=\; \boldsymbol{\sigma}^\text{tr}_{n+1} - \mathbb{C}^e : \Delta\boldsymbol{\varepsilon}^p \;=\; \boldsymbol{\sigma}^\text{tr}_{n+1} - 2\mu\,\Delta\boldsymbol{\varepsilon}^p, \] where the second form again uses that \(\Delta\boldsymbol{\varepsilon}^p\) is deviatoric, so \(\text{tr}(\boldsymbol{\sigma}_{n+1}) = \text{tr}(\boldsymbol{\sigma}^\text{tr}_{n+1})\) (the hydrostatic part is preserved by the corrector).
22.8 Consistent (Algorithmic) Tangent — Outline
The consistent tangent \(\mathbb{C}^\text{alg}_{n+1} = d\boldsymbol{\sigma}_{n+1}/d\boldsymbol{\varepsilon}_{n+1}\) for this combined-hardening model is notably more involved than the linear-hardening J2 case handled in L08. Two new complications appear:
- The power-law isotropic term \(H_\text{iso}\,(\bar{\varepsilon}^p)^m\) makes \(\Delta\gamma\) the solution of a non-linear equation, so \(d(\Delta\gamma)/d(\Delta\boldsymbol{\varepsilon})\) must be obtained by implicit differentiation of \(R(\Delta\gamma) = 0\).
- The kinematic term adds a contribution \(H_\text{kin}\,\Delta\boldsymbol{\varepsilon}^p\) to the effective-stress update, which itself depends on \(\Delta\gamma\).
Option A — analytical tangent via implicit differentiation. Differentiate every equation in Steps 3–5 with respect to \(\Delta\boldsymbol{\varepsilon}\) and eliminate intermediate differentials using \(R'(\Delta\gamma)\) derived above. The resulting \(\mathbb{C}^\text{alg}_{n+1}\) retains the same \((\kappa\,\mathbf{I}\otimes\mathbf{I},\ \mathbb{P}_\text{dev},\ \mathbf{n}\otimes\mathbf{n})\) structure as the J2-linear case in L08, but the scalar coefficients \(\theta_1, \theta_2\) are modified by \(H_\text{kin}\) and by the effective hardening slope \(\tfrac{2}{3}\,H_\text{iso}\,m\,(\bar{\varepsilon}^p_{n+1})^{m-1}\). The normal \(\mathbf{n}\) is evaluated at \(\mathbf{s}^\text{tr}_\xi\) (the return is radial in effective-stress space).
Option B — numerical tangent. Perturb each independent component of \(\boldsymbol{\varepsilon}_{n+1}\) in turn, rerun the stress update, and form \(\mathbb{C}^\text{alg}_{n+1}\) by finite differences. More expensive (six stress updates per Gauss point in 3D) but robust; useful during model development and for cross-checking analytical derivations.
Either way, the computed \(\mathbb{C}^\text{alg}_{n+1}\) must be used (not the continuum \(\mathbb{C}^{ep}\)) for the global Newton–Raphson assembly in FEM — otherwise quadratic convergence is lost. See the “Consistent (Algorithmic) Tangent” section of L08 for the underlying argument.
22.9 Suggested Self-Study Exercises
- Recover linear isotropic hardening. Set \(H_\text{kin} = 0\) and \(m = 1\) and show that \(R(\Delta\gamma) = 0\) reduces to the closed form \(\Delta\gamma = f^\text{tr} / (2\mu + \tfrac{2}{3}H_\text{iso})\) used in L08.
- Recover pure kinematic hardening. Set \(H_\text{iso} = 0\) and show that \(R\) becomes linear in \(\Delta\gamma\); derive the closed form.
- Implement the local Newton solve for \(m = 1/2\). Pick representative values (\(\mu = 80\) GPa, \(\sigma_{y0} = 250\) MPa, \(H_\text{iso} = 500\) MPa, \(H_\text{kin} = 1\) GPa, \(\bar{\varepsilon}^p_n = 0.01\)) and a trial state with \(f^\text{tr}/\sigma_{y0} \in \{0.01,\ 0.1,\ 1\}\). Plot \(R(\Delta\gamma)\) and report the number of Newton iterations to reach \(|R| < 10^{-10}\).
- Derive the analytical consistent tangent. Follow Option A above to obtain \(\mathbb{C}^\text{alg}_{n+1}\) in closed form. Compare against the numerical tangent of Option B (e.g., central differences with perturbation \(10^{-6}\)) and quantify the component-wise error.
- Extend to Armstrong–Frederick kinematic hardening. Replace the Prager rule with \[ \dot{\boldsymbol{\alpha}} \;=\; \tfrac{2}{3}\,C\,\dot{\boldsymbol{\varepsilon}}^p - \gamma_\text{AF}\,\boldsymbol{\alpha}\,\dot{\bar{\varepsilon}}^p. \] What changes in Steps 3 and 4? Does the return remain radial in effective-stress space? What form does \(R(\Delta\gamma)\) take now?