Elastoplastic Analysis for Circular Tunnel Based onModified Lade Criterion considering Strain Softening and Dilatancy

Modified Lade criterion can not only describe the strength properties of many kinds of rocks well but also has simple and practical parameters. Although the elastoplastic solution of circular tunnel has been extensively investigated, the method based onmodified Lade criterion considering the effect of the intermediate principal stress, strain-softening behavior, and dilatancy has not yet been studied. In this paper, a new numerical procedure based on modified Lade criterion is proposed to calculate the elastoplastic solutions for surrounding rock of the circular tunnel. *e comparisons of stress, displacement, and plastic zone radius are carried out between the presented method and published literatures under axisymmetric and nonaxisymmetric original in situ stress field. Finally, a series of parametric analyses are executed and discussed. It can be concluded that the lateral pressure coefficient, λ, influences both the size of plastic zone and the development direction. *e plastic zone radius shows a negative power function change with increasing critical deviatoric plastic strain and increases slightly with increasing dilation angle, ψ.


Introduction
Elastoplastic analysis of circular tunnels is one of the most fundamental problems in the fields of tunneling, mining, and geotechnical engineering. After tunnel excavation, the original state of stress equilibrium is disrupted, which causes redistribution of the stress of surrounding rock, elastoplastic deformation, or even failure at the same time. To ensure the safety of tunnel excavation, the stress and displacement solutions around a tunnel are required in the design stage of a tunnel. Based on Mohr-Coulomb (M-C) strength criterion, Fenner [1] first derived the analytical solution of circular tunnel surrounding rock, which was then modified by Kastner considering the nonhydrostatic stress field in 1949 [2]. For a long time, M-C and Hoek-Brown (H-B) failure criteria have been employed to carry out the elastoplastic analysis for the circular tunnel. Ogawa and Lo [3] investigated the influences of dilatancy and failure criteria (M-C and H-B criteria) on stress and displacement for a circular tunnel. Brown et al. [4] presented two new elastoplastic solutions for circular axisymmetric tunnel problem with H-B criterion. Sharan [5][6][7] took various rock mass behaviors into account to derive elastoplastic solutions with H-B and generalized H-B criterion. Alonso et al. [8] presented a one-dimensional numerical solution for circular tunnels with strain-softening surrounding rock. Lee and Pietruszczak [9] proposed a new numerical method accommodating both the M-C and H-B criteria to execute the elastoplastic analysis for circular opening with strain-softening rock mass. Based on the idea of the simplification of strain-softening process, Wang et al. [10] treated the strainsoftening process as a series of stress drop and plastic flow with increasing plastic strain and proposed a new approach to calculate the elastoplastic solutions with M-C and H-B criteria. According to similar ideas, Zhang et al. [11] assumed the surrounding rock properties were uniform in a pretty small region and presented a multistep brittle-plastic model. Cui et al. [12] introduced a variable critical plastic softening parameter to carry the elastoplastic analysis of a circular opening in strain-softening rock mass.
Although various exact solutions or numerical solutions have been proposed in the above works, the influence of intermediate principal stress was ignored in the strength criteria, which has significant effect on the mechanical behaviors of surrounding rock [13][14][15]. erefore, the yield criteria considering intermediate principal stress were gradually adopted to investigate the elastoplastic solutions of tunnel surrounding rock and achieved significant results, such as Drucker-Prager (D-P) criterion [16,17], the unified strength theory [18][19][20][21][22], and Lade-Duncan (L-D) criterion [23]. However, the assumption of axisymmetric in situ stress field is still introduced to calculate the elastoplastic solution of rock mass, which is not in line with the actual stress state under most circumstances [24]. In addition, the parameters of these strength criteria are difficult to be obtained. Modified Lade criterion (MLC) was proposed by Ewy [25] based on original Lade's criterion (LC) [26,27], which is a simplified version of LC and is compelled to agree with M-C strength criterion. e MLC can not only describe the impact of intermediate principal stress on mechanical properties of surrounding rock but also have simple parameters. Its parameters S a and η can be calculated as a function of the well-known c and φ from the M-C criterion. In addition, based on the comparison result of various strength criteria for different types of rock limited by polyaxial experimental data, the MLC shows the better agreement with the test data [28,29]. erefore, in this study, a new numerical procedure based on the MLC is proposed to carry out the elastoplastic analysis of surrounding rock under nonaxisymmetric original stress filed.

