Definition
==================
This article uses notation defined in [Plasticity and the Return-Map Algorithm](../ReturnMap).  It also uses ideas and notation described in the [Tensile yield article](../Tensile)

Denote the so-called "cohesion" by $$C$$, the "friction angle" by $$\phi$$, and the "dilation angle" by $$\psi$$.  These may may remain unchanged during plastic deformation ("perfect plasticity"), or may change with a [Hardening or Softening law](../HardSoft).

Often in rock mechanics, rock strength is written in terms of the UCS, which is $$2C\cos\phi/(1-\sin\phi)$$.

The yield function written in terms of the principal stresses, $$\sigma_{I} \geq \sigma_{II} \geq \sigma_{III}$$, is
$$f(\phi) = \mbox{\frac{1}{2}}(\sigma_{I} - \sigma_{III}) + \mbox{\frac{1}{2}}(\sigma_{I} + \sigma_{III})\sin\phi - C\cos\phi$$
For later purposes I have written this as a function of the friction angle, $$\phi$$.  Using identities written in the [Tensile yield article](../Tensile), the yield function may be re-written as
$$f(\phi) = \sigma_{m}\sin\phi + \bar{\sigma}K(\theta, \phi) - c\cos\phi$$
where
$$K(\theta, \phi) = \cos\theta - \frac{1}{\sqrt{3}}\sin\phi\sin\theta \ .$$

The flow potential is
$$r_{ij} = \frac{\partial f(\psi)}{\partial \sigma_{ij}} \ .$$
Note the appearance of the dilation angle $$\psi$$ in this formula, meaning that this plasticity is nonassociative (unless $$\psi = \phi$$).
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".

Geometrical description and some useful identities
==================

In principal stress space, the yield surface is a "hexagonal pyramid".  Sometimes this is called a "hexagonal 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.

[image:57 align:right]
The Mohr-Coulomb yield surface in principal stress space is a hexagonal pyramid.  The inadmissible region lies outside this pyramid.  In this figure, the principal stresses are labelled as sI, sII and sIII.  Also shown is the mean stress, $$\sigma_{m} = (\sigma_{I} + \sigma_{II} + \sigma_{III})/3$$, which is labelled by sM.  The axis of the pyramid lies along this direction.  Also shown is the octahedral plane, which is perpendicular for $$\sigma_{m}$$.  The radial coordinate on this plane is $$\bar{\sigma}$$, which is labelled as sBar in this figure.

[image:58 align:right]
The tensile yield surface on the octahedral plane (with $$\sigma_{m}$$ constant) and the meridional plane (with Lode angle $$\theta$$ constant). On the octahedral plane the yield surface appears as a hexagon: it is a regular hexagon for friction angle $$\phi=0$$, as shown in green, and becomes more triangular for larger $$\phi$$.  Also shown is the smoothing described below. The invariants are $$\sigma_{m}$$ (labelled as sM), $$\bar{\sigma}$$ (labelled as sBar), and the Lode angle, $$\theta$$ is labelled as th.

Smoothing the yield surface
==================

### Introduction

The geometrical analysis has shown that the yield surface has a "tip", and the appearance of the maximum and minimum principal stresss 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:

1. 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} > C\cot\phi$$ would be _defined_ to flow to $$\sigma_{I} = \sigma_{II} = \sigma_{III} = C\cot\phi$$.  The advantage of this is that it is computationally very efficient for simple linear yield functions like the one considered here.  However, I have not yet implemented this approach.

2. Use 6 planar surfaces and MultiSurface plasticity.  When used alone, this is usually faster than the smoothed version, below.  However, when used in conjunction with other plasticity model (e.g. J2, or WeakPlane), it may not be much faster.  Also, because the flow directions change so dramatically around the "tip" and "edges" of the pyramid, yield-function tolerances typically have to be set higher than the smoothed version.  This plasticity model has been implemented as a TensorMechanicsPlasticityMohrCoulombMulti Material class in MOOSE.

2. "Smooth" (alter) the yield surface so the potential function gives the desired flow direction.  The same problem is also encountered in [Tensile plasticity](../Tensile), and a very similar approach is used here to smooth the pyramid's tip and edges.

### Smoothing the pyramid's tip

The yield function is smoothed using the hyperbolic technique described in [Tensile plasticity](../Tensile):
$$f = \sigma_{m}\sin\phi + \sqrt{\epsilon^{2} + \bar{\sigma}^{2}K^{2}} - C\cos\phi$$
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.  The MOOSE code uses the identical form,
$$f(\phi) = \sigma_{m}\sin\phi + \sqrt{\epsilon^{2} + \mbox{\frac{1}{4}}(\sigma_{I} - \sigma_{III} + (\sigma_{I} + \sigma_{III} - 2\sigma_{m})\sin\phi)^{2}} - C\cos\phi \ ,$$
in the case where edge smoothing is not used ($$\epsilon$$ is a user-entered value).  This is because eigenvalues and their derivatives are computationally cheap to calculate via a BLAS function, and are tidy compared with computations involving the Lode angle.

