A Simple Solution and Dilatancy Characterization for a Circular Tunnel Excavated in the Nonlinear Strain-Softening and Nonlinear Dilatancy Rock Mass

,


Introduction
Te elastoplastic analysis of underground excavation or opening under axisymmetric conditions has been a classical mechanical issue [1][2][3], and a considerable number of studies have been performed on rock stress and strain determination [4,5]. Tese achievements are utilized to develop ground reaction curves (GRCs) for tunnels with diferent shapes, excavated in rock masses of diferent qualities with varying stress states [6,7], which are further used for tunnel design and response analysis [8,9].
Te stability or elastic-plastic state of rock masses are primarily descripted by various strength criteria [10], while the majority of the elastoplastic analysis was performed on the circular tunnels excavated in the rock materials compatible with the two-dimensional Mohr-Coulomb (M-C) or Hoek-Brown (H-B) criterion [11,12]. However, numerous fndings from in situ stress measurements demonstrate that deep openings are usually dominated by the pronounced three-dimensional unequal stress conditions [13][14][15]. Also, most studies are conducted on the basis of the consideration of the maximum and minimum principal stresses, where the efect of the intermediate principal stress efect on mechanical properties of the rock mass has been neglected. Tis is detrimental to the design of the tunnel support system [16,17].
Additionally, other achievements obtained from the experiments [18], numerical simulations [19], and theoretical analyses [20] emphasize the role of the intermediate principal stress in the strength potential and self-stabilizing efect of the rock masses. For example, Kabwe [21] achieved the consideration of intermediate principal stress efect in the proposed numerical solution by the inclusion of the stress Lode parameter in the circumscribed Drucker-Prager (D-P) yield criterion and further suggested that yield criteria accounting for the intermediate principal stress efect should be employed in analytical solutions to gain an accurate ground reaction estimation. Some other researchers also noted that the rock mass strength depends entirely on the intermediate principal stress, which acts as the confning stress [18,22]. Tese also refect the problem of conservativeness of tunnel design based on the two-dimensional theory [23]. Terefore, a strength criterion with a nonlinear strength envelope is indispensable for accurate descriptions of rock mechanical behaviors.
Te three-dimensional USC proposed by Yu [24] for geotechnical materials can address the limitations that traditional yield criteria fail to consider the intermediate principal stress efect. Tis criterion was derived from the twin-shear element unifed strength theory and takes into account the efects of all the three principal stresses [25]. Tereafter, the criterion is adopted for elastoplastic analysis of a tunnel excavated in elastic-perfectly-plastic (EPP) rock mass, in which the intermediate principal stress was assumed to act along the tunnel axis [26]. Zhang et al. [27] proposed solutions for radial displacement and support force of a tunnel built in the USC rock mass. Subsequently, the USC yield criterion has been successfully employed in studies, e.g., [28,29]. As for a yield criterion considering all the three principal stresses, the USC has proved to be an alternative. In the follow-up studies, the double slip angle phenomenon occurs in twin-shear element USC for some specifc stress states [30]. To overcome this unfavorable situation, Hu and Yu [30] further proposed the triple-shear element USC (TS-USC), which takes into account the efect of all the principal shear stresses on the yield failure. It provides an efective solution for mechanical response analysis of geotechnical materials [31,32].
In rock mechanics and geomechanics, precise descriptions of rock mass behaviors remain the most challenging task, especially for the mechanical properties in the postfailure stage, as it signifcantly afects the stability of the rock mass [33]. Although there are lots of analytical solutions or approaches to characterize the elastic-plastic response of the rock masses, most of them are based on the EPP [17,34] or elastic-brittle (EB) [35,36] rock masses. Nevertheless, a large number of underground excavation or opening involves the soft rock strata, where these average quality rock masses are characterized by the strain-softening [37,38] and dilatancy behaviors [39,40]. To analyze the strain-softening behaviors of the rock masses in a tunnel, Lee and Pietruszczak [41] and Park et al. [42] proposed methods to divide the potential plastic zone into concentric rings. Based on these methods, many research studies have been performed on the deformation and failure of the rock masses [43][44][45]; however, the dilation angle in these solutions was simplifed as a constant or linear parameter. Dilatancy is closely associated with the process of deformation and failure of the rock masses [46,47], and thus, the dilation angle should not be oversimplifed in the elastoplastic analysis.
Yi et al. [48] also stated that strain-softening and dilatancy behaviors in rock masses bring a signifcant synergistic efect, and considerable errors would be induced if neglecting the strain softening and dilatancy during a deep tunnel analysis. Terefore, it is more appropriate to investigate rock displacement laws and tunnel stability with considering the strainsoftening and dilatancy behaviors of rock masses.
As a result of this review, it turns out that the intermediate principal stress efect, strain softening, and dilatancy all exert notable infuence on the tunnel stability and should be evaluated for accurately computing the stresses and displacements around the excavated tunnels. Tis study mainly aims to develop a simple numerical solution for displacements and stresses determination of the tunnels surrounded by the rock mass exhibiting the strain-softening and dilatancy behaviors. Te present paper is organized as follows. Te second section states the derivation framework. In this section, the TS-USC, the nonlinear strain-softening model, and the nonlinear dilatancy model are applied to descript the intermediate principal stress efect as well as the postfailure behaviors of the rock masses. In the third section, with reference to the fnite diference method proposed by Lee and Pietruszczak [41], the numerical solutions for the stresses and displacements are derived for the strainsoftening and dilatancy rock masses. In the fourth section, some published rigorous algorithms, as well as the illustrative example, are referred and compared to confrm the validity of the proposed solution. In the ffth section, a parametric analysis is conducted to evaluate the infuence of some typical parameters on the dilatancy response of rock masses, and some discusses and suggestions are then presented. Te study has been technically concluded in the last section.

