# A Return Mapping Algorithm for Nonlinear Yield Criteria with the Equivalent Mohr–Coulomb Strength Parameters

Authors: Chao Hu; Fei Ye; Wenxi Fu; Zhuo Chen
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

This paper proposes a modified return mapping algorithm for a series of nonlinear yield criteria. The algorithm is established in the principal stress space and ignores the effect of the intermediate principal stress. Three stress return schemes are derived in this paper: return to the yield surface, return to the curve, and return to the apex point. The conditions used for determining the correct stress return type are also constructed. After the proposed algorithm is programmed in the finite element software, we merely need the equivalent Mohr–Coulomb (M-C) strength parameters, the derivatives of their functions, and the tensile strength of these nonlinear yield criteria. In addition, the Hoek–Brown (H-B) yield criterion is taken as an example to validate the proposed method. The results show that the updated stresses and the final principal stresses obtained by the proposed method are in good agreement with those obtained by other methods. Furthermore, the proposed method is more suitable for the associated plastic-flow rule.


---

## Body

## 1. Introduction

The stress update is a vitally important ingredient to return the predictor stress to the yield surface in the elastoplastic finite element computation. In recent decades, numerous stress update algorithms have been proposed. These algorithms can be mainly classified into three categories: explicit integration algorithms [1–3], return mapping algorithms [4–6], and exact stress integration algorithms [7, 8]. Traditionally, the return mapping algorithm, also called the predictor-corrector scheme, is constructed in the six-dimensional stress space and is usually cumbersome since it requires that the second derivative of the plastic potential and nondifferentiable yield surface be smoothed. However, these disadvantages are well overcome by constructing return mapping algorithms in the principal stress space [9–12]. The geometric features of the yield surfaces in the principal stress space can be graphically visualized and are easily and exactly used for determining the return position of updated stress, thereby avoiding the smoothing of nondifferentiable yield surfaces, as other methods involve [12–18].At present, scholars have proposed many return mapping algorithms in the principal stress space for different yield criteria. Clausen et al. [10, 11] proposed a kind of implicit integration algorithm for the linear yield criteria with multiple yield planes built upon the works of Pankaj and Bićanić [19]. Lin and Li [20] presented a return mapping algorithm for unified strength theory; this algorithm is suitable not only for the unified strength theory but also for other criteria (e.g., the Mohr–Coulomb (M-C) yield criterion and the Tresca yield criterion) due to the flexibility of unified strength theory. Fang et al. [21] proposed an implicit numerical integration algorithm for the unified yield criterion in the principal stress space; this algorithm is more useful and convenient than the algorithm proposed by Lin and Li [20], as linear isotropic hardening and nonassociated plasticity can be considered. For generalized Hoek–Brown (H-B) plasticity, Clausen et al. [13] proposed a simple stress-update algorithm in the principal stress space. Sørensen et al. [22] developed a model for the nonassociated H-B material, in which strain hardening can be considered. Moreover, Karaoulanis [23] proposed a return mapping algorithm for nonsmooth multisurface plasticity based on a spectral representation of stresses and strains, and the algorithm can be used for the M-C and Drucker–Prager yield criteria. Many alternative return mapping algorithms can be applied to stress updating in the finite element computation; however, there is still a troublesome problem. For any yield criterion (e.g., the H-B yield criterion [24–26], Christensen yield criterion [27–31], or the parabolic-type yield criterion [32, 33]) that is not built-in finite element software, a corresponding return mapping algorithm is needed when the yield criterion is adopted in the elastoplastic finite element computation. Therefore, we need to code the algorithm first to implement the yield criterion in the finite element software when the yield criterion is not built-in, which is usually a difficult and time-consuming task.In this paper, a modified return mapping algorithm is proposed for a series of nonlinear yield criteria. This algorithm is established in the principal stress space, and the effect of the intermediate principal stress has not been taken into account. In the derivation of the algorithm, isotropic linear elasticity and perfect plasticity are assumed, and the equivalent strength parameters are adopted. Essentially, the equivalent strength parameters are scalar functions of the principal stresses [34]. Once the principal stresses are determined, there should be an intersection line or point defined by the primary yield surface of the M-C yield criterion and that of the nonlinear yield criteria in the principal stress space. Figures 1 and 2 show an example of the H-B yield criterion. If the predictor stress is returned to the line or the point, the updated stress obtained by different algorithms should be identical. When the modified return mapping algorithm code has been programmed in the finite element software, we merely need the equivalent M-C strength parameters, the derivatives of their functions, and the tensile strength of these nonlinear yield criteria, and then, these criteria can be quickly implemented. In addition, the H-B yield criterion is taken as an example to validate the proposed method.Figure 1
The envelope of the H-B yield criterion in the Mohr diagram and the parameters.Figure 2
The yield surface of H-B and M-C yield criterion obeyingσ1≥σ2≥σ3 in the principal stress space.


