# 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.

# 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})