COMPARING THE HOEK-BROWN AND MOHR-COULOMB FAILURE CRITERIA IN FEM ANALYSIS

. This paper revisits the issue of a potential substitutions of the Hoek-Brown failure model by the standard Mohr-Coulomb model in the stability analysis of rock masses. The derivation of equivalent shear strength parameters of the Mohr-Coulomb proposed by Hoek et al. [1] is addressed with emphases on the suitable range of stresses for which the equivalence of the two failure criteria applies. To that end, a simple numerical analysis of the oedometric test is carried out. It is seen that a correct choice of the upper limit of the minimum compressive principal stress is crucial for the Mohr-Coulomb model to provide predictions comparable to the Hoek-Brown model. This issue is addressed next in the light of the solution of slope stability problem. All the presented results were derived with the help of the GEO5 FEM ﬁnite element software [2].


Introduction
Application of the Mohr-Coulomb (MC) model is a common approach when addressing the stability of soil slopes. Although for a rock mass the Hoek-Brown (HB) failure criterion [1,[3][4][5] proved more appropriate, the use of the Mohr-Coulomb model in the stability assessment of rock slopes is still quite popular in practice. This can be attributed to the fact that the factor of safety can easily be defined in terms of the basic shear strength parameters such as the cohesion c and the angle of internal friction ϕ. This is supported by the possibility of estimating these parameters through simple empirical formulas on the basis of HB model parameters [1]. This, however, requires the knowledge of the range of stresses for which the linear MC model provides results comparable to those predicted by the nonlinear HB model. So it is perhaps more straightforward to proceed in the footsteps of Benz et al. [6] who introduced a simple modification to the original HB model, which brings the two models to the same footing and allows for a standard strength reduction approach to estimate the factor of safety in the light of the Hoek-Brown model. Both issues are examined in this paper.

Hoek-Brown model
The Hoek-Brown model was first introduced by Hoek [3] in application to the analysis of an underground excavation of intact rock masses. Later in [4] this criterion was extended to weak rock masses by adjusting the material constants according to a geological quality. The geological strength index (GSI) is used as a qualitative classification number. By implementing the influence of prior excavation process, the current form of the HB model covers the whole range of rock masses. The generalized Hoek-Brown failure criterion is defined in terms of principal stresses as where σ 1 and σ 3 are the maximum and minimum effective principal stresses, respectively 1 . The material constants m b , s and a follow from empirical equations as where the disturbance factor D account for prior excavations and ranges from 0 (undisturbed rock) to 1 (disturbed by production blasting). The uniaxial compressive strength σ ci and the model parameter m i are considered for an intact rock. The rock mass stiffness is assumed constant and independent of the current stresses. Following [1], the modulus of elasticity is estimated by " 10 ((GSI−10)/40) .
Equation (5) is called when σ ci ≤ 100 MPa and suggests a rapid decrease of the modulus with a decreasing uniaxial compressive strength. Equation (6) is then used when σ ci > 100 MPa. The above equations were exploited in this paper regardless of their updated version presented in [5].

Analogy with the Mohr-Coulomb criterion
In the finite element analysis the failure criterion has a meaning of a yield function bounding the elastic response. The GEO5 FEM program assumes the stressstrain relation to be elastic perfectly plastic. The Hoek-Brown yield function then follows from Eq. (1) in the form When presented in the principal stress space, the Hoek-Brown yield function is similar to the MC model independent of the intermediate principal stress σ 2 . In case of plastic yielding, the stress return mapping is driven by a non-associated flow rule with the plastic potential function given by where ψ is the dilation angle and for the assumed perfect plasticity it remains constant. Similar to MC model, the HB model plots as irregular hexagon, but unlike the MC model it experiences a nonlinear variation in the deviatoric plane, see also the projection into σ 1 -σ 3 plane displayed in Fig. 1.

Equivalent shear strength parameters
The simplicity of the MC model and its prevailing popularity among design engineers promote its use in applications concerning the stability issues of a rock mass. The similarity of the MC and HB models was examined by Hoek et al. in [1,4] in order to derive the shear strength parameters ϕ, c controlling the MC failure criterion as, compare with Eq. (1), where σ cm is the uniaxial compressive strength of the rock mass and k is the slope of the MC plot in the σ 3 and σ 1 principle stress space, see Fig. 2. Effective values of the shear strength parameters c, ϕ are then given by Comparing Eqs. (1) and (9) it is seen that no mathematical relationship between the two models can be determined. To that end, an iterative approach balancing the area above and below the MC model, see Fig. 2, was developed in [1] to get (1), and σ 1n = |σ 1,max |/σ ci depends on the upper limit of the confining stress σ 1,max seen in Fig. 2. It is this parameter we accord our attention to in the next section.

Simple oedometer test
To investigate the influence of the choice of the upper stress limit σ 1,max we consider a simple oedometric test. Assuming elasticity we get vol. 30/2021 Hoek  where σ 3 < 0 is the applied vertical compressive stress.
As mentioned in the previous section, arriving at a comparable response predicted by the MC and HB models using the equivalent shear strength parameters strongly depends on the adopted range of stresses bounded by σ t and σ 1,max . While σ t depends on the HB model parameters, the choice of σ 1,max is governed by specific loading conditions. Thus an incorrect choice of this parameter may lead to erroneous predictions of the MC model.
To see this, we begin with a value σ 1n = 0.25 recommended in [4] for triaxial simulations. Assuming the values of the HB model provided in Table 1 we get from Eqs. (11) and (12) Substituting Eq. (13) into the Hoek-Brown yield function (7) gives the limit value of the vertical stress for the first plastic yielding as σ z = −16.165 MPa.
Moving beyond elasticity calls for numerical simulations. To that end, we adopted the GEO5 FEM software and set up the most simple geometrical model of an oedometer under plane strain conditions consisting of only two constant strain 3-node triangular elements. Boundary and loading conditions of the model are evident in Fig. 3. The model was loaded by vertical compression such as to pass the elasticity limit point.  This simulation clearly demonstrates the influence of the chosen σ 1n ratio on the degree of equivalence. Unlike the triaxial test 2 , where the stress state is known, the evolution of principal stresses in oedometer is controlled by the yield function, i.e. by the strength parameters of the rock mass. Therefore, the analytical determination of the limit stress value σ 1max is not trivial.
The consequence of an incorrect choice of σ 1n = 0.25 appears in Fig. 4(a)   The corresponding evolution of stresses is shown in Fig. 5(a). A reasonable choice of σ 1n ratio is further supported in Figs. 5(b) and 6(a) comparing the yield surfaces for individual settings. The actual stress path corresponding to both models appears in 6(b) suggesting a close match for the expected stress range.
It is worth mentioning that the stress state generated either in oedometer or triaxial apparatus is constant within the sample. However, in real applications a non-homogeneous stress state arises. One thus has to be extremely careful when specifying the expected range of stresses if adopting just one pair of equivalent shear strength parameters, essentially independent of the actual stress state. Nevertheless, for some specific applications simple relations for the determination σ 1,max have been proposed in [1]. One such a relationship for the case of slope stability analysis is examined in the next section.

Slope stability analysis
The slope stability analysis typically aims at providing the factor of safety. For the family of elastic perfectly plastic MC models the factor of safety is defined as the ratio of the current shear strength parameter, e.g. cohesion, to the minimum one for which the convergence of the nonlinear problem still exists. This approach is thus not directly applicable to the HB model. However, the authors in [6] introduced an elegant way to overcome this obstacle. So the analysis presented next proceeds in their footsteps.

Hoek-Brown model in stability analysis
The approach outlined in [6] builds upon a gradual reduction of the strength of a rock mass by modifying the yield function as vol. 30/2021 Hoek where η is the reduction factor. A different reduction factor is introduced to reduce the shear strength parameters ϕ, c as where ϕ c and c c are the actual strength parameters of the subsoil material and γ ϕ and γ c are their corresponding reduction coefficients. A typical setting in the stability analysis assumes γ ϕ = γ c = γ. The objective now is to determine the relationship between the two reduction factors η and γ. This is achieved by constructing a tangent to the Hoek-Brown yield function at the current state of stress sitting on the yield function and comparing this slope to the slope of the MC yield surface to get Notice that η is a function of the current stress state and must be evaluated for every element of the mesh in every loading step separately. For complete derivation of Eq. (17) we refer the interested reader to [6].

Estimating σ 1,max in stability analysis
The application of the MC model requires setting the value of σ 1,max to estimate the shear strength parameters. By performing a large parametric study on a variety of slope geometries and rock mass properties using Bishop's limit state analysis, Hoek et al. [1] proposed the following relationship for σ 1,max where H is the height of the slope, γ is the unit weight of the rock mass and σ cm is the rock mass strength defined by Eq. (10) The rock mass strength expressed using parameters of the Hoek-Brown model is then given by substituting Eqs. (11) and (12) into Eq. (20). It holds

Practical use of the methods
To test reliability of Eq. (19), we performed the slope stability analysis of a homogeneous rock slope employing both constitutive models. The geometry of the computational model of a slope with the height of H = 10 m and inclination of 35.5 • is shown in Fig. 7 together with the finite element mesh. The domain was discretized using either linear 3-node (T3) or quadratic 6-node (T6) triangular elements. An average element edge size was about 0.5 m. The material data listed in the upper part of Table 2 were taken from [6]. The model parameters presented in the lower part of this table were determined from Eqs. (2) -(6). The analysis was carried out using standard shear strength parameters reduction approach implemented in GEO5 FEM. Adopting    In the case of HB model the parameter η was calculated for the prescribed reduction factor γ using Eq. (17). The results are summarized in Table 3 for both types of approximations. The result acquired by employing the Plaxis 2D [8] software is also presented for the sake of comparison. Note that it was derived for the 15-node triangular element and a slightly different finite element mesh. As expected, the T3 element delivers with the same mesh refinement the most optimistic predictions. The results provided by the T6 element, the default setting in GEO5 FEM, are comparable to that of Plaxis 2D. Also, comparing the HB and MC predictions in terms of both the factor of safety and graphical representation of the slip surfaces in Fig. 8 supports Eq. (19) in providing a suitable estimate of σ 1,max .
It should, however, be pointed out that we considered a particularly simple case of a homogeneous slope. Thus if possible, the use of the HB model for the slope stability analysis as proposed in [6] is preferable.

Conclusions
The principal objective of this contribution was to highlight the importance of a correct choice of the stress ranges in the definition of equivalent shear strength parameters when replacing the Hoek-Brown failure model with the Mohr-Coulomb model. This was demonstrated by running a numerical analysis of a simple oedometer test for two different estimates of σ 1n ratio. While the first choice assuming simply the value recommended for triaxial simulations resulted in vastly different predictions when compared to a direct analysis with the HB model, the second choice based on a numerical estimate of the value of σ 1,max led to comparable predictions.
Next, the slope stability analysis was performed taking advantage of the procedure proposed in [6] so that standard approach based on a gradual reduction of strength parameters was possible to apply to both models. Although for the simple case of a homogeneous slope the estimate of σ 1,max proposed by Hoek et al. [1] delivers satisfactory results with the MC model, the approach based on a direct application of the HB model is preferable, particularly in the case of more complex subsoil conditions.