The Elastoplastic Solutions of Deep Buried Roadway Based on the Generalized 3D Hoek-Brown Strength Criterion considering Strain-Softening Properties

State Key Laboratory of Coal Resources and Safe Mining, School of Mines, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China Key Laboratory of Deep Coal Resource Mining, Ministry of Education of China, School of Mines, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China Institute of Mining Engineering and Geology, Xinjiang Institute of Engineering, Urumqi 830091, China Rizhao Natural Resources and Planning Bureau, Shandong, Rizhao 276800, China School of Science, Yangzhou Polytechnic Institute, Yangzhou 225127, China


Introduction
The deformation of rock mass is related to geological conditions, in situ stress, support technologies, and excavation disturbance. The prediction of rock mass deformation behavior is very useful for deep coal mining [1]. At great depth with strong deep structural activity or high deformability, rock squeezing may occur, and the repair rate of large deformed roadways can reach 70% [2,3]. In the deep coal mining, due to the influences of "three high and one disturbance" (high geo-stress, high ground temperature, high permeable pressure, and the disturbance of mining activities), the unique stress characteristic and plastic yielding were engendered by the mechanics of roadway [4][5][6].
For deep underground engineering, the elasto-plastic [7][8][9][10], elasto-brittle-plastic [11][12][13], and strain-softening models [14][15][16] were established to finish off the issue of tunnel based on elasto-plastic approaches [17][18][19] in the past period. For the first two methods mentioned above, analytical solutions are available accounting to a simple point of view. However, experimental and field observation show that the deformation and failure mechanism of rock mass in the postpeak stage of rock was revealed by the strainsoftening model [19][20][21][22][23][24]. Hence, lots of scholars have studied the strain-softening model in detail, especially for Brown et al. [17], Carranza-Torres [18], Alonso [7], and Lee and Pietruszczak [25]. There are two pivotal parameters, which are the critical softening parameter η * and the dilatancy angle ψ, to research the features of strain-softening of deep underground engineering. It is worth noting that the confining stress acts on the above two parameters according to the existing findings [7,8,20,26,27]. In addition, the dilatancy angle was only deemed to a constant to research the features of strain-softening [7,8,25,26]. Nevertheless, the variation of the critical softening parameter η * and the dilatancy angle ψ should be seemed reasonable [20,27,28].
The above scholars mainly study the related problems of circular tunnel, but for deep underground engineering, there are similar problems in roadway in coal mine. Above all, five kinds of cross-section shape of roadway are commonly used in mining design, including rectangular, straight wall arch, trapezoidal, round, and oval [29]. The straight-wall semicircular arch roadway is widely used as the cross-section shape of underground roadway in coal mine because of its fast forming and high stability [30]. The impact of straight wall arch roadway on the surrounding rock was equal to the outer round with the same diameter of the circular roadway section [31]. In conclusion, for the convenience of analysis and calculation, the typical arched cross-section shape of the coal mine roadway with straight wall is equivalent to the outer circular cross-section roadway with the same diameter for the following specific analysis (as shown in Figure 1).
What is more, the elasto-plastic, elasto-brittle-plastic, and strain-softening models have been built based on Mohr-Coulomb (M-C) and Hoek-Brown (H-B) criteria [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][32][33][34][35][36]. The Mohr-Coulomb criterion is linear which is relatively convenient to use, but the criterion does not reflect the nonlinear destruction characteristics of rock bodies. Since the Hoek-Brown strength criterion has been widely used in rock engineering, its prediction of rock behavior is sufficient and can be easily applied to a series of rock engineering problems. However, although many evidences show that in many cases, the intermediate principal stress does affect the strength of the rock, the H-B strength standard does not consider the influence of the intermediate principal stress. In order to overcome this shortcoming, the criterion proposed by Zhu [37] and Zhang [38] can predict the same strength as the original H-B strength criterion under triaxial compression and extension and is regarded as a true 3D version of the original H-B strength criterion. 3D Hoek-Brown strength criterion (GZZ strength criterion) was the strength criterion suggested by International Society for Rock Mechanics (ISRM). Therefore, the GZZ strength criterion can be used as a failure criterion for constructing rock quality constitutive models and can be implemented in 3D finite element (FE) codes to perform 3D numerical analysis of rock engineering problems. The main objective of this study is to compare results considering the effect of a variable η * and ψ on the strain-softening behavior of rock masses in equivalent circular roadways to select the correct combinative models under the GZZ strength criterion.

