Potential flow is assumed, so

$$\begin{array}{}\mathrm{(3)}\phantom{\rule{3.0em}{0ex}}& d{\mathit{\epsilon}}^{pl}=\frac{d{\overline{\epsilon}}^{pl}}{g}\frac{\partial G}{\partial \mathit{\sigma}},\end{array}$$

where g can be written as

$$g=\frac{1}{c}\mathit{\sigma}:\frac{\partial G}{\partial \mathit{\sigma}}$$

and G is the flow potential, chosen as a hyperbolic function in the meridional stress plane and a smooth elliptic function in the deviatoric stress plane:

$$\begin{array}{}\mathrm{(4)}\phantom{\rule{3.0em}{0ex}}& G=\sqrt{{\left({\u03f5c|}_{0}\mathrm{tan}\psi \right)}^{2}+{\left({R}_{mw}q\right)}^{2}}-p\mathrm{tan}\psi ,\end{array}$$

where $\psi \left(\theta ,{f}^{\alpha}\right)$ is the dilation angle measured in the p–${R}_{mw}q$ plane at high confining pressure; ${c|}_{0}={c|}_{{\overline{\epsilon}}^{pl}=0}$ is the initial cohesion yield stress; and $\u03f5$ is a parameter, referred to as the eccentricity, that defines the rate at which the function approaches the asymptote (the flow potential tends to a straight line as the eccentricity tends to zero). This flow potential, which is continuous and smooth in the meridional stress plane, ensures that the flow direction is defined uniquely in this plane. The function asymptotically approaches a linear flow potential at high confining pressure stress and intersects the hydrostatic pressure axis at 90°. A family of hyperbolic potentials in the meridional stress plane is shown in Figure 4.

Figure 4. Family of hyperbolic flow potentials in the meridional plane.
The flow potential is also continuous and smooth in the deviatoric stress plane (the $\mathrm{\Pi}$-plane); we adopt the deviatoric elliptic function used by Menétrey and Willam (1995):

$${R}_{mw}\left(\mathrm{\Theta},e\right)=\frac{4\left(1-{e}^{2}\right){\mathrm{cos}}^{2}\mathrm{\Theta}+{\left(2e-1\right)}^{2}}{2\left(1-{e}^{2}\right)\mathrm{cos}\mathrm{\Theta}+\left(2e-1\right)\sqrt{4\left(1-{e}^{2}\right){\mathrm{cos}}^{2}\mathrm{\Theta}+5{e}^{2}-4e}}{R}_{mc}\left(\frac{\pi}{3},\varphi \right),$$

where $\mathrm{\Theta}$ is the deviatoric polar angle defined previously, ${R}_{mc}\left(\frac{\pi}{3},\varphi \right)=\left(3-\mathrm{sin}\varphi \right)/6\mathrm{cos}\varphi $, and e is a parameter that describes the “out-of-roundedness” of the deviatoric section in terms of the ratio between the shear stress along the extension meridian ($\mathrm{\Theta}=0$) and the shear stress along the compression meridian ($\mathrm{\Theta}=\frac{\pi}{3}$). The elliptic function has the value ${R}_{mw}(\mathrm{\Theta}=0,e)={R}_{mc}(\frac{\pi}{3},\varphi )/e$ along the extension meridian and has the value ${R}_{mw}(\mathrm{\Theta}=\frac{\pi}{3},e)={R}_{mc}(\frac{\pi}{3},\varphi )$ along the compression meridian; this ensures that the flow potential matches the yield surface at the triaxial compression and extension in the deviatoric plane provided that e is defined appropriately (see further discussion later). Although the elliptic function is defined only in the sector $0\le \mathrm{\Theta}\le \pi /3$, the polar radius ${R}_{mw}\left(\mathrm{\Theta},e\right)$ extends to all polar directions $0\le \mathrm{\Theta}\le 2\pi $ using the three-fold symmetry shown in Figure 5.

Figure 5. Menétrey-Willam flow potential in the deviatoric plane.
By default, the out-of-roundedness parameter, e, is dependent on the friction angle $\varphi $; it is calculated by matching the flow potential to the yield surface in both triaxial tension and compression in the deviatoric plane:

$$e=\frac{3-\mathrm{sin}\varphi}{3+\mathrm{sin}\varphi}.$$

Alternatively, e can also be considered to be an independent material parameter; in this case the user can provide its value directly. Convexity and smoothness of the elliptic function requires that $1/2<e\le 1$. The upper limit, $e=1$ (or $\varphi =$ 0°), leads to ${R}_{mw}(\mathrm{\Theta},e=1)={R}_{mc}(\frac{\pi}{3},\varphi )$, which describes the Mises circle in the deviatoric plane. The lower limit, $e=1/2$ (or $\varphi =$ 90°), leads to ${R}_{mw}(\mathrm{\Theta},e=1/2)=2{R}_{mc}(\frac{\pi}{3},\varphi )\mathrm{cos}\mathrm{\Theta}$ and would describe the Rankine triangle in the deviatoric plane (this limiting case is not permitted within the Mohr-Coulomb model described here).

Flow in the meridional stress plane can be close to associated when the angle of friction, $\varphi $, and the angle of dilation, $\psi $, are equal and the eccentricity parameter, $\u03f5$, is very small; however, flow in this plane is, in general, nonassociated. Flow in the deviatoric stress plane is always nonassociated. Therefore, the use of this Mohr-Coulomb model generally requires the solution of nonsymmetric equations.