Strain-Softening Analyses of a Circular Bore under the Inﬂuence of Axial Stress

: Strain-softening analyses were performed around a circular bore in a Mohr–Coulomb rock mass subjected to a hydrostatic stress ﬁeld in cross section and out-of-plane stress along the axis of the bore. Numerical procedures that simplify the strain-softening process in a step manner were employed, and on the basis of the theoretical solutions of the elastic–brittle–plastic (EBP) medium, the strain-softening results of the displacements, stresses and the plastic zones around the circular bore were obtained. The numerical solution was validated based on the fact that the strain-softening process became EBP when the softening slope was very steep and elastic-perfectly plastic (EP) when the softening slope was near zero. The results illustrated that the stresses and displacements in the rock mass surrounding the bore were affected by axial stress and that a proper consideration of out-of-plane stress was necessary. Moreover, the presented results can be used for the veriﬁcation of numerical codes.


Introduction
Many analytical solutions are presented while the excavation of circular and spherical cavities is considered in pre-stressed fields. Different yield criteria (including Mohr-Coulomb (M-C) and Hoek-Brown (H-B) criteria) are employed, and the corresponding EP, EBP and strain-softening constitutive relations are taken into account.
For elasto plastic and elastic-brittle-plastic constitutive relations, Brown et al. [1], Detournay [2], Ogawa and Lo [3], Wang [4], Carranza-Torres [5], Carranza-Torres and Fairhurst [6], Yu [7], Park and Kim [8] and Sharan [9] presented analytical and semianalytical solutions for predicting the stresses and displacements around the cavities. Park and Kim [8] commented on this kind of problem and derived closed-form solutions on the basis of M-C and H-B yield criteria.The displacements obtained by different assumptions were compared. Recently, Frikhaand Bouassida [10] presented an analytical model dealing with an expanded cylindrical cavity in an elasto plastic medium with a variable plastic potential of flow. Tran et al. [11] analysed the earth pressure around the cylindrical shafts in soft ground by means of experimental and numerical method and discussed the shaft-soil interaction.
With consideration of the strain-softening behavior of geomaterial, Alonso et al. [12], Guan et al. [13], Park et al. [14], Lee and Pietruszczak [15] and Wang et al. [16] presented different forms of solutions for the circular bore in M-C and H-B rock and ground reaction or response curve (GRC). Their results were in good agreement with each other. Zhou and Randolph [17] studied the shear bands ahead of the advancing cylindrical and spherical penetrometers. Chen et al. [18] predicted the evolution of the plastic region around a wellbore in Drucker-Prager rocks.
Normally, the problem is treated as an axisymmetric plane strain to derive analytical or semi-analytical solutions. The geostress is presumed to be hydrostatic and the geomaterial isotropic. The out-of-plane stress is assumed to be the second principle stress with σ z = ν(σ x + σ y ) .
Reed [19] studied the abovementioned axisymmetric plane strain problem in which the initial stress was of equal magnitude in three axes, and he found that the axial stress did not remain the intermediate one The effect of the out-of-plane stress on the cavity wall displacement was significant for geomaterial with a large drop in strength at yield.
Lu et al. [20] analyzed the circular bore with different magnitudes of the out-of-plane stress. The numerical analyses show that the results are much influenced by the out-ofplane stress.
Wang et al. [21] presented the theoretical solution for predicting stresses and displacements around a circular bore in M-C rock mass, subjected to a hydrostatic stress field in the cross section and out-of-plane stress along the axis of the bore. The EBP relationship with a non-associated flow rule was used in the analysis. The comparison of the results illustrated that the effect of the out-of-plane stress on the displacements and plastic zones cannot be ignored.
As an extension of the brittle-plastic solution of the circular bore [21], the strainsoftening analysis is performed and the corresponding solutions were obtained. Numerical schemes were presented and verified, and the results obtained by different constitutive relationships were compared.

Problem Definition
A long circular bore with radius r 0 excavated in a homogeneous, infinite isotropic rock mass can be considered to be an axisymmetric problem if the far-field stress is uniform as shown in Figure 1. During bore excavation, the internal pressure (p in ) is gradually reduced, and displacement and stresses are only the functions of radius r when gravity is ignored. In the initial stage, the medium is elastic and the solution is Since the problem is axisymmetric, σ z , σ θ and σ r are the three principle stresses. It should be noted that compression stress is taken to be positive in this paper. When the internal pressure continues to decrease, a plastic zone appears around the bore.
If the M-C criterion is employed, it is written as

