Mohr-Coulomb return mapping: the corners and the apex

Implementing Mohr-Coulomb is not the same problem as writing it down. Written down it is a line of algebra that every geotechnical engineer knows. Implemented, it is a pyramid with six faces, six edges and a point, and the algebra applies to the faces only.

This article works one stress increment through the standard single-surface return and shows what it returns. The result is not an approximation to the right answer; it is inadmissible.

Every number in this article is recomputed from its inputs by a script that runs on each build of this site.

The surfacePermalink to “The surface”#

With compression positive and $\sigma_1 \ge \sigma_2 \ge \sigma_3$, the Mohr-Coulomb yield function is

$$f = \frac{\sigma_1 - \sigma_3}{2} - \frac{\sigma_1 + \sigma_3}{2}\sin\varphi - c\cos\varphi$$

which is exactly the tangency condition on a Mohr circle, rearranged. Note what it does not contain: $\sigma_2$. The criterion is indifferent to the intermediate principal stress, and that indifference is the source of the corners — in principal stress space the surface is a hexagonal pyramid whose axis is the hydrostatic line, and the six sextants of the deviatoric plane correspond to the six possible orderings of the three principal stresses.

The plastic potential has the same form with the dilatancy angle in place of the friction angle:

$$g = \frac{\sigma_1 - \sigma_3}{2} - \frac{\sigma_1 + \sigma_3}{2}\sin\psi$$

Take $E = 20000$ kPa, $\nu = 0.30$, $\varphi' = 30°$, $c' = 10$ kPa and $\psi = 0$, giving $G = 7692.3077$ kPa and $\lambda = 11538.4615$ kPa.

The elastic predictorPermalink to “The elastic predictor”#

Start from an isotropic state of 100 kPa. There

$$f = 0 - 50 - 8.660254 = -58.660 \ \text{kPa}$$

which is comfortably elastic. Apply a strain increment of $(0.010,\ -0.005,\ -0.005)$ — constant volume, pure shear in the principal directions. With zero volumetric strain the Lamé term drops out and

$$\Delta\sigma = 2G\,\Delta\varepsilon = (153.8462,\ -76.9231,\ -76.9231)\ \text{kPa}$$ $$\sigma^{\text{trial}} = (253.8462,\ 23.0769,\ 23.0769)\ \text{kPa}$$ $$f(\sigma^{\text{trial}}) = 115.3846 - 69.2308 - 8.660254 = 37.494\ \text{kPa} > 0$$

The trial state is outside the surface, so a plastic correction is required.

The return, done the way it is written downPermalink to “The return, done the way it is written down”#

The return is along the elastic tangent applied to the flow direction:

$$\sigma^{n+1} = \sigma^{\text{trial}} - \Delta\lambda\, \mathbf{D}\,\frac{\partial g}{\partial\sigma}, \qquad \Delta\lambda = \frac{f(\sigma^{\text{trial}})} {\dfrac{\partial f}{\partial\sigma}^{\!\mathsf T}\mathbf{D}\dfrac{\partial g}{\partial\sigma}}$$

With $\psi = 0$ the flow direction is $\partial g/\partial\sigma = (0.5,\ 0,\ -0.5)$, whose trace is zero, so the elastic operator reduces to

$$\mathbf{D}\frac{\partial g}{\partial\sigma} = 2G\,(0.5,\ 0,\ -0.5) = (7692.3077,\ 0,\ -7692.3077)$$

and with $\partial f/\partial\sigma = (0.25,\ 0,\ -0.75)$ the denominator is 7692.3077, giving

$$\Delta\lambda = \frac{37.494}{7692.3077} = 0.00487$$ $$\sigma^{n+1} = (216.3526,\ 23.0769,\ 60.5705)\ \text{kPa}$$

Check the yield function on that state and it is satisfied to machine precision. The return worked. Every step above is correct.

What is wrong with the answerPermalink to “What is wrong with the answer”#

$$\sigma_2 = 23.0769 < \sigma_3 = 60.5705$$

The derivation assumed $\sigma_1 \ge \sigma_2 \ge \sigma_3$, and the state it produced does not satisfy it. The intermediate principal stress has ended up 37.494 kPa below the minor one — which, not coincidentally, is exactly the trial value of the yield function, because $\sigma_3$ was raised by $\Delta\lambda \cdot 2G / 2$ while $\sigma_2$ was left alone.

The stress point has been returned to the plane containing one face of the pyramid, but it has left the sextant in which that face is the surface. Outside its sextant, that plane is not the yield surface — it is an extension of it into a region governed by a different face. The point satisfies the equation of a surface that does not exist there.

A program that stops here does not crash and does not warn. It returns a converged increment, and the error is a stress state that is wrong in a way no residual check will find, because the residual it checks is the one it just satisfied.

The three regionsPermalink to “The three regions”#

The correct treatment recognises that a return can land in three kinds of place.

On a face. The single-surface return applies, the ordering is preserved, and the result is admissible. This is the common case and it is the one the textbook covers.

On an edge. Two faces meet, and the returned state must satisfy both yield functions at once. The return is to the intersection line, with two plastic multipliers determined together. In principal stress space the two edge families are $\sigma_1 = \sigma_2$ (triaxial extension) and $\sigma_2 = \sigma_3$ (triaxial compression), and the case above belongs to the second.

At the apex. The pyramid closes to a point on the hydrostatic axis, at

$$p_{\text{apex}} = c\cot\varphi = \frac{10}{\tan 30°} = 17.321\ \text{kPa}$$

of hydrostatic tension. A trial state whose return would carry it beyond that point has nowhere on the surface to go, and the only admissible answer is the apex itself — all three principal stresses equal, all deviatoric stress removed in one step. The gradient of $f$ is undefined there, so no formula involving $\partial f/\partial\sigma$ can produce it.

Why this is worth handling exactlyPermalink to “Why this is worth handling exactly”#

The alternative in common use is to round off the corners — replacing the hexagon with a smooth surface that approximates it, usually a hyperbolic or Drucker-Prager-like fit. That makes every return a face return and removes the problem entirely.

It also changes the material. A smoothed surface is a different yield criterion, with different strengths in the directions where the smoothing bites hardest, and those are the triaxial extension and compression states that laboratory testing actually measures. The error is small in the middle of a face and largest exactly where the calibration data lives.

Handling the regions exactly costs a branch and, for the edge case, a two-surface solve. That is a bounded amount of work, done once, in return for a criterion that means what the laboratory reported. It is also the kind of thing that has to be verified at the stress point and again in a boundary value problem, because a model verified at the material point is not yet a model that works — and the second test is where these faults are found.

One more consequence worth stating: the flow rule above used $\psi = 0$, and that choice is not a detail either. It sets both the return direction and the volume change, and what it costs is a separate question with a larger answer than most people expect.