Basic Assumption.
Due to the effects of geological structure and depth, it is unreasonable to assume that the in situ stress of tunnel surrounding rock is an axisymmetric stress field. erefore, the elastoplastic analysis for the circular tunnel is based on the following assumptions: (1) e initial vertical in situ stress is P 0 , and the horizontal in situ stress is λP 0 , where λ is the lateral pressure coefficient (2) e circular tunnel is deep-buried and infinite, and the plane strain assumption is satisfied (3) Rock mass shows initial elasticity, isotropy, and strain-softening behavior obeying MLC

Solutions in the Elastic
Zone. According to the basic assumption, the principle of stress superposition of structural mechanics can be employed to carry the mechanical analysis in the elastic zone [30,31]. erefore, the original nonaxisymmetric in situ stress field (Model A) can be decomposed into the sum of uniform in situ stress field (Model B) with λP 0 and vertical stress field (Model C) with (1-λ)P 0 ( Figure 1).

Stress Solutions.
e radial and tangential stress (σ Br and σ Bθ ) in the elastic zone of Model B with uniform stress λP 0 are given as where P i is the support pressure, λ is lateral pressure coefficient, R 0 is the radius of the circular tunnel, and r is radial distance from any point to the center of the tunnel. e radial, tangential, and shear stress (σ Cr , σ Cθ , and τ Crθ ) in the elastic zone of Model C with vertical stress (1-λ) P 0 are given as 2 Advances in Civil Engineering where θ represents the angle from the horizontal direction.
From equations (1) and (2), the total stress components in the elastic zone of the Model A can be obtained as

Displacement Solutions.
According to Hooke's law, the radial displacement (u) of any point in the surrounding rock can be expressed as where and E is Young's modulus of surrounding rock and v is Poisson's ratio.

Description of the MLC.
Due to no limitation of the experimental data and the assumption of failure envelope shape, the LC can reflect well the strength properties of most of the rocks tested under various stress conditions [32]. However, there is no direct relationship between the parameters from LC and c and φ from the most commonly used M-C criterion, which limits the application of LC at some extent. For the purpose of more applicable, the parameters of MLC are determined directly by c and φ from M-C criterion. erefore, the MLC not only can take Advances in Civil Engineering the effect of intermediate principal stress into account but also is easier to be applied in engineering, which is expressed as where I 1 ′ and I 3 ′ are the appropriate first and third stress tensor invariants, σ 1 , σ 2 , and σ 3 are the principal stresses, S a and η are material constants, which can be determined by the cohesion, c, and internal friction angle, φ, from the M-C criterion, and P p represents the pore pressure.
In order to take the dilatancy into account at the plastic stage, the relationship of the principal stress and dilation angle is employed as follows [33]: where ψ is the dilation angle. Substituting equation (9) into equation (7), the MLC can be reformulated as

Plastic Potential Function and Softening
Parameter. e plastic potential function can be written as where c p represents the deviatoric plastic strain and k(c p ) indicates the coefficient of dilation. According to the plastic flow law, the relation between radial and tangential plastic strain increments can be expressed as Considering the strength parameters as a function of deviatoric plastic strain c p , it is reasonable to assume the bilinear functions of c p can reflect the change of strength parameters as where ω indicates one of the c, φ, and ψ, ω p and ω r are the peak value and residual value of any one of the strength parameters (c, φ, and ψ), and c * p represents the critical deviatoric plastic strain.

Stress and Displacement Solution in the Plastic Zone.
e stress in the plastic zone is related to the properties of the tunnel and surrounding rock, not to the original in situ stress field. When the stress in the plastic zone of circular tunnel is in equilibrium, it should satisfy the axisymmetric static equilibrium equation as Based on equation (3), the stress components of the boundary between elastic and plastic zone (r � R p ) can be given as where R p is the radius of plastic area and σ Rp represents the radial stress of the elastic-plastic boundary. Substituting equation (16) into equation (10), there is