Critical Pressure
In consideration of axial stress, there will be 3 cases. In the first case, the axial stress is the major principle stress at the bore wall, i.e., σ z > σ θ > σ r . The critical pressure is This means that the plastic zone begins to develop in the surrounding rock if the internal pressure p in ≤ p c1 .
As the internal pressure is reduced, σ z decreases and σ θ increases. At a certain stage, σ θ = σ z at the bore wall. Thereafter, the plastic region is as shown in Figure 2. In the inner part of the plastic region, σ θ = σ z > σ r , whereas in the outer part σ θ > σ z > σ r . As the inner pressure decreases, strain-softening continues until residual strength is reached. An interface (indicated by radius R s ) between strain-softening and the residual strength zone will appear. Evolution of plastic type I (a) plastic zone with σ z > σ θ > σ r ; (b) inner plastic zone with σ z = σ θ > σ r and outer plastic zone with σ z > σ θ > σ r .
In this case, when We know that the out-of-plane stress is the intermediate principle stress when q 2 > q > q 1 .
When the internal pressure p in ≤ p c2 and continues to decrease, the plastic region is as shown in Figure 3. (a) plastic zone with σ θ > σ z > σ r ; (b) inner plastic zone with σ θ = σ z > σ r and outer plastic zone with σ θ > σ z > σ r .
In the third case, the critical internal pressure is When the internal pressure p in ≤ p c3 and continues to decrease, the plastic region is as shown in Figure 4. . Evolution of plastic type III. (a) plastic zone with σ θ > σ r > σ z ; (b) inner plastic zone with σ θ > σ r = σ z and outer plastic zone with σ θ > σ r > σ z . Figure 5 shows a typical stress-strain curve for strain-softening geomaterial. It consists of a linear part 'OP', strain-softening part 'PB' and residual part 'BC'. During the strain-softening stage, the strength parameters decrease as the softening parameter increases. Here, the shear plastic strain η = ε p 1 − ε p 3 is used as the softening parameter. The strength parameters are equal to their residual values as η ≥ η * . (η * is the limit of shear plastic strain). To avoid an egative slope, the strain-softening part is divided into manysegments such as 'PA 2 ', 'A 2 A 4 ', ' A 4 A 6 ' and 'A 6 B' shown in Figure 5. Each segment is simplified [16]. Equation (7) describes the parameter evolution:

Solution Procedures
As p in is lower than critical pressure (p c ), the plastic region will develop in the surrounding rock mass. Assuming that the inner pressure reduces gradually from p c to p in , the pressure decrease will be divided into n steps and denoted by ∆p k (k = 1, n). After the n steps of the pressure decrease, the internal pressure p in(n) = p c − n ∑ k=1 ∆p k .
If the inner pressure changes from p c to p in(1) = p c − ∆p 1 , a plastic annulus with thickness r 1 −r 0 appears as shown in Figure 6a.  (1) to p in (2) ; (c) from p in(j−1) to p in(j) .
If the pressure decreases further from p in (1) to p in(2) = p c − ∆p 1 − ∆p 2 , plastic deformation will continue in the first plastic annulus (with thickness r 1 −r 0 ). Meanwhile, a new plastic annulus with thickness r 2 −r 1 will appears as shown in Figure 6b.
By analogy, when the pressure decreases to p in(j) = p c − j ∑ k=1 ∆p k , there will be j annuli in the surrounding rock mass. The jth annulus, adjacent to the elastic zone and with inner radius r j−1 and outer radius r j , is illustrated in Figure 6c. Plastic deformation will go on in the other j−1 plastic annuli until the residual strength is reached during the pressure decrease.
In each of the plastic annuli, incremental stress, displacement, et al. will be obtained, and the solution process will be discussed in detail later. Now, the ith load increment (or the ith pressure decrease ∆p i ) is to be considered, I which ∆σ r(i) and ∆σ θ(i) are the corresponding incremental radial and tangential stresses respectively, and σ r(i−1) and σ θ(i−1) are the stress components at the end of the i−1th load increment.
The equilibrium equation for the incremental stresses is, but Equation (8) can be rewritten as