## 2. Equivalent Mohr–Coulomb Strength Parameters

The M-C yield criterion in the principal space is expressed as follows:(1)f=kσ1−σ3−n=0,where k=1−sinϕ/1+sinϕ, n=2ccosϕ/1+sinϕ, σ1, σ2, and σ3 denote the major principal stress, the intermediate principal stress, and the minor principal stress, respectively, and σ1≥σ2≥σ3. Note that compressive stresses are considered to be positive. The plastic potential is chosen as follows:(2)g=mσ1−σ3−p=0,where m=1−sinψ/1+sinψ and p=2κcosψ/1+sinψ. The plastic potential resembles the shape of the yield surface. ψ and κ are set to be identical to φ and c, respectively, for associated material behaviour and otherwise for nonassociated material behaviour.A general methodology to determine the equivalent M-C strength parameters was proposed by Balmer [34]. The equivalent parameters can be given by(3)ϕi=sin−1∂σ1/∂σ3−1∂σ1/∂σ3+1,(4)ci=σ1−σ3∂σ1∂σ3−σ3+σ1−σ3∂σ1/∂σ3tanϕi.Analogously, the equivalent parameters of plastic potentials can be obtained by the approach proposed by Balmer.


## 3. Return Mapping in the Principal Stress Space

In the previous sections, the approach for obtaining the equivalent M-C strength parameters and the equivalent parameters of the plastic potentials has been presented. In this section, the return mapping algorithm is established based on the principal stress space. As all manipulations of the algorithm are carried out in the principal stress space [10, 11], the components of the vectors are expressed here via the principal stress, and unless otherwise stated, the last three components of these vectors are equal to zero and are omitted.When the predictor stress,σB=σ1Bσ2Bσ3BT, is located outside of the yield surface of M-C yield criterion, it must be brought back to the yield surface. Essentially, there are three types of stress returns [10, 11, 20]:(1)
Return to the yield surface(2)
Return to the curve(3)
Return to the apex


### 3.1. Return to the Yield Surface

The crucial procedure to bringσB back to the yield surface is to evaluate the plastic corrector, Δσp.Δσp was derived by Crisfield [35] and is written in the principal stress space as follows:(5)Δσp=fσBDbaTDb=fσBrpwithrp=DbaTDb,where rp is the direction of the plastic corrector, fσB represents the magnitude of the plastic corrector, a and b are the gradients of the yield surface and the plastic potential, respectively, and D is a matrix of elastic constants for linearly elastic isotropic materials. a, b, and D are given by(6)a=∂F∂σ=∂kϕi,ci∂σ1σ1+kϕi,ci−∂nϕi,ci∂σ10∂kϕi,ci∂σ3σ1−1−∂nϕi,ci∂σ3T,b=∂G∂σ=∂mψi,κi∂σ1σ1+mψi,κi−∂pψi,κi∂σ10∂mψi,κi∂σ3σ1−1−∂pψi,κi∂σ3T,D=E1+μ1−2μ1−μμμμ1−μμμμ1−μ,where E is Young’s modulus and μ is Poisson’s ratio.OnceΔσp is derived, the updated stress, σC=σ1Cσ2Cσ3CT, can be given by the following formula:(7)σC=σB−Δσp.As equation (7) is a nonlinear equation, an iterative method must be utilized to find σC. Because the slope of a line that connects σC with σB and the slope of the direction of rp must be identical [9], the following equation can be obtained:(8)hσC=σ3B−σ3Cσ1B−σ1C−μ∂mψi,κi/∂σ1σ1+mψi,κi−∂pψi,κi/∂σ1+1−μ∂mψi,κi/∂σ3σ1−1−∂pψi,κi/∂σ31−μ∂mψi,κi/∂σ1σ1+mψi,κi−∂pψi,κi/∂σ1+μ∂mψi,κi/∂σ3σ1−1−∂pψi,κi/∂σ3=0.Incorporating equation (1), an iterative scheme can be constructed on the basis of the Newton–Raphson method:(9)σi+1C=σiC−H′σiC−1HσiC,where σC can be rewritten as σC=σ1Cσ3CT because the updated intermediate principal stress need not be considered; H=hfT, where H′ is the Jacobian matrix, σiC is the current value at the iteration step, i, and σi+1C, which is calculated by equation (9) in the next value at the iteration step, i+1. Iterations are performed until(10)H<ε,where ε is the tolerance. If a σC satisfying equation (10) is found, then the updated intermediate principal stress, σ2C, can be found by inserting σ1C and σ3C into equation (7).The iterative procedure of equation (9) must start with an initial value, σ0C. Usually, a simple method can generate a suitable initial value, i.e., σ0C=σB−Δ, where ∆ is a very small positive number to ensure that the σ1B−σ1C term in equation (8) does not vanish. In addition, it should be noted that the σiC in the iteration process should be greater than σB. In the case in which the algorithm does not converge, a new corrected σi+1C can be estimated by(11)σi+1C=σB−σiC−σB2.