Problem Defnition.
Te problem is defned in Figure 1, where a circular tunnel with a radius of R 0 is subjected to hydrostatic stress feld σ 0 . Te support structure provides a uniformly distributed radial support force p i at the excavation surface, and the radii of plastic and residual zones in the rock mass are R p and R s , respectively. Te tangential and radial stresses at the elastic-plastic interface of the rock mass are σ θp and σ rp , respectively. When p i is less than p i * (p i * is the support force at the interface of the softening and residual zone), the softening zone continues to develop and then yields the residual zone. Te rock strength in the softening zone gradually decreases to the residual value. Te tangential and radial stresses at the softening and residual zone interface are σ θs and σ rs , respectively.

Triple-Shear Element Unifed Strength Criterion.
Te TS-USC comprehensively captures the basic strength properties of rocks, such as the tensile and compressive inequalities and the intermediate principal stress efect. Te double slip angle problem observed in twin-shear element USC for some specifc stress states is also overcome. It is expressed as follows [30]: where σ 1 , σ 2, and σ 3 are the maximum, intermediate, and minimum principal stresses; σ t � (2c cos φ)/(1 + sin φ), α � (1 − sin φ)/(1 + sin φ); c and φ are the cohesion and the angle of internal friction; and b is the intermediate principal stress infuence parameter. Te values of b that satisfes the externally convex type of plane π-limit line is as follows [24]: σ 2 satisfes the following conditions [26]: Combination of equations (1) and (3) gives the expression for the TS-USC under plane strain condition as follows:

Nonlinear Strain-Softening Behavior.
Neither the EPP nor the EB model is utilized to well describe strain-softening behaviors exhibited in average quality rock mass during tunnelling. Terefore, the softening parameter is introduced to characterize this behavior as follows [40]: where η is softening parameter, and ε p θ and ε p r are tangential plastic strain and radial plastic strain, respectively.
When performing the theoretical model derivation, the tangential stress σ θ and radial stress σ r in the rock mass are assumed to be equal to the σ 1 and σ 3 , respectively. By combining equation (4), the TS-USC under plane strain state that considering the η is given as follows: Alonso et al. [49] stated that a circular opening can be modeled as a gradual stress reduction from the initial state σ 0 to support force p i . Figure 2 gives the stress path of a circular opening. η * is the critical softening parameter of the rock mass, and it is an indicator to determine the transition of the rock mass from the softening state to the residual state [49,50]. When η � 0, the rock mass is in the elastic state, which corresponds to the initial ground stress Advances in Civil Engineering state at point L. Te radial stress at the elastic-plastic interface (point M) is σ rp , which can be solved by the elastic solution. Te rock mass in the zone MN is in softening state, that is, 0 <η < η * . Ten, when η ≥ η * in the zone NO, the rock is in the residual state.
Introducing the softening factor δ, an adopted nonlinear softening model is expressed as follows [51]: where ω p and ω r represent the peak and residual values of rock mass strength parameters c and φ, respectively. For the δ, it is varying from − 1 to 1. A linear softening model is obtained when δ � 0 and a nonlinear softening model for other δ values, while the maximum upward and downward curvatures are obtained when δ � 1 and − 1, respectively [51]. Based on the incremental theory of plasticity, the plastic potential function can be expressed as follows [52]: where K ψ is the dilatancy coefcient, and it can be obtained from the following equation [41]: where ψ is the dilatancy angle.

