Mohr-Coulomb plasticity.

The theory and algorithms used to code the classical Mohr-Coulomb plasticity model are primarily taken from [100] and [21] and the references therein.

The Mohr-Coulomb plasticity is non-associative with a piecewise linear yield surface $f$ an a piecewise linear plastic potential $g$. The basic idea is that the shear stress at which plastic flow occurs increases with increasing pressure, similar to contact sliding. This is illustrated in Figure 140 showing the shear stress $\tau$ versus normal stress $\sigma_n$. $\sigma_1 \ge \sigma_2 \ge \sigma_3$ are the principal stresses.

Figure 140: Mohr-Coulomb yield surface

The yield surface satisfies:

$$\tau=c-\sigma_n \tan \varphi,$$ (388)

where $c$ is the cohesion. This is equivalent to (cf. Figure 140)

$$\sigma_1 \left( \frac{1+\sin \varphi }{1 - \sin \varphi } \right) -\sigma_3 -2 c \frac{\cos \varphi }{1 - \sin \varphi }=0,$$ (389)

or

$$k \sigma_1 - \sigma_3 - 2 c \sqrt{k} =0,$$ (390)

where

$$k=\frac{1+\sin \varphi }{1 - \sin \varphi }.$$ (391)

$\varphi$ is called the friction angle. Similarly, the plastic potential is defined by:

$$g = m \sigma_1 - \sigma_3,$$ (392)

where

$$m=\frac{1+\sin \psi }{1 - \sin \psi }.$$ (393)

$\psi$ is called the dilation angle. It describes to what extent the volume of the material changes due to shear motion. For a positive value of $\psi$ the volume increases, as for dense sand. For a negative value the volume decreases, as for loose sand. In the latter case the grains fit better due to the motion. You can easily illustrate this by pouring suger in a bowl. Shaking the bowl the volume will decrease and more sugar will fit.

As mentioned, the yield surface is piecewise linear, so the gradient of the surface is not continuous. This complicates matters. Now, looking at the yield surface and plastic potential it can be observed that it contains only principal stresses. Assuming the material to be elastically isotropic, this allows us to perform all operations in principal stress space, thereby reducing the tensors to vectors. This was first proposed in [20].

Replacing $\sigma_1, \sigma_2, \sigma_3$ by $x,y,z$ for simplicity, the yield surface as defined in Equation (390) defines in three-dimensional space one of six faces of a irregular pyramid with central axis in $(1,1,1)$ direction. The other faces correspond to sectors in which the order of $\sigma_1, \sigma_2, \sigma_3$ is different.

Figure 141: Section plane $\perp$ (1,1,1) through (0,0,0)

Figure 141 shown a view on a plane orthogonal to the $(1,1,1)$ axis and through the origin. The six sectors are bordered by planes going through the x-, y- and z- axes and containing the $(1,1,1)$-axis. The cross section of this plane with the yield surface is an irregular hexagon (dashed line). Since the principal stresses can always be rearranged such that $x \ge y \ge z$ it is sufficient to look at that sector (labeled 1) and the neighboring ones (labeled 6 and 2).

Figure 142: Pressure distribution in the lid-driven cavity

Figure 142 shows the yield surface viewed from the apex and looking in the direction of $(-1,-1,-1)$. In the sectors of interest the normal to the surface $\boldsymbol{a}$ is shown. For sector 1 this is immediately clear from Equation (390) which can be written as $(k,0,-1) \cdot (\sigma_1,\sigma_2,\sigma_3) -2c \sqrt{k} =0$. A vector along the interection lines of the pyramid faces is obtained by taking the cross product of the normal of the neighboring faces,e.g. $\boldsymbol{r_6} = \boldsymbol{a_6} \times \boldsymbol{a_1}$. The apex of the pyramid is obtained by substituting $\sigma_1=\sigma_3$ in Equation (390): it yields

$$s_a := \sigma_1=\sigma_2=\sigma_3=\frac{2c \sqrt{k} }{k-1} = \frac{c}{\tan \varphi }.$$ (394)

Therefore, the equation of the intersection line between section 6 and 1 satisfied $(s_a,s_a,s_a)+\lambda (1,k,k)$. Its intersection with the plane $x+y+z=0$ leads to

$$\lambda=\frac{3 s_a}{1+2k},$$ (395)

yielding the location of point A in Figure 141. For point B one gets

$$\lambda=\frac{3 s_a}{2+k},$$ (396)

