Many books provide detailed introductions to plasticity. Two such books are:
JC Simo and TJR Hughes "Computational Inelasticity" Springer 1997. This book has the advantages that it is available online to subscribing institutions, and also describes some aspects of multi-surface plasticity which is discussed below.
I Doghri "Mechanics of Deformable Solids. Linear and nonlinear, analytical and computational aspects" Springer 2000. Provides some extra insights into computational aspects of plasticity
The return-mapping algorithm is central to simulations involving plasticity, taking much of the computational time and lines of code. It is discussed in detail in the aforementioned references. The point of the algorithm is: at each quadrature point given the stress and other internal parameters, and a specified strain increment, determine the final values of stress and the internal parameters.
If the specified strain increment causes an elastic response, then the problem is trivially solved by performing an elastic calculation to find the final stress. This often happens, for instance, if the strain is decreasing from an initially high value (so the strain increment is negative), for this usually corresponds to letting an initially-stressed system relax back to equilibrium.
However, if the specified strain increment causes a plastic response, the return-mapping algorithm is used to determine the final stress (and internal parameters). This involves forming a so-called "trial stress" that is actually inadmissible (impossible to obtain in a real physical situation), and using a Newton-Raphson procedure to let this trial stress "flow" back to the admissible region. The sections below describe this process mathematically.
This is the most general case of plasticity, and in paragraphs below the more common specializations to single-surface, associative and perfect plasticity are also mentioned
Denote the components of the stress tensor by $$$\sigma_{ij}$$$, with $$$i=0,\ldots, d-1$$$ in $$$d$$$ dimensions. (One can also use $$$i=1,\ldots, d$$$ if preferred). Denote the $$$Q$$$ internal variables by $$$q_{a}$$$, with $$$a=1,\ldots, Q$$$. Physically these are "hardening parameters" that allow the yield surface to evolve with time. For instance, in J2 plasticity described below, the equivalent plastic strain is a hardening parameter. In multi-surface plasticity, there are $$$N$$$ yield functions, $$$f_{\alpha}$$$, ($$$\alpha = 1,\ldots,N$$$) that describe the admissible region: $$$$E_{\sigma,\ q} = \{ (\sigma, q)\ |\ f_{\alpha}(\sigma, q) \leq 0 \ \ \ \forall \alpha\}$$$$ It is assumed this describes a convex region. This assumption is pretty trivial for the cases below when $$$N=1$$$, but for the multi-surface case this should be explicitly checked. It is also assumed that the $$$f_{\alpha}$$$ are smooth functions, which is actually not true for many common yield functions (such as Mohr-Coulomb), so in such cases a smoothing technique must be used to modify the yield function, or some rules created which govern the behavior around the non-smooth regions. However, the $$$f_{\alpha}$$$ may intersect one another in a non-smooth way.
Denote the components of the total strain tensor by $$$\epsilon_{ij}$$$. It is related to the deformation $$$u$$$, via $$$$\epsilon_{ij} = \frac{1}{2}(\nabla_{i}u_{j} + \nabla_{j}u_{i}) \ ,$$$$ (at least for non-Cosserat and small strains). As explained in the above references, it may be split into its elastic and plastic parts: $$$$\epsilon_{ij} = \epsilon^{\mathrm{e}}_{ij} + \epsilon^{\mathrm{p}}_{ij}$$$$ where $$$\sigma_{ij} = E_{ijkl}\epsilon^{\mathrm{e}}_{ij}$$$, and $$$E_{ijkl}$$$ is the elasticity tensor. Upon plastic loading, where the total strain $$$\epsilon$$$ is prescribed, our goal is to find the plastic part, $$$\epsilon^{\mathrm{p}}$$$ as well as the hardening parameters.
An increment is a change of a parameter under a change in "time". I have put "time" in quotes because often in solid mechanics a sequence of quasi-static approximations are used to approximate a sequence of loading scenarios, so the flow of true time is rather obscured. Let us just introduce the notation $$$$\dot{p} \equiv p - p_{\mathrm{old}} \ ,$$$$ for any variable, $$$p$$$. Here $$$p_{\mathrm{old}}$$$ is the variable's value at the end of the previous "time" step, and we would typically like to determine $$$p$$$, which is the variable's value at the end of the current "time" step.
To this end, the following quantities may be introduced:
$$$r^{\alpha}_{ij}(\sigma, q)$$$ are $$$\alpha=1,\ldots,N$$$ "flow potentials".
$$$h^{\alpha}_{a}(\sigma, q)$$$ are $$$\alpha=1,\ldots,N$$$ "hardening potentials".
Then, given an increment in displacement, $$$\dot{u}$$$, which naturally implies an increment in total strain, $$$\dot{\epsilon}_{ij}$$$, the increment in plastic strain is defined to be $$$$\dot{\epsilon}^{\mathrm{p}}_{ij} = \sum_{\alpha=1}^{N}\gamma^{\alpha}r^{\alpha}_{ij} \ . \label{plastic.strain.flow.eqn}$$$$ This is the flow rule for plastic strain. Similarly, the increment in the internal variables is $$$$\dot{q}_{a} = -\sum_{\alpha=1}^{N}\gamma^{\alpha}h^{\alpha}_{a} \ .$$$$ These are the flow rules for the internal parameters.
In these flow rules, the so-called "consistency parameters" $$$\gamma^{\alpha}$$$ ($$$\alpha = 1,\ldots,N$$$) have been introduced, and they obey the Kuhn-Tucker and consistency relationships $$$\gamma^{\alpha} \geq 0 \ \ \mbox{and}\ \ \gamma^{\alpha} f_{\alpha} = 0 \ \ \mbox{and, if$$$f_{\alpha}=0$$$, then }\ \ \gamma^{\alpha} \dot{f}_{\alpha} = 0 \ .$$$ In these expressions there is no sum over $$$\alpha$$$.
At the "old" time, the stress and internal variables are always within the admissible region, or on its boundary. The strain increment is then applied, and a so-called "trial stress" $$$$\sigma^{\mathrm{trial}}_{ij} = E_{ijkl}(\epsilon_{ij} - \epsilon_{ij}^{\mathrm{p,\ old}}) \ ,$$$$ may be formed. Note that if there were no increment in plastic strain for this "time interval" ($$$\epsilon^{\mathrm{p}} = \epsilon^{\mathrm{p,\ old}}$$$) then this trial stress would be the correct stress.
Assume that $$$\dot{q}_{a}=0$$$ for all $$$a$$$. Suppose it is found that $$$$f_{\alpha}(\sigma^{\mathrm{trial}}, q^{\mathrm{old}}) \leq 0$$$$ for all $$$\alpha$$$. Then clearly $$$\sigma = \sigma^{\mathrm{trial}}$$$, and $$$\gamma^{\alpha}=0$$$ is a solution to the Kuhn-Tucker and consistency equations, and the flow rules also hold. This is the trivial case of "plastic unloading" where plastic strain does not increase.
In this section, I will first present the basic theory according to Simo and Hughes, and in a section below I will discuss enhancements of the algorithm to cover certain pathological cases that they choose to ignore.
If assuming $$$\dot{q}_{a}=0$$$ and using $$$q_{a}$$$ and the trial stress it is found that $$$$f_{\alpha}(\sigma^{\mathrm{trial}}, q^{\mathrm{old}}) > 0$$$$ for any $$$\alpha$$$, then the trial stress is inadmissible and so some algorithm must be used to determine the correct values for $$$\sigma$$$ and $$$q$$$. Let us summarize the conditions that must be met (i.e., the equations that must be solved). They are:
The Kuhn-Tucker constraints: $$$\gamma^{\alpha} \geq 0$$$ and $$$f_{\alpha}(\sigma, q) \leq 0$$$ and $$$\gamma^{\alpha}f_{\alpha}(\sigma, q) = 0$$$ (no sum on $$$\alpha$$$). These mean that either: (a) $$$\gamma^{\alpha}=0$$$ and $$$f_{\alpha} \leq 0$$$; or (b) $$$\gamma^{\alpha}>0$$$ and $$$f_{\alpha}=0$$$.
The consistency constraint: if $$$f_{\alpha}(\sigma^{\mathrm{old}}, q^{\mathrm{old}})=0$$$ then $$$\gamma^{\alpha}\dot{f}_{\alpha}(\sigma, q) = 0$$$ (no sum on $$$\alpha$$$).
The plastic flow rule: $$$\dot{\epsilon}^{\mathrm{p}}_{ij} = \sum_{\alpha=1}^{N}\gamma^{\alpha}r^{\alpha}_{ij}(\sigma, q)$$$. It is convenient to write the LHS of this equation as $$$\dot{\epsilon}^{\mathrm{p}} = \epsilon^{\mathrm{p}} - \epsilon^{\mathrm{p,\ old}} = E^{-1}(\sigma^{\mathrm{trial}} - \sigma)$$$.
The internal flow rules: $$$\dot{q}_{a} = -\sum_{\alpha=1}^{N}\gamma^{\alpha}h^{\alpha}_{a}(\sigma, q)$$$.
These equations are solved for $$$\sigma$$$, $$$q_{a}$$$, and $$$\gamma^{\alpha}$$$ using the Newton-Raphson technique, often employing a line-search, and techniques borrowed from linear-programming to handle the inequalities in Kuhn-Tucker constraints. The initial values are $$$\sigma = \sigma^{\mathrm{trial}}$$$ and $$$q_{a} = q_{a}^{\mathrm{old}}$$$ and $$$\gamma^{\alpha}=0$$$. Further discussion is in subsequent sections.
Suppose that $$$f_{\alpha}(\sigma^{\mathrm{old}}, q^{\mathrm{old}}) = 0$$$ for all $$$\alpha$$$. Also, suppose that $$$f_{\alpha}(\sigma^{\mathrm{trial}}, q^{\mathrm{old}}) > 0$$$ for all $$$\alpha$$$. Then the following equations must be solved using the Newton-Raphson approach:
$$$f_{\alpha}(\sigma, q) = 0$$$
$$$\sigma^{\mathrm{trial}}_{ij} - \sigma_{ij} = E_{ijkl}\sum_{\alpha = 1}^{N}\gamma^{\alpha}r^{\alpha}_{kl}(\sigma, q)$$$
$$$q_{a} - q_{a}^{\mathrm{old}} = -\sum_{\alpha=1}^{N}\gamma^{\alpha}h^{\alpha}_{a}(\sigma, q)$$$
The initial values for the Newton-Raphson procedure are $$$\sigma = \sigma^{\mathrm{trial}}$$$ and $$$q_{a} = q_{a}^{\mathrm{old}}$$$ and $$$\gamma^{\alpha}=0$$$. Even in this simple case, a problem becomes apparent. For what if, after finding the solution, $$$\gamma^{\alpha}<0$$$ for some $$$\alpha$$$? This violates the Kuhn-Tucker conditions. This is discussed in the next section.
In MOOSE, multi-surface plasticity is implemented through the FiniteStrainMultiPlasticity Material. This assumes that there is exactly one internal parameter per single-surface plasticity model, and that the functions for single-surface plasticity model depend only on its internal parameter. (Actually, the case of zero internal parameters is also trivially covered when the yield function, etc, only depends on the stress. But there must not be 2 or more internal parameters per single-surface plasticity model.) In this case $$$$Q = N$$$$ and $$$$\begin{array}{rcl} f_{\alpha} & = & f_{\alpha}(\sigma, q_{\alpha}) \\ r^{\alpha}_{ij} & = & r^{\alpha}_{ij}(\sigma, q_{\alpha}) \\ h^{\alpha}_{a} & = & \delta^{\alpha}_{a}h^{\alpha}(\sigma, q_{\alpha}) \end{array}$$$$ where there is no sum on $$$\alpha$$$.
The Newton-Raphson procedure attempts to solve three types of equation:
Denoted by "$$${\mathtt{f}}$$$" in the code. $$$f_{\alpha} = 0$$$ for all active $$$\alpha$$$, up to a tolerance specified by the $$${\mathtt{\mbox{yield_function_tolerance}}}$$$ of the plastic model
Denoted by "$$${\mathtt{epp}}$$$" in the code. $$$0 = -E^{-1}_{ijkl}(\sigma_{kl}^{\mathrm{trial}} - \sigma_{kl}) + \sum_{\mathrm{active}\ \alpha}\gamma^{\alpha}r^{\alpha}_{ij}(\sigma, q)$$$, up to a tolerance specified by $$${\mathtt{\mbox{ep_plastic_tolerance}}}$$$. In the code there is a variable $$${\mathtt{delta\_dp}}$$$ that is shorthand for $$$E^{-1}_{ijkl}(\sigma_{kl}^{\mathrm{trial}} - \sigma_{kl})$$$, or $$$\dot{\epsilon}^{\mathrm{p}}_{ij}$$$.
Denoted by "$$${\mathtt{ic}}$$$" in the code. $$$0 = q_{\alpha} - q_{\alpha}^{\mathrm{old}} + \gamma_{\alpha}h^{\alpha}_{\alpha}$$$, up to a tolerance specified by $$${\mathtt{\mbox{internal_constraint_tolerance}}}$$$ of the plastic model. There is no sum on $$$\alpha$$$ in this expression.
In addition to these constraints, the Kuhn-Tucker and consistency conditions must also be satisfied. If the above constraints are satisfied, then these last conditions amount to: if $$$f_{\alpha}=0$$$ (up to a tolerance), then $$$\gamma^{\alpha}\geq 0$$$; otherwise $$$\gamma^{\alpha}=0$$$.
These conditions, and the line-search scheme are described in detail below.
The FiniteStrainMultiPlasticity Material may also be used for single-surfaces. However, if there are greater than one internal parameter, FiniteStrainPlasticBase must be used instead. This Material is not described in this article.
A paper has been written on this subject.