RETURN MAPPING SCHEME FOR THE HOEK-BROWN MODEL WITH TENSION CUT-OFF

. 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.


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 function 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 making 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 criterion may also represent a less or more significant reduction of the tensile strength.

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 m i provided by series of triaxial testing, • geological strength index (GSI) used as a classification number of the geological quality of rock mass estimated from in situ geological observation, • disturbance factor D describing the rate of disturbance by prior excavations or blasting, also estimated from in situ geological observation.
The material constants m b , s and a govern the shape of the yield surface and are given by empirical equations as

Hoek-Brown yield function
Assuming a standard order of principle stresses σ 1 > σ 2 > σ 3 , the generalized yield function is defined as where σ 1 and σ 3 are the maximal and minimal effective principal stresses, respectively. Unlike typical geotechnical notation, the compressive stresses are assumed negative and σ 1 has a meaning of confining 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 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 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. 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 Note that formulation (8) is not suitable for the complex return mapping procedure, as the derivative of the yield surface along the triaxial compression and triaxial extension edges (Lode's angle θ = ±30 • , respectively) is singular, see [5] with reference to the MC model.

Plastic potential
Plastic potential functions are formulated on the basis of the MC yield criterion as where ψ is well known dilation angle. Constant rate of dilation is assumed in the present formulation. Vectors normal to the plastic potential surfaces governing the direction of the plastic yielding are therefore independent of the current stresses and take the form More advanced formulations of the potential functions are also available in literature, see, e.g., Carranza-Torres and Fairhurst [6] who proposed the plastic potential function based on the mobilized value of the dilation angle ψ mob depending on current stress state. However, such formulations go beyond the present scope.

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 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 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 vol. 40/2023 Hoek-Brown plasticity model with tension cut-off stress space 1 . To avoid failing the return process we consider an upper limit of the tensile strength 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

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 plasticity. 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 relation to be linearly elastic perfectly plastic for the HB model. The total strain increment decomposes into elastic and plastic parts as The increment of stress relates to the elastic part of the total strain and follows the linear relation For non-associated plasticity the plastic part of the strain increment is defined as 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 predictor and plastic corrector steps. The former step gives the trial stresses in the form 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 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.

Single yield surface return strategy
Implicit return mapping is applied using the Newton-Raphson method in the form where df 13 is the derivative of the yield function associated with the main sector and ∆λ is the plastic strain increment. The updated stress is provided by where the direction of the plastic strain increment n 1 g remains constant during the return mapping, recall Equations (12)-(14). Writing the direction of principal plastic corrector as gives Substituting Equations (29)-(31) into Equation (4) gives the yield function f 13 in the form and the corresponding derivative w.r.t. ∆λ as Return to the Rankine yield surface T 1 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 The tensile yield surface defined in terms of the plastic multiplier is given by and the derivative w.r.t. ∆λ is therefore constant 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 = D el n 1 g (non-associated flow rule). The procedure then continues with the second step bringing the stresses back onto the HB yield surface 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 expected: This is a one-step return along the direction given by σ pl = D el n 1 t (associated flow rule).

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 compression for the determination of the multi-surface plasticity return, two yield functions, f 13 and f 23 , are active and the iterative scheme takes the form where the Jacobian matrix is composed of the partial derivatives of the yield functions with respect to two plastic multipliers ∆λ 1 and ∆λ 2 . Therefore The vector of updated principle stresses for the current iterative step is given by Vectors providing the orientation of the plastic corrector are defined as and individual components of the stress vector for j = 1 → 3 Expressing the yield function in terms of two plastic multipliers provides the yield function f 13 as and its partial derivatives in the form 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 f 13 and f 12 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 T 1 and T 2 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.

Complex return mapping scheme
Here, only the return from an inadmissible trial stress space is discussed. As already explained in the previous sections, a direct return of the trial stress for the condition σ tr 1 > σ t onto the HB yield surface fails as the HB yield function is undefined. The return 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 function 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 follows, 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 T 1 and f 13 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. After this attempt, violation of the HB condition may eventually occur. If it does, the HB criterion is substituted with the linear function f sub , see Equation (46), which satisfies T 1 condition (σ 1 = σ t ), and the trial stress is shifted to the intersection of two planes using multi-surface plasticity.
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 including all potential scenarios is then summarized in Figure 4.

Conclusions
A detailed description of the return mapping algorithm for the case of inadmissible trial stresses, here defined as the stress space for which σ tr 1 > σ t , is presented. To the authors knowledge the number of contributions 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 tension 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 surface 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.