σ z : 1st Principle Stress
When σ z > σ θ > σ r , the yield criterion is formulated as Equation (10) can be rewritten as and the condition of compatibility is taken into account, Assuming that the plastic potential uses the same formula as the yield criterion, the incremental radial and axial plastic strain are written as The incremental plastic part of tangential strain is where β i = 1+sin ψ i 1−sin ψ i and ψ i is the dilation angle corresponding to the ith load increment. The strains are usually divided into elastic and plastic parts: With ∆ε p θ(i) = 0 and ∆ε z(i) = 0 considered, we have The incremental elastic strains should obey Hooke's law, and they are in the following form By substituting Equation (22) into Equations (20) and (21), and then into Equation (23), the following equation is obtained: Equations (9), (11) and (23) will be used to get ∆σ z(i) , ∆σ θ(i) and ∆σ r(i) . Substitution of Equations (9) and (10) into Equation (23) produces a second-order differential equation for incremental radial stress.
Solving thedifferential equation and the general solution is In Equation (25), F i and G i are determined by boundary conditions. To get F i and G i , Equation (26) is used: From Equations (9) and (11), Equation (26) can be written as By substituting Equation (25) into Equation (27), Equation (28) is obtained: Equations (25) and (28) will be used to obtain incremental stress and displacement in the plastic region.
For the ith load increment, the total number of the plastic annulus is i. In each plastic annulus, F i , G i and the outer radius (r i = R) of the ith plastic annulus need to be solved. Therefore, total number of unknown variables is 2i + 1.
The radial stress and displacement are continuous between the adjacent plastic annuli, so there are 2(I − 1) equations. In addition, at In total, there are 2i + 1 equations.
So far, we have the same number of equations as that of unknown variables. The abovementioned equations form a set, and the unknown variables are solved by iterative methods.
As the pressure p in decreases further, σ z decreases and σ θ increases. Once σ θ = σ z , the plastic area is divided into two.
For the outer plastic zone with σ z > σ θ > σ r , the solution for radial stress and displacement increments was discussed above and can be obtained from Equations (25) and (28).
In the inner plastic zone, σ z = σ θ > σ r , and the yield criteria governing the inner plastic zone are For the ith load increment, Equation (30) is written as The incremental stress components (∆σ θ(i) and ∆σ r(i) ) can be obtained by the equilibrium Equation (8), yield criterion (31), and ∆σ z(i) = ∆σ θ(i) . The solution processes for the incremental stress components are the same as those in Wang et al. [16].
Because σ z = σ θ in the inner plastic zone, the stresses lie on the intersection of the two yield surfaces (Equations (29) and (30)). The correct flow rule is obtained by summing the contributions from the two plastic potentials [19].
The increments of plastic strain are formulated as By using Equations (12) and (17)- (19), the differential equation for incremental displacement in the inner plastic zone is given by By means of Equations (22) and (31), we have where At present, the right hand side is known in Equation (35). After ∆u (i) is obtained, the incremental displacement at the interface of the inner and outer plastic zones becomes known and can be used as boundary condition for solving Equation (35). Thus, ∆u will be easily determined by numerical integration.
So far, solution processes have been presented for the two plastic zones with σ z > σ θ > σ r and σ z = σ θ > σ r .

σ z : 2nd Principle Stress
In the plastic zone with σ θ > σ z > σ r , the plastic region governed by the M-C criterion is The solution for incremental stress (∆σ θ(i) and ∆σ r(i) ) and displacement (∆u (i) ) is exactly the same as those in Wang et al. [16]. The incremental stress along the axis is ∆σ z(i) = ν(∆σ θ(i) + ∆σ r(i) ).
For the area with σ θ = σ z > σ r , incremental stresses and displacement were obtained in the last sub-section.
where (45) is in the same form as in Equation (25) with H = − g 3 g 2 ,