Nonlinear Dilatancy Behavior.
Te ψ in rocks is not the constant, as it decreases with the increased confning pressure [53]. Alejano and Alonso [50] proposed a model for determining the peak of dilation angle by ftting the experimental data.
where σ c is the unconfned compression strength.
For a deep circular tunnel, Detournay [54] noted that K ψ decays exponentially from the peak value K p ψ and proposed a ftted expression as follows: By combining equations (9) and (11), ψ is obtained as follows: Substituting equations (10) into (12), the nonlinear dilatancy model is obtained as follows:

Solution for the Elastic Rock Mass.
Te radial stress σ e r , tangential stress σ e θ , and radial displacement u e r in the elastic zone are expressed as follows [55,56]: Elastic zone

Sofening zone
Residual zone Figure 2: Excavation stress path and softening parameter region conversion [49].

Advances in Civil Engineering
where E is the modulus of elasticity, and μ is Poisson's ratio.

Solution for the Plastic Rock Mass.
Te equilibrium equation is as follows [42]: Te radial and circumferential strains of the rock mass under small strain conditions are as follows [41]: Te stress at the elastic-plastic interface can be obtained as follows: Substituting equations (17) into (6) yields the radial stress at the elastic-plastic interface as follows: Te fnite diference method is employed to derive stress and strain of the rock mass. Figure 3 gives a schematic diagram of radial stress evolution of the rock mass. Te plastic zone is divided into n circles, where σ θ(0) and σ r(0) are tangential and radial stresses at the elastic-plastic interface, and σ θ(n) and σ r(n) are the tangential and radial stresses at the excavation surface. Te inner and outer boundaries of the ith ring correspond to the radii r (i− 1) and r (i) , respectively, and the stress components corresponding to the inner and outer boundaries are σ θ(i− 1) , σ r(i− 1) and σ θ(i) , σ rt(i) , respectively.
As shown in Figure 3, the stress increment in each ring of plastic zone is as follows: Te radial stress in each ring is calculated by Substituting equations (20) into (6), the tangential stress of each ring can be given as follows: Each ring corresponds to diferent rock strength parameters c and φ can be calculated by the following equation [51]: Te equation (15) in fnite diference form can be expressed as follows: Separating r (i− 1) and r (i) , the ratio of the radii of the inner and outer boundaries in the i-th ring expressed in terms of stress is obtained as follows: Te elastic strain increment of each ring in the plane strain condition can be derived from the corresponding stress increment of each ring [57].
where ∆ε e r(i) and ∆ε e θ(i) are radial and tangential elastic strain increments, respectively.
Te incremental η in the i-th ring is expressed as follows: It can be obtained from the nonassociated fow rule as follows [57]: where ∆ε p r(i) and ∆ε p θ(i) are radial and tangential plastic strain increments; and K ψ(i) is the corresponding dilatancy coefcient of the i-th ring.
Te total strain in each ring is the total of elastic and plastic strains, then the individual strain components in the i-th ring are expressed as follows: By combining equations (27) and (28), the following relationship between the radial strain ε r(i) and the tangential strain ε θ(i) in the i-th ring is obtained as follows:

Advances in Civil Engineering
Furthermore, and ε r(i) and ε θ(i) in the i-th ring can be obtained as follows: Te compatibility equation of deformation is as follows [58]: Te fnite diference form of equation (31) in the i-th ring is as follows: Te ratio of the inner and outer boundary radii of the i-th ring expressed in terms of strain is as follows : Substituting equations (30) into (33) yields Te plastic strain increment for the i-th ring is as follows: Combining equations (25), (34), and (35), the expression for η of the i-th ring is obtained as follows:

Plastic Radius.
Te boundary conditions in the solution are as follows: Combining equation (14) yields the residual radius as follows: Superimposing n equation (33) yields Te radius of plastic zone is as follows: Furthermore, the radius of the i-th ring is as follows: 3.4. Displacement in the Plastic Zone. Substituting equations (35) and (41) into (16), the radial displacement of the i-th ring can be obtained as follows:

Verification Examples and Comparison
Te results in the literature [41,49] are referred to verify the soundness of the solution proposed in this paper. Te parameters used for validation are as follows: σ 0 , R 0 , E, and μ are 20 MPa, 3 m, 10 GPa and 0.25, respectively. φ p and φ r are 30°and 22°, while c p and c r are 1.0 MPa and 0.7 MPa, respectively. η * is fxed as 0.008. It is worth noting that the TS-USC of the proposed solution is converted to the M-C yield criterion (by setting b to 0), for consistency with that in the literature [41]. A comparison of the results obtained from the two solutions is displayed in Figure 4, which verifed the validity of the solution proposed in this study.
Also, other interesting attempts are performed to solve the tunnel unloading problem and investigate the infuence of the postpeak drop modulus on elastoplastic behaviors of the rock mass [59][60][61]. Here, the results reported in [60,61] were chosen as representative cases for comparison. According to Alejano et al. [60], three strain-softening models with considering the constant/variable drop modulus (M) and the constant/variable dilatancy (ψ) were proposed for the purpose of characterizing the postfailure behaviors of the poor, average, and good quality rock masses. Te results in terms of p i and rock displacement u 0 , with a reference to the average quality rock mass type, are depicted in Figure 5 gained by the proposed model are all larger than those determined by the abovementioned three models. It delivers a clear fact that the deliberate characterization of the postfailure behaviors is essential for tunnel stability analysis. For comparison purposes, the computed rock displacements for the unloaded case by the González-Cao et al. [61] together with results obtained with the proposed model under two diferent b values are represented in Figure 5(b). It can be seen that u 0 values under the same p i levels determined by the proposed model is also slightly larger than those calculated by the referred model. Obviously, the intermediate principal stress efect contributes to the selfstabilization of the rock mass. Tis can be confrmed by the decrease of u 0 values when b increases from 0 to 1. Te rock displacement diferences between the models decrease with the increased GSI values (which means the improvement of the rock mass quality). It is notable that, for diferent quality rock masses, some errors might arise if stability analysis is not based on the suitable models, particularly for the poor to average quality rock masses. Terefore, to obtain a reliable rock mass response, it is crucial to develop reliable models to descript rock behaviors.

Parametric Analysis on Dilatancy Behavior of the Rock Mass in the Plastic Zone
Parametric analysis on dilatancy behaviors of the rock mass in the plastic zone is conducted in this section. Te relevant parameters utilized for investigation are referred from the literature [60]. Te details are listed in Table 1.

Infuence of b.
When b is the variable, ψ under diferent p i , η * , and δ show similar curve characteristics; so, a representative scenery of p i � 0, η * � 0.006, and δ � 0.5 are selected to analyze the variation feature of ψ in the rock mass, as shown in Figure 6. Te red circle in the fgure is the demarcation point of the softening zone and the residual zone of the rock mass. It can be observed that, for diferent b, ψ tends to decrease from deep rock mass to the excavation surface, and the maximum ψ value is found at the elasticplastic interface. ψ in the softening zone is signifcantly larger than that in the residual zone, indicating that dilatancy efect in the softening zone is more notable than that in the residual zone. Te reduction rate of ψ in the softening zone is also considerably higher than that in the residual zone. It is evident that intermediate principal stress only diminishes the plastic radius but also determines the progression pattern of ψ with depth variations of the rock mass. Te magnitude of ψ at the elastic-plastic interface is dependent on b. When b � 0, ψ at elastic-plastic interface is 5.85°; ψ at the interface increases with the increased b; and ψ acquires the maximum value of 7.47°when b � 1. As the intermediate principal stress efect increases (i.e., the increasing b values), the peak values of ψ increases; also, the dilatancy efect of the rock mass is more signifcant when the intermediate principal stress is considered than when it is not considered.