which is farther away from the origin since $k \ge 1$. The normal to the plastic potential surface is labeled $\boldsymbol{b}$ and amounts for sector 1 to $(m,0,-1)$. Similar expressions apply to the other sectors.

The governing equations of an elastic material with Mohr-Coulomb plasticity read:

(i)

Elastic stress-strain relations

$$\boldsymbol{ \sigma} = \boldsymbol{D} \cdot \boldsymbol{\epsilon}^e$$ (397)

(ii)

Internal variable relationship

$$c(\epsilon^{\text{peq}})$$ (398)

(iii)

Yield surface

$$f(\boldsymbol{ \sigma},c)=0 \Leftrightarrow \boldsymbol{a} \cdot \boldsymbol{\sigma }- 2 c \sqrt{k}=0$$ (399)

(iv)

Evolution equations

$$\dot{ \boldsymbol{ \epsilon}}^p= \dot{ \gamma} \frac{ \partial g ... ...mbol{ \sigma},c)}{\partial \boldsymbol{ \sigma}} = \dot{ \gamma} \boldsymbol{b}$$ (400)

(v)

Kuhn-Tucker equations

$$\dot{ \gamma} \ge 0, \quad f(\boldsymbol{ \sigma},c) \le 0, \quad \dot{ \gamma} f(\boldsymbol{ \sigma},c) =0$$ (401)

(vi)

Consistency condition

$$\dot{ \gamma} \dot{ f}(\boldsymbol{ \sigma},c)=0.$$ (402)

Notice that in the above equatons $\boldsymbol{\sigma }$ and $\boldsymbol{\epsilon}$ are vectors in principal space.

In order to explain the numerical procedure it is assumed that all variables are known at the end of increment $n$ (corresponding to time $t$) and that their values are sought at the end of increment $n+1$ (corresponding to time $t + \Delta t$). The input to the algorithm is a change on total strain $\boldsymbol{\Delta \epsilon }$.

At first it is assumed that now plasticity occurs in the new increment, i.e.:

$$\boldsymbol{ \epsilon} _{n+1} ^p = \boldsymbol{ \epsilon}_n^p$$ (403)

$$c_{n+1} = c_{n}$$ (404)

$$\gamma _{n+1} = \gamma_n$$ (405)

$$\boldsymbol{ \sigma}_{n+1} = \boldsymbol{D} \cdot (\boldsymbol{\epsilon}^e_n + \boldsymbol{\Delta \epsilon })$$ (406)

If

$$\boldsymbol{a} \cdot \boldsymbol{\sigma }_{n+1} - 2 c_{n+1} \sqrt{k} \le 0,$$ (407)

the solution is found. If not, $\Delta \gamma_{n+1} := \gamma_{n+1}-\gamma_n > 0$ and due to the Kuhn-Tucker equations $f(\boldsymbol{\sigma },c)=0$ is required. The stress at $n+1$ can be written as:

$$\boldsymbol{\sigma }_{n+1} = \boldsymbol{D} \cdot \boldsymbol{\epsilon }^e_{n+1}$$ (408)

$$= \boldsymbol{D} \cdot (\boldsymbol{\epsilon }_{n+1} - \boldsymbol{\epsilon }_{n+1}^p)$$ (409)

$$= \boldsymbol{D} \cdot (\boldsymbol{\epsilon }_{n+1} - \boldsymbol{\epsilon }_{n}^p+ \boldsymbol{\epsilon }_{n}^p- \boldsymbol{\epsilon }_{n+1}^p)$$ (410)

$$= \boldsymbol{\sigma }_{n+1}^ {\text{trial}} - \boldsymbol{D} \cdot \boldsymbol{\Delta \epsilon }_{n+1}^p$$ (411)

$$= \boldsymbol{\sigma }_{n+1}^ {\text{trial}} - \Delta \gamma \boldsymbol{D} \cdot \boldsymbol{b}.$$ (412)

In the last equation the evolution equation was used (flow rule). Defining $\boldsymbol{s}:= \boldsymbol{D} \cdot \boldsymbol{b}$, satisfying the yield surface now requires:

$$\boldsymbol{a} \cdot ( \boldsymbol{\sigma }_{n+1}^ {\text{trial}}... ...\gamma_{n+1} \; \boldsymbol{s}) -2 c(\epsilon _{n+1}^{\text{peq}}) \sqrt{k} =0.$$ (413)

Notice that $\boldsymbol{s}$ points in the direction of the stress correction w.r.t. the trial stress. In order to solve this equation a relationship between $\epsilon _{n+1}^{\text{peq}}$ and $\Delta \gamma$ is needed. Now, $\dot{\epsilon} _{n+1}^{\text{peq}}$ is defined as:

