4. Metal plasticity at finite strain

When a metal component undergoes large deformation — a crash box folding, a sheet-metal hem being formed, a bolt being torqued into yield — the small-strain plasticity theory used in most introductory FEA courses is no longer adequate. Strains exceed a few percent, the material rotates substantially, and the difference between the undeformed and deformed configurations can no longer be ignored. This chapter develops the finite-strain metal plasticity framework used in explicit and implicit solvers: the von Mises yield criterion, isotropic and kinematic hardening, the radial return mapping algorithm, and the objective integration of stress that keeps the constitutive response independent of rigid-body rotation.

Kinematics: Splitting Deformation into Elastic and Plastic Parts

At finite strain, the multiplicative decomposition is the standard starting point. The deformation gradient is split as

F = F_e · F_p

where F_p maps the reference configuration to a stress-free intermediate configuration and F_e maps that intermediate configuration to the current configuration. The velocity gradient in the current configuration is l = Ḟ · F⁻¹, and its symmetric part d = sym(l) is the rate of deformation. The additive split d = d_e + d_p holds in the current configuration, which is what makes the return-mapping algorithm tractable.

For metals, plastic flow is assumed isochoric (volume-preserving), so tr(d_p) = 0. This is a direct consequence of the pressure-insensitivity of the von Mises yield surface. The elastic part of the deformation is typically small even when total strains are large, which justifies a hypoelastic stress–strain relation of the form

σ̊ = C : d_e = C : (d − d_p)

where σ̊ is an objective stress rate and C is the fourth-order elastic stiffness tensor. The choice of objective rate is not cosmetic — it determines whether the integration scheme is truly frame-indifferent.

The von Mises Yield Criterion and Hardening

The von Mises yield function is expressed in terms of the deviatoric stress s = σ − (1/3)tr(σ)I:

f = √(3/2 · s : s) − σ_y(ε̄_p) = σ_eq − σ_y(ε̄_p)

Here σ_eq is the equivalent (von Mises) stress and σ_y is the current yield stress. The yield surface is a cylinder in principal stress space, centered on the hydrostatic axis. Plastic flow is normal to the yield surface (associated flow), giving the plastic flow rule

d_p = γ̇ · ∂f/∂σ = γ̇ · (3/2) · s / σ_eq

with γ̇ ≥ 0 the plastic multiplier. The equivalent plastic strain rate is ε̄̇_p = γ̇, which makes the hardening law directly expressible in terms of a single scalar internal variable.

Isotropic hardening expands the yield surface uniformly. A common form is the Swift power law:

σ_y(ε̄_p) = K · (ε₀ + ε̄_p)^n

or the simpler linear form σ_y = σ_y0 + H·ε̄_p. Isotropic hardening captures the Bauschinger effect poorly — it predicts that yield strength increases equally in tension and compression after plastic flow.

Kinematic hardening translates the yield surface instead of expanding it, using a back-stress tensor α:

f = √(3/2 · (s − α) : (s − α)) − σ_y0

The evolution of the back stress follows the Armstrong–Frederick rule:

α̇ = (2/3)·C·d_p − γ·α·ε̄̇_p

The first term is linear Prager hardening; the second is a dynamic recovery term that saturates the back stress. This combination reproduces the Bauschinger effect and is essential for cyclic loading, springback prediction, and low-cycle fatigue. In practice, combined isotropic–kinematic hardening is used: the yield surface both expands and translates, with the isotropic part capturing monotonic hardening and the kinematic part capturing load-path reversal.

Radial Return Mapping

The return-mapping algorithm is the workhorse of implicit finite-strain plasticity. Given the state at step n and an increment of deformation, the algorithm proceeds in two stages.

Step 1 — Elastic predictor. Assume the increment is purely elastic and compute a trial stress:

σ_trial = σ_n + C : Δε
s_trial = dev(σ_trial)
σ_eq_trial = √(3/2 · s_trial : s_trial)

If σ_eq_trial ≤ σ_y(ε̄_p,n), the step is elastic and the trial state is accepted. Otherwise, plastic flow has occurred and a return mapping is required.