5.2.
Infuence of η * . Te thresholds of η * is generally in range of 0.001-0.01 [62]. Terefore, η * are set as 0.001, 0.002, 0.004, 0.006, 0.008, and 0.01, respectively, in the analysis. Te results of ψ with varying η * for the selected scenery of p i � 0,  Figure 7. It confrms that ψ at the elastic-plastic interface is constant (all values are 6.59°) for diferent η * , indicating that the magnitude of ψ is not infuenced by η * . Similarly, ψ decreases gradually from deep rock mass to the excavation surface, but the decay rate varies signifcantly under diferent η * . Te decay rate of ψ increases sharply with the decrease of η * . Te smaller the η * , the less the range of the softening zone, which causing the residual zone to occupy a larger proportion. For the rock mass with the η * value of 0.001, it clearly observes a drop-down decrease in ψ from the peak, refecting brittle failure caused by large plastic of the rock mass at the elasticplastic interface.

5.3.
Infuence of δ. δ controls postpeak softening behavior of the rock mass, and δ values of − 1, − 0.5, 0, 0.5, and 1 are selected to analyze evolutionary characteristics of ψ with diferent softening behaviors. Te results of ψ for p i � 0, η * � 0.006, and b � 0.4 are displayed in Figure 8. δ of the rock mass shows an efect on the size of the plastic radius, which gradually increases from 35 m to 37.8 m as the δ goes from − 1 to 1. Te fxed value of 6.59°for ψ at the elastic-plastic interface indicates that δ does not infuence the magnitude of ψ. For diferent δ, the variation pattern of ψ in the residual zone is basically the same in the radial direction, and they all exhibit a two-stage approximately linear decrease. In the softening zone, the variability of ψ is signifcant. Te rate of reduction of ψ, compared to the linear softening (δ � 0), is more drastic for depressed nonlinear strain-softening behavior (δ � 0.5 and 1). ψ in the softening zone under convex nonlinear strain-softening behavior falls of rapidly as it approaches the residual zone.

Infuence of p i .
Te infuence of varying p i on ψ, for the scenery of η * � 0.006, b � 0.4, and δ � 0.5, are displayed in Figure 9. ψ tends to decrease from deep rock mass to the excavation surface, while ψ at the elastic-plastic interface is constant, i.e., 6.59°. ψ at residual-softening interfaces under diferent p i are also the same, indicating that support force does not control the magnitude of the ψ values in plastic zone of rock mass for the given η * , b, and δ values.

Conclusion
(1) Te TS-USC is introduced to establish an axisymmetric elastoplastic fnite diference solution of a deep circular opening, which considers the intermediate principal stress, nonlinear dilatancy, and nonlinear strain-softening behaviors of the rock mass. Te solution introduces b and δ to characterize the action degree of intermediate principal stress and strain-softening behavior of the rock mass. Te variations of the introduced parameters can conveniently refect the variability of derivation results controlled by diferent yield criteria, dilatancy, and strain-softening behaviors, which are of wider theoretical and practical application. (2) Parametric analysis revealed that the maximum value of ψ under diferent b occurs at the elasticplastic interface. Te values and reduction rates of ψ in the softening zone are signifcantly larger than those in the residual zone. Te peak value of ψ increases with the increased b, indicating that intermediate principal stress mainly afects the peak value of ψ, and the dilatancy efect of rock mass is more pronounced when considering the intermediate principal stress than when it is not. (3) Te magnitude of ψ at the elastic-plastic interface is not afected by η * , but the decay rate of ψ due to the change of η * varies signifcantly. Te decay rate of ψ dramatically increases as η * decreases. Te smaller the η * , the smaller the softening zone and the more the residual zone. Te rock strength decreases rapidly, resulting in a drop-type reduction of the ψ from the peak when η * � 0.001, while causing brittle failure with large plastic deformation. (4) Te softening behavior of the rock mass afects the magnitude of the plastic radius but does not afect values of ψ at the elastic-plastic interface. Te variation pattern of ψ in the residual zone is similar for diferent δ, but the diference in variation of ψ in the softening zone is obvious. (5) For a given intermediate principal stress state, critical softening parameter and softening behavior of the rock mass (refecting by parameters of b, η * , and δ), the value and development form of the ψ in the plastic zone is independent of the support force.

Data Availability
Te data used to support the fndings of this study are included within the article.

Conflicts of Interest
Te authors declare that there are no conficts of interest.  Advances in Civil Engineering