### 3.2. Return to the Curve

IfσC is located on the curve defined by the intersection between two yield planes, then the line connecting the updated stress with the predictor stress is perpendicular to the vertical direction of the plane defined by the direction of the plastic corrector of two adjacent yield surfaces. Therefore, if the predictor stress is returned to the curve (satisfying σ1=σ2), defined by the intersection of the primary yield plane and the adjacent yield surface satisfying σ2≥σ1≥σ3 in the principal stress space, the following equation can be obtained:(12)ℓ=rp×r1pTσC−σB=0,where ‘‘×’’represents the cross product between rP and r1p and r1p is the plastic corrector of the adjacent yield surface. Note that the gradient and the plastic potential of the adjacent yield surface can be easily obtained by interchanging the components of a and b in equation (6). Since σC is located on both yield surfaces, it must also satisfy the following equations:(13)fσC=0,f1σC=0,where f1 is the adjacent yield surface. Obviously, σC cannot be immediately determined by simultaneous solution of equations (12) and (13) because these equations are nonlinear. To find σC, which is analogous to equation (9), the following iterative scheme can be established:(14)σi+1C=σiC−M′σiC−1MσiC,where M=ℓfT and M′ is the Jacobian matrix. Iterations are performed until M<ε, and then, σ2C can be found by σ2C=σ1C.The approach for obtaining an initial value in equation (14) is the same as that of the return to the yield surface case. In addition, σiC in the iteration process should also be greater than σB. Otherwise, the new σi+1C is estimated by equation (11).IfσB should be brought back into the curve (satisfying σ2=σ3), defined by the intersection of the primary yield plane and the adjacent yield plane satisfying σ1≥σ3≥σ2, the method of finding σC is completely analogous to the return of the predictor stress to the curve (satisfying σ1=σ2). Therefore, it is not necessary to describe the method again here.


### 3.3. Return to the Apex

IfσB returns to the apex point, its spatial position can be directly determined because the point can be determined by the intersection of two or more curves, and a calculation is not needed [9, 10]. σC is given by(15)σC=σt,where σt=σtσtσtT. The tensile strength σt can easily be obtained by setting σ1=σ3=σt in the criterion and resolving the equation.


## 3.1. Return to the Yield Surface

