Worked example

Damage mechanics in a UMAT

Damage models fail in ways plasticity models do not: the tangent loses positive definiteness by design, and the answer depends on the mesh unless you do something about it.

What ddsdde does for this model today: the schema holds it, the verification harness runs the tests that apply, the Jacobian check works on any material point you can express, and the diagnosis reads the job. Code generation currently emits J2 plasticity only — a damage model still means writing the subroutine yourself. This page is the guide for doing that.

A damage model is a stiffness that decreases. Everything difficult about implementing one follows from that single fact.

Choosing the damage variable

The usual form degrades the effective stress: sigma = (1 - D) C : eps_e, with D running from zero to one. Three decisions come with it.

Scalar or tensor? A scalar D is isotropic damage and is much easier. A tensor damage variable represents direction-dependent degradation and must be rotated with the material in finite strain, which a scalar must not.
What drives it? Equivalent plastic strain, an energy release rate, a strain measure. Whatever it is, it has to be a state variable too, so damage typically costs at least two STATEV slots rather than one.
Where does it saturate? A model that lets D reach exactly one leaves an element with no stiffness at all, and the assembled matrix becomes singular. A residual stiffness, or a cap slightly below one, is not a fudge — it is a modelling decision that should be recorded as such.

Irreversibility, and how it gets lost

Damage does not heal. That sounds trivial to enforce and is the most common bug in this class of model, because the natural way to write the update does not enforce it:

The wrong way and the right way

C     Wrong: D follows the driving variable up and back down again.
      D = DAMAGE_FUNCTION(KAPPA)

C     Right: D can only increase.
      DNEW = DAMAGE_FUNCTION(KAPPA)
      D    = MAX(DOLD, DNEW)

The symptom of getting this wrong is a material that recovers stiffness on unloading. Under monotonic loading it is invisible; under any cyclic loading it produces hysteresis loops that grow rather than stabilise.

The same applies to the history variable that drives damage. If kappa is the largest strain measure ever reached, it needs the same MAX. Take it against the value from the start of the increment, not from the last equilibrium iteration.

The tangent, which is now unsymmetric and indefinite

Differentiating sigma = (1 - D) C : eps_e gives two terms:

dsigma/deps = (1 - D) C  -  (C : eps_e) (x) dD/deps

The second term is the one that causes trouble. Three consequences follow, and all three are worth knowing before the first run rather than after it.

It is generally unsymmetric. Unless dD/deps happens to be parallel to C : eps_e, the tangent is not symmetric, and Abaqus will symmetrise it silently unless the step carries UNSYMM=YES.
It can be indefinite. Once softening begins the tangent loses positive definiteness, and negative eigenvalue warnings are expected rather than a sign of a bug. That is a real problem for diagnosis: the message that usually means a wrong tangent now means the model is working.
It is only defined where damage is growing. During elastic unloading dD/deps is zero and the tangent is (1 - D) C. Getting this branch wrong makes unloading increments cost far more iterations than they should.

Mesh dependence, and what to do about it

A local softening model localises into a band one element wide. Refine the mesh and the band gets narrower, the energy dissipated goes down, and the predicted failure load falls without limit. The result is not converging on anything.

This is a property of the mathematics, not a bug in your subroutine, and no amount of care in the implementation fixes it. The usual remedies:

Approach	What it does	What it costs
Characteristic length (`CELENT`)	Scale the softening branch by the element size so the dissipated energy per unit area is mesh independent	Cheap and local. Sensitive to element shape and to the loading direction relative to the mesh
Viscous regularisation	Add rate dependence so the localisation band has a width set by the viscosity and the loading rate	One extra parameter with no physical meaning, and results that depend on the time step
Non-local or gradient formulations	Make damage depend on an averaged quantity over a length scale	Substantially harder to implement, and generally needs more than a UMAT

CELENT is supplied to your UMAT for exactly this purpose. Using it makes the model mesh-size aware; it does not make it mesh independent, and the distinction is worth stating in the report.

What to verify

The general checks all still apply. These are the ones specific to damage:

Damage never decreases. Load, unload, reload. D at the end of the cycle must be at least what it was at the peak.
Unloading is elastic at the degraded stiffness. Unload from a damaged state and confirm the slope is (1 - D) times the original.
Dissipation is non-negative in every increment. A damage model that creates energy has a sign error, and it shows up here before it shows up anywhere else.
The tangent still matches finite differences. Both in the damaging branch and in the unloading branch, separately. The unloading branch is the one people get wrong.
Two mesh densities. If the answers differ by a lot, your regularisation is not working, and that is a modelling finding rather than an implementation bug.