The consistent tangent operator
==================

In any timestep, MOOSE solves for the new values of the displacement variables.  In any nonlinear iteration of this timestep, it proposes some new values and checks whether the residual is zero when using these proposed values.  If not, then the Jacobian is used to propose updated values for the displacement that will (hopefully!) reduce the residual closer to zero.  The Jacobian encodes the change in the residual with respect to changes in the displacement variables.

To calculate the Jacobian, MOOSE needs to know how the stress changes with respect to changes in the displacement variables.  To this end, define the consistent tangent operator, $$H_{ijkl}$$, through
$$\delta \sigma_{ij} = H_{ijkl} \delta \epsilon_{kl} \ .$$
In this formula $$\delta \epsilon$$ is an arbitrary change in the total strain (which occurs because the displacements are changed) and $$\delta \sigma$$ is the consequential change in the stress.  $$H_{ijkl}$$ enters MOOSE's Jacobian matrix.

In a purely elastic situation $$H_{ijkl}=E_{ijkl}$$ (the elasticity tensor) but the plastic situation is more complicated.   Consider the plastic case from now on.  In this case the proposed values of displacements for the current timestep resulted in an inadmissible stress, $$\sigma^{\mathrm{trial}}$$, that was brought back to the yield surface via the return-map process.  A small change in these proposed values of displacement (what we need to calculate the Jacobian) will also result in an inadmissible stress that is slightly different from $$\sigma^{\mathrm{trial}}$$ that will then be brought back to a slightly different position on the yield surface.  This new return process will produce a slightly different consistency parameters (plastic multipliers, $$\gamma$$) and final internal variables $$q$$.

General multisurface plasticity
==================

The following uses notation from the [Return Map article](../ReturnMap).  Use the symbol $$\delta$$ to denote a change in a quantity due to the small change in displacement variables.  Then
$$\delta \sigma_{ij} = H_{ijkl} \delta \epsilon_{kl} = E_{ijkl} (\delta \epsilon_{ij} - \delta \epsilon_{ij}^{\mathrm{plastic}}) \ .$$
But the stress after return-mapping from the changed displacement must still lie on the yield surface, so
$$0 = \delta f^{\alpha} = \frac{\partial f^{\alpha}}{\partial \sigma_{ij}} \delta\sigma_{ij} + \frac{\partial f^{\alpha}}{\partial q}\delta q \ .$$
Also, the normality condition $$\epsilon^{\mathrm{plastic}} = \gamma^{\alpha} r^{\alpha}$$ must be satisfied so that
$$\delta \epsilon^{\mathrm{plastic}} = \delta\gamma^{\alpha} r^{\alpha} + \gamma^{\alpha} \left( \frac{\partial r^{\alpha}}{\partial \sigma}\delta \sigma + \frac{\partial r^{\alpha}}{\partial q_a}\delta q_{a} \right)$$
Finally, the change of the internal parameters reads
$$0 = \delta q_{a} + \delta\gamma^{\alpha} h^{\alpha}_{a} + \gamma^{\alpha} \left( \frac{\partial h^{\alpha}_{a}}{\partial \sigma}\delta \sigma + \frac{\partial h^{\alpha}_{a}}{\partial q_b}\delta q_{b} \right)$$

Using straightforward but rather tedious algebra, these equations may be solved for $$H_{ijkl}$$.  The result is implemented in MOOSE.  However, the result is not particularly illuminating and I don't write it here.  The steps are quite similar to the Drucker-Prager case below.

Example of Drucker-Prager plasticity
==================
Consider just the case of a single yield surface described by
$$f = \sqrt{J_{2} + \epsilon^{2}} + B \mbox{Tr} \sigma - A$$
This is Drucker-Prager plasticity, with the cone's tip smoothed in a hyperbolic fashion.  Suppose there is just one internal parameter with hardening potential $$h=-1$$.  The above expressions contain $$\delta\gamma$$, $$\delta q$$ and $$\delta\epsilon_{ij}^{\mathrm{plastic}}$$, which must all be eliminated to yield an expression involving only $$\delta\sigma$$ and $$\delta\epsilon$$ from which $$H$$ can be deduced using
$$\delta \sigma_{ij} = H_{ijkl}\delta \epsilon_{kl}$$
This process is now performed explicitly.

