Integrating a law over a step

A nonlinear law is not a formula that turns a strain into a stress: the stress also depends on the history, carried by internal variables that evolve. A finite-element code, or a test machine, asks the law a discrete question: the gradients were at the beginning of a time step of length , the internal variables were ; the gradients are at its end; what are the fluxes, the internal variables, and how do the fluxes vary with ? This page derives the answer the code computes, and why it is exact.

The objects of a step

A law has three kinds of variables, which the code stores as flat vectors (in sans-serif in this page):

	meaning	example (poroelastoplasticity)
	the gradients: the primal fields at the point	strain and pore pressure
	the fluxes conjugate to them	stress and change of fluid content
	the internal variables: what the point remembers	plastic strain , cumulated plastic strain

A symmetric second-order tensor enters these vectors by its six Kelvin–Mandel components . The factor makes the basis orthonormal: the double contraction   is the dot product of the two vectors, and a fourth-order tensor with the minor symmetries acts on them as a   matrix. Nothing in a law has to remember a factor 2 on the shears.

The residual route

Integrating the evolution laws of the internal variables over the step by an implicit scheme (backward Euler in all the laws shipped) gives as many equations as internal variables, written as a residual

For von Mises plasticity, for instance,   and the residual collects the discrete normality rule and the consistency condition (see J2Plasticity):

This is the standard form of the implicit integration of inelastic laws (Simo and Hughes, 1998). The residual is solved by Newton's method, with the Jacobian obtained by automatic differentiation of the residual the user wrote. Convergence is judged on the size of the Newton correction, which has the units of whatever the units of the residual.

The elastic predictor

A rate-independent law has two regimes. If the stress computed with the internal variables frozen at (the trial state) is admissible, the step is elastic and  ; only otherwise is the residual solved. A law gives this test as its elastic function; a damage law gives its threshold the same way. On an elastic step the internal variables do not depend on , which the tangent below uses.

The consistent tangent

A global Newton method on the structure needs the derivative of the fluxes at the end of the step with respect to the gradients at the end of the step, through the integration algorithm: is itself a function of , defined implicitly by  . Differentiating  ,

and the chain rule on gives the tangent

On an elastic step the second term vanishes. This is the tangent consistent with the integration algorithm, the one that keeps the quadratic convergence of the structural Newton method (Simo and Taylor, 1985); the tangent of the continuous law (the "continuum" tangent) does not. The four partial derivatives are again computed by automatic differentiation, so a law never contains a derivative written by hand.

For a law of several fields the formula is unchanged: is the full matrix of the derivatives of every flux with respect to every gradient. Its blocks are the coupling terms a coupled finite-element formulation assembles, , , , for poroelasticity.

The potential route

Many laws derive from two potentials, in the framework of generalized standard materials (Germain et al., 1983): a free energy of the gradients and of the internal variables , which gives the fluxes,  , and a convex dissipation potential , which gives the evolution of the internal variables through  .

Over a step, replace the rate by the difference quotient   and consider the incremental energy

Its stationarity in ,   , is exactly the backward Euler discretization of the evolution law, so the internal variables at the end of the step minimize when the potentials are convex. This variational form of the update is due to Ortiz and Stainier (1999). In the code it becomes a residual like any other,

with its Jacobian the Hessian of : the potential route reuses the whole residual machinery, Newton's method and the consistent tangent included, and automatic differentiation provides the first and second derivatives of the two potentials.

Each flux may carry a sign: a poroelastic free energy written in the strain and the pore pressure gives the stress as but the change of fluid content as (see BiotPoroelasticity).

Derivatives with respect to the parameters, exact

Identifying a law means differentiating its response with respect to its parameters (moduli, yield stress, viscosity). The internal variables depend on them through  , and by the same argument as for the tangent,

The code obtains this without forming it. It first solves   on the values alone, which gives a root . It then takes one Newton step in dual numbers, from that root:

The value of is zero, and its derivative with respect to is , since carries none. The derivative of is therefore exactly the formula above. Newton's method doubles the order of agreement at each step, so such steps give every derivative up to order  . This is what nested differentiation needs: one step for a gradient, two for a Hessian.

No tolerance is ever applied to a derivative, and no finite difference is used anywhere. The same holds for a dual strain or a dual time step, which is how the tests check the consistent tangent: the Jacobian of the whole step, computed with dual numbers carried through the predictor and the Newton iterations, equals .

Driving a point: mixed control

A material test prescribes some components by their gradient (an axial strain) and others by their flux (free lateral faces: zero lateral stress). Let mark the components driven by their gradient. At each time the driver looks for the free gradients such that the fluxes reach their targets , by Newton's method on

The consistent tangent makes this iteration converge quadratically; with a linear hardening law a single correction is exact. A step whose local or global iteration fails is retried as two half steps, recursively, which is also what a finite-element driver does.