Elasto-Perfect Plasticity with Mohr-Coulomb Failure

This is a continuation of the previous post on Mohr-Coulomb Failure. The references and nomenclature there apply to this post as well.

With the Mohr-Coulomb failure criterion defined, the elasto-plastic constitutive matrix can be developed. Since pure plasticity is assumed, thus once failure is reached and F=0 , the stresses can “rearrange” themselves to remain at the failure surface but cannot move from it unless the strains are reduced to the point where the stress state is within the failure surface and F<0. This is easier to see if we look at the diagram above. We should also note that, for this presentation, the plane strain/axisymmetric case will be considered; the derivation would proceed a little differently for other cases such as plane stress or three-dimensional stress.

To do this, both Equations 8 and 10 from the previous post are differentiated. Considering that

$\sigma=\left[\begin{array}{c} \sigma_{x}\\ \sigma_{y}\\ \sigma_{z}\\ \tau_{xy} \end{array}\right]$ (1)

for the two-dimensional model, then

$\frac{\partial F}{\partial\sigma}=a_{F}=C_{1}a_{1}+C_{2}a_{2}+C_{3}a_{3}$ (2)

where

$C_{1}=\frac{sin\phi}{3}$
$C_{2}=cos\theta\left[\left(1+tan\theta tan3\theta\right)+\frac{sin\phi\left(tan3\theta-tan\theta\right)}{\sqrt{3}}\right]$ (3)
$C_{3}=\frac{\sqrt{3}sin\theta+cos\theta sin\phi}{2J'_{2}cos3\theta}$

and

$a_{1}=\left[\begin{array}{c} 1\\ 1\\ 1\\ 0 \end{array}\right]$
$a_{2}=\frac{1}{2\sqrt{J'_{2}}}\left[\begin{array}{c} \sigma'_{x}\\ \sigma'_{y}\\ \sigma'_{z}\\ 2\tau_{xy} \end{array}\right]$ (4)
$a_{3}=\left[\begin{array}{c} \sigma'_{y}\sigma'_{z}+\frac{J'_{2}}{3}\\ \sigma'_{x}\sigma'_{z}+\frac{J'_{2}}{3}\\ \sigma'_{y}\sigma'_{x}-\tau_{xy}^{2}+\frac{J'_{2}}{3}\\ -2\sigma'_{z}\tau_{xy} \end{array}\right]$

Inspection of Equation 3 will show that C₂ and C₃ are undefined for values of θ =30°. This is the well known “corner” problem which has occupied the literature for many years, as summarized by Abbo et.al. (2011). For this study the “corner cutting” of Owen and Hinton (1980) is used, and

$C_{1}=\frac{sin\phi}{3}$
$C_{2}=\frac{1}{2}\left[\sqrt{3}-\frac{sin\phi}{\sqrt{3}}\right],\:\theta\approx\pm30^{\circ}$ (5)
$C_{3}=0$

Although the accuracy of this “corner cutting” can be improved with methods such as those of Abbo et.al. (2011), they come at more computational expense.

In like fashion for the plastic potential function,

$\frac{\partial Q}{\partial\sigma}=a_{Q}=C_{4}a_{1}+C_{5}a_{2}+C_{6}a_{3}$ (6)

where

$C_{4}=\frac{sin\psi}{3}\\C_{5}=cos\theta\left[\left(1+tan\theta tan3\theta\right)+\frac{sin\psi\left(tan3\theta-tan\theta\right)}{\sqrt{3}}\right]\\C_{6}=\frac{\sqrt{3}sin\theta+cos\theta sin\psi}{2J'_{2}cos3\theta}$ (7)

and by extension

$C_{4}=\frac{sin\psi}{3}\\C_{5}=\frac{1}{2}\left[\sqrt{3}-\frac{sin\psi}{\sqrt{3}}\right],\:\theta\approx\pm30^{\circ}\\C_{6}=0$ (8)

and the vectors a₁ , a₂ , a₃ are the same as with Equation 2.

At this point we need to pause and consider how we plan to present our result. It’s certain possible, for example, to multiply through Equations 2, 3, and 4 and use this result in subsequent derivations. This approach suffers from two problems: the element expressions in the matrices and vectors become very complicated, and the computation cost is increased by the repetitive calculations. For some good advice on this subject–for this and other problems–we would suggest that you review this.

Now the elasto-plastic constitutive matrix can be considered. Here the difference between path indepedence and path dependence becomes significant. In the elastic case, we can use the elastic consitutive matrix D^e and, knowing the strains, compute the stresses, and perform the inverse as well. With plasticity, it is necessary to compute the strains and stresses incrementally, moving from one stress and strain state to the next and adding the results to what has gone before. We do this with both elastic, plastic and combined strains and stresses for computational consistency. The elastic results will be the same except for the additional accumulation of numerical error due to the larger number of computations.

For any strain increment partially or totally beyond the yield surface, that strain increment will contain both elastic and plastic portions (Griffiths and Willson (1986)), or

$d\epsilon=d\epsilon^{e}+d\epsilon^{p}$ (9)

Strain increments occur normal to the plastic potential surface, thus

$d\epsilon^{p}=\lambda_{p}a_{Q}$ (10)

For elastic materials, the incremental stress-strain relationship is simply

$d\sigma=D^{e}d\epsilon^{e}$ (11)

where, for plane strain and axisymmetric problems, as before,

$D^{e}=\left[\begin{array}{cccc} {\frac{\left(1-\nu\right)\lambda}{\nu}} & \lambda & \lambda & 0\\ \noalign{\medskip}\lambda & {\frac{\left(1-\nu\right)\lambda}{\nu}} & \lambda & 0\\ \noalign{\medskip}\lambda & \lambda & {\frac{\left(1-\nu\right)\lambda}{\nu}} & 0\\ \noalign{\medskip}0 & 0 & 0 & G \end{array}\right]$ (12)

and also, for axisymmetric problems only,

$\epsilon^{e}=\left[\begin{array}{c} \epsilon_{r}\\ \epsilon_{z}\\ \epsilon_{\theta}\\ \epsilon_{rz} \end{array}\right]$ (13)

For a perfectly elasto-plastic material, i.e., one without hardening or softening, stress changes will take place only during elastic action. Thus in these cases Equations 9, 10 and 11 can be rearranged to yield

$d\sigma=D^{e}\left(d\epsilon-\lambda_{p}a_{Q}\right)$ (14)

Since stresses on the failure surface can “rearrange” themselves without leaving same surface,

$a_{F}^{t}d\sigma=0$ (15)

Combining and rearranging Equations 9 through 15 yields the following relationships:

$\lambda_{p}=\frac{a_{F}^{t}D^{e}d\epsilon}{a_{F}^{t}D^{e}a_{Q}}$ (16)

$d\sigma=\left(D^{e}-\frac{D^{e}a_{Q}a_{F}^{t}D^{e}}{a_{F}^{t}D^{e}a_{Q}}\right)d\epsilon$ (17)

$d\sigma=\left(D^{e}-D^{p}\right)d\epsilon$ (18)

where

$D^{p}=\frac{D^{e}a_{Q}a_{F}^{t}D^{e}}{a_{F}^{t}D^{e}a_{Q}}$ (19)

and finally

$D^{ep}=D^{e}-D^{p}$ (20)

We need to break this down for greater understanding. Let us first define

$a_{\beta}=\frac{1}{a_{F}^{t}D^{e}a_{Q}}$ (21)

This is a scalar quantity. (Admittedly sometimes it is difficult to follow matrix equations because they combine scalar and vector/matrix quantities.) We can thus rewrite Equation 19

$D^{p}=a_{\beta}\left(D^{e}a_{Q}a_{F}^{t}\right)D^{e}$ (22)

and Equation 20 can be likewise rewritten as

$D^{ep}=\left(I-a_{\beta}\left(D^{e}a_{Q}a_{F}^{t}\right)\right)D^{e}$ (23)

where $I =$ the identity matrix.

If plasticity has not yet been initiated, the quantity $a_{\beta}\left(D^{e}a_{Q}a_{F}^{t}\right)$ (which is a 4 x 4 matrix) is null, and the elasticity matrix governs. Once it is initiated the magnitude of $a_{\beta}\left(D^{e}a_{Q}a_{F}^{t}\right)$ begins to increase and D^ep begins to decrease, which indicates a ‘softening’ of the material. This is what we would expect with plasticity. The whole topic of matrix magnitude is not a straightforward one and is usually viewed in relationship to matrix norms, which are discussed in Gourdin and Boumharat (1989).

In any case

$d\sigma=D^{ep}d\epsilon=\left(D^{e}-D^{p}\right)d\epsilon=\left(I-a_{\beta}\left(D^{e}a_{Q}a_{F}^{t}\right)\right)D^{e}d\epsilon$ (24)

The first thing to note about Equation 20 is that, if same is used to reconstruct a tangent stiffness matrix K_T in a true Newton stepping scheme, same K_T will not be symmetric if φ ≠ ψ. For many engineering materials, and especially those where both of these quantities are zero, this is not an issue; the flow rule is associated, D^ep is symmetric and K_T will be also. “Regular” engineering materials are used for, say, the steel or concrete portions of a system, but for many geotechincal applications interest in plastic deformation of these components is limited. It is also not an issue with purely cohesive soils. For situations with cohesionless soils, this is not the case; generally φ ≫ ψ and a non-associated flow rule is necessary. This aspect will be important in many decisions regarding the structure and types of schemes used in the model. An overview of the values of the dilitancy angle ψ can be found in Warrington (2019).

Turning to the inclusion of plasticity itself, it is certainly possible to compute the plastic stresses from Equation 20, and also possible to explicitly derive Equation 19 (Griffiths and Willson (1986)). However, it is not always optimal to do so, either from the standpoint of a workable algorithm or from a computational efficiency standpoint. To compute the final stress state in a load or time step where failure takes place, some type of iteration or multiple steps are required. Potts and Zdravkovic (1999) state that there are two basic types of algorithms to accomplish this: substepping (such as Sloan (1987)) and return (Ortiz and Simo (1986)). For STADYN a return algorithm was chosen, and implemented as follows:

For a load or time step, the estimated incremental strain was computed.
The estimated incremental stress was computed, based on the assumption that the strains were still in the elastic region.
The resulting incremental stresses were added to the stresses at the beginning of the step. If these were elastic, then plasticity was not considered and the incremental stresses were added to the original ones. If they were not, then the plasticity routine was invoked.
The plasticity routine began by computing a_F , a_Q , F , θ and a_β.
The plasticity constant λ_p was determined by a modification to Equation 16, namely
$\lambda_p = F \alpha_\beta$ (25)
The incremental strains in the return step were computed by the equation (see Equation 14)
$d\epsilon_{r}=\lambda a_{Q}$ (26)
The return incremental stresses were thus
$d\sigma_{r}=D^{e}d\epsilon_{r}$ (27)
Both of these were subtracted from the current estimated strain and stress, thus
$d\epsilon_{ep}^{n+1}=d\epsilon_{ep}^{n}-d\epsilon_{r}\\d\sigma_{ep}^{n+1}=d\sigma_{ep}^{n}-d\sigma_{r}$ (28)
The current stress state was checked against the previous stress state. If the norm of the vector difference of the two stress states was within the convergence tolerance, the iteration was stopped and the computed elasto-plastic stresses and strains were accepted. If not, the cycle was repeated.

Potts and Zdravkovic (1999) criticize the return algorithms because they obtain a result based on information and computations in illegal stress space. However, overall the experience in this study is that the return algorithm worked well, with most of the convergence to the new failure surface stress state taking place in the first iteration and the rest refinement steps. A fair way of differentiating between the two is that, with substepping, one starts with the existing stress state and works one’s way to the failure surface, while the return algorithm starts by overshooting the failure surface and then coming back to it. The goal in both cases is to return to the failure surface, and in principle the result should be the same. Another way of differentiating between the two is that substepping is, with its avoidance of illegal stresses, more of an engineering type of approach to the problem, while the return algorithm is a more strictly mathematical method of arriving at a solution.