Step 2 — Plastic corrector. The consistency condition f(σ_n+1, ε̄_p,n+1) = 0 must be satisfied. For von Mises plasticity with associated flow, the return direction is radial in the deviatoric plane, giving the closed-form update:

Δγ = (σ_eq_trial − σ_y(ε̄_p,n)) / (3G + H)
s_n+1 = s_trial − 2G·Δγ·(3/2)·s_trial/σ_eq_trial
ε̄_p,n+1 = ε̄_p,n + Δγ

where G is the shear modulus and H is the hardening modulus dσ_y/dε̄_p. For nonlinear hardening, H depends on Δγ, so the scalar equation for Δγ is solved by Newton iteration. The consistent tangent modulus — the derivative of σ_n+1 with respect to Δε — must be derived from the algorithmic update, not the continuum elastoplastic tangent, to preserve quadratic convergence in implicit solvers.

Objective Integration and the Jaumann Rate

At finite strain, the stress update must be objective — the computed stress must not change when a rigid-body rotation is superimposed on the deformation. The most common approach in commercial codes is the Jaumann rate with corotational integration:

σ̊_J = σ̇ − w·σ + σ·w

where w = skew(l) is the spin tensor. The Jaumann rate is objective but has a known deficiency: it predicts spurious shear stress under large simple shear. Alternatives include the Green–Naghdi rate (using the rotation from the polar decomposition F = R·U) and the Truesdell rate. For most metal forming and crash applications with moderate rotations per step, the Jaumann rate with small time increments is adequate. For severe shear localization or high-speed machining, the Green–Naghdi rate is preferred.

A practical implementation detail: the stress must be rotated from the previous configuration to the current one before the return mapping, then rotated back if a corotational frame is used. In explicit codes, the rotation is often applied incrementally using the incremental rotation tensor ΔR = (I − 0.5·Δt·w)⁻¹·(I + 0.5·Δt·w), which is second-order accurate and unconditionally stable for the rotation update.

Engineering Example: Uniaxial Tension to 15% Strain

Consider a 6061-T6 aluminum tensile specimen modeled with one 8-node brick element under displacement control. The material parameters are representative of published data: E = 68.9 GPa, ν = 0.33, σ_y0 = 276 MPa, isotropic Swift hardening with K = 410 MPa, ε₀ = 0.002, n = 0.12. The specimen is pulled to 15% engineering strain in 100 equal increments.

At each increment, the solver computes the trial stress, checks the yield condition, and if yielded, solves for Δγ by Newton iteration. The expected output is a true stress–true strain curve that rises from 276 MPa at yield to approximately 340–350 MPa at 15% strain. The equivalent plastic strain at the end of the analysis should be approximately 0.13–0.14, slightly less than the total strain because of elastic unloading.

Typical failure symptoms and their causes:

• Stress oscillates or diverges in implicit analysis: the consistent tangent modulus was not derived from the algorithmic update, or the return mapping did not converge within the iteration limit. Check the Newton loop tolerance and the tangent implementation.
• Stress is too high at large strain: the objective rate is missing or incorrect, causing the stress to be integrated in a fixed frame. Verify that the spin tensor is included in the stress update.
• Plastic strain is negative or exceeds total strain: the yield condition check is using the wrong equivalent stress definition, or the hardening modulus has the wrong sign. Verify that σ_eq = √(3/2·s:s) and that H > 0 for hardening.
• No plastic flow occurs despite stress above yield: the trial stress is being compared to the wrong yield stress (e.g., using the initial yield instead of the current hardened value). Ensure the hardening variable is updated before the yield check.

Summary

Finite-strain metal plasticity rests on four pillars: the multiplicative decomposition of deformation, the von Mises yield criterion with isotropic and/or kinematic hardening, the radial return mapping algorithm for stress integration, and an objective stress rate to handle large rotations. The return mapping is exact for linear isotropic hardening and requires local Newton iteration for nonlinear or kinematic hardening. The consistent tangent modulus is essential for quadratic convergence in implicit solvers. When implementing or debugging a finite-strain plasticity model, always validate against the single-element uniaxial tension analytical solution first, then check objectivity by superimposing a rigid-body rotation and confirming that the stress tensor rotates with the material rather than remaining fixed in space.