The crucial procedure to bringσB back to the yield surface is to evaluate the plastic corrector, Δσp.Δσp was derived by Crisfield [35] and is written in the principal stress space as follows:(5)Δσp=fσBDbaTDb=fσBrpwithrp=DbaTDb,where rp is the direction of the plastic corrector, fσB represents the magnitude of the plastic corrector, a and b are the gradients of the yield surface and the plastic potential, respectively, and D is a matrix of elastic constants for linearly elastic isotropic materials. a, b, and D are given by(6)a=∂F∂σ=∂kϕi,ci∂σ1σ1+kϕi,ci−∂nϕi,ci∂σ10∂kϕi,ci∂σ3σ1−1−∂nϕi,ci∂σ3T,b=∂G∂σ=∂mψi,κi∂σ1σ1+mψi,κi−∂pψi,κi∂σ10∂mψi,κi∂σ3σ1−1−∂pψi,κi∂σ3T,D=E1+μ1−2μ1−μμμμ1−μμμμ1−μ,where E is Young’s modulus and μ is Poisson’s ratio.OnceΔσp is derived, the updated stress, σC=σ1Cσ2Cσ3CT, can be given by the following formula:(7)σC=σB−Δσp.As equation (7) is a nonlinear equation, an iterative method must be utilized to find σC. Because the slope of a line that connects σC with σB and the slope of the direction of rp must be identical [9], the following equation can be obtained:(8)hσC=σ3B−σ3Cσ1B−σ1C−μ∂mψi,κi/∂σ1σ1+mψi,κi−∂pψi,κi/∂σ1+1−μ∂mψi,κi/∂σ3σ1−1−∂pψi,κi/∂σ31−μ∂mψi,κi/∂σ1σ1+mψi,κi−∂pψi,κi/∂σ1+μ∂mψi,κi/∂σ3σ1−1−∂pψi,κi/∂σ3=0.Incorporating equation (1), an iterative scheme can be constructed on the basis of the Newton–Raphson method:(9)σi+1C=σiC−H′σiC−1HσiC,where σC can be rewritten as σC=σ1Cσ3CT because the updated intermediate principal stress need not be considered; H=hfT, where H′ is the Jacobian matrix, σiC is the current value at the iteration step, i, and σi+1C, which is calculated by equation (9) in the next value at the iteration step, i+1. Iterations are performed until(10)H<ε,where ε is the tolerance. If a σC satisfying equation (10) is found, then the updated intermediate principal stress, σ2C, can be found by inserting σ1C and σ3C into equation (7).The iterative procedure of equation (9) must start with an initial value, σ0C. Usually, a simple method can generate a suitable initial value, i.e., σ0C=σB−Δ, where ∆ is a very small positive number to ensure that the σ1B−σ1C term in equation (8) does not vanish. In addition, it should be noted that the σiC in the iteration process should be greater than σB. In the case in which the algorithm does not converge, a new corrected σi+1C can be estimated by(11)σi+1C=σB−σiC−σB2.


## 3.2. Return to the Curve

IfσC is located on the curve defined by the intersection between two yield planes, then the line connecting the updated stress with the predictor stress is perpendicular to the vertical direction of the plane defined by the direction of the plastic corrector of two adjacent yield surfaces. Therefore, if the predictor stress is returned to the curve (satisfying σ1=σ2), defined by the intersection of the primary yield plane and the adjacent yield surface satisfying σ2≥σ1≥σ3 in the principal stress space, the following equation can be obtained:(12)ℓ=rp×r1pTσC−σB=0,where ‘‘×’’represents the cross product between rP and r1p and r1p is the plastic corrector of the adjacent yield surface. Note that the gradient and the plastic potential of the adjacent yield surface can be easily obtained by interchanging the components of a and b in equation (6). Since σC is located on both yield surfaces, it must also satisfy the following equations:(13)fσC=0,f1σC=0,where f1 is the adjacent yield surface. Obviously, σC cannot be immediately determined by simultaneous solution of equations (12) and (13) because these equations are nonlinear. To find σC, which is analogous to equation (9), the following iterative scheme can be established:(14)σi+1C=σiC−M′σiC−1MσiC,where M=ℓfT and M′ is the Jacobian matrix. Iterations are performed until M<ε, and then, σ2C can be found by σ2C=σ1C.The approach for obtaining an initial value in equation (14) is the same as that of the return to the yield surface case. In addition, σiC in the iteration process should also be greater than σB. Otherwise, the new σi+1C is estimated by equation (11).IfσB should be brought back into the curve (satisfying σ2=σ3), defined by the intersection of the primary yield plane and the adjacent yield plane satisfying σ1≥σ3≥σ2, the method of finding σC is completely analogous to the return of the predictor stress to the curve (satisfying σ1=σ2). Therefore, it is not necessary to describe the method again here.


## 3.3. Return to the Apex

IfσB returns to the apex point, its spatial position can be directly determined because the point can be determined by the intersection of two or more curves, and a calculation is not needed [9, 10]. σC is given by(15)σC=σt,where σt=σtσtσtT. The tensile strength σt can easily be obtained by setting σ1=σ3=σt in the criterion and resolving the equation.


## 4. Determination of the Correct Stress Return