### Smoothing the pyramid's edges

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\phi\sin\theta & \ \ \mbox{otherwise} \end{array} \right.$$
Here $$\theta_{T}$$ is a user-input value determining where on the octahedral plane the smoothing should start.   This parameter must satisfy $$\theta_{\mathrm{lower}} < \theta_{T} \leq \pi/6$$.  The lower bound, $$\theta_{\mathrm{lower}}>0$$, exists through demanding convexity, and is a complicated function of $$\phi$$, which is written in Equation (C.18) of the aforementioned paper.  The MOOSE code throws an exception if the user inputs a value which violates the bound.  (Usually this will not be the case as for most materials $$\phi < 45$$deg, so that $$\theta_{\mathrm{lower}} < 7$$deg, and users would usually choose $$\theta_{T} = 25$$deg.)  The constants $$A$$, $$B$$ and $$C$$ depend on $$\theta_{T}$$ and $$\phi$$ and are rather lengthy: they are written in Equations (15), (16) and (17) of the aforementioned reference.

The test suite
==================

1. Exception tests.  The nontrivial tests are that $$\theta_{T} \geq \theta_{\mathrm{lower}}$$.

2. small\_deform\_1.  An identical stretch is applied to a single unit element in the $$x$$, $$y$$ and $$z$$ directions.  With $$C=10$$, $$\phi=60$$degrees, and tip\_smoother = 4, then $$\sigma_{m} = (10\cos 60 - 4)/\sin 60 = 1.1547$$.  MOOSE yields the correct answer.

3. small\_deform\_2.  Repeated stretches in the each direction are applied to a single unit element, with $$\sigma_{m}=0$$.  This maps out the yield surface in the octahedral plane for zero mean stress.