$$\dot{\epsilon}_{n+1}^{\text{peq}} = \sqrt{\frac{2}{3} } \Vert \dot {\boldsymbol{\epsilon } }_{n+1}^p \Vert,$$ (414)

from which:

$$\epsilon_{n+1}^{\text{peq}} = \int_{0}^{t+\Delta t} \sqrt{\frac{2}{3} } \Vert \dot {\boldsymbol{\epsilon } }^p \Vert dt$$ (415)

$$= \int_{0}^{t} \sqrt{\frac{2}{3} } \Vert \dot {\boldsymbol{\epsilon... ...t+\Delta t} \sqrt{\frac{2}{3} } \Vert \dot {\boldsymbol{\epsilon } }^p \Vert dt$$ (416)

$$= \epsilon_n^{\text{peq}} + \sqrt{\frac{2}{3} } \Vert {\boldsymbol{\Delta \epsilon } }_{n+1}^p \Vert$$ (417)

$$= \epsilon_n^{\text{peq}} + \sqrt{\frac{2}{3} } \Delta \gamma_{n+1} \; \Vert \boldsymbol{b} \Vert$$ (418)

$$= \epsilon_n^{\text{peq}} + \sqrt{\frac{2}{3} } \Delta \gamma_{n+1} \; \sqrt{1+m^2}$$ (419)

Defining $k^*:=2 \sqrt{k}$ and $m^*:= \sqrt{2(1+m^2)/3}$ the yield equation now amounts to:

$$f=\boldsymbol{a} \cdot ( \boldsymbol{\sigma }_{n+1}^ {\text{trial... ...; \boldsymbol{s}) - k^* c(\epsilon_n^{\text{peq}}+ \Delta \gamma_{n+1} m^* )=0.$$ (420)

This is a nonlinear equation in $\Delta \gamma_{n+1}$. Using the Newton-Raphson procedure the first derivative is needed, which yields:

$$\frac{\partial f}{ \partial \Delta \gamma_{n+1}} = -\boldsymbol{a... ...oldsymbol{s} - k^*m^* \frac{\partial c}{\partial \epsilon_{n+1}^{\text{peq}}} .$$ (421)

Starting with an initial guess $\Delta \gamma_{n+1}^0=0$, one arrives at $\Delta \gamma_{n+1}^{(k+1)} = \Delta \gamma_{n+1}^{(k)} + \Delta\Delta \gamma_{n+1}^{(k)}$ by solving the equation:

$$f(\Delta \gamma_{n+1}^{(k)} + \Delta\Delta \gamma_{n+1}^{(k)}) \a... ...{n+1}^{(k)}) + f'(\Delta \gamma_{n+1}^{(k)}) \Delta\Delta \gamma_{n+1}^{(k)}=0,$$ (422)

or

$$\left [ - \boldsymbol{a} \cdot \boldsymbol{s} - k^*m^* \left . \f... ...a_{n+1}^{(k)} \boldsymbol{s} \right ] - k^* c(\epsilon_{n+1}^{\text{peq},(k)}),$$ (423)

where

$$\epsilon_{n+1}^{\text{peq},(k)}=\epsilon_{n}^{\text{peq}}+m^* \Delta \gamma _{n+1}^{(k)}.$$ (424)

Figure 143: Limits of the return regions for Mohr-Coulomb

This works as long a the return is to a face of the yield surface. Since the return vector for sector 1 is $\boldsymbol{s_1}=\boldsymbol{D} \cdot \boldsymbol{b_1}$ one can graphically plot the region in three-dimensional principal stress space which is returned to the yield face in sector 1 (Figure 143). It is the space for which

$$(\boldsymbol{\sigma }_{n+1}^ {\text{trial}}-\boldsymbol{s_a} ) \c... ...al}}-\boldsymbol{s_a} ) \cdot (\boldsymbol{s_1} \times \boldsymbol{r_6}) \le 0.$$ (425)

This is called region I. The space within sector 1 in between region I and sector 2 and underneath the apex is mapped onto the intersection line of face 1 and face 2 of the yield surface (characterized by the equation $\boldsymbol{s_a} +\lambda \boldsymbol{r_1}, \lambda \in \mathbb{R}$). It is characterized by the equations

$$(\boldsymbol{\sigma }_{n+1}^ {\text{trial}}-\boldsymbol{s_a} ) \c... ...al}}-\boldsymbol{s_a} ) \cdot (\boldsymbol{s_1} \times \boldsymbol{s_2}) \le 0.$$ (426)

