Acta Polytechnica CTU Proceedings https://doi.org/10.14311/APP.2023.40.0117 Acta Polytechnica CTU Proceedings 40:117–122, 2023 © 2023 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague RETURN MAPPING SCHEME FOR THE HOEK-BROWN MODEL WITH TENSION CUT-OFF Tereza Žalská∗, Michal Šejnoha Czech Technical University in Prague, Faculty of Civil Engineering, Thákurova 7, 166 29 Prague 6, Czech Republic ∗ corresponding author: tereza.zalska@fsv.cvut.cz Abstract. The present paper revisits a stress update procedure for the Hoek-Brown plasticity model enhanced by a Rankine type of tension cut-off failure criterion. Limiting tensile stresses not only supports the behavior of weak rock masses with a very low tensile strength but in combination with the original Hoek-Brown model improves robustness of the stress update procedure, i.e. the stress return mapping algorithm. Herein, the primary focus is on the stress return from a space of inadmissible trial stresses for which neither the Hoek-Brown yield surface nor its derivative can be evaluated. All potential scenarios are thoroughly discussed both in the framework of single- and multi-surface plasticity. The presented procedures were implemented and verified with the help of the Geo5 FEM software. Keywords: Hoek-Brown criterion, Rankine criterion, nonlinear yield function, rock mechanics, finite element method, tension cut-off. 1. Introduction In geological engineering, shear strength properties of the rock mass are mostly expressed by the parameters entering the formulation of the Hoek-Brown (HB) plasticity model. These can be transformed into the pair of basic shear strength parameters, i.e, cohesion c and friction angle φ, exploiting similarity between the Hoek-Brown and Mohr-Coulomb (MC) yield functions. Then the numerical analysis can be simply governed by the MC model. However, formulation of the HB model is based on the assumption of a nonlinear development of the rock mass strength with increasing loading. On the other hand, equivalent shear strength parameters derived by Hoek and Brown in [1] are independent of the current stress state. Therefore, this approach requires a proper estimation of the stress range for obtaining comparable results, see [2]. In light of this, a direct application of the HB model appears more suitable, which also promoted a recent implementation of the HB model into Geo5 FEM software [3]. In a numerical solution, violation of the yield func- tion is followed by a correction step bringing the stresses back to the yield surface. Formulation of a suitable stress return algorithm typically calls for the determination of the first derivative of a given yield function. Point out that the HB yield function is not defined for the maximum trial principal stress exceeding the tensile strength of the rock mass mak- ing the standard return mapping inapplicable. To overcome this obstacle, the present paper introduces a stress update algorithm that incorporates the Rankine failure criterion in conjunction with the HB model when moving from inadmissible to admissible stress space. Within this concept, the Rankine failure cri- terion may also represent a less or more significant reduction of the tensile strength. 2. Governing equations Application of the generalized Hoek-Brown criterion covers the whole range of rock masses from an intact one to a highly disturbed weak rock mass. The input parameters of the HB model are [4]: • uniaxial compressive strength σci and Hoek-Brown constant mi provided by series of triaxial testing, • geological strength index (GSI) used as a classifica- tion number of the geological quality of rock mass estimated from in situ geological observation, • disturbance factor D describing the rate of distur- bance by prior excavations or blasting, also esti- mated from in situ geological observation. The material constants mb, s and a govern the shape of the yield surface and are given by empirical equations as mb = mi exp GSI − 100 28 − 14D , (1) s = exp GSI − 100 9 − 3D , (2) a = 1 2 + 1 6 ( exp −GSI 15 − exp −20 3 ) . (3) 2.1. Hoek-Brown yield function Assuming a standard order of principle stresses σ1 > σ2 > σ3, the generalized yield function is defined as f13(σ1,σ3) = σ1 − σ3 − σci ( s − mb σ1 σci )a , (4) where σ1 and σ3 are the maximal and minimal ef- fective principal stresses, respectively. Unlike typical geotechnical notation, the compressive stresses are assumed negative and σ1 has a meaning of confining 117 https://doi.org/10.14311/APP.2023.40.0117 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en Tereza Žalská, Michal Šejnoha Acta Polytechnica CTU Proceedings pressure. Equation (4) represents the main sector of the yield surface in the principle stress space. For specific cases of triaxial compression (Lode’s angle θ = 30◦) and triaxial extension (θ = −30◦) two more sectors need to be considered. These are defined as f23(σ2,σ3) = σ2 − σ3 − σci ( s − mb σ2 σci )a , (5) f12(σ1,σ2) = σ1 − σ2 − σci ( s − mb σ1 σci )a . (6) In the principal stress space, the HB model plots as a irregular hexagonal pyramid with curved edges and surfaces, see Figure 1. Assuming the uniaxial tensile strength be equal to the biaxial one (σ3 = σ1 = σt) yields the tensile strength of the rock mass adopting Equation (4) as σt = sσci mb . (7) For the intact rock it holds GSI = 100 which gives a = 0.5. As the minimum value of GSI is equal to 0, the exponent a ranges from 0.5 to approximately 0.666. It is obvious that setting σmax > σt generates singularity in Eqautions (4)–(6). The yield condition is not defined beyond the apex and hence it cannot be evaluated. σm σ2 σ3 σ1 Figure 1. Hoek-Brown yield surface plotted in prin- cipal stress space. Expressing the yield function in terms of stress invariants, i.e., the mean stress σm, the deviatoric stress measure J and the Lode’s angle θ gives fHB (σm,J,θ) = (2J cos θ) 1/a − sσ1/aci + + ( cos θ − sin θ √ 3 ) mbJσ 1/a−1 ci + mbσmσ 1/a−1 ci . (8) Note that formulation (8) is not suitable for the com- plex return mapping procedure, as the derivative of the yield surface along the triaxial compression and triaxial extension edges (Lode’s angle θ = ±30◦, re- spectively) is singular, see [5] with reference to the MC model. 2.2. Plastic potential Plastic potential functions are formulated on the basis of the MC yield criterion as g13(σ1,σ3) = σ1 1 + sin ψ 1 − sin ψ − σ3, (9) g23(σ2,σ3) = σ2 1 + sin ψ 1 − sin ψ − σ3, (10) g12(σ1,σ2) = σ1 1 + sin ψ 1 − sin ψ − σ2, (11) where ψ is well known dilation angle. Constant rate of dilation is assumed in the present formulation. Vec- tors normal to the plastic potential surfaces governing the direction of the plastic yielding are therefore inde- pendent of the current stresses and take the form n1g = { 1 + sin ψ 1 − sin ψ ; 0; −1 } T, (12) n2g = { 0; 1 + sin ψ 1 − sin ψ ; −1 } T, (13) n3g = { 1 + sin ψ 1 − sin ψ ; −1; 0 } T. (14) More advanced formulations of the potential func- tions are also available in literature, see, e.g., Carranza- Torres and Fairhurst [6] who proposed the plastic po- tential function based on the mobilized value of the dilation angle ψmob depending on current stress state. However, such formulations go beyond the present scope. 2.3. Rankine yield function Yield surfaces of the Rankine type are often used as the additional tension cut-off surfaces. The model assumes the material to be isotropic and unlike the HB model, the stress return mapping is in general driven by an associated flow rule. The Rankine yield surfaces are defined as T1(σ1,σt) = σ1 − σt, (15) T2(σ2,σt) = σ2 − σt, (16) T3(σ3,σt) = σ3 − σt, (17) where σt ≤ σt is the prescribed tensile strength. The vectors normal to individual sectors of the yield surface are for the assumed associated flow rule given by n1t = {1; 0; 0} T, (18) n2t = {0; 1; 0} T, (19) n3t = {0; 0; 1} T. (20) Note that apart from a physically accepted reduction of a tensile stress σt, recall Equation (7), the return to the Rankine yield surface will be adopted as an intermediate step when an attempt to return to the HB yield function is made from an inadmissible trial 118 vol. 40/2023 Hoek-Brown plasticity model with tension cut-off stress space1. To avoid failing the return process we consider an upper limit of the tensile strength σmaxt = σt − σ∗ where σ∗ = 1 kPa is assumed in Geo5 FEM. For the sake of completeness we present the Rankine yield function, which plots as a regular triangular pyramid in the principal stress space, also in terms of stress invariants as fR(σm,J,θ) = 2 √ 3J cos ( θ + π 6 ) + + 3σm − 3σt. (21) 3. Stress update algorithm In the presented approach, the rock mass is treated as isotropic and no hardening is considered. Stress update algorithm respects the standard theory of plas- ticity. Elastic region of the strain is bounded by the yield surface. In case of plastic yielding both the yield and consistency condition must be satisfied. The Geo5 FEM software assumes the stress-strain rela- tion to be linearly elastic perfectly plastic for the HB model. The total strain increment decomposes into elastic and plastic parts as dε = dεel + dεpl. (22) The increment of stress relates to the elastic part of the total strain and follows the linear relation dσ = Del( dε − dεpl). (23) For non-associated plasticity the plastic part of the strain increment is defined as dεpl = ∆λ ∂g(σ) ∂σ , (24) with ∆λ being positive for plastic yielding, otherwise it is equal to 0. An associated plastic flow rule assumes g(σ) = f(σ). As already mentioned, a non-associated flow rule is applied when returning to the HB yield surface, while associated plasticity is adopted for the Rankine yield surface. Derivation of the plastic strain increment follows standard two step procedure employing elastic predic- tor and plastic corrector steps. The former step gives the trial stresses in the form σtr = σj−1 + Del dε. (25) Once the yield condition is violated, i.e., f(σtr ) > 0, the return mapping is carried out to bring the trial stresses back to the yield surface. Depending on the orientation of the elastic guess and the direction of the plastic strain increment, reordering of the principal stresses may occur, i.e. σ1 > σ2 > σ3 may no longer apply. As will be graphically shown in the next section, reordering of principal stresses requires formulation 1The trial stress space is assumed inadmissible whenever condition σmax > σt occurs regardless of whether or not the HB model is defined. of the return mapping algorithm in the framework of multi-surface plasticity. This is because the updated stress, although theoretically located on the yield surface, gives Lode’s angle which falls outside the limits −30◦ ≤ θ ≤ 30◦. In such a case, the correct stress lies on the triaxial edge, where two sectors of the yield surface intersect so both corresponding surfaces become active in return mapping. 3.1. Single yield surface return strategy Implicit return mapping is applied using the Newton- Raphson method in the form (∆λ)i+1 = (∆λ)i − (f13)i ( df13)i , (26) where df13 is the derivative of the yield function associated with the main sector and ∆λ is the plastic strain increment. The updated stress is provided by σi+1 = σtr − (∆λ)i+1Deln1g, (27) where the direction of the plastic strain increment n1g remains constant during the return mapping, re- call Equations (12)–(14). Writing the direction of principal plastic corrector as Deln1g = σ pl = { σ pl 1 ,σ pl 2 ,σ pl 3 } T (28) gives σ1(∆λ) = σtr1 − ∆λσ pl 1 , (29) σ2(∆λ) = σtr2 − ∆λσ pl 2 , (30) σ3(∆λ) = σtr3 − ∆λσ pl 3 . (31) Substituting Equations (29)–(31) into Equation (4) gives the yield function f13 in the form f13(∆λ) = (σtr1 − σ tr 3 ) − ∆λ(σ pl 1 − σ pl 3 )− − σci [ s − mb σci ( σtr1 − ∆λσ pl 1 )]a (32) and the corresponding derivative w.r.t. ∆λ as df13(∆λ) = σpl3 − σ pl 1 − − ambσ pl 1 [ s − mb σci ( σtr1 − ∆λσ pl 1 )]a−1 . (33) Return to the Rankine yield surface T1 follows the same minimization procedure. Owing to the linear form of the yield function, the correct increment of the plastic strain is found in one iteration step. Therefore ∆λ = − T1 dT1 . (34) The tensile yield surface defined in terms of the plastic multiplier is given by T1 = σtr1 − ∆λσ pl 1 − σt (35) 119 Tereza Žalská, Michal Šejnoha Acta Polytechnica CTU Proceedings σ1 σ2 σ3 0 50 100 150 200 0 20 40 60 80 100 f HB f R f HBf R -σm (kPa) J ( k P a ) f R f HB -σ1 (kPa) -σ 3 ( k P a ) σm Figure 2. Non-associated return to the HB surface from trial stress violating T1 in two consecutive steps (ψ = 0): first - from trial stress (red point) onto the Rankine yield surface (orange point), second - from Rankine surface to HB surface (green point) plotted in the deviatoric plane (left), σm − J plane (middle), σ1 − σ3 plane (right). and the derivative w.r.t. ∆λ is therefore constant dT1 = −σpl1 . (36) Returning from the trial stress space violating the Rankine failure criterion considers two scenarios: • Correct return to the HB yield surface is expected: This is a two-step return where returning to the Rankine yield surface is performed first assuming the direction of plastic corrector given by the HB model, i.e., σpl = Deln1g (non-associated flow rule). The procedure then continues with the second step bringing the stresses back onto the HB yield sur- face along the same direction. The total plastic strain increment is then the sum of plastic strain increments generated by both steps. See Figure 2. • Correct return to the Rankine yield surface is ex- pected: This is a one-step return along the direction given by σpl = Deln1t (associated flow rule). 3.2. Triaxial lines return strategy Two active sectors of the HB yield surface apply to both the state of triaxial compression and triaxial extension. Considering triaxial stress state in com- pression for the determination of the multi-surface plasticity return, two yield functions, f13 and f23, are active and the iterative scheme takes the form{ ∆λ1 ∆λ2 }i+1 = { ∆λ1 ∆λ2 }i −(H−1)i { f13 f23 }i , (37) where the Jacobian matrix is composed of the partial derivatives of the yield functions with respect to two plastic multipliers ∆λ1 and ∆λ2. Therefore H =   ∂f13 ∂(∆λ1) ∂f13 ∂(∆λ2) ∂f23 ∂(∆λ1) ∂f23 ∂(∆λ2)   . (38) The vector of updated principle stresses for the current iterative step is given by σi+1 = σtr − (∆λ1)i+1Deln1g − (∆λ2) i+1Deln2g. (39) Vectors providing the orientation of the plastic correc- tor are defined as σ pl 1 = D eln1g = { σ pl 11; σ pl 12; σ pl 13 }T , (40) σ pl 2 = D eln2g = { σ pl 21; σ pl 22; σ pl 23 }T (41) and individual components of the stress vector for j = 1 → 3 σj (∆λ1, ∆λ2) = σtrj − ∆λ1σ pl 1j − ∆λ2σ pl 2j. (42) Expressing the yield function in terms of two plastic multipliers provides the yield function f13 as f13(∆λ1, ∆λ2) = = (σtr1 − σ tr 3 ) − ∆λ1(σ pl 11 − σ pl 13) − ∆λ2(σ pl 21 − σ pl 23)− − σci [ s − mb σci ( σtr1 − ∆λ1σ pl 11 − ∆λ2σ pl 21 )]a (43) and its partial derivatives in the form H11 = σ pl 13 − σ pl 11 − ambσ pl 11[ s − mb σci ( σtr1 − ∆λ1σ pl 11 − ∆λ2σ pl 21 )]a−1 , (44) H12 = σ pl 23 − σ pl 21 − ambσ pl 21[ s − mb σci ( σtr1 − ∆λ1σ pl 11 − ∆λ2σ pl 21 )]a−1 . (45) The other yield function and components of the H matrix would be expressed similarly and the same approach is applied to the state of triaxial extension with sectors f13 and f12 being simultaneously active. Attention should be also paid to the triaxial corner of the tensile yield surface. In this case, the proposed approach is adopted for the intersection of the yield functions T1 and T2 and, similarly to a single surface plasticity with one iteration step needed. Again, both non-associated and associated plasticity is considered in return mapping scheme depending on whether the stress is expected to be brought back to HB model or directly to the Rankine yield surface, recall the last paragraph in the previous section. As an example of the former case, see Figure 3. 120 vol. 40/2023 Hoek-Brown plasticity model with tension cut-off Figure 3. Demonstration of the return to the triaxial corners of the HB yield surface from trial stress vio- lating tensile yield condition (ψ = 0 is assumed): σtr violates T1 condition (TC corner) σtr ′ violates both T1, T2 conditions (TC corner) σtr∗ violates T1 condition (TE corner). 4. Complex return mapping scheme Here, only the return from an inadmissible trial stress space is discussed. As already explained in the pre- vious sections, a direct return of the trial stress for the condition σtr1 > σt onto the HB yield surface fails as the HB yield function is undefined. The re- turn mapping therefore begins with evaluating the tensile condition (15) employing the reduced tensile strength σt ≤ σt − σ∗. Violating this condition calls a non-associated return onto the Rankine yield func- tion controlled by the plastic potential formulated for the HB model. If the stress point violates the HB criterion, a non-associated return from the point on the Rankine surface onto the HB yield surface fol- lows, see Figure 2. Otherwise, associated plasticity is adopted for bringing the original trial stress back onto the Rankine tensile yield surface, recall the last paragraph in Section 3.1. Special attention should be paid to the intersection of T1 and f13 surfaces with reference to Figure 5. The proposed algorithm first predicts that the return mapping will continue to the HB surface after being preceded by a non-associated return to the tensile yield surface. If no violation of the HB criterion is detected after the first step, the associated return of the trial stress to the tensile surface is called. Af- ter this attempt, violation of the HB condition may eventually occur. If it does, the HB criterion is substi- tuted with the linear function fsub, see Equation (46), which satisfies T1 condition (σ1 = σt), and the trial stress is shifted to the intersection of two planes using multi-surface plasticity. fsub(σ3) = σt − σ3 − σci ( s − mb σt σci )a (46) It is appropriate to point out that reordering of the principle stresses, and thus detection of the triaxial corner return, may occur during the described steps. In such a case, the above two-step procedure would also be performed in the framework of multi-surface plasticity. A graphical representation is evident in Figure 3. The complete return mapping scheme in- cluding all potential scenarios is then summarized in Figure 4. Figure 4. Full return mapping scheme including all potential scenarios. 121 Tereza Žalská, Michal Šejnoha Acta Polytechnica CTU Proceedings σ1 σ2 σ3 σ1 σ2 σ3 σ1 σ2 σ3 10 0 10 20 0 5 10 15 20 25 10 0 10 20 0 5 10 15 20 25 10 0 10 20 30 0 5 10 15 20 25 f HB f R -σm (kPa) J ( k P a ) J ( k P a ) -σm (kPa) -σm (kPa) J ( k P a ) Delng1 D eln t 1 Figure 5. Sequence of the return attempts leading to correct stress lying on the intersection of T1 and f13 (ψ = 0) plotted in the deviatoric planes (top) and σm − J planes (bottom). 5. Conclusions A detailed description of the return mapping algorithm for the case of inadmissible trial stresses, here defined as the stress space for which σtr1 > σt, is presented. To the authors knowledge the number of contribu- tions addressing this issue is limited. In [7] a smooth variant of the HB model is presented and the tension cut-off condition is introduced via an allowable tensile mean stress instead of actual tensile strength. In [8, 9] the authors also deal with the issue of tensile stresses. However, the approach is not explained in all details and its effectivity is not obvious. The present approach exploits simplicity of the Rankine yield surface, which is often used as a ten- sion cut-off criterion. The return of inadmissible trial stresses to the main sector or triaxial edges of the HB yield surface assumes two consecutive steps. First, the stress is brought back to the Rankine yield surface defined by the reduced tensile strength. Second, the implicit iterative procedure for shifting the stress onto the HB surface is applied. Both steps assume the same orientation of the plastic corrector. It is shown, that the prediction of all possible situations, which can lead to exceeding the tensile strength and failing the solution, results in the complex algorithm further complicated by the assumption of associated plasticity for tensile yielding. This assumption becomes active once the two-step return procedure onto the HB sur- face fails and the trial stress is then brought back directly to the Rankine yield surface. This algorithm was implemented into Geo5 FEM software and thoroughly tested. Acknowledgements The support provided by the SGS project No. SGS22/030/OHK1/1T/11 and by the GAČR project No. 22-12178S is gratefully acknowledged. References [1] E. Hoek, C. Carranza-Torres, B. Corkum. Hoek-Brown failure criterion – 2002 edition. In Proceedings of the 5th North American Rock Mechanics Symposium – NARMS-TAC, pp. 267–273. 2002. [2] T. Poklopová, M. Šejnoha. Comparing the Hoek-Brown and Mohr-Coulomb failure criteria in FEM analysis. In Acta Polytechnica CTU Proceedings, vol. 30, pp. 69–75. 2021. https://doi.org/10.14311/APP.2021.30.0069 [3] Fine-Ltd. Geo5 FEM, 2022. [2022-09-13]. http://www.fine.cz/ [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. https://doi.org/10.1016/S1365-1609(97)80069-X [5] J. Huang, D. V. Griffiths. Observations on return mapping algorithms for piecewise linear yield criteria. International Journal of Geomechanics 8(4):253–265, 2008. https: //doi.org/10.1061/(ASCE)1532-3641(2008)8:4(253) [6] C. Carranza-Torres, C. Fairhurst. The elasto-plastic response of underground excavations in rock masses that satisfy the Hoek-Brown failure criterion. International Journal of Rock Mechanics and Mining Sciences 36(6):777–809, 1999. https://doi.org/10.1016/S0148-9062(99)00047-9 [7] R. G. Wan. Implicit integration algorithm for Hoek-Brown elastic-plastic model. Computers and Geotechnics 14(3):149–177, 1992. https://doi.org/10.1016/0266-352X(92)90031-N [8] J. Clausen. Efficient Non-Linear Finite Element Implementation of Elasto-Plasticity for Geotechnical Problems. Ph.D. thesis, Aalborg Universitet, 2007. [9] J. Clausen, L. Damkilde. An exact implementation of the Hoek-Brown criterion for elasto-plastic finite element calculations. International Journal of Rock Mechanics and Mining Sciences 45(6):831–847, 2008. https://doi.org/10.1016/j.ijrmms.2007.10.004 122 https://doi.org/10.14311/APP.2021.30.0069 http://www.fine.cz/ https://doi.org/10.1016/S1365-1609(97)80069-X https://doi.org/10.1061/(ASCE)1532-3641(2008)8:4(253) https://doi.org/10.1061/(ASCE)1532-3641(2008)8:4(253) https://doi.org/10.1016/S0148-9062(99)00047-9 https://doi.org/10.1016/0266-352X(92)90031-N https://doi.org/10.1016/j.ijrmms.2007.10.004 Acta Polytechnica CTU Proceedings 40:117–122, 2023 1 Introduction 2 Governing equations 2.1 Hoek-Brown yield function 2.2 Plastic potential 2.3 Rankine yield function 3 Stress update algorithm 3.1 Single yield surface return strategy 3.2 Triaxial lines return strategy 4 Complex return mapping scheme 5 Conclusions Acknowledgements References