Description of Problem
For the convenience of analysis and calculation, the typical arched cross-section shape of the coal mine roadway with straight wall is equivalent to the outer circular cross-section roadway with the same diameter for the following specific analysis.
Based on the above analysis of roadway, for the sake of further study on the properties of surrounding rock stress, the following assumptions are required: (1) The roadway is equivalent circular (2) The rock mass is characterised as isotropic, continuous, infinite, and initially elastic, and its mechanical properties obey the GZZ strength criterion As shown in Figure 2, the excavation radius of the equivalent circular roadway is R 0 . Here, the hydrostatic stress field σ 0 is applied to the whole area before excavation. At the areas where the magnitude of the internal support strength p i is no larger than the critical value p c , a plastic zone is assumed to be formed around the equivalent circular roadway. When p i < p c * * , p c * * is deemed to the pivotal support pressures; the equivalent circular roadway was got into the plastic softening zone. When p i < p c * , p c * is deemed to the pivotal support pressures, the equivalent circular roadway was got into the plastic residual zone. In the case of elastic-plastic behavior, 2 Geofluids the explicit expression of plastic radius R p can be derived [9,11]. Additionally, when consider the strain-softening behavior, the plastic zone can be separated by softening and residual zones through the application of an interface with radius R s (as displayed in Figure 2). At the boundary of the elasto-plastic zones, σ θp and σ rp are the tangential stress and the radial stress, respectively. Additionally, at the boundary of the softening residual zones, σ θs and σ rs are the tangential and the radial stress. According to the stress equilibrium at the boundary of the elasto-plastic and the softeningresidual zones, p c * * and p c * are equal to σ rp and σ rs , while the R p and R s are the radii of plastic and plastic residual zone.
2.1. GZZ Strength Criterion. The standard Hoek Brown strength criterion can be used to calculate the threedimensional stress state [33]. But the effects of intermediate principal stress is not considered in the initial form, which does influence the rock strength [39,40]. Therefore, it is believed that the 3D Hoek-Brown strength criterion is more suitable for describing yield state of rock mass under 3D stress states. As is shown in Figure 3, the 3D Hoek-Brown strength criterion is established by Zhang and Zhu [37]. It is also the strength criterion suggested by International Society for Rock Mechanics (ISRM) and can be expressed as follows.
In Equation (1), m b , a, and s are all empirical parameters reflecting rock mass characteristics. a is rock mass characteristic parameter. m b is the empirical parameter of the rock, which can be determined by experiment or rock types. s reflects the degree of rock fragmentation. σ c is the uniaxial compression strength. τ oct is the octahedral shear stress, and σ m is the average value of main stress. τ oct and σ m can be calculated by the Equations (2) and (3), respectively.
Where m b = m i and s = 1 for intact rock, these quantities can be calculated by the geological strength index GSI and a factor D, where the value of D varies from 0 to 1. GSI is a factor that is related to the degree of disturbance caused by blast damage and stress relaxation, which varies between 10 and 100. m i is determined by mechanical properties of rock mass. The values of m i varies from 2 to 32, which is related to the compositions and mineralogical properties of rock mass.