and is called region II of sector 1. Similary for the space in between region I and sector 6 and underneath the apex satisfying

$$(\boldsymbol{\sigma }_{n+1}^ {\text{trial}}-\boldsymbol{s_a} ) \c... ...al}}-\boldsymbol{s_a} ) \cdot (\boldsymbol{s_6} \times \boldsymbol{s_1}) \le 0.$$ (427)

This is region III. Finally, the region above the apex corresponding to

$$(\boldsymbol{\sigma }_{n+1}^ {\text{trial}}-\boldsymbol{s_a} ) \c... ...ial}}-\boldsymbol{s_a} ) \cdot (\boldsymbol{s_6} \times \boldsymbol{s_1}) \ge 0$$ (428)

is returned to the apex and is called region IV.

Returning the trial stress to the intersection line of two yield faces implies that the final state has to satisfy both face equations. The equivalent plastic strain now satisfies:

$$\epsilon_{n+1}^{\text{peq}} = \epsilon_n^{\text{peq}} + m^*( \Delta \gamma_{1,{n+1}} + \Delta \gamma_{2,{n+1}})$$ (429)

and the yield surface equations amount to:

$$\boldsymbol{a_1} \cdot ( \boldsymbol{\sigma }_{n+1}^ {\text{trial... ...\Delta \gamma_{2,n+1} \; \boldsymbol{s_2}) - k^* c(\epsilon_{n+1}^{\text{peq}}) =$$ 0 (430)

$$\boldsymbol{a_2} \cdot ( \boldsymbol{\sigma }_{n+1}^ {\text{trial... ...\Delta \gamma_{2,n+1} \; \boldsymbol{s_2}) - k^* c(\epsilon_{n+1}^{\text{peq}}) = 0 .$$ (431)

Applying Newton-Raphson now leads to:

$$[A] \cdot \begin{Bmatrix}\Delta \Delta \gamma_1 \\ \Delta \Delta \gamma_2 \end{Bmatrix} _{n+1}^{(k)} = \{B\},$$ (432)

where

$$[A]= \begin{bmatrix}- \boldsymbol{a_1} \cdot \boldsymbol{s_1} - k... ...on^{\text{peq}} } \right \vert _{\epsilon_{n+1}^{\text{peq},(k)}} \end{bmatrix}$$ (433)

and

$$\{B\}= \begin{Bmatrix}\boldsymbol{a_1} \cdot \left [ \boldsymbol{... ...oldsymbol{s_2} \right ] - k^* c(\epsilon_{n+1}^{\text{peq},(k)}) \end{Bmatrix}.$$ (434)

For a return to the apex the yield equations of sector 1, 2 and 6 have to be satisfied yielding 3 equations in 3 unknowns.

To calculate the consistent tangent the change of stress due to a change in total strain is needed. The change of stress in the present increment amounts to:

$$d(\boldsymbol{\Delta \sigma }) = d(\boldsymbol{D} \cdot \boldsymbol{\Delta \epsilon}^e)$$ (435)

$$= d(\boldsymbol{\boldsymbol{D} \cdot \Delta \epsilon}-\boldsymbol{D} \cdot \boldsymbol{\Delta \epsilon}^p)$$ (436)

$$= d(\boldsymbol{D} \cdot \boldsymbol{\Delta \epsilon}-\boldsymbol{D} \cdot \frac{\partial g}{\partial \boldsymbol{\sigma } } \Delta \gamma )$$ (437)

$$= \boldsymbol{D} \cdot d \boldsymbol{\Delta \epsilon } - \boldsymbo... ...mbol{D} \cdot \frac{\partial g}{\partial \boldsymbol{\sigma } } d \Delta \gamma$$ (438)

$$= \boldsymbol{D} \cdot d \boldsymbol{\Delta \epsilon } - \boldsymbol{s} d \Delta \gamma$$ (439)

Now, a relationship between $d \Delta \gamma$ and $d \boldsymbol{\Delta \epsilon }$ is needed. This is obtained from the yield equation:

$$f(\boldsymbol{\sigma }_n + \boldsymbol{\Delta \sigma } ,c(\gamma + \Delta \gamma))=0,$$ (440)

from which

$$\left . \frac{\partial f}{\partial \boldsymbol{\sigma } } \right ... ...c}{\partial \gamma} \right \vert _ {\gamma + \Delta \gamma} d \Delta \gamma =0.$$ (441)

Substituting the previously obtained expression for $d(\boldsymbol{\Delta \sigma })$ yields (dropping the bound of the derivatives to simplify the notation; the derivatives are to be taken at the end of the present increment):

$$\frac{\partial f}{\partial \boldsymbol{\sigma } } \cdot \boldsymb... ...c{\partial f}{\partial c} \frac{\partial c}{\partial \gamma} d \Delta \gamma =0$$ (442)

from which

$$d \Delta \gamma = \frac{\boldsymbol{a} \cdot \boldsymbol{D} \cdot... ... \boldsymbol{s} + k^*m^* \frac{\partial c}{\partial \epsilon ^{\text{peq}}} } .$$ (443)

Finally one obtains for the elastoplastic consistent tangent

$$\boldsymbol{D}^{ep} = \boldsymbol{D} - \frac{\boldsymbol{s} \otim... ...t \boldsymbol{s} + k^*m^* \frac{\partial c}{\partial \epsilon ^{\text{peq}}}} .$$ (444)

The derivation of the elastoplastic consistent tangent for a return to the intersection between face 1 and face 2 is similar. Now, a change in stress increment amounts to:

$$d \boldsymbol{\Delta \sigma } = \boldsymbol{D} \cdot d \boldsymbo... ...1} d \Delta \gamma_1 - \boldsymbol{D} \cdot \boldsymbol{b_2} d \Delta \gamma_2.$$ (445)

A relationship between the change in $\Delta \gamma$ and the change in $\boldsymbol{\Delta \epsilon }$ can again be obtained by considering the adjacent yield faces (i=1,2):

$$f_i(\boldsymbol{\sigma}_n + \boldsymbol{\Delta \sigma }, c(\gamma + \Delta \gamma_1 + \Delta \gamma_2))=0,$$ (446)

from which

$$\left . \frac{\partial f_i}{\partial \boldsymbol{\sigma } } \righ... ... \right \vert _ {\gamma + \Delta \gamma_1+\Delta \gamma_2} d \Delta \gamma_1 =0$$ (447)

or

$$\begin{bmatrix}\boldsymbol{a_1} \cdot \boldsymbol{s_1} + k^*m^* \... ...{a_2} \cdot \boldsymbol{D} \cdot d\boldsymbol{\Delta \epsilon } \end{Bmatrix} .$$ (448)

Denoting the inverse of the left hand matrix by $[B]$ the solution of the above system can be written as:

$$d \Delta \gamma_1 = B_{11} \boldsymbol{a_1} \cdot \boldsymbol{D} \cdot d\boldsymbol{\... ...{12} \boldsymbol{a_2} \cdot \boldsymbol{D} \cdot d\boldsymbol{\Delta \epsilon }$$ (449)

$$d \Delta \gamma_2 = B_{12} \boldsymbol{a_1} \cdot \boldsymbol{D} \cdot d\boldsymbol{\... ...22} \boldsymbol{a_2} \cdot \boldsymbol{D} \cdot d\boldsymbol{\Delta \epsilon },$$ (450)

leading to:

$$\boldsymbol{D}^{ep} = \boldsymbol{D} - B_{11} \boldsymbol{s_1} \o... ...bol{D} - B_{11} \boldsymbol{s_2} \otimes \boldsymbol{a_2} \cdot \boldsymbol{D}.$$ (451)

The derivation for the return to the apex is similar (now three equations have to be satisfied).

The tangent derived here only applies to the normal directions. It has to be complemented by a shear part $\hat{T}$ in the form [20]:

$$\boldsymbol{C}= \begin{bmatrix}\boldsymbol{D}^{ep} & 0 \\ 0 & \boldsymbol{\hat{T}} \end{bmatrix}$$ (452)

where

$$\boldsymbol{\hat{T}} = \begin{bmatrix}\frac{\sigma_{1,n+1}-\sigma... ...1}}{\sigma_{2,n+1}^{\text{trial}}-\sigma_{3,n+1} ^{\text{trial}}} \end{bmatrix}$$ (453)

Finally, the tangent matrix $\boldsymbol{C}$ has to be transformed back into the global system yielding $\boldsymbol{C'}$:

$$\boldsymbol{C'} = \boldsymbol{A}^T \cdot \boldsymbol{C} \cdot \boldsymbol{A},$$ (454)

where $\boldsymbol{A}$ is the transformation matrix in equation (A5) in [20]. Notice that equation (A6) in [20] is wrong and has to be replaced by the equation above. Similarly the stress is to be transformed into the global system by:

$$\boldsymbol{\sigma '} = \boldsymbol{A}^T \cdot \boldsymbol{\sigma }.$$ (455)