Worked example

J2 plasticity, end to end

The model every implementer starts with, worked through completely — including the numbers each check produces, so you can compare against your own.

Rate-independent von Mises plasticity with associated flow, integrated by radial return. It is the benchmark because a mechanician can judge it by eye, and because it exercises every part that is hard: a local solve, a state variable, a branch, and a tangent that differs from the elastic one.

The model

Elastic predictor, yield function, hardening, consistency, flow, return

sigma_trial = sigma_n + C_e : delta_eps

f = sqrt(3/2 s:s) - sigma_y(eps_p)         von Mises yield

sigma_y = sigma_y0 + H eps_p               linear isotropic hardening

q_trial - 3 mu dgamma - sigma_y(eps_p + dgamma) = 0     consistency

n = (3 / (2 q_trial)) s_trial              flow direction, ||n||^2 = 3/2

sigma = sigma_trial - 2 mu dgamma n         radial return

For linear hardening the consistency condition solves in closed form, dgamma = (q_trial - sigma_y) / (3 mu + H). For a saturating law such as Voce it does not. The increment then needs a local Newton iteration, which is the more interesting case and the one worth building for.

The consistent tangent, derived rather than quoted

Parameterise the return by the equivalent plastic strain increment and differentiate with respect to the strain increment. Because n is deviatoric, C_e : n = 2 mu n, so:

d q_trial / d eps = 2 mu n
d dgamma  / d eps = 2 mu n / (3 mu + H)
d n       / d eps = (3 mu / q_trial) I_dev - (2 mu / q_trial) n (x) n

C = K 1(x)1 + 2 mu theta I_dev - (4 mu / 3) thetabar n(x)n

theta    = 1 - 3 mu dgamma / q_trial
thetabar = 3 mu / (3 mu + H) - (1 - theta)

Setting dgamma = 0 gives theta = 1 and recovers the continuum tangent. In Voigt form with engineering shear the shear diagonal is mu theta, not twice it. The n(x)n block needs no extra factor: the doubling in the contraction cancels the halving in gamma = 2 eps.

The subroutine

This is the generated implementation, verbatim. The Newton loop uses a guarded update rather than an exit. Fixed-form Fortran has no loop exit, and a GOTO in a subroutine other people must modify is a trade not worth making. Once converged the body does nothing.

The local solve, with the hardening slope from the declared law

        DEQPL = 0.D0
        TOLER = 1.D-8 * SYIELD
        DO KNEWT = 1, 30
          SYIELD = SY0 + QINF * (1.D0 - EXP(-(BSAT * (EQPLAS + DEQPL))))
C         Slope at the end of the increment, differentiated from the law above.
          HARDP = QINF * (EXP(-(BSAT * (EQPLAS + DEQPL))) * BSAT)
          RESID = SEQ - 3.D0 * XMU * DEQPL - SYIELD
          IF (ABS(RESID) .GT. TOLER) THEN
            DEQPL = DEQPL + RESID / (3.D0 * XMU + HARDP)
          END IF
        END DO

HARDP is the derivative of the hardening law evaluated at the end of the increment. Using the initial slope instead is the most common way a nonlinear hardening implementation loses its quadratic convergence while still giving the right answer.

What the checks produce

Numbers from the shipped benchmark, so you have something to compare your own against. Anything below about 1e-6 is the numerical method rather than the model.

Material Jacobian against a central-difference reference, ten states
State	Relative Frobenius error
Uniaxial extension, elastic	4.0e-12
Hydrostatic, elastic	8.5e-12
Uniaxial, first plastic increment	1.6e-11
Uniaxial, developed plasticity	3.5e-11
Elastic unloading after plastic flow	9.2e-13
Simple shear, plastic	4.6e-12
Combined tension and shear, plastic	4.1e-11
All three shear components, plastic	3.2e-12

Convergence classification: quadratic, iteration matrix spectral radius 7.1e-9. On the verification side, sixteen tests pass and two do not apply. The rotation test is skipped because the model declares small-strain kinematics. The tangent-modulus test is skipped for the saturating law, because a Voce law has no single closed-form tangent modulus to check against.

Deliberately breaking it

The instructive part. Each of these is a real mistake, and what separates them is not the size of the error but its behaviour.

Injected fault	Jacobian check	Physics tests
Elastic tangent returned after yield	FAIL — matrix unchanged across yield	All pass
Continuum instead of algorithmic tangent	FAIL — error vanishes with the increment, slope 1.12	All pass
Shear slots in Explicit order	FAIL at one probe out of ten	All pass
Shear diagonal written as 2G	FAIL — halving it reconstructs the reference	All pass
Sign error in the plastic correction	FAIL — the discrepancy is a single rank-one term	All pass
PROPS read in the wrong order	PASS	FAIL — yield begins at 2000 rather than 250
State never written back to STATEV	PASS	FAIL — an inelastic increment leaves STATEV untouched
Plastic strain accumulated with a sign	PASS	FAIL — peak stress falls between cycles

The single-element deck

One C3D8, three symmetry planes, uniaxial stress

*ELEMENT, TYPE=C3D8, ELSET=ONE
1, 1, 2, 3, 4, 5, 6, 7, 8
*SOLID SECTION, ELSET=ONE, MATERIAL=UMATMAT
*MATERIAL, NAME=UMATMAT
*USER MATERIAL, CONSTANTS=5
200000., 0.3, 250., 180., 12.,
*DEPVAR
1,
*BOUNDARY
XSYM, 1
YSYM, 2
ZSYM, 3
*STEP, NLGEOM=NO, INC=1000
*STATIC
0.05, 1.0, 1.E-8, 0.05

A unit cube means the prescribed displacement is the strain and the reaction is the stress. Compare S11 at each increment against the expected values, then read the message file. Three to four equilibrium iterations an increment is a consistent tangent; eight to twelve is an approximate one.