Advances in Civil Engineering
where According to equation (17), σ Rp can be obtained numerically by using appropriate algorithm such as Cardano formula. en, σ θb can be calculated by equation (16).
As can be seen from equation (3), the total stress components (σ Ar , σ Aθ , τ Aθr ) in the elastic zone is function of θ, which means that the stress components (σ rb , σ Rp , and σ θb ) of the boundary between elastic and plastic zone (r � R p ) and the radius of plastic area (R p ) are also function of θ and is not a constant. Actually, equations (16) and (17) can also explain this. erefore, in order to calculate the stress and displacement in the plastic zone, the plastic area can be divided into n (n ≥ 500) concentric annuli with the same radial stress under different θ and σ rb , as shown in Figure 2 [9]. e normalized radius, ρ i , is used to represent the inner radius of i th annulus and can be expressed as e normalized radii of the elastic-plastic boundary can be depicted by ρ 0 , and there is erefore, the radial stress increment of i th annulus can be given by By substituting equation (10) into equation (15), there is where Based on equation (15), the tangential stress σ θ(i) of i th annulus can be obtained by To the plane strain problem, the elastic strain increment of i th annulus can be written by Hooke's law as follows: e total strains in the plastic zone can be resolved into elastic and plastic strain, and there is while the compatibility equation must be satisfied as Based on equations (26) and (27), the compatibility equation can be rewritten as Combining with equations (13), (25), (24), and (28), Δε p θ(i) can be obtained by the approximating differential formation, and there is where Based on equation (26), the total strain of i th annulus in the plastic zone can be expressed as So, the radial displacement u (i) of i th annulus can be calculated from the normalized radial displacement, and there is Because R p is not a constant, different R p can be obtained based on different θ with equations (19)- (30). For example, when θ � 0°(marked as θ |θ � 0°) , R p can be Advances in Civil Engineering calculated by the proposed method and marked as R p|θ � 0°. In order to calculate R p |θ � 0°, the plastic area can be divided into n (n ≥ 500) concentric annuli with the same radial stress (Δσ r|θ�0°) , and then, R p |θ � 0°c an be obtained by equations (19)- (30), as shown in Figure 3(a). en, R p|θ � 1°, R p|θ � 2°, . . ., R p|θ � 90°c an be calculated with the same method. e plastic zones were plotted in Figures 3(b)-3(d) to make it easy to understand when θ � 20°, θ � 40°, and θ � 90°. e flowchart for the sequence of calculations is shown in Figure 4.

Verification under Axisymmetric Stress Field.
In order to verify the proposed numerical procedure based on MLC, σ 2 � σ 3 and ψ � −90°in equation (9) are assumed to ignore the effect of intermediate stress principal, and then, the results from this study with λ � 1 can be compared with the M-C criterion. e dataset was taken from [9] as input data: R 0 � 5.0 m, P 0 � 3 MPa, P i � 0 MPa, E � 10 GPa, v � 0.2, φ p � 30°, φ r � 26°, c p � 0.5 MPa, and c r � 0.2 MPa. Figure 5 presents the comparisons of stress and displacement from degenerated MLC and M-C criterion. It can be seen that radial stress, tangential stress, and radial displacement calculated from the proposed method are in good agreement with the published literature.

Verification under Nonaxisymmetric Stress Field.
Because the original in situ stress field is not ideal uniform stress field, it is significant to consider the effects of lateral pressure coefficient, λ. To verify the results under nonaxisymmetric stress field, the parameters used from [31] are R 0 � 2.965 m, P 0 � 10 MPa, 3 MPa, and λ � 1.0, 1.5, and 2.0, respectively. σ 2 � σ 3 and ψ � −90°are still used in equation (9) to ignore the effect of intermediate principal.
e results obtained by this study and published literatures are shown in Table 1. From Table 1, both the results from this study are slightly smaller than [31,34] under different λ. When λ � 1.5 and 2.0, the results of this study with θ � 90°are a litter bigger than [35]. In general, the results from proposed numerical procedure are basically consistent with the literatures.

