DOI:10.14311/APP.2021.30.0069 Acta Polytechnica CTU Proceedings 30:69–75, 2021 © Czech Technical University in Prague, 2021 available online at http://ojs.cvut.cz/ojs/index.php/app COMPARING THE HOEK-BROWN AND MOHR-COULOMB FAILURE CRITERIA IN FEM ANALYSIS Tereza Poklopová∗, Veronika Pavelcová, Michal Šejnoha Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Prague 6, Czech Republic ∗ corresponding author: tereza.poklopova@fsv.cvut.cz Abstract. 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 finite element software [2]. Keywords: Equivalent shear strength, finite element method, Hoek-Brown criterion, nonlinear yield function, rock mechanics. 1. 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–5] proved more appropri- ate, the use of the Mohr-Coulomb model in the stabil- ity 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 sup- ported by the possibility of estimating these param- eters 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 per- haps more straightforward to proceed in the footsteps of Benz et al. [6] who introduced a simple modifica- tion 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. 2. Hoek-Brown model The Hoek-Brown model was first introduced by Hoek [3] in application to the analysis of an under- ground excavation of intact rock masses. Later in [4] this criterion was extended to weak rock masses by adjusting the material constants according to a geo- logical quality. The geological strength index (GSI) is used as a qualitative classification number. By imple- menting 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 σ3 = σ1 − σci ! mb −σ1 σci + s "a , (1) where σ1 and σ3 are the maximum and minimum ef- fective principal stresses, respectively1. The material constants mb, s and a follow from empirical equations as mb = mi exp GSI − 100 28 − 14D , (2) s = exp GSI − 100 9 − 3D , (3) a = 1 2 + 1 6 ! exp −GSI 15 − exp −20 3 " , (4) where the disturbance factor D account for prior ex- cavations and ranges from 0 (undisturbed rock) to 1 (disturbed by production blasting). The uniaxial compressive strength σci and the model parameter mi are considered for an intact rock. The rock mass stiffness is assumed constant and independent of the current stresses. Following [1], the modulus of elas- ticity is estimated by Em(GPa) = ! 1 − D 2 " # σci 100 10((GSI−10)/40), (5) Em(GPa) = ! 1 − D 2 " 10((GSI−10)/40). (6) 1Unlike typical geotechnical notation, the tensile stresses are assumed positive. 69 http://dx.doi.org/10.14311/APP.2021.30.0069 http://ojs.cvut.cz/ojs/index.php/app T. Poklopová, V. Pavelcová, M. Šejnoha Acta Polytechnica CTU Proceedings Equation (5) is called when σci ≤ 100 MPa and sug- gests a rapid decrease of the modulus with a decreas- ing uniaxial compressive strength. Equation (6) is then used when σci > 100 MPa. The above equa- tions were exploited in this paper regardless of their updated version presented in [5]. 3. Analogy with the Mohr-Coulomb criterion In the finite element analysis the failure criterion has a meaning of a yield function bounding the elastic re- sponse. The GEO5 FEM program assumes the stress- strain relation to be elastic perfectly plastic. The Hoek-Brown yield function then follows from Eq. (1) in the form f (σ1, σ3) = σ1 − σ3 − σci ! mb −σ1 σci + s "a . (7) 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 re- turn mapping is driven by a non-associated flow rule with the plastic potential function given by g(σ1, σ3) = 1 + sin ψ 1 − sin ψ σ1 − σ3, (8) 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 vari- ation in the deviatoric plane, see also the projection into σ1-σ3 plane displayed in Fig. 1. Figure 1. Hoek-Brown yield function plotted in plane of maximum and minimum principal stresses. 3.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), σ3 = −σcm + kσ1, (9) 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 c = σcm 2 √ k , sin ϕ = k − 1 k + 1 . (10) Comparing Eqs. (1) and (9) it is seen that no math- ematical relationship between the two models can be determined. To that end, an iterative approach bal- ancing the area above and below the MC model, see Fig. 2, was developed in [1] to get c = σci[(1 + 2a)s + (1 − a) mbσ1n](s + mbσ1n)a−1 A $ 1 + (6amb(s + mbσ1n)a−1)/A , (11) A = (1 + a)(2 + a), sin ϕ = % 6amb(s + mbσ1n)a−1 2(1 + a)(2 + a) + 6amb(s + mbσ1n)a−1 & , (12) for the range of stresses of σt > σ1 > σ1,max, where the bi-axial tensile stress σt = σcis/mb is found by setting σ1 = σ3 = σt in Eq. (1), and σ1n = |σ1,max|/σci depends on the upper limit of the confining stress σ1,max seen in Fig. 2. It is this param- eter we accord our attention to in the next section. Figure 2. Nonlinear Hoek-Brown criterion and equivalent Mohr-Coulomb model in σ1-σ3 plane. 4. 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 σ1 = σ2 = ν 1 − ν σz , σ3 = σz , (13) 70 vol. 30/2021 Hoek-Brown and Mohr-Coulomb failure criteria Parameter Symbol Value Unit Poisson’s ratio ν 0.3 [−] Uniaxial compressive strength σci 20.0 [MPa] Hoek-Brown constant mi 8.0 [−] Geological Strength index GSI 30.0 [−] Disturbance factor D 0.0 [−] Reduced Hoek-Brown constant mb 0.656680 [−] Parameter s 0.000419 [−] Parameter a 0.522340 [−] Modulus of elasticity Em 1414.20 [MPa] Table 1. Input and calculated parameters of rock mass. Figure 3. Geometry, boundary and loading condi- tions of FEM model of oedometer. 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 parame- ters 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 gov- erned by specific loading conditions. Thus an incor- rect choice of this parameter may lead to erroneous predictions of the MC model. To see this, we begin with a value σ1n = 0.25 rec- ommended in [4] for triaxial simulations. Assuming the values of the HB model provided in Table 1 we get from Eqs. (11) and (12) c = 649.0 kPa, ϕ = 22.8◦. Substituting Eq. (13) into the Hoek-Brown yield func- tion (7) ν 1 − ν σz − σz − σci ! mb ν 1 − ν −σz σci + s "a = 0, (14) gives the limit value of the vertical stress for the first plastic yielding as σz = −16.165 MPa. Moving beyond elasticity calls for numerical simu- lations. 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 triangu- lar 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. (a) elastic stress path (b) Figure 4. Case 1 σ1n = 0.25: a) evolution of stresses based on HB and MC models, b) intersections of elas- tic stress path with yield surfaces. This simulation clearly demonstrates the influence of the chosen σ1n ratio on the degree of equivalence. Unlike the triaxial test2, where the stress state is known, the evolution of principal stresses in oedome- ter 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) plotting the evolution of stresses in terms of invariant stress measures, where σm is the mean stress and J is the square root of the second invariant of the deviatoric stress. Clearly, the HB model predicts the deviation from the elastic 2See for example [7] for more details on the HB model ver- ification using triaxial simulations. 71 T. Poklopová, V. Pavelcová, M. Šejnoha Acta Polytechnica CTU Proceedings (a) elastic stress path (b) Figure 5. Case 2 σ1n = 0.79: a) evolution of stresses based on HB and MC models, b) intersections of elas- tic stress path with yield surfaces. response precisely at the elastic limit point for the ap- plied load σz = −16.165 MPa, while the MC model still suggests an elastic behavior. The inadequacy of the MC model setting is better seen in Fig. 4(b) plot- ting the elastic stress-path with the onset of yielding identified by the plus sign for both yield surfaces. Ev- idently, the assumed range of stresses for which the equivalence holds, recall Fig. 2, is outside the gener- ated stress state controlled by the MC model. The new estimate of σ1n ratio was obtained numer- ically adopting the HB model to get the maximum lateral principle stress σ1,max = −15.79 MPa for the maximum applied load σz = σ3 = 30 MPa. This gives σ1n = 0.79 and consequently new estimates of strength parameters c, ϕ c = 1345.5 kPa, ϕ = 15.6◦. The corresponding evolution of stresses is shown in Fig. 5(a). A reasonable choice of σ1n ratio is fur- ther 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 gen- erated either in oedometer or triaxial apparatus is constant within the sample. However, in real appli- cations 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 elastic stress path (a) (numerical solution) (b) Figure 6. a) Influence of estimated stress scale, b) evolution of stresses based on HB and MC models. 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 anal- ysis is examined in the next section. 5. Slope stability analysis The slope stability analysis typically aims at provid- ing the factor of safety. For the family of elastic per- fectly 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 con- vergence of the nonlinear problem still exists. This approach is thus not directly applicable to the HB model. However, the authors in [6] introduced an el- egant way to overcome this obstacle. So the analysis presented next proceeds in their footsteps. 5.1. 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 f (σ1, σ3) = σ1 − σ3 − σci η ! mb −σ1 σci + s "a , (15) 72 vol. 30/2021 Hoek-Brown and Mohr-Coulomb failure criteria Parameter Symbol Value Unit Unit weight γ 25.0 ' kN/m3 ( Poisson’s ratio ν 0.3 [−] Uniaxial compressive strength σci 30.0 [MPa] Hoek-Brown constant mi 2.0 [−] Geological Strength Index GSI 5.0 [−] Disturbance factor D 0.0 [−] Reduced Hoek-Brown constant mb 0.067225 [−] Parameter s 0.000026 [−] Parameter a 0.019210 [−] Modulus of elasticity Em 410.73 [MPa] Table 2. Input and calculated parameters of rock mass. where η is the reduction factor. A different reduc- tion factor is introduced to reduce the shear strength parameters ϕ, c as tan ϕd = tan ϕc γϕ , cd = cc γc , (16) where ϕc and cc are the actual strength parameters of the subsoil material and γϕ and γc are their corre- sponding reduction coefficients. A typical setting in the stability analysis assumes γϕ = γc = γ. The ob- jective 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 η = 1 2 ) ** +γ , 2 − f̃ ′ - .///01 + , 1 γ2 − 1 - 1 −f̃ ′ 22 1 2 − f̃ ′ 22 + f̃ ′ 3 44 5 , (17) where f̃ ′ = ∂f̃ (σ1) ∂σ1 = −mba ! mb −σ1 σci + s "a−1 . (18) 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 deriva- tion of Eq. (17) we refer the interested reader to [6]. 5.2. Estimating σ1,max in stability analysis The application of the MC model requires setting the value of σ1,max to estimate the shear strength param- eters. 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 |σ1max| = 0.72σcm ! σcm γH "−0.91 , (19) 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) σcm = 2c cos ϕ 1 − sin ϕ , (20) 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 σcm = σci (mb + 4s − a(mb − 8s))(mb/4 + s)a−1 2(1 + a)(2 + a) . (21) 5.3. Practical use of the methods To test reliability of Eq. (19), we performed the slope stability analysis of a homogeneous rock slope em- ploying 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 do- main 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 ma- terial 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). Figure 7. Geometry of FEM model, boundary con- ditions and finite element mesh. The analysis was carried out using stan- dard shear strength parameters reduction ap- proach implemented in GEO5 FEM. Adopting 73 T. Poklopová, V. Pavelcová, M. Šejnoha Acta Polytechnica CTU Proceedings FEM tool GEO5 FEM Plaxis 2019 Material model HB HB MC MC HB Element type T3 T6 T3 T6 T15 Safety factor 1.44 1.37 1.42 1.33 1.35 Table 3. Results of slope stability analysis. (a) (b) (c) (d) Figure 8. Potential slip surfaces represented by the localized equivalent deviatoric plastic strain: a) HB criterion, T3 elements, b) HB criterion, T6 elements, c) equivalent MC criterion, T3 elements, d) equivalent MC criterion, T6 elements. Eqs. (19), (11) and (12) yields |σ1max| = 189 kPa, c = 20 kPa, ϕ = 21◦. In the case of HB model the parameter η was cal- culated for the prescribed reduction factor γ using Eq. (17). The results are summarized in Table 3 for both types of approximations. The result acquired by em- ploying 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 dif- ferent finite element mesh. As expected, the T3 el- ement 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, compar- ing 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 con- sidered 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 prefer- able. 6. 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 di- rect 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 homoge- neous 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. Acknowledgements The support provided by the SGS project No. SGS20/038/OHK1/1T/11 is gratefully acknowl- edged. References [1] E. Hoek, C. Carranza-Torres, B. Corkum. Hoek- Brown failure criterion - 2002 edition. Proceedings of the 5th North American symposium - NARMS-TAC 2002. [2] Fine-Ltd. GEO5 FEM, user manual. http://www.fine.cz/, 2020. 74 vol. 30/2021 Hoek-Brown and Mohr-Coulomb failure criteria [3] E. Hoek, E. T. Brown. Underground excavations in rock. London: Institution of Mining and Metallurgy, 1980. [4] E. Hoek, E. T. Brown. Practical estimates of rock mass strength. International Journal of Rock Mechanics and Mining Sciences 34.8:1165–1186, 1997. [5] E. Hoek, E. T. Brown. The Hoek-Brown failure criterion and GSI - 2018 edition. Journal of Rock Mechanics and Geotechnical Engineering 2018. [6] T. Benz, R. Schwab, R. Kauther, P. A. Vermeer. A Hoek-Brown Criterion with instrinsic material strength factorization. International Journal of Rock Mechanics and Mining Sciences 45:210–222, 2008. [7] T. Poklopová. Prediction of the response of geotechnical structure using Hoek-Brown model - implementation, verification. Master’s thesis, Czech Technical University in Prague, Faculty of Civil Engineering, 2020. In Czech. [8] P. Ltd. PLAXIS manual. http://www.plaxis.nl/. 75