This article uses notation defined in Plasticity and the Return-Map Algorithm
Denote the so-called "tensile strength" of the material by $$$T$$$. It may be constant, or governed by a Hardening or Softening law.
The yield function is $$$$f = \mbox{max}(\sigma_{I},\ \sigma_{II},\ \sigma_{III}) - T \ .$$$$ Here $$$\sigma_{I}$$$, $$$\sigma_{II}$$$ and $$$\sigma_{III}$$$ are the principal stresses (the eigenvalues of $$$\sigma$$$). Therefore this type of plasticity models situations where the material fails under tensile load. This plasticity is associative, so the flow potential is $$$$r_{ij} = \frac{\partial f}{\partial \sigma_{ij}} \ .$$$$ The hardening law is defined through $$$$h = -1 \ .$$$$ so that $$$q=q_{\mathrm{old}}+\gamma$$$, so that $$$q$$$ is the accumulated sum of $$$\gamma$$$ over all "time steps".
In principal stress space, the yield surface is a "triangular pyramid". Sometimes this is called a "triangular cone", or simply a "cone" in the literature, since many other types of yield functions also map out conical-like regions of admissible stress space.
The triangular pyramid is shown in the first Figure, where the shaded region shows the admissible region. In that figure the principal stresses are labelled sI, sII and sIII. These stresses form an orthogonal basis.
Another frequently-used basis involves the three stress invariants, which form a cylindrical-like basis in stress space, and is shown in the second Figure. The "mean stress" $$$\sigma_{m} = \frac{1}{3}\mbox{Tr}\,\sigma$$$, (labelled as sM in the figure) is the axial direction. The "deviatoric shear stress", $$$\bar{\sigma} = \sqrt{\frac{1}{2}s_{ij}s_{ij}}$$$, where the "deviatoric stress" $$$s_{ij} = \sigma_{ij} - \sigma_{m}\delta_{ij}$$$, measures along the radial direction, and is labelled sBar in the Figure. Finally, the "Lode angle" is an angular coordinate, $$$\theta$$$ defined by $$$\sin 3\theta = -\frac{3\sqrt{3}J_{3}}{2\bar{\sigma}^{3}}$$$ with $$$J_{3} = \mbox{det}s$$$, and $$$-\pi/6 \leq \theta \leq \pi/6$$$. These final bounds mean that only the region $$$\sigma_{I} \geq \sigma_{II} \geq \sigma_{III}$$$ is being considered, but the other regions are clearly related.
There are two commonly-used planes related to this latter basis. The "octahedral plane" has $$$\sigma_{m}$$$ constant. The "meridional plane" has $$$\theta$$$ constant. The tensile yield surface is shown on these planes in the third Figure. On the octahedral plane, the yield surface appears as a triangle, with vertices at $$$\theta = 30$$$deg. The smoothing described below, smooths these vertices (which are edges of the triangular pyramid). On the meridional plane the yield surface appears as a line with slope dependent on $$$\theta$$$, and the triangular pyramid has its tip at $$$\sigma_{m}=T$$$. The smoothing described below rounds this tip.
Finally, there are some useful identities relating these two bases. From now on, assume that the principal stresses are ordered in descending order, that is $$$$\sigma_{I} \geq \sigma_{II} \geq \sigma_{III} \ .$$$$ The formulae presented below only hold for this ordering $$$$\begin{array}{rcl} \sigma_{m} & = & \frac{1}{3}(\sigma_{I} + \sigma_{II} + \sigma_{III}) \\ \sigma_{I} & = & \frac{2}{\sqrt{3}}\bar{\sigma}\sin\left(\theta + \frac{2\pi}{3}\right) + \sigma_{m} \\ \sigma_{II} & = & \frac{2}{\sqrt{3}}\bar{\sigma}\sin\theta + \sigma_{m} \\ \sigma_{III} & = & \frac{2}{\sqrt{3}}\bar{\sigma}\sin\left(\theta - \frac{2\pi}{3}\right) + \sigma_{m} \\ \sigma_{I} - \sigma_{m} & = & \bar{\sigma}\left(\cos\theta - \frac{1}{\sqrt{3}}\sin\theta\right) \\ \end{array}$$$$
Although the yield function is particularly simple, the geometrical analysis has shown that the yield surface has a "tip", and the appearance of the maximum principal stress in the yield function has resulted in "edges" on the yield surface. In principal there is nothing wrong with this, but it means that the flow direction defined by $$$r_{ij} = \partial f/\partial \sigma_{ij}$$$ is not defined everywhere in stress space. There are three usual ways of getting around this problems:
Decompose stress space into a number of regions, and define different flow rules in those regions. For instance, a stress configuration with $$$\sigma_{I} = \sigma_{II} = \sigma_{III} > T$$$ would be defined to flow to $$$\sigma_{I} = \sigma_{II} = \sigma_{III} = T$$$. The advantage of this is that it is computationally very efficient for simple linear yield functions like the one considered here.
Use multi-surface plasticity, with three planar surfaces. This has been implemented in MOOSE as the TensorMechanicsPlasticTensileMulti Material model. It is usually faster than the smoothed version (below)when used alone, but when used in conjunction with other plasticity models (Such as J2 or WeakPlane), it is not a lot faster.
"Smooth" (alter) the yield surface so the potential function gives the desired flow direction. The same problem is also encountered in Weak Plane Shear, and a very similar approach is used here to smooth the pyramid's tip, and a well-established method is used to smooth the pyramid's edges.
Using the identities written above, write $$$f = \sigma_{I} - T = \sigma_{m} + \bar{\sigma}K - T$$$, with $$$K(\theta) = \cos\theta - \frac{1}{\sqrt{3}}\sin\theta$$$. The advantage of doing this is that $$$\bar{\sigma}$$$ is the "radial coordinate", in particular $$$\bar{\sigma}>0$$$, and now the yield function looks very similar to the weak-plane shear yield function. The yield function is smoothed using the same hyperbolic technique: $$$$f = \sigma_{m} + \sqrt{\epsilon^{2} + \bar{\sigma}^{2}K^{2}} - T$$$$ with some small parameter $$$\epsilon$$$. This smoothing is discussed comprehensively in OC Zienkiewicz and GN Prande "Some useful forms of isotropic yield surfaces for soil and rock mechanics" (1977), In: G Gudehus (Ed.) "Finite Elements in Geomechanics" Wiley, Chichester, pp 179-190. Because eigenvalues are computationally cheap to calculate via a BLAS function, the MOOSE code uses the identical form, $$$$f = \sigma_{m} + \sqrt{\epsilon^{2} + (\sigma_{I} - \sigma_{m})^{2}} - T \ ,$$$$ in the case where edge smoothing is not used ($$$\epsilon$$$ is a user-entered value). (As an aside: It is interesting reading the old papers where it appears that the $$$(\sigma_{m}, \bar{\sigma})$$$ version of the yield function was often used for computational efficiency. It is quite cute that finding the eigenvalues of a 3x3 symmetric matrix was deemed expensive!)
The edges of the pyramid also need smoothing. This is achieved in MOOSE using a technique described in AJ Abbo, AV Lyamin, SW Sloan, JP Hambleton "A C2 continuous approximation to the Mohr-Coulomb yield surface" International Journal of Solids and Structures 48 (2011) 3001-3010. The function $$$K$$$ mentioned above is replaced with $$$$K(\theta) = \left\{ \begin{array}{ll} A + B\sin 3\theta + C \sin^{2}3\theta & \ \ \mbox{for}\ \ \theta > \theta_{T} \\ \cos\theta - \frac{1}{\sqrt{3}}\sin\theta & \ \ \mbox{otherwise} \end{array} \right.$$$$ Here $$$0 \leq \theta_{T} \leq \pi/6$$$ is a user-input value determining where on the octahedral plane the smoothing should start. This formulation is slightly different than the version in the aforementioned reference which has $$$|\theta|>\theta_{T}$$$, since the smoothing is not needed for $$$\theta < -\theta_{T}$$$. The constants $$$A$$$, $$$B$$$ and $$$C$$$ depend on $$$\theta_{T}$$$ alone and are rather lengthy: they are written in Equations (15), (16) and (17) of the aforementioned reference.
small_deform_1. A single unit element is stretched in the $$$z$$$ direction only, and the stress is observed to return to $$$\sigma_{zz} = T$$$. No smoothing of the pyramid's tip is necessary ($$$\epsilon = 0$$$), and edge smoothing is also not needed, since the return occurs along the direction $$$\sigma_{II} = \sigma_{III} = 0$$$. This tests that the stress returns to the correct point in the most basic of situations.
small_deform_2. A single unit element is stretched in all directions. In this case smoothing of the pyramid's tip must be used. The parameters chosen are $$$T=1$$$ and $$$\epsilon = 0.5$$$. The stress returns to the yield surface correctly at point $$$\sigma_{xx} = \sigma_{yy} = \sigma_{zz} = T - \epsilon$$$. This checks that the stress returns to the pyramid's "tip" (which is rounded) correctly.
small_deform_3. A single element is stretched by some amount $$$E$$$ in the $$$z$$$ and $$$x$$$ directions. This causes the return direction to be along the hyper surface $$$\sigma_{I} = \sigma_{II}$$$. This checks that the stress returns to the pyramid's edge (which is rounded) correctly. The tensile strength is set to $$$T=1$$$, the tip smoothing is unnecessary ($$$\epsilon=0$$$), and the edge smoothing parameter is $$$\theta_{T} = 25$$$degrees. This latter choice means Abbo et al's $$$A + B + C = 0.609965$$$. The trial stress has $$$\sigma_{zz}^{\mathrm{trial}} = \sigma_{xx}^{\mathrm{trial}}$$$, so that $$$\sigma_{m}^{\mathrm{trial}} = 2\sigma_{zz}^{\mathrm{trial}}/3$$$ and $$$\bar{\sigma}^{\mathrm{trial}} = \sigma_{zz}^{\mathrm{trial}}/\sqrt{3}$$$. If this were on the yield surface then $$$(2/3)\sigma_{zz}^{\mathrm{trial}} + K\sigma_{zz}^{\mathrm{trial}}/\sqrt{3} - 1 = 0$$$, meaning $$$\sigma_{zz}^{\mathrm{trial}} = 0.9815$$$. The applied deformation, $$$E$$$ is chosen to marginally exceed this stress, and the return-map algorithm returns to the correct point on the yield surface.
small_deform_4. Similar to small_deform_3, but with a larger deformation, showing that the stress returns to the pyramid's "edge" correctly.
small_deform_5. Repeated stretches are applied to a single unit element in the $$$z$$$ and $$$x$$$ directions. The return always occurs along the direction $$$\sigma_{I} = \sigma_{II}$$$, so edge smoothing is necessary. The repeated stretches allows the yield surface to be mapped out in the meridional plane for Lode angle $$$\theta = 30$$$degrees. The result agrees perfectly with expectations, as shown in the Figure below.
small_deform_6. Repeated stretches are applied to a single unit element in the $$$z$$$ direction only. The return always occurs along the direction $$$\sigma_{II} = \sigma_{III}$$$, so edge smoothing is irrelevant. The repeated stretches allow the yield surface to be mapped out in the meridional plane for Lode angle $$$\theta = -30$$$degrees. The result agrees perfectly with expectations, as shown in the Figure below.
small_deform_7. Similar to small_deform_6, but with 'cap' smoothing. The result agrees perfectly with expectations, as shown in the Figure below.
many_deforms.i This performs 1000 random large deformations of a unit element and checks that the stress returns to the yield surface each time.
small_deform_hard1. This checks the internal-parameter evolution in the case of hardening/softening. A single element is stretched with strain $$$$\epsilon_{zz} = t/10^{6} \ ,$$$$ The Young's modulus is chosen to be 20 MPa, initial tensile strength $$$T=10$$$ Pa, tensile residual=10Pa, tensile strength rate=0. Then the trial stress will be $$$$\sigma^{\mathrm{trial}}_{zz} = 20\,\mbox{Pa} \ .$$$$ This stress is inadmissible, so it must return to $$$\sigma_{zz} = T = 10\,\mbox{Pa}$$$. Moose produces the correct answer.
small_deform_hard2. A single element is stretched in x and z directions with strain $$$$\epsilon_{zz} = t/10^{6} \ ,$$$$ The Young's modulus is chosen to be 20 MPa, initial tensile strength $$$T=10$$$ Pa, tensile residual=10Pa, tensile strength rate=1e6. Then the trial stress will be $$$$\sigma^{\mathrm{trial}}_{xx} = \sigma^{\mathrm{trial}}_{zz} = 20\,\mbox{Pa} \ .$$$$ Therefore, the trial stress is inadmissible, and the following equations must be solved The trial stress has $$$\sigma_{zz}^{\mathrm{trial}} = \sigma_{xx}^{\mathrm{trial}}$$$, so that $$$\sigma_{m}^{\mathrm{trial}} = 2\sigma_{zz}^{\mathrm{trial}}/3$$$ and $$$\bar{\sigma}^{\mathrm{trial}} = \sigma_{zz}^{\mathrm{trial}}/\sqrt{3}$$$. If this were on the yield surface then $$$(2/3)\sigma_{zz}^{\mathrm{trial}} + K\sigma_{zz}^{\mathrm{trial}}/\sqrt{3} - 1 = 0$$$, meaning $$$\sigma_{zz}^{\mathrm{trial}} = 0.9815$$$.
small_deform_hard3. Repeated stretches are made in the $$$z$$$ direction to a single unit element. Cubic hardening is used with $$$T_{0}=1$$$, $$$T_{\mathrm{res}}=0.5$$$, and $$$l_{T}=0.00001$$$. The result is as expected, and shown in the figure below.