Often in geological settings, rocks are deposited in layers, and between the layers there are joints, which are usually weak, and which we call "weak planes" here. Upon applying stresses/strains, these joints often fail by opening (Weak-plane tensile failure), or by shearing (Weak-plane shear failure). Here we describe Weak-plane shear plasticity, which is non-associative and includes hardening/softening.
This article uses notation defined in Plasticity and the Return-Map Algorithm
Define a vector, $$$n$$$, which is normal to the weak plane. The user may choose two options:
Usually the first choice would be made.
Denote the so-called "cohesion" by $$$C$$$, the "friction angle" by $$$\phi$$$, and the "dilation angle" by $$$\psi$$$. These may be constant, or governed by a Hardening or Softening law
At the start of each "time-step" (increment), rotate the quad-point's reference frame so that $$$n$$$ aligns with the $$$z$$$ axis. Denote the rotated stress by $$$\tilde{\sigma}$$$. Then the shear yield function is $$$$f(\phi) = \sqrt{\tilde{\sigma}_{zx}^{2} + \tilde{\sigma}_{zy}^{2}} + \tilde{\sigma}_{zz}\tan\phi - C \ ,$$$$ and the flow potential is $$$$r_{ij} = \frac{\partial f(\psi)}{\partial \tilde{\sigma}_{ij}} \ .$$$$ Notice the appearance of the dilation angle $$$\psi$$$ in this formula. Associative plasticity is obtained for $$$\psi=\phi$$$. The hardening law is defined via $$$$h = -1 \ ,$$$$ so that $$$q = q^{\mathrm{old}} + \gamma$$$, ie, $$$q$$$ is the accumulated sum of $$$\gamma$$$ over all "time steps".
The yield surface is shown in the figure below. Concentrating on just the $$$\sigma_{zx}$$$, $$$\sigma_{zy}$$$ and $$$\sigma_{zz}$$$ components of stress, it is a cone. Often the yield surface is represented in 2D as a line by just plotting the radial component, $$$$\tau = \sqrt{\tilde{\sigma}_{zx}^{2} + \tilde{\sigma}_{zy}^{2}} \ ,$$$$ which is called the shear stress. Such a diagram is also shown in the figure below. Finally, the figure also shows a typical return-map. The angle of return in stress space is related to $$$\psi$$$. It is not equal to $$$\psi$$$ exactly, since the elasticity tensor, $$$E_{ijkl}$$$ is also involved in determining the direction. However, if $$$\psi=0$$$, the return is directly "downwards" in this figure (ie, in a radial direction).
There are a number of subtleties that need to be discussed:
There is a problem at the tip of the cone. One of the first assumptions of plasticity is that the yield surface is smooth, and that is not the case here. There are several ways around this, and the MOOSE implementation uses a smoothing routine. Define a "smoother" parameter, $$$a$$$, with dimensions of stress. Then the yield function used by MOOSE is $$$$f(\phi) = \sqrt{\tilde{\sigma}_{xz}^{2} + \tilde{\sigma}_{yz}^{2} + a^2} + \tilde{\sigma}_{zz}\tan\phi - C \ ,$$$$ (and still $$$r = \partial f(\psi)/\partial \tilde{\sigma}$$$). The cone's tip gets rounded, and $$$\tilde{\sigma}_{zz} = (C-a)/\tan\phi$$$ at zero shear stress. This is shown in a figure below. In practical situations the choice of $$$a$$$ is largely irrelevant, as the weak-plane shear yield surface is used in conjunction with a weak-plane tensile yield surface (multi-surface plasticity) with tensile strength substantially less than $$$(C-a)/\tan\phi$$$. I suggest choosing $$$a = 0.1C$$$ (unless $$$C=0$$$, of course). There are two types of smoothing currently implemented in MOOSE, and these are activated by setting the 'tip_scheme' parameter. These are: (i) 'hyperbolic' smoothing. Here $$$a$$$ is a constant entered by the user. (ii) 'cap' smoothing. Here $$$$a^2 = \epsilon^2 + p(\tilde{\sigma}_{zz} - \sigma_{\mathrm{start}})^2$$$$ The function $$$p$$$ is $$$$p(x) = \left\{ \begin{array}{ll} 0 & \ \ \ \mbox{for}\ \ \ x\leq 0 \\ x(1-\exp(-rx)) & \ \ \ \mbox{for}\ \ \ x>0 \end{array} \right.$$$$ The three user-entered parameters are $$$\epsilon$$$, $$$\sigma_{\mathrm{start}}$$$, and $$$r$$$. This type of smoothing degenerates to the hyperbolic case for $$$\tilde{\sigma}_{zz}\leq \sigma_{\mathrm{start}}$$$, but for large $$$\sigma_{zz}$$$ (compared with $$$1/r$$$), it produces a hemispherical cap. This is advantageous in the return-map algorithm if the trial stresses have large $$$\sigma_{zz}$$$. Therefore I recommend using this smoothing unless there are physical reasons against it.
As with any plasticity, an arbitrary elasticity tensor may cause problems. One of the fundamental assumptions of the return-mapping procedure is that $$$$\frac{\partial f}{\partial \sigma_{ij}}E_{ijkl}r_{kl} > 0 \ .$$$$ Physically, $$$-E_{ijkl}r_{kl}$$$ is the return direction in stress space, so physically this inequality translates as the condition that $$$f$$$ should decrease in that direction. Obvious problems arise with the returning concept described here if that is not the case. In this case, the conditions are $$$$0 = E_{xzzz} = E_{yzzz}$$$$ and $$$$\left( \begin{array} E E_{xzxz} & E_{xzyz} \\ E_{xzyz} & E_{yzyz} \end{array} \right)$$$$ is not negative definite. These conditions hold for all materials that I'm familiar with.
The current MOOSE implementation explicitly only considers the case when the stress is symmetric. Therefore, the yield function is written in the code as $$$$f = \sqrt{\frac{1}{4}(\tilde{\sigma}_{xz}+\tilde{\sigma}_{xz})^{2} + \frac{1}{4}(\tilde{\sigma}_{yz}+\tilde{\sigma}_{zy})^{2} + a^2} + \tilde{\sigma}_{zz}\tan\phi - C \ ,$$$$ For the case of Cosserat mechanics when $$$\sigma$$$ is not symmetric, this will need to be altered.
A number of tests have been implemented demonstrating that the MOOSE code correctly implements the return-map algorithm for weak-plane shear plasticity, both for small strains and large strains. These are:
small_deform1. Denote the so-called "normal stress", $$$\tilde{\sigma}_{zz}$$$ by $$$N$$$ for ease of notation. The parameters in this test are $$$C=1$$$, $$$\tan\phi=1/2$$$, and perfect plasticity is used, so that the cone's tip lies at $$$N=2$$$. Poisson's ratio is zero, and Lame $$$\mu = 10^6$$$. Apply a deformation so the trial stress has $$$\tau^{\mathrm{trial}}=10$$$ and $$$N^{\mathrm{trial}}=2$$$. The return-map must then solve $$$f=0$$$, and $$$$\tau = \tau^{\mathrm{trial}} - \frac{1}{2}\mu\gamma$$$$ and $$$$N = N^{\mathrm{trial}} - \mu\tan(\psi)\gamma$$$$ Finally, the dilation angle is chosen to be $$$\tan\psi = 2/18$$$. Then the returned value should be $$$\tau=1$$$ and $$$N=0$$$.
small_deform2. A sequence of boundary displacements are applied to a single element. MOOSE uses the return-map algorithm to map out the weak-plane shear envelope (perfect plasticity is used). The result is shown to the right, and clearly MOOSE returns correctly to the yield surface.
small_deform3. Ten-thousand random small deformations are applied to a single element using perfect plasticity. It is checked that MOOSE returns to the yield surface each time.
small_deform4. A sequence of boundary displacements are applied to a single element, and 'cap' smoothing is used. MOOSE uses the return-map algorithm to map out the weak-plane shear envelope (perfect plasticity is used). The result is shown to the right, and clearly MOOSE returns correctly to the yield surface.
large_deform1. The normal direction is chosen $$$n=(0,0,1)$$$, and it rotates with the mesh. A single element is rotated through 90 degrees about the $$$x$$$ axis, and then stretched in $$$z$$$ direction. No plastic deformation is observed, in agreement with expectations.
large_deform2. The normal direction is chosen $$$n=(0,0,1)$$$, and it rotates with the mesh. A single element is rotated through 45 degrees about the $$$x$$$ axis, and then stretched in the $$$y=z$$$ direction. This should create a pure tensile load ($$$\tau=0$$$), which should return to the yield surface. The parameters chosen are $$$C=10^6$$$, $$$\tan\phi=1$$$, and $$$a=5\times 10^5$$$ (smoothing parameter), so the tip of the smoothed cone is at $$$5\times 10^5$$$ (perfect plasticity is used). The result of the stretch should be $$$$\sigma_{yy} = \sigma_{yz} = \sigma_{zz} = 2.5\times 10^5 \ ,$$$$ and this is indeed the result obtained by MOOSE.
large_deform3. Ten-thousand random large deformations are applied to a single element with $$$n$$$ rotating with the mesh. Perfect plasticity is used. It is checked that MOOSE returns to the yield surface each time.
large_deform4. Ten-thousand random large deformations are applied to a single element with $$$n$$$ rotating with the mesh. Cap smoothing is used. Perfect plasticity is used. It is checked that MOOSE returns to the yield surface each time.
small_deform_harden1. Repeated stretches in the $$$z$$$ direction are applied to a single element (with $$$n=(0,0,1)$$$). Exponential hardening is used: $$$C_{0}=1000$$$, $$$C_{\mathrm{res}}=2000$$$, and $$$\zeta_{C}=40000$$$, and the smoothing is done with $$$a=500$$$. The successive return-maps to the yield surface induce hardening of the cohesion, which agrees with the expected result as shown in the Figure.
small_deform_harden2. Apply repeated stretches and shears to map out the yield function in a similar way to small_deform2, but this time use exponential hardening/softening. The following parameters are used $$$C_{0}=1000$$$, $$$C_{\mathrm{res}}=700$$$, and $$$\psi_{0}=5$$$, $$$\psi_{\mathrm{res}}=1$$$, and $$$\phi_{0}=45$$$, and $$$\phi_{\mathrm{res}}=30$$$. Smoothing is done with $$$a=500$$$. The Figure shows the initial yield surface, and the fully developed yield surface (ie when $$$\zeta q\gg 1$$$).
small_deform_harden3. This is similar to small_deform_harden1, but with exponential hardening occurring through an evolving friction angle: $$$\phi_{0}=45$$$, $$$\phi_{\mathrm{res}}=30$$$, $$$\zeta_{\phi}=40000$$$. The results are shown in the Figure.
small_deform_harden4. The same as small_deform_harden1, but with cubic hardening.
large_deform_harden3. The idea is the same as large_deform3, but with exponential hardening/softening: $$$C_{0}=1000$$$, $$$C_{\mathrm{res}}=0$$$, $$$\phi_{0}=45$$$, $$$\phi_{\mathrm{res}}=10$$$, $$$\psi_{0}=5$$$, $$$\psi_{\mathrm{res}}=2$$$.