The stress areas are separated by boundary planes,P1 = 0 and P2 = 0, at the apex point:(16)P1σ=rp×r1pTσ−σt=0,P2σ=r2p×rpTσ−σt=0,where r2p is the plastic corrector of the adjacent yield plane satisfying σ1≥σ3≥σ2. Here, r1p and r2p are evaluated at the apex point. If σB is located outside the yield plane (i.e., fσB>0), an apex return is valid when(17)P1σB≥0,P2σB≥0.In the case of some associated plasticity,ψ equals 90° and κ reaches infinity at the apex point, such that r1p and r2p calculated at the apex tend towards infinity. In this case, the vectors r1p and r2p can be evaluated as follows:(18)rp=1−μμμT,r1p=μ1−μμT,r2p=1−μμμT.IfσB should not be returned to the apex point (i.e., the conditions in equation (17) are not satisfied), a surface return first is executed as outlined in Section 3.1. The return is suitable if the components of σC satisfy the following:(1)
Ifσ1C≥σ2C≥σ3C, σB returns to the surface(2)
Ifσ2C≥σ1C, σB returns to the curve (satisfying σ1=σ2) through the procedure outlined in Section 3.2(3)
Ifσ2C≤σ3C, σB returns to the curve (satisfying σ2=σ3) through the procedure outlined in Section 3.2With the stress return schemes and the conditions for choosing the proper return schemes at hand, it is necessary to determine a consistent matrix. The general method for calculating the consistent matrix in the principal stress space has been derived by Clausen et al. [9, 10]. The details, therefore, need not be given here.


## 5. Application Example

In the following, the proposed method is applied to the H-B yield criterion as an application example. Considering that the proposed method needsσt and the equivalent M-C strength parameters of the H-B yield criterion, the equivalent M-C strength parameters can be obtained by inserting the H-B yield criterion [26] into equations (3) and (4). The equivalent parameters are given by(19)ϕi=arcsinαmbs+mbσ3/σciα−12+αmbs+mbσ3/σciα−1,ci=s+mbσ3/σciα−1σcis+mbσ3/σci−αmbσ321+αmbs+mbσ3/σciα−1,where σci is the uniaxial compressive strength and s, mb, and α are strength parameters of the H–B yield criterion. σt can easily be obtained by setting σ1=σ3=σt in the H-B yield criterion [26]:(20)σt=−sσcimb.To assess the accuracy and validity of the proposed method, three stress points are selected to implement the analysis of stress return. Points A (1119504.4, 30000, and 10000), B (1119504.4, 1119504.4, and 10000), and C (1119504.4, 10000, and 10000) are located on the yield surface, on the yield curve satisfyingσ1=σ2 and on the yield curve satisfying σ2=σ3, respectively. The magnitude of the elastic trial stress increment is given as 30000, 20000, and 10000. Notably, this increment ensures an elastic trial stress outside of the yield surface. The adopted parameters of the H-B yield criterion are σci = 2 MPa, s = 0.3, mb = 1.55, and α = 0.5. In addition, the deformation parameters E = 50 × 109 Pa and μ = 0.2 are adopted. The different parameters mbg = 0 and αg = 1 are utilized for the nonassociated case. For comparison purposes, the approach proposed by Clausen and Damkilde [9] is applied to the analysis. The relative error of the updated stress obtained by the two methods is given by(21)relative error=σ∗−σσ∗×100%,where σ∗ is the principal stress obtained by Clausen’s method and σ is the principal stress obtained by the proposed method. The return values calculated by the two methods are listed in Tables 1 and 2. As shown in these tables, the results obtained by the two methods are in good agreement. The very small errors are mainly caused by the round-off errors in the numerical calculation, which are usually unavoidable.Table 1
The values of the updated stress obtained by the proposed method and Clausen’s method and the relative errors (associated case).
The predictor stressThe proposed methodClausen’s methodRelative error (%)Updated stressPoint A1149504.41149049.51149049.505000050385.350385.302000022381.222381.201179504.41178556.11178556.107000070781.470781.403000034855.334855.30Point B1149504.41150206.11150206.101139504.41150206.11150206.102000022868.122868.101179504.41180902.61180902.601159504.41180902.61180902.603000035851.935851.90Point C1149504.41159399.71159399.703000026744.526744.502000026744.526744.501179504.41199245.21199245.205000043665.743665.703000043665.743665.70Table 2
The values of the updated stress obtained by the proposed method and Clausen’s method and the relative errors (nonassociated case).
The predictor stressThe proposed methodClausen’s methodRelative error (%)Updated stressPoint A1149504.41147693.81147693.805000050000.050000.002000021810.621810.601179504.41175813.21175813.207000070000.070000.003000033691.233691.20Point B1149504.41152139.61152139.601139504.41152139.61152139.602000023682.423682.401179504.41184630.81184630.801159504.41184630.81184630.803000037436.837436.80Point C1149504.41156701.51156701.503000025605.725605.702000025605.725605.701179504.41193827.61193827.605000041353.641353.603000041353.641353.60Moreover, to further validate the proposed method, the elastoplastic finite element analysis of a practical problem is carried out based on the H-B criterion. The excavation of a circular tunnel under a plane strain assumption is taken as an example. The dimensions and the boundary conditions of the tunnel and the adopted finite element meshes are shown in Figures3 and 4, respectively. To obtain the initial stress field, the outer boundary is subjected to a hydrostatic pressure of q = 100 MPa. The material is assumed to be weightless. The elements used by the test are four-node plane bilinear quadrilateral elements, and a total of 2,025 elements are used, as shown in Figure 4. The displacements perpendicular to the left and lower boundaries are restricted. In the H-B criterion, σci is set to 200 MPa, and the other parameters are identical to those mentioned above. First, the initial stress is generated; then, the excavation is carried out.Figure 3
Geometry and boundary conditions for the circular tunnel excavation example.Figure 4
Example of the element mesh with 2,025 elements.Figure5 shows the relative errors of the principal stresses obtained by Clausen’s method and the proposed method. As shown in this figure, the results obtained by the proposed method are consistent with those obtained by Clausen’s method. Moreover, the relative errors indicate that there is little difference between the cases of associated plasticity and nonassociated plasticity. In view of these results, the proposed method is more suitable for the associated plasticity than the nonassociated plasticity.Figure 5
Relative errors of the final principal stresses: (a) associated case and (b) nonassociated case.
(a)(b)