Influence of Lateral Pressure Coefficient, λ.
e lateral pressure coefficient, λ, is used to represent the inhomogeneity of original in situ stress field. λ < 1 means the vertical stress is greater than the horizontal stress, while λ > 1 is on the contrary. λ � 1 indicates the original in situ stress is uniform stress field. In order to investigate the effect of λ on the radius of plastic zone, the datasets from Section 3.2 was taken as input data with λ � 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.5, and 3.0, respectively. e plastic zone and relationship between λ and R mp /R 0 are shown in Figure 6, where R mp represents the maximum plastic radius.
From Figure 6, when λ < 1, R mp /R 0 decreases gradually with the increase of λ. e plastic zone is mainly located on both sides of the circular tunnel with λ � 0.25. With increasing λ, the shape of plastic zone is sliding into an ellipse with the major axis in the horizontal direction. When λ � 1.0 with the uniform stress field, a circular plastic zone is developed. When λ > 1, R mp /R 0 presents positive logarithm relevant to λ. e elliptical plastic zone appears again, while the major axis turns into the vertical direction. e plastic zone develops mainly along the top and bottom when λ � 3.0. erefore, λ affects not only the size of plastic zone but also the development direction, which should be reasonably taken into account in design and construction of the tunnel.

Influence of Critical Deviatoric Plastic Strain
represents the start point of residual behavior of surrounding rock, which can indicate the softening degree to a certain extent. e surrounding rock enters the residual state quickly with no softening behavior when c * p is very small. To study the influence of c * p on the radius of the plastic zone, the datasets from Section 3.2 was taken as input data with λ � 1.0, and c * p � 0.005, 0.01, 0.05, 0.1, 0.2, and 0.3, respectively. Relationship between R p /R 0 and c * p is plotted in Figure 7. Figure 7 shows that, with the increase of c * p , R p /R 0 decreases by means of a negative power function. Especially, Elastic-plastic boundary Figure 2: Normalized plastic zone analysis model. 6 Advances in Civil Engineering R p /R 0 drops from 8.73 to 3.2 rapidly when c * p raises from 0.005 to 0.05. en, with the increase of c * p continuously, the rate of decrease gradually reduces. In this way, it is pretty significant to determine the correct c * p for the calculation of the plastic zone.

Influence of the Dilation Angle, ψ.
To study the influences of the dilation angle, ψ, on the plastic zone, the datasets from Section 3.2 was still taken as input data with λ � 1.0, and ψ � 5°, 10°, 15°, 20°, 25°, and 30°, respectively. Relationship between R p /R 0 , radial displacement, and different ψ is shown in Figure 8.
From Figure 8(a), with the increase of ψ, R p /R 0 is approximately linear rise, whereas the slope of regression equation is 0.004, which means the gradient of increase is very small. When ψ increases from 5°to 30°, R p /R 0 is just 4.9% more. In this way, ψ has no significant effect on the radius of plastic zone. Figure 8(b) plots the relationship between radial displacement and ψ, where the radial displacement increases exponentially with increasing ψ. R p /R 0 increases more than 3 times when ψ raises from 5°to 30°.
at means the dilation angle has an important role to affect the radial displacement, which are coincident with [3]. Hence, the correctness of the presented method is also verified to a certain extent.
<90 o ?   Advances in Civil Engineering  Advances in Civil Engineering 9

Conclusions
Based on the relationship between the principal stresses and ψ, the MLC was first introduced to consider the effects of the intermediate principal stress and dilatancy on the elastoplastic behaviors of surrounding rock under the nonaxisymmetric stress field. A new numerical procedure was then proposed and verified. At last, a series of parametric analyses were performed. e following conclusions can be obtained: (1) Due to the convenience of parameters, the presented method can be efficiently employed to comprehensively carry out the elastoplastic analysis for the circular tunnel under both axisymmetric and nonaxisymmetric stress field. (2) e λ influences not only the size of plastic zone but also the developed direction. When λ < 1, the plastic zone grows along the horizontal direction, while develops vertically when λ < 1. erefore, it is pretty significant to consider the inhomogeneity of the original in situ stress field. (3) It is important to determine accurately the value of c * p based on experiments, which has great influence on the plastic zone radius, especially when it is smaller. In addition, the effect of ψ on the plastic zone radius is much greater than that on radial displacement.

Data Availability
e data used to verify the proposed numerical procedure are on pages 594-597 of [9] (https://www.sciencedirect.com/ science/article/abs/pii/S0886779807001125).   10 Advances in Civil Engineering