The solution for Equation
. F i and G i in Equation (25) will be determined by boundary conditions. By means of (22), (43) and (44), complex manipulations produce with 1) )/T. In a manner similar to the solution process in Section 3.3, F i and G i can be determined for all plastic annuli, and thenincremental radial stress and displacement can be obtained by Equations (25) and (46).
In the plastic region, both σ r and σ θ decrease when the pressure p in decreases further. Once σ r = σ z , there will be 2 plastic zones.
For the outer plastic zone with σ θ > σ r > σ z , the solution for radial stress and displacement increments is discussed above.
For the inner plastic zone with σ θ > σ r = σ z , the solution for radial and tangential stresses is the same as those in Wang et al. [16]. The stress in the axis of the bore is σ z = σ r .
The solution for displacement in the inner plastic zone is presented as follows: In the inner plastic zone, the yield criteria are The flow rule is obtained by summing the contributions from the two plastic potentials: In the ith load increment, By using Equations (12) and (17)- (19), the differential equation is in the same form as in Equation (35) By means of Equations (22) and (39), we have where After ∆u (i) is obtained in the outer plastic zone, the incremental displacement at the interface of the inner and outer plastic zones becomes known and used as a boundary condition for solving Equation (35). Here, it should be noted that f (r) in Equation (35) is formulated with Equation (53) in this case. Then, ∆u will be easily determined by numerical integration.
So far, solution processes have been discussed for the two plastic zones with σ θ > σ r > σ z and σ θ > σ r = σ z .

Convergence of Numerical Solutions
While the numerical analysis is carried out, the internal pressure decrease from p in = p c to p in = 0 is divided n times and each pressure decrease is ∆p i = (p c − p in )/n. First, the influence of load step n on the numerical results is investigated.
A summary of the parameters is given in Table 1. Several values of n are used to demonstrate the convergence of the numerical solutions as shown in Figure 7.

Numerical Analysis
When the parameters presented in Table 1 are applied, we learn that the out-of-plane stress is the intermediate principle stress if q = 30 MPa and that it is the major principle stress if q = 60 MPa.

Axial Stress: q = 30 MPa
For, q = 30 MPa, p c = 5.77 MPa, when the internal pressure is less than p c , the plastic zone will develop in which σ θ > σ z > σ r . While the internal pressure continues to decrease, two kinds of plastic zones form at certain stages. After the support pressure is completely released, the distribution of the stresses and displacement in the plastic zones is displayed in Figure 8. Both plastic zones are observed. In the outer plastic zone, σ θ > σ z > σ r and the normalized plastic radius R/r 0 = 1.754. In the inner plastic zone, σ θ = σ z > σ r . The normalized interface radius (R/r 0 ) is 1.6. The dimensionless radial displacement at the bore wall is 4.35. To understand the differences among various constitutive models, the distribution of the stresses and displacement in the plastic zones, which are obtained by means of the EBP and EP models, are shown in Figures 9 and 10.  A contrastive analysis of the results presented in Figures 8-10 shows that the dimensionless radial displacement at the bore wall, the normalized plastic radius and the normalized interface radius obtained by the strain-softening model fall between the EBP and EP models.
When the geomaterial behaves in an EBP manner, the plastic region in the surrounding rock mass and the inward displacement at the bore wall are the largest. For EP media, the plastic region and inward displacement are the smallest.
Two types of plastic zones (σ θ = σ z > σ r and σ θ > σ z > σ r ) around the bore appear for the 3 models. Figure 11 shows the evolution of the radii of the elastic-plastic and softening-residual interfaces, denoted by R and R s , for different strain-softening slope. Figure 12 shows the corresponding GRCs. Figure 11. Evolution of (a) plastic radii R and (b) softening radii R s (the radii between softening and residual zones) in the strain-softening media with different softening slopes when q = 30 MPa. It should be noted that η* is a parameter related to plastic strain and that it controls the strain-softening slope, i.e., a different η* means a different strain-softening slope. The larger the η*, the smaller the slope of the softening. When η* is large enough, the slope of the softening is zero, and strain-softening model becomes the EP one. While it becomes very small, the slope of the softening is infinite (i.e., the peak strength changes into residual strength abruptly), and strain-softening model becomes EBP one.
In Figures 11 and 12, the results are very close to that of the EBP model when η* = 0.002. It can be inferred that the results are close to those of the EP model when η* > 0.016. Figure 11b shows that part of the plastic region is in the residual strength stage and the other part is in the strain-softening stage when η* = 0.002, 0.004 and 0.008. The strain-softening zone area becomes small as η* decreases.

Out-of-Plane Stress: q = 60 MPa
Here the geometric and mechanical parameters presented in Table 1 is used except that q = 60 MPa.
While the internal pressure is less than p c (=6.5 MPa), the plastic zone first develops when σ z > σ θ > σ r . Only one plastic zone with σ θ = σ z > σ r is formed in this case as the support pressure is completely released. The distribution of the stresses and displacement in the plastic zones is shown in Figure 13. The normalized plastic radius R/r 0 = 1.85 and the dimensionless radial displacement at the bore wall is 5.59.   A comparison of the results presented in Figures 13-15 indicates that the dimensionless radial displacement at the bore wall and the normalized plastic radius obtained by the strain-softening model fall between those of the EBP and EP models.
The normalized plastic radii are 1.88, 1.85 and 1.31, and the dimensionless radial displacements are 5.6, 5.59 and 1.82, corresponding to the EBP model, strain-softening model and EP model, respectively. For these 3 models, a similar phenomenon is observed concerning the magnitude of the normalized plastic radius and the dimensionless radial displacement at the bore wall.
Two types of plastic zone (σ θ = σ z > σ r and σ θ > σ z > σ r ) around the bore appear only for the EP model. Figure 16 shows the evolution of the radii of the elastic-plastic and softening-residual interfaces for different strain-softening slopes. Figure 17 shows the corresponding GRCs.
In Figures 16 and 17, the results are very close to those of the EBP model when η* = 0.004. However, we have the same curves as those obtained by EBP model if η* = 0.002. This is a little different from the results shown in Figures 11 and 12 when q = 30 MPa and η* = 0.002. It looks as though the influence of the out-of-plane stress is great. Figure 16b shows that part of the plastic region is in the residual strength stage and the other part is in the strain-softening stage when η* = 0.004 and 0.008. Figure 16. Evolution of (a) plastic radii R and (b) softening radii R s (the radii between softening and residual zones) in the strain-softening media with different softening slopes when q = 60 MPa. Numerical results obtained from strain-softening analysis are summarized in Table 2 for q = 30 and 60 MPa. Obviously, both plastic radii and the inward displacement at the bore wall increase when the out-of-plane stress ascends from q = 30 to 60 MPa.

Influences of the Dilation Angle on Numerical Results
In plasticity theory, either an associative flow rule or a non-associative flow rule is used to derive the incremental plastic strain. When the associative flow rule is used, the plastic potential function is exactly the same as the yield function and the dilation angle equals the internal friction angle in the M-C criterion, whereas, in the non-associative flow rule, the plastic potential function is different from the yield function.
It is an observable fact that dilatancy is dependent on shear plastic strain [2,22]. Moreover, the dilation angle is not equal to the internal friction angle. According to Hoek and Brown [22], the following relationships are recommended, and the dilation angle and the friction angle are in the range ψ = ϕ/4 for good-quality rock and ψ = 0 for poor-quality rock. This means that the dilation angle is much smaller than the internal friction angle and that the non-associative flow rule is more appropriate for elasto plastic analysis. Figure 18 shows GRC in the strain-softening media with η* = 0.004 and q = 30 MPa. Different evolution laws of dilation angle are considered.
First, the associative flow rule was employed, and ψp = ϕp = 50.9 • and ψr = ϕr = 39 • . The dimensionless radial displacement increased as the internal pressure decreased. The dimensionless radial displacement was 31.8 when the internal pressure was completely released. When ψp = 1/2, ϕp = 25.45 • , ψr = 1/2, and ϕr = 19.5 • , the flow rule was non-associative. The dimensionless radial displacement was 8.3 when the internal pressure was zero. When ψp = 15.5 • and ψr = 7.5 • , the dimensionless radial displacement is 5.4. GRC was also presented in Figure 18 for comparison. In this case, the flow rule was also non-associative. Figure 19 shows GRC in the strain-softening media with η* = 0.004 and q = 60 MPa, and similar analyses were carried out. When the associative flow rule was used, the dimensionless radial displacement was 59.5. When the non-associative flow rule was employed, the dimensionless radial displacement was 11.7 for ψp = 1/2 ϕp and ψr = 1/2 ϕr, and 7.2 for ψp = 15.5 • and ψr = 7.5 • . The dimensionless radial displacements are summarized in Table 3 for different outof-plane stresses. The results indicated that the displacements obtained by the associative flow rule were larger than those by the non-associative flow rule. When the dilation angle was much smaller than the internal friction angle, the difference was huge. When the out-of-plane stress increases, the difference between the displacements obtained by associative flow rule and those by non-associative flow rule also increases.

Discussions of the Numerical Analysis
In last subsection, a set of parameters was used to study the influence of the out-ofplane stress on bore stability. The plastic radius and the distribution of the stresses and displacements were presented and compared for q = 30 and 60 MPa. Figures 8 and 13 indicate that the axial in situ stress affects displacement, plastic radius and type of plastic zone in the rock mass surrounding the bore in strain-softening analysis. For the same 'η*', a comparison of Figures 11,12,16 and 17 show that the displacement and plastic radius are larger when q = 60 than when q = 30 MPa.
Different evolutionary laws of the dilation angle were considered in the strainsoftening analysis, and the calculation results were presented in Table 3. In the first case, the dilation angle (ψ) and the internal friction angle (ϕ) were equal and had the same evolution law in the post-failure stage, and the displacements at the bore wall were the largest. When the dilation angle was half of the internal friction angle and both evolved in a similar manner, the displacements at the bore wall were much smaller than the displacements in case 1.Case 3 showed displacements obtained by means of the parameters in Table 1 and were the smallest.
When q = 30 MPa, the displacement in case 1 was about 4 times that of case 2, and about 6 times that of case 3. When the out-of-plane stress increased to q = 60 MPa, the displacement in case 1 was about 5 times that of case 2, and about 8 times that of case 3. The differences were obvious for various evolution laws of the dilation angle. Moreover, the ascending out-of-plane stress made the difference increase.
Normally, most engineering rock masses display neither in an elastic-perfectly plastic nor in an EBP way. They behave in a strain-softening manner as shown in Figure 5. During the strain-softening process, a progressive loss of strength occurs and the parameters (cohesion c, internal friction angle ϕ and dilation angle ψ) will vary from their peak values to residual ones. The EBP model overestimates the deformation and the plastic zone of the geomaterials. However, the EP model underestimates the deformation and the plastic region of geomaterials. If the strain-softening analysis is carried out, the numerical results are more appropriate for practical use. When deep bores are excavated in hard rock masses, high geostress occurs. In this situation, the influence of axial stress cannot be ignored when the problem is treated as a normal plane strain.

Concluding Remarks
This paper discussed an approach for obtaining the solutions of the axi-symmetric bore problem considering out-of-plane stress. Based on EBP analysis [21], incremental formulations for strain-softening solution of the problem were presented.
The emphasis was put on the strain-softening characteristics of the surrounding M-C rock masses, which were treated as an infinite isotropic media subject to a uniform stress field (σ0) in the cross section of the bore and out-of-plane stress (q) in the axis. The influences of the out-of-plane stress on displacement, stresses and plastic radius were analyzed.
Strain-softening analyses were carried out for both q = 30 and 60 MPa. Numerical results showed that the plastic radius and displacement at the bore wall obtained by q = 60 MPa were bigger than those for q = 30 MPa. The distribution of the displacement, stresses and plastic zone were affected by the magnitude of the out-of-plane stress.
Different softening slopes (i.e., different critical shear plastic strains) were studied, and it was shown that the strain-softening solution converged to the brittle-plastic solution when the critical shear plastic strain was small enough and that converged to the EP solution when the critical shear plastic strain was large enough. Meanwhile, the out-ofplane stress influenced the numerical results.
Different evolutionary laws of dilatancy in the post-failure stage were considered, and the corresponding results showed that the influence of out-of-plane stress on the bore deformation increased when dilation angle increased. Normally, the bore deformation is overestimated when the associative flow rule is used.
In reality, most rocks encountered in geotechnical engineering are neither ideally brittle nor ideally plastic; they behave in strain-softening manner. Results obtained by strainsoftening analysis might be more appropriate for evaluating the stability of underground structures.