## 6. Discussion

In this paper, a return mapping algorithm suitable for a series of nonlinear yield criteria ignoring the effect of the intermediate principal stress is proposed. When the proposed algorithm coded in the finite element software is utilized to a yield criterion ignoring the effect of the intermediate principal stress, it merely requires the equivalent strength parameters of these criteria, the derivatives of their functions, and the tensile strength because the equivalent strength parameters are combined during algorithm derivation. The application example also demonstrates its capability of exactly executing stress updates. Compared with other return mapping algorithms in the principal stress space, a yield criterion can be more easily implemented in a numerical simulation with the finite element method when the algorithm is utilized to a yield criterion due to the merit of the proposed algorithm.Usually, the equivalent M-C strength parameters of a nonlinear yield criterion can be obtained by the Balmer method [34]. In fact, the Balmer method is suitable for nonlinear criteria whose expression can be written as σ1=fσ3. In this case, an explicit expression of the equivalent parameters can be given. For other yield criteria that are not written as σ1=fσ3; however, there is not an explicit expression of the equivalent parameters, causing the application of the proposed algorithm to these criteria to be very inconvenient. Moreover, it should be noted that some yield criteria have such a complicated expression of the equivalent parameters that the derivatives of the functions of these equivalent parameters will become more complicated. To solve this problem, some commercial mathematical software programs (e.g., Mathematica and MATLAB) can be used to calculate the equivalent parameters and the derivatives of their functions.The equivalent M-C strength parameter is essentially a function of the principal stress, which gives us a hint. Under the framework of the proposed algorithm, we need to obtain only the functions of the strength parameters concerning the principal stress through the test, and we do not need to know which soil or rock mass material obeys which yield criterion, which will sufficiently improve the flexibility, and which is convenient to numerically analyse a special engineering problem.Finally, it is important to note that because this algorithm is based on the implicit return mapping algorithm in the principal stress space for the M-C yield criterion; the proposed algorithm is suitable for nonlinear yield criteria that ignore the effect of the intermediate principal stress.


## 7. Conclusions

(1)
Combining the equivalent M-C strength parameters, three stress return schemes (return to the yield surface, return to the curve, and return to the apex point) are derived. The conditions for choosing the proper stress return scheme are constructed, and a modified return mapping algorithm is thus established for the nonlinear yield criteria that ignore the intermediate principal stress.(2)
The stress return results based on the H-B criterion show that the updated stresses obtained by the proposed algorithm are in good agreement with those obtained by Clausen’s method.(3)
The results of the numerical simulation of the excavation of a circular tunnel illustrate that the final principal stresses obtained by the proposed method are also in good agreement with those obtained by Clausen’s method. The comparison between the associated and the nonassociated plasticity cases shows that the proposed algorithm is more suitable for the associated plastic-flow rule.



---

*Source: 8823988-2020-11-24.xml*