Strain-Softening
Behaviour. The strain softening is usually a transitional yield process or plastic potential, and its behaviour is dominated by a softening parameter η which can be described as follows: where ε 1 p and ε 3 p represent the major and minor principal  3 Geofluids plastic strain, which are equal to the tangential ε θ p and radial strains ε r p . η is the plastic shear strain. As displayed in Figure 4, the mechanical properties including the strength and deformation parameters is different when the deviatoric plastic strain change. This relationship can be described by a bilinear function with the application softening parameter η as follows: where ω is one of the strength parameters and η * represent the critical plastic softening parameter as shown in Figure 4. When η reaches critical plastic deviator strain, it begins to enter the residual state. Provided that the experiment data are available, the η * may be assigned for each parameter. Here, for the purpose of the clarity, the single value η * is assigned for each medium. The superscripts "p" and "r" represent the peak and residual strength. As displayed in Equation (6), at the stage of the plastic softening, the strength parameters has the negative linearly relationship with the η * . Once the parameter increases to the value higher than the critical value η * , the strength parameters remain stable. The peak and residual values of m b , s, and a can be obtained through the application of the geological strength index (GSI) [23]. The residual value of GSI is able to be obtained from the peak value of GSI [21,22].
The nonassociated flow regulation is described as follows: where K ψ is the dilatancy coefficient, K ψ = ð1 + sin ψÞ/ð1 − sin ψÞ, and ψ is the angle of dilatancy, which are different with the change of the confining stress and mechanical properties of the rock [20].
2.3. Variable Path of σ r . As described in the previous studies [7], the process of roadway excavation is simulated as the process of roof confining pressure decreasing gradually from σ 0 to p i . Figure 5 shows the excavation stress path and yield criterion from the original hydrostatic stress to the state of the final support. Point A is defined as the state of the initial hydrostatic. At the stage of the excavation process, σ r decreases gradually to the state of the elasto-plastic boundary (point B) where the confining stress σ rp (referring to p * * ) can be calculated by the elastic solution. Then, the status of the rock mass reaches the plastic softening area. During the process between point B and point C, since the strength parameters varies at the area of plastic softening, the GZZ function change with η (0 < η < η * ). Once the value of η reaches η * , σ r decreases to the values of σ rs which refers to p * . At this point, the state of the rock mass is transferred from the softening to the residual. During the process between point C and point D, the residual failure criterion can be employed with the application of the parameter of residual strength ω r . Then, if p i is lower than the σ rs , σ r decreases from σ rs at the boundary of the softening-residual to p i on the surface of the excavation, which is imposed by the supporting systems. The results indicate with the application of the support systems, the σ r deceases from σ rp to p i at the plastic zones.
For the soft rock masses after the excavation of the equivalent circular roadway, at the plastic zone, the values of the σ r decrease from the boundary of the elasto-plastic to the of boundary excavation along the path of stress that displayed in Figure 5. The stress evolution in the equivalent circular roadway is different from that in the compression test. In the compression test, the confining pressure stress remains unchanged, while the axial stress increases until the occurrence of the failure. Therefore, if not conder the effects of the variable confining stress, the accuracy and practicability of the assumption that the critical plastic parameter η * and dilatancy angle ψ remain unchanged should be taken with caution. Due to the variable property of σ r in the plastic zone of strain softening behavior, the parameters related to the confining stress such as the η * and ψ are also variable.

The Dilatancy Angle
Model. The experimental results [26] show the angle of dilatancy angle is directly related to the failure of rock mass. Therefore, the dilatancy model is of great significance in the calculation of mine excavation displacement. In essence, in rock excavation engineering, in order to better reflect the mechanism evolution in the process of rock failure and study the potential dilatancy process in the process of rock plastic deformation, it is necessary to establish an appropriate dilatancy model.
Previous studies have described the dilatancy angle of rock mass is not constant, but decreases with the increase of confining pressure [42][43][44][45]. In addition, the dilatancy decreases gradually until it become constant at the residual strength during the failure process. If the plastic deformation is large enough, the dilatancy will gradually decrease to zero. The maximum values of the dilatancy angle appear near the peak strength at the elastoplastic boundary of plastic zone. Then, the angle of dilatancy decreases in the strain softening stage until it become stable at the residual strength.