4. small\_deform\_3.  Repeated stretches are applied to a single unit element in the $$x$$ and $$z$$ directions, and smaller stretches along the $$y$$ directions, so that $$\sigma_{I} = \sigma_{II}$$, which means the Lode angle is 30deg.  (The Poisson's ratio is set to zero for ease of calculation.)  This allows the yield surface in the meridional plane to be mapped out, and the result is shown in the Figure below.

5.  small\_deform\_4.  Repeated stretches are applied to a single unit element in the $$z$$ direction, and smaller stretches in the $$x$$ and $$y$$ directions, so that $$\sigma_{II} = \sigma_{III}$$, which means that the Lode angle is $$-30\,\mbox{deg}$$.  This allows the yield surface in the meridional plane to be mapped out, and the result is shown in the Figure below.

6. many\_deforms.  1000 random large deformations are applied to a single unit element, and it is checked that the algorithm returns to the yield surface.

7. small\_deform\_hard1.  Identical stretches are applied to a single unit element in the $$x$$, $$y$$, and $$z$$ directions.  The following parameters are chosen: $$C_{0}=10$$, $$C_{\mathrm{res}}=2$$, $$\zeta_{C}=10000$$, $$\phi=60$$deg, tip\_smoother $$=4$$.  MOOSE returns correctly to $$\sigma_{m} = (C\cos\phi - 4)/\sin\phi$$, and this allows the hardening relationship for $$C$$ to be checked, and the results are shown in the Figure below.

7. small\_deform\_hard_cubic.  This is identical to small\_deform\_hard1, except that cubic hardening is chosen, and $$l_{c}=0.0003$$.  The results are shown in the Figure below.

8. small\_deform\_hard2.  Identical stretches are applied to a single unit element in the $$x$$, $$y$$, and $$z$$ directions.  The following parameters are chosen: $$\phi_{0}=60$$deg, $$\phi_{\mathrm{res}}=10$$deg, and $$\zeta_{\phi}=5000$$.  With tip\_smoother denoted by $$T$$, algorithm should return to $$\sigma_{m} = (C\cos\phi - T)/\sin\phi$$, or when $$T=C$$, then $$\phi = 2\pi n - 2\arctan(\sigma_{m}/C)$$.  This allows the hardening relationship for $$\phi$$ to be checked, and the results are shown in the Figure below.

9. small\_deform\_hard3.  Repeated stretches in the $$x$$ and $$z$$ directions are applied to a single unit element, along with smaller stretches in the $$y$$ direction, so that $$\sigma_{I}=\sigma_{II}$$.  This means lode angle $$=30$$deg.  With large $$\zeta$$ values for the friction angle and cohesion, the hardened yield surface can be mapped out in the meridional plane.  The result is in accordance with expectations as shown in the Figure below.

10. uni_axial1.  A single unit element is considered.   The $$x$$ displacements are fixed to zero on the $$x=x_{\mathrm{min}}$$ boundary.  The $$y$$ displacements are fixed to zero on the $$y=y_{\mathrm{min}}$$ boundary.  The $$z$$ displacements are fixed to zero on the $$z=z_{\mathrm{min}}$$ boundary.  The $$z=z_{\mathrm{max}}$$ boundary is moved slowly towards the $$z_{\mathrm{min}}$$ boundary.  Therefore, this simulates a uniaxial compression test, with only nonzero stress component being $$\sigma_{zz}$$, which is negative.  This means $$\sigma_{m}=\sigma_{zz}/3$$, and $$\bar{\sigma}=|\sigma_{zz}|/\sqrt{3}$$.  Also, since $$\sigma_{zz}<0$$, the return-map occurs on the hypersurface $$\sigma_{I}=\sigma_{II}$$, that is, the Lode angle is $$\theta=30$$deg.  The material fails plastically when $$\frac{1}{3}\sigma_{zz}\sin\phi + \frac{1}{\sqrt{3}}K|\sigma_{zz}| - C\cos\phi = 0$$.  Using hardening for the friction angle allows this relationship to be checked.  The parameters chosen are: Young's modulus 100GPa, Poisson's ration 0.3; $$C=10$$MPa; $$\phi_{0}=0$$; $$\phi_{\mathrm{res}}=40$$deg; $$\zeta_{\phi}=10^4$$; and $$\theta_{T}=25$$deg.  The MOOSE results agree perfectly with expectations, as shown in the Figure below.

11. uni_axial2.  A classic uniaxial compression test is performed to observe shear-band formation.  A 3D sample is supported with roller boundary conditions on 3 faces, and is squashed from the top with an applied deformation.  This means that the only nonzero stress will be $$\sigma_{yy}$$, and that will be negative.  The mesh contains a notch that initiates the shear band.  It may be shown for this situation that the sample will fail along a plane whose normal is at angle $$\pi/4 + \phi/2$$ to the $$y$$ direction.  This is a difficult test to perform numerically due to well-known mesh dependencies, and results of such tests are rarely shown in the literature.  For instance, a mesh composed of cubical elements typically causes the shear band to form at a 45 degree angle rather than the desired $$\pi/4 + \phi/2$$ angle.  Therefore, I use a carefully constructed mesh that contains elements oriented in a radial pattern around the notch.  The mesh is shown in the Figures below.  The results for friction angle $$\phi = 2$$deg are shown in the Figure, and the shear-band orientation agrees well with the expected result.  This is a "heavy" test so is not run automatically every time the test suite is run.

12. uni_axial3.  This is an identical test to uni_axial2, except that friction angle $$\phi=40$$deg is chosen.  The results are shown in the Figure, and the shear-band orientation agrees well with the expected result.  This is a "heavy" test so is not run automatically every time the test suite is run.

[image:47 align:right]
Yield surface on the octahedral plane.  Smoothing is evident for $$|\theta| > \theta_{T} = 25$$deg.  Results from small\_deform\_2.

[image:48 align:right]
Yield surface on the meridional plane, for $$\theta=\pm 30$$deg.  Tip smoothing is evident.  Results from small\_deform\_3, and small\_deform\_4.

[image:49 align:right]
Hardening of the cohesion with exponential hardening.  Results from small\_deform\_hard1.

[image:66 align:right]
Hardening of the cohesion with cubic hardening.  Results from small\_deform\_hard\_cubic.

[image:50 align:right]
Hardening of the friction angle with exponential hardening.  Results from small\_deform\_hard2.

[image:51 align:right]
Comparison of the initial yield surface, and the fully-hardened yield surface on the meridional plane.  Results from small\_deform\_hard3.

[image:52 align:right]
Compressive stress at yield as a function of friction angle for an unconfined uniaxial compressive test.  Results from uni_axial1.

[image:55 align:right]
Mesh used in the uni_axial2 and uni_axial3 tests.  Roller boundary conditions are applied as shown, as well as on the $$z=z_{\mathrm{max}}$$ face, and the top surface is moved slowly downwards.  The notch initiaties fracture.  The mesh is specially designed to study the angle of the resulting shear band.

[image:56 align:right]
Close-up view of the notch for the mesh used in the uni_axial2 and uni_axial3 tests.

[image:53 align:right]
Shear band formation in a uniaxial compression test.  Friction angle = 2 degrees.  Results from uni_axial2.

[image:54 align:right]
Shear band formation in a uniaxial compression test.  Friction angle = 40 degrees.  Results from uni_axial3