Specialising to a single-surface with one hardening parameter and $$h=-1$$, the above expressions read
$$\delta \sigma_{ij} = H_{ijkl}\delta \epsilon_{kl} = E_{ijkl} (\delta \epsilon_{ij} - \delta \epsilon_{ij}^{\mathrm{plastic}})$$
and
$$-\frac{\partial f^{\alpha}}{\partial q}\delta \gamma = \frac{\partial f^{\alpha}}{\partial \sigma_{ij}} \delta\sigma_{ij}$$
Combining these two gives
$$-\frac{\partial f^{\alpha}}{\partial q}\delta \gamma = \frac{\partial f^{\alpha}}{\partial \sigma_{ij}} \delta\sigma_{ij} = \frac{\partial f^{\alpha}}{\partial \sigma_{ij}} E_{ijkl} (\delta \epsilon_{ij} - \delta \epsilon_{ij}^{\mathrm{plastic}})$$
For a single-surface with one internal parameter, the normality condition reads
$$\delta \epsilon^{\mathrm{plastic}} = \delta\gamma^{\alpha} r + \gamma \left( \frac{\partial r}{\partial \sigma}\delta \sigma + \frac{\partial r}{\partial q}\delta \gamma \right)$$
and using this in the preceding equation yields
$$-\frac{\partial f^{\alpha}}{\partial q}\delta \gamma = \frac{\partial f}{\partial \sigma_{ij}} E_{ijkl} \left[\delta \epsilon_{kl} - \delta\gamma^{\alpha} r_{kl} - \gamma \left( \frac{\partial r_{kl}}{\partial \sigma_{mn}}\delta \sigma_{mn} - \frac{\partial r_{kl}}{\partial q}\delta \gamma \right) \right]$$
Defining
$$h = -\frac{\partial f}{\partial q} + \frac{\partial f}{\partial \sigma_{ij}}E_{ijkl} \left(r_{kl} - \gamma \frac{\partial r_{kl}}{\partial q} \right) \ ,$$
and solving the above expression for $$\delta \gamma$$ gives
$$h\delta\gamma = \frac{\partial f}{\partial \sigma_{ij}} E_{ijkl} \left(\delta \epsilon_{kl} - \gamma \frac{\partial r_{kl}}{\partial \sigma_{mn}}\delta \sigma_{mn} \right)$$
Substituting this into the expression for $$\delta \epsilon^{\mathrm{plastic}}$$ gives
$$\delta \epsilon_{ij}^{\mathrm{plastic}} = \frac{1}{h} \frac{\partial f}{\partial \sigma_{pq}}E_{pqkl} \left(\delta \epsilon_{kl} - \gamma \frac{\partial r_{kl}}{\partial \sigma_{mn}}\delta \sigma_{mn} \right) \left( r_{ij} + \gamma \frac{\partial r_{ij}}{\partial q} \right) + \gamma \frac{\partial r_{ij}}{\partial \sigma_{kl}} \delta \sigma_{kl}$$
This may be substituted into the expression $$\delta \sigma_{ij} = E_{ijkl} (\delta \epsilon_{ij} - \delta \epsilon_{ij}^{\mathrm{plastic}})$$ to finally yield an expression that contains only $$\delta \sigma$$ and $$\delta \epsilon$$, and thus $$H_{ijkl}$$ can be derived.

Drucker Prager plasticity is often used with isotropic elasticity.  In this case the above expressions can be greatly simplified.  Suppose the elasticity tensor is isotropic:
$$E_{ijkl} = \lambda \delta_{ij}\delta_{kl} + \mu (\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})$$