Geofluids
In order to improve the understanding and calculation formula of rock dilatancy, many attempts have been made. Alejano and Alonso [20] described a feasible function of peak dilatancy angle of rock according to the equation of peak dilatancy angle of joint that obtained by Barton and Bandis [46]. The proposed function for the peak dilatancy angle of rock mass is described as follows: where ψ p represents the angle of the peak dilatancy, φ p represent the angle of the peak friction, and σ 3 represents the confining stress which refers to σ r in underground roadways. Based on the GZZ strength criterion, the peak friction angle variation with confining pressure can be obtained by using the derivative of peak failure criterion function to confining pressure [20]. In order to get the friction angle, the derivative of peak failure criterion function to confining pressure σ 3 can be described as follows: When combining Equations (10) and (11), the function of peak friction angle is obtained as follows: While incorporating Equation (12) into Equation (9), the function of peak dilatancy angle is reobtained as follows: From the equation, it can be seen that ψ p is directly related to σ 3 as well as the mechanical properties of the rock.
In addition, based on the function of peak dilatancy angle as well as the results obtained by Detournay's [27], the coefficient of the decayed dilatancy is described as follows: where K ψ p is the original dilatancy coefficient.
Equations (9)- (14) indicate that the dilatancy coefficient K ψ changes with the change of η, η * , and σ 3 , and the coefficient decays from the initial peak K ψ p , which also depends on σ 3 as well as the mechanical properties of the rock. Generally, the dilatancy model is nonlinear in the plastic region. Due to the complex dilatancy characteristics of rock mass, Hoek and Brown recommend that the dilatancy angles are ϕ/4, ϕ/8, and 0 are related to the rock mass with good quality, fair quality, and poor quality [19]. According to the recommendation, a more reasonable constant dilatancy model is proposed by Alejano et al. which can be applied to rock mass with different quality [8]. The equation is described as follows: where φ is the friction angle. GSI p is the peak of GSI. Given the strength parameters of GZZ strength criterion, the optimal curve that is suitable for the GZZ strength criterion with

Elastic zone
Softening plastic zone Residual plastic zone 0 r2 r1

Geofluids
confining pressure from 0 to σ rp can be carried out and thus obtained the value of φ [19,23].

The Variable Critical Plastic
Parameter. Through a large number of triaxial compression tests under different confining pressures, the deformation capacity of rock mass after failure is studied. Considering the strain softening characteristics, the strength parameters of rock mass are decreased (as shown in Figure 6). Figure 6 displays the geometry of the strain softening model with confining stress for the estimation of η * , which can be calculated based on the equation described as follows: where E represent Young's modulus; M represent the drop modulus of the postpeak behaviour of rock mass; K ψ is the coefficient of dilatancy. Previous indoor and field test results show that the drop modulus M is directly related to the quality of rock mass quality [19,22] as well as the confining pressure stress [24,48]. Additionally, Alonso et al. [49] propose the function related to M as follows: 3. The Existing Model's Extension

Stress and Displacement Properties in the Elastic Zone.
Based on the assumption of plane strain axial symmetry as well as the elastoplastic mechanics, the mechanical properties of the surrounding rock such as the stress-strain and displacement in elastic zone (r ≥ R p ) is obtained using the function described as follows: In Equation (18), r is the distance between the calculated points and the median of the equivalent circular roadway. R 0 is the radius of the equivalent circular laneway, and σ 0 is the original rock stress. E and ν represent the Young's modulus and Poisson ratio, respectively.

Stress and Deformation Properties in the Plastic Zone.
Based on the plane strain axial symmetry assumption for the equivalent circular roadway, the equilibrium equation is described as follows: In Equation (19), σ r and σ θ represent radial and circumferential stress, respectively.
For the deep equivalent circular roadway, based on Hook's regulation, the tangential strain can be estimated by the stress conditions of the rock mass with the function provided as follows: Since the variable strength parameters and the variables ψ and η * caused by the change of σ 3 in the plastic softening area cannot be solved in closed form, especially for the rock mass that fails in the nonlinear failure patterns. Previous researchers use the methods of finite difference that is developed by Brown et al. [17] to identify the solution to the problem of strain softening in plastic zone [25,50]. In the proposed method, the plastic zone is separated into a finite number of concentric rings, which is shown in Figure 7. In the process of the plastic deformation, the tangential strain ε θ gradually increases while the radial strain increases in a nonflow manner with ψ. As shown in Figure 7, the plastic zone is separated into a finite number of the concentric rings as regards the radial stress [25]. The characteristic of adjacent rings is that the stress increment is very small.
The plastic zone is assumed to be formed by the concentric annuli n which is displayed in Figure 7. The ring of the number of i is defined by two circles of normalized radii ρ ði-1Þ = r ði-1Þ /R p and ρ ðiÞ = r ðiÞ /R p . It should be noted that to satisfy the equilibrium conditions, the thickness of each  6 Geofluids annulus is automatically defined in the numerical process, and its thickness is not equal. For the convenience of calculation, all radius divided by the radius R P of plastic zone are normalized (as shown in Equation (21)).
Then, the radius of the inner and outer boundary of the number i ring are ρ ðiÞ and ρ ði−1Þ , respectively. At the outermost side of the plastic zone (which is the interface of the elastoplastic zone), ρ ð0Þ is equal to 1. At the tunnel wall, ρ ðnÞ = R 0 /R p . It should be noted that the difference between the radial stresses on the inner and outer sides of each ring ðΔσ r = σ rðiÞ − σ rði−1Þ Þ is defined as a constant value. In the study, the equivalent increment of confining pressure is applied for each annulus. This increment is described as follows: According to the calculation of increment of confining pressure using Equation (22), the redial pressure of the ith annulus is obtained using the function described as follows: Based on the equilibrium differential function (Equation (10)), when n is large enough, the tangential stress σ θðiÞ in the plastic zone is estimated with the equation described as follows: If the number of annuli n is sufficiently large enough, the equilibrium equation can be related to the normalized radius ρ = r/R p . For the convenience of calculation, all radius divided by the radius R P of plastic zone are normalized. Then, the radius of the inner and outer boundary of the number i ring is ρ ðiÞ and ρ ði−1Þ , respectively. At the outermost side of the plastic zone (the interface of the elasto-plastic zone) ρ ð0Þ = 1. At the roadway wall, ρ ðnÞ = R 0 /R p . Therefore, the stress and strains at the boundary of the elasto-plastic can be obtained based on Brown et al. [17], which is described as follows: Once the stress components are obtained, the strain parameters can be calculated according to the physical equation. The strain parameters including the elastic and plastic parts can be described as follows: where ε rðiÞ , ε θðiÞ , ε rði−1Þ , and ε θði−1Þ represent the radial strains and tangential strains at the ith and ði − 1Þth radius. Δε e rðiÞ and Δε e θðiÞ represent the increments of the elastic radial and tangential strain at the ith radius. Δε p rðiÞ and Δε p θðiÞ represent the increment of the plastic radial and tangential strain at the i th radius, respectively.
According to the Hooke's Law, the increments of the elastic strain can be calculated by the increments of the stress using the function as follows: where Δσ θðiÞ represent the increment of the tangential stress, Δσ θðiÞ = σ θðiÞ − σ θði−1Þ . The increments of the plastic strain obey the nonassociated flow regulation as follows: When both K ψ and η * are related to σ rðiÞ , K ψðiÞ can be calculated using the equation described as follows:

Geofluids
Combining Equations (27) and (29), the following equation can be obtained: Then, the finite difference form of the equilibrium function (Equation (19)) can be expressed as follows: By combining Equation (18) and Equation (29), the rela-tionship between r ðiÞ and r ði−1Þ can be obtained as follows: For the purpose of solving the problem of the strain components, Equation (20) can be rewritten as follows: When combining Equations (28), (31), and (34), the radial displacement u ðiÞ is obtained using the function described as:

Calibration Examples
In order to verify and confirm the accuracy of the proposed methodologies, the results obtained from the proposed method are compared to that obtained using other methods. The verification of the proposed methods in this study is separated into three parts. The mechanical properties of rock mass used in the study are assumed to meet the strength criterion with strain-softening behavior. The first part focuses on showing the advantages of three-dimensional H-B strength criterion (GZZ), compared with two-dimensional H-B strength criterion (narrow H-B criterion and generalized H-B criterion). On the other hand, it focuses on the comparison of the current three-dimensional D-P criteria to show its correctness and rationality of the theoretical derivation.
Firstly, the differences between the three-dimensional H-B criterion and the two-dimensional H-B criterion (narrow H-B criterion and broad H-B criterion) are discussed. A series of data are selected from the publication by Sharan [51], which are R 0 = 5 m, σ 0 = 30 MPa, p i = 5 MPa, E = 5:5 GPa, ν = 0:25, σ c = 30 MPa, m p = 1:7, s p = 0:0039, m r = 1:0, s r = 0:0, a = 0:5, ψ = 30°, n = 500. The comparison between the three-dimensional H-B criterion and the twodimensional narrow H-B criterion is carried out, and the specific change rule is shown in Figure 8. In addition, when a ≠ 0:5, the difference between the three-dimensional H-B criterion and the two-dimensional generalized H-B criterion is compared and verified under the condition that the above data are kept unchanged. The specific change rule is shown in Figure 9.
It could be shown from Figures 8 and 9 with the increase of the ratio value of radius at a certain point to the radius of the equivalent circular roadways increases, the radial stress increases approximately linearly, while the circumferential stress increases approximately linearly first and then decreases gradually. As shown in Figures 8(a) and 9(a), the radial stress of GZZ criterion is generally lower than that of 2D H-B criterion obtained by Sharan [51], but the linear growth rate of GZZ criterion is higher, which highlights the role of intermediate main stress. As shown in Figures 8(b) and 9(b), the tangential stress of GZZ criterion is generally lower than that of 2D H-B criterion obtained by Sharan [51]. On one hand, the turning point of GZZ criterion is relatively later than that of 2D H-B criterion. On the other hand, the subsequent decreasing trend range of 2D H-B criterion is very small, which also shows the role of intermediate main stress. As shown in Figure 10, the threedimensional H-B criterion and the three-dimensional D-P criterion were compared. Both of the calculation results produce the plastic zone, as shown in the dotted line part in the figure. The three-dimensional solution enters the state of elastic stress at the shallow depth of roadway, which is beneficial to the stability of roadway. This is mainly due to the favorable mechanical action of the intermediate main stress. With the decrease of the plastic zone's radius, the stress gradient in plastic zone and the stress value at the same depth are increased. In the plastic zone, the intermediate main stress increases with the increase of surrounding rock depth and converges to the original rock stress at the elastic-plastic interface.
In order to better verify the theoretical derivation, the calculation results based on this derivation are compared with the existing research results. Li et al. [52] deduced the analytical solution of circular tunnel based on GZZ strength criterion, but it did not consider the strain softening characteristics of surrounding rock. If ω p in Equation (6) is equal to ω r , the solution in this paper can be reduced to the solution without considering the strain softening characteristics of surrounding rock. In order to study the accuracy of the theoretical derivation, this solution is compared with the solution provided by Li et al. [52]. Figure 11 shows the comparison of radial and circumferential stress curves of surrounding rock under different uniaxial compressive strength. Table 1 shows the comparison of plastic zone radius of surrounding rock under different uniaxial compressive strength. The uniaxial compressive strengths are 20, 30, 50, and 70 MPa, respectively. It can be seen that the calculation results in this paper almost coincide with the results given by Li et al. [52], and the radius of plastic zone is the same. Therefore, the correctness and rationality of the theoretical derivation are verified.

Discussion
Three-dimensional nonlinear GZZ strength criterion considered the influence of intermediate main stress. In this part, According to the analysis of Xia et al. [53], the softening law of the surrounding rock strength parameters with strain involved in the GZZ guidelines is described as follows. The value a is increased with the increase of the tangential strain. When the tangential strain is between 0 and 0.001, the value a remains at 0.5; then, with the increase of strain, the value a increases linearly to 0.6, and the tangential strain is 0.003. After that, the value a remains constant with the increase of the tangential strain; the s value is increased with the increase of the tangential strain. When the tangential strain is between 0 and 0.001, the value s remains at 0.0004. With the increase of strain, the s value decreases linearly until it is 0, and the tangential strain is 0.003. After that, the s value remains constant with the increase of the tangential strain; the m b value is increased with the increase of the tangential strain. When the tangential strain is between 0 and 0.001, the m b value remains at 0.5. Then, with the increase of strain, the m b value decreases linearly to 0.25, and the tangential strain value is 0.003. After that, the m b value remains constant with the    10 Geofluids increase of the tangential strain; the ψ value is increased with the increase of the tangential strain. The softening law is that when the tangential strain is between 0 and 0.001, the ψ value is kept at 15°. With the increase of strain, the ψ value decreases linearly to 5°, and the tangential strain is 0.003. After that, the ψ value remains constant with the increase of the tangential strain.
In the sensitivity analysis, it is determined that the relevant parameters of the deep equivalent circular roadway are the original rock stress is p 0 = 15 MPa, the excavation radius of the equivalent circular roadway is R 0 = 5 m, E = 20 Gpa, μ = 0:35, the strength parameter m i = 6, the surrounding rock uniaxial compressive strength σ c is 20 MPa, the initial value of the geological strength index GSI p is 50, and the residual values GSI r are 10, 20, 30, 40, and 50, respectively [53]. When GSI r is 50, surrounding rock softening is not considered which is as the control group. Assuming the support reaction p i is 1.5 MPa, the critical plastic deviator strain η * of surrounding rock is 0.5%.

Impact of Rock Dilation.
Dilatancy plays an important role in the deformation of the rock boundary of roadway excavation [54]. Previous studies described that the expansion of surrounding rock of the roadway has little impact on the stress distribution [55,56]. However, the yield process of surrounding rock is accompanied by the expansion behavior. The volume of plastic zone keeping constant during the yield process is not consistent with the actual situation of rock mass [57]. In this paper, different dilatancy angles are introduced into the proposed theoretical calculation equation, and the influences of expansion on stress distribution are studied.
In addition, the main focuses of this paper are on the stability analysis of the roadway excavation process, and the rock mass involved is basically weak rock formations. Previous studies have divided rocks into soft rocks and hard rocks [58]. Since the mechanical characteristics of each rock are similar, the mechanical parameters of mudstone samples are selected to calculate the difference of stress and radius of plastic zone under different dilation angles. Compared with the change of main stress with dilation angle in Table 2, the sensitivity of stress distribution to dilation angle is weak. Considering the three-dimensional case of intermediate main stress, the dilation angle has little effect on the plastic zone of mudstone. It also verifies the conclusion that the dilation angle has little effect on the plastic zone stress [55,56]. The results indicate the influence of the intermediate main stress and the strength potential of rock materials. The self-bearing capacity of materials is also considered to make sure that the analysis of the stability of the surrounding rock is reliable. Figure 12 shows the distribution curve of displacement in plastic zone along the depth of surrounding rock under different dilatancy angles. It can be seen from Figure 12 that the displacement in the plastic zone of soft rock increases with the increase of dilation angle, which is consistent with the engineering practice. For soft rock, within 3 times of excavation radius, the expansion is very significant. When the dilatancy angle is large, the distribution of displacement in plastic zone along the depth of surrounding rock presents prominent nonlinear characteristics. The surrounding rock around the inner wall of the roadway is disturbed by the roadway excavation, and the strain energy stored in rock mass is also released along the disturbed area. As a result, the displacement distribution in plastic zone increases nonlinearly with the decrease of surrounding rock depth. The displacement of the inner wall of the tunnel develops the most. Therefore, for soft rocks such as mudstone, the plastic zone's displacement was affected by the increasing dilation angle, so attentions are needed to be paid to the influence of rock dilation in the calculation of stress distribution. Figure 13 shows the stress distribution of surrounding rock under different softening degree. For the radial stress (σ r ), under the same support pressure (p i ), the higher the softening degree of surrounding rock, the smaller the radial stress. However, due to the influence of the size of the plastic zone, the distribution of the tangential stress σ θ in the elastic and plastic zone is not the same. In the inner wall of the roadway, the higher the softening degree of the surrounding rock, the smaller the tangential stress. It can be seen from Figure 13(b) that the softening degree of surrounding rock has a great influence on the scope of plastic area. The higher the softening degree, the greater the plastic zone of surrounding rock. However, regardless of the degree of softening, the maximum value of tangential stress is the same under various conditions. In other words, under the same support pressure (p i ), the strain softening of surrounding rock can change the stress distribution of igure 12: Effect of dilation angle on plastic zone of roadway surrounding rock displacement of the soft rock. 11 Geofluids surrounding rock, but cannot change the values of maximum stress. For the surrounding rock with high softening degree, the stress changes slowly due to the large radius of plastic zone. However, regardless of the degree of softening, the stress of surrounding rock at the boundary of plastic zone is constant.

Effect of Softening Degree.
In terms of the strain of surrounding rock, it can be seen from Figure 13 that the strain softening of surrounding rock has a great influence on its strain, and the strain of surrounding rock in plastic zone is particularly sensitive to softening. Under the same support pressure (p i ), when GSI res are 40, 30, and 20, the radial strain at the inner wall of the tunnel is 1.2, 1.6, and 3.4 times of that without strain softening, respectively. When GSI res is 10, the radial strain at the inner wall of the tunnel is 18.6 times larger than that under the condition of without strain softening. The change regulation of tangential strain is similar to that of radial strain. It can be seen that after the surrounding rock softens, although the stress in the plastic zone is relatively smaller, its strain increases exponentially, resulting in a large increase in the deformation of the surrounding rock. Figure 14 shows the characteristic curve of surrounding rock under different softening degree. As can be seen from Figure 14, when the support pressure is large, the softening effect can be ignored due to the small deformation of roadway. However, with the decrease of support pressure, the deformation of surrounding rock changes significantly with the degree of softening. Under the same support pressure, the deformation of surrounding rock is several times or even dozens of times. For example, when the support pressure (p i ) is 1.5 MPa, the displacement of the inner wall of the roadway is 78 mm without considering the softening of surrounding rock. When GSI res are 40, 30, and 20, the radial displacement at the inner wall of the roadway is 110 mm, 180 mm, and 500 mm, respectively. When GSI res is less than 20, the displacement of the inner wall of the roadway is more than a few meters, which means the roadway lose its stability. In the same way, under the condition of the same deformation of roadway in the inner wall of the roadways, the support pressure required for the surrounding rock with different softening degree varies several times.

Conclusions
In this paper, based on GZZ strength criterion considering the impact of intermediate main stress, an equivalent circular roadway solution method which can reflect the strainsoftening and dilatancy properties of surrounding rock is proposed. The stress, strain, and deformation of surrounding rock after the excavation of underground roadways are obtained by numerical method. The influence of the intermediate main stress on the mechanical behavior of the plastic  12 Geofluids zone of the surrounding rock and the relationship between the mechanical properties of the surrounding rock and the softening parameters are studied. Then, the failure mechanism of roadway with soft rock is analyzed, and the main conclusions are as follows: (1) Compared with the traditional H-B strength criterion, the value of surrounding rock stress calculated by GZZ strength criterion considering the effect of intermediate main stress is smaller, while the radius of plastic zone, softening zone, and surrounding rock strain are larger. Therefore, it is more conservative to use the GZZ strength criterion for calculation (2) Under the same support pressure, the strain-softening of surrounding rock can change the stress distribution, but cannot change the value of the maximum stress. In other words, regardless of the degree of strain-softening, the stress of surrounding rock at the boundary of plastic zone is constant. But for the surrounding rock with high strain-softening degree, because of its large plastic zone radius, the surrounding rock stress near the inner wall of the tunnel is relatively small, and the displacement is large. When the support pressure is small, the deformation of surrounding rock changes significantly with the degree of strain-softening. Under the same support pressure, the deformation of surrounding rock is several times or even dozens of times.
Similarly, under the condition of the same deformation of the surrounding rock in the inner wall of the roadways, the required supporting reaction force may increase several times or even dozens of times considering the soften of the surrounding rock (3) The dilatancy properties of rock have significant effects on the displacement of the plastic zone of the surrounding rock mass. Under the effects of dilatancy, the displacement of the plastic zone is obviously larger than that without considering the effects of dilatancy. The displacements of the plastic zone of the rock mass increase with the increase of the dilatancy angle (4) The stress of surrounding rock is redistributed due to the excavation and unloading of roadway. Shear slip failure of rock structure occurs under high ground stress, and the self-supporting capacity of the surrounding rock is greatly reduced. The deformation of surrounding rock of roadway increases rapidly when the support effect is enough. It is the main reason of large deformation and failure of deep underground roadway with soft rock. Therefore, in the design and calculation of roadway supporting system, the strain softening characteristics of surrounding rock should be properly considered to avoid large deformation and failure of roadway

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare no conflict of interest.