# Elastoplastic Analysis for Circular Tunnel Based on Modified Lade Criterion considering Strain Softening and Dilatancy

Authors: Yansheng Deng; Jianxiu Wang; Baoping Zou
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

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 on modified 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. The 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. The plastic zone radius shows a negative power function change with increasing critical deviatoric plastic strain and increases slightly with increasing dilation angle, ψ.


---

## Body

## 1. 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–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 strain-softening 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–15]. Therefore, 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–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. The MLC can not only describe the impact of intermediate principal stress on mechanical properties of surrounding rock but also have simple parameters. Its parameters Sa 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]. Therefore, 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.


## 2. Problem Description

### 2.1. 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. Therefore, the elastoplastic analysis for the circular tunnel is based on the following assumptions:(1)
The initial vertical in situ stress isP0, and the horizontal in situ stress is λP0, where λ is the lateral pressure coefficient(2)
The 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


### 2.2. 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]. Therefore, the original nonaxisymmetric in situ stress field (Model A) can be decomposed into the sum of uniform in situ stress field (Model B) with λP0 and vertical stress field (Model C) with (1-λ)P0 (Figure 1).Figure 1
The mechanical model of the circular tunnel under nonaxisymmetric in situ stress field.


#### 2.2.1. Stress Solutions

The radial and tangential stress (σBr and σBθ) in the elastic zone of Model B with uniform stress λP0 are given as(1)σBr=λP01−PiR02r2,σBθ=λP01+PiR02r2,where Pi is the support pressure, λ is lateral pressure coefficient, R0 is the radius of the circular tunnel, and r is radial distance from any point to the center of the tunnel.The radial, tangential, and shear stress (σCr, σCθ, and τCrθ) in the elastic zone of Model C with vertical stress (1-λ)P0 are given as(2)σCr=12P01−λ1−R02r2−12P01−λ1+3R04r4−4R02r2cos2θ,σCθ=12P01−λ1+R02r2+12P01−λ1+3R04r4cos2θ,τCrθ=−12P01−λ1+2R02r2−3R04r4sin2θ,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(3)σAr=12P01+λ1−R02r2−12P01−λ1+3R04r4−4R02r2cos2θ+PiR02r2,σAθ=12P01+λ1+R02r2+12P01−λ1+3R04r4cos2θ−PiR02r2,τArθ=−12P01−λ1+2R02r2−3R04r4sin2θ.


#### 2.2.2. Displacement Solutions

According to Hooke’s law, the radial displacement (u) of any point in the surrounding rock can be expressed as(4)u=εr=1−v2ErΔσθ−v1−vΔσr,where(5)Δσθ=121+λP0R02r2+121−λP03R04r4cos2θ−PiR02r2,Δσr=−121+λP0R02r2−121−λP03R04r4−4R02r2cos2θ+PiR02r2,and E is Young’s modulus of surrounding rock and v is Poisson’s ratio.Substituting equation (5) into equation (4), there is(6)u=1+v2EP0r1+λR02r2+1−λ3R04r4cos2θ−2PiR02P0r2−v1−λ4R02r2cos2θ.


### 2.3. Solutions in the Plastic Zone

#### 2.3.1. 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. Therefore, the MLC not only can take the effect of intermediate principal stress into account but also is easier to be applied in engineering, which is expressed as(7)I1′3I3′=27+η,(8)I1′=σ1+Sa−Pp+σ2+Sa−Pp+σ3+Sa−Pp,I3′=σ1+Sa−Ppσ2+Sa−Ppσ3+Sa−Pp,η=4tan2φ9−7sinφ1−sinφ,Sa=ctanφ,where I1′ and I3′ are the appropriate first and third stress tensor invariants, σ1, σ2, and σ3 are the principal stresses, Sa and η are material constants, which can be determined by the cohesion, c, and internal friction angle, φ, from the M-C criterion, and Pp 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]:(9)σ2=1+sinψσ1+1−sinψσ32,where ψ is the dilation angle.Substituting equation (9) into equation (7), the MLC can be reformulated as(10)3+sinψ/2σ1+3−sinψ/2σ3+3Sa3σ1+Sa1+sinψσ1+1−sinψσ3/2+Saσ3+Sa=27+η.


#### 2.3.2. Plastic Potential Function and Softening Parameter

The plastic potential function can be written as(11)Gσθ,σr,γp=σθ−kγpσr,(12)kγp=1+sinψγp1−sinψγp,where γp represents the deviatoric plastic strain and k(γ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(13)dεrp=−kγpdεθp.Considering the strength parameters as a function of deviatoric plastic strainγp, it is reasonable to assume the bilinear functions of γp can reflect the change of strength parameters as(14)ωγp=ωp−ωp−ωrγpγp∗,0<γp<γp∗,ωr,γp≥γp∗,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 γp∗ represents the critical deviatoric plastic strain.


#### 2.3.3. Stress and Displacement Solution in the Plastic Zone

The 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(15)rdσrdr+σr−σθ=0.Based on equation (3), the stress components of the boundary between elastic and plastic zone (r = Rp) can be given as(16)σrb=σRp,σθb=P01+λ+2P01−λcos2θ−σRp,where Rp 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(17)−sin3ψ−BsinψσRp3+−3sin2ψ3Sa+3+sinψ2+3BAsinψ2+BA2+BSaσRp2+−3sinψ3Sa+A3+sinψ22+BASa+Sa2−A2sinψ−BASa−BA22σRp2+3Sa+A3+sinψ23−ABSaA1+sinψ+2Sa2+BSa2=0,where(18)A=P01+λ+2P01−λcos2θ,B=27+η.According to equation (17), σRp can be obtained numerically by using appropriate algorithm such as Cardano formula. Then, σθ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 = Rp) and the radius of plastic area (Rp) are also function of θ and is not a constant. Actually, equations (16) and (17) can also explain this. Therefore, 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].Figure 2
Normalized plastic zone analysis model.The normalized radius,ρi, is used to represent the inner radius of ith annulus and can be expressed as(19)ρi=riRp.The normalized radii of the elastic-plastic boundary can be depicted byρ0, and there is(20)ρ0=RpRp=1.Therefore, the radial stress increment ofith annulus can be given by(21)Δσr=Pi−σRPn.By substituting equation (10) into equation (15), there is(22)2+M+Nσr+1+Mrdσr/dr+3Sa3σr+rdσr/dr+SaMσr+Mrdσr/dr+Nσr+Saσr+Sa27−η=0,where(23)M=1+sinψ2,N=1−sinψ2.Based on equation (15), the tangential stress σθ(i) of ith annulus can be obtained by(24)σθi=σri−σri−1ρi−ρi−1ρi+σri.To the plane strain problem, the elastic strain increment ofith annulus can be written by Hooke’s law as follows:(25)ΔεrieΔεθie=1+νE1−ν−ν−ν1−νΔσriΔσθi.The total strains in the plastic zone can be resolved into elastic and plastic strain, and there is(26)εrεθ=εreεθe+εrpεθp,while the compatibility equation must be satisfied as(27)dεθdr+εθ−εrr=0.Based on equations (26) and (27), the compatibility equation can be rewritten as(28)dεθpdρ+εθp−εrpρ=−dεθedρ−εθe−εreρ.Combining with equations (13), (25), (24), and (28), Δεθip can be obtained by the approximating differential formation, and there is(29)Δεθip=−Δεθieρ¯i−Δρiεθie−εrie−Δρiεθi−1p−εri−1pρ¯i+1+kγiΔρi,where(30)ρ¯i=ρi+ρi−12,(31)γi=γi−1+Δεθip−Δεrip.Based on equation (26), the total strain of ith annulus in the plastic zone can be expressed as(32)εriεθi=εri−1εθi−1+ΔεrieΔεθie+ΔεripΔεθip.So, the radial displacementu(i) of ith annulus can be calculated from the normalized radial displacement, and there is(33)ui=εθiρiRp.BecauseRp is not a constant, different Rp can be obtained based on different θ with equations (19)–(30). For example, when θ = 0°(marked as θ|θ = 0°), Rp can be calculated by the proposed method and marked as Rp|θ = 0°. In order to calculate Rp|θ = 0°, the plastic area can be divided into n (n ≥ 500) concentric annuli with the same radial stress (Δσr|θ=0°), and then, Rp|θ = 0° can be obtained by equations (19)–(30), as shown in Figure 3(a). Then, Rp|θ = 1°, Rp|θ = 2°, …, Rp|θ = 90° can be calculated with the same method. The plastic zones were plotted in Figures 3(b)–3(d) to make it easy to understand when θ = 20°, θ = 40°, and θ = 90°. The flowchart for the sequence of calculations is shown in Figure 4.Figure 3
Rp with different θ: (a) θ = 0°,(b) θ = 20°, (c) θ = 40°, and (d) θ = 90°.
(a)(b)(c)(d)Figure 4
Flowchart for the sequence of calculations.


## 2.1. 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. Therefore, the elastoplastic analysis for the circular tunnel is based on the following assumptions:(1)
The initial vertical in situ stress isP0, and the horizontal in situ stress is λP0, where λ is the lateral pressure coefficient(2)
The 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


## 2.2. 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]. Therefore, the original nonaxisymmetric in situ stress field (Model A) can be decomposed into the sum of uniform in situ stress field (Model B) with λP0 and vertical stress field (Model C) with (1-λ)P0 (Figure 1).Figure 1
The mechanical model of the circular tunnel under nonaxisymmetric in situ stress field.


### 2.2.1. Stress Solutions

The radial and tangential stress (σBr and σBθ) in the elastic zone of Model B with uniform stress λP0 are given as(1)σBr=λP01−PiR02r2,σBθ=λP01+PiR02r2,where Pi is the support pressure, λ is lateral pressure coefficient, R0 is the radius of the circular tunnel, and r is radial distance from any point to the center of the tunnel.The radial, tangential, and shear stress (σCr, σCθ, and τCrθ) in the elastic zone of Model C with vertical stress (1-λ)P0 are given as(2)σCr=12P01−λ1−R02r2−12P01−λ1+3R04r4−4R02r2cos2θ,σCθ=12P01−λ1+R02r2+12P01−λ1+3R04r4cos2θ,τCrθ=−12P01−λ1+2R02r2−3R04r4sin2θ,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(3)σAr=12P01+λ1−R02r2−12P01−λ1+3R04r4−4R02r2cos2θ+PiR02r2,σAθ=12P01+λ1+R02r2+12P01−λ1+3R04r4cos2θ−PiR02r2,τArθ=−12P01−λ1+2R02r2−3R04r4sin2θ.


### 2.2.2. Displacement Solutions

According to Hooke’s law, the radial displacement (u) of any point in the surrounding rock can be expressed as(4)u=εr=1−v2ErΔσθ−v1−vΔσr,where(5)Δσθ=121+λP0R02r2+121−λP03R04r4cos2θ−PiR02r2,Δσr=−121+λP0R02r2−121−λP03R04r4−4R02r2cos2θ+PiR02r2,and E is Young’s modulus of surrounding rock and v is Poisson’s ratio.Substituting equation (5) into equation (4), there is(6)u=1+v2EP0r1+λR02r2+1−λ3R04r4cos2θ−2PiR02P0r2−v1−λ4R02r2cos2θ.


## 2.2.1. Stress Solutions

The radial and tangential stress (σBr and σBθ) in the elastic zone of Model B with uniform stress λP0 are given as(1)σBr=λP01−PiR02r2,σBθ=λP01+PiR02r2,where Pi is the support pressure, λ is lateral pressure coefficient, R0 is the radius of the circular tunnel, and r is radial distance from any point to the center of the tunnel.The radial, tangential, and shear stress (σCr, σCθ, and τCrθ) in the elastic zone of Model C with vertical stress (1-λ)P0 are given as(2)σCr=12P01−λ1−R02r2−12P01−λ1+3R04r4−4R02r2cos2θ,σCθ=12P01−λ1+R02r2+12P01−λ1+3R04r4cos2θ,τCrθ=−12P01−λ1+2R02r2−3R04r4sin2θ,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(3)σAr=12P01+λ1−R02r2−12P01−λ1+3R04r4−4R02r2cos2θ+PiR02r2,σAθ=12P01+λ1+R02r2+12P01−λ1+3R04r4cos2θ−PiR02r2,τArθ=−12P01−λ1+2R02r2−3R04r4sin2θ.


## 2.2.2. Displacement Solutions

According to Hooke’s law, the radial displacement (u) of any point in the surrounding rock can be expressed as(4)u=εr=1−v2ErΔσθ−v1−vΔσr,where(5)Δσθ=121+λP0R02r2+121−λP03R04r4cos2θ−PiR02r2,Δσr=−121+λP0R02r2−121−λP03R04r4−4R02r2cos2θ+PiR02r2,and E is Young’s modulus of surrounding rock and v is Poisson’s ratio.Substituting equation (5) into equation (4), there is(6)u=1+v2EP0r1+λR02r2+1−λ3R04r4cos2θ−2PiR02P0r2−v1−λ4R02r2cos2θ.


## 2.3. Solutions in the Plastic Zone

### 2.3.1. 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. Therefore, the MLC not only can take the effect of intermediate principal stress into account but also is easier to be applied in engineering, which is expressed as(7)I1′3I3′=27+η,(8)I1′=σ1+Sa−Pp+σ2+Sa−Pp+σ3+Sa−Pp,I3′=σ1+Sa−Ppσ2+Sa−Ppσ3+Sa−Pp,η=4tan2φ9−7sinφ1−sinφ,Sa=ctanφ,where I1′ and I3′ are the appropriate first and third stress tensor invariants, σ1, σ2, and σ3 are the principal stresses, Sa and η are material constants, which can be determined by the cohesion, c, and internal friction angle, φ, from the M-C criterion, and Pp 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]:(9)σ2=1+sinψσ1+1−sinψσ32,where ψ is the dilation angle.Substituting equation (9) into equation (7), the MLC can be reformulated as(10)3+sinψ/2σ1+3−sinψ/2σ3+3Sa3σ1+Sa1+sinψσ1+1−sinψσ3/2+Saσ3+Sa=27+η.


### 2.3.2. Plastic Potential Function and Softening Parameter

The plastic potential function can be written as(11)Gσθ,σr,γp=σθ−kγpσr,(12)kγp=1+sinψγp1−sinψγp,where γp represents the deviatoric plastic strain and k(γ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(13)dεrp=−kγpdεθp.Considering the strength parameters as a function of deviatoric plastic strainγp, it is reasonable to assume the bilinear functions of γp can reflect the change of strength parameters as(14)ωγp=ωp−ωp−ωrγpγp∗,0<γp<γp∗,ωr,γp≥γp∗,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 γp∗ represents the critical deviatoric plastic strain.


### 2.3.3. Stress and Displacement Solution in the Plastic Zone

The 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(15)rdσrdr+σr−σθ=0.Based on equation (3), the stress components of the boundary between elastic and plastic zone (r = Rp) can be given as(16)σrb=σRp,σθb=P01+λ+2P01−λcos2θ−σRp,where Rp 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(17)−sin3ψ−BsinψσRp3+−3sin2ψ3Sa+3+sinψ2+3BAsinψ2+BA2+BSaσRp2+−3sinψ3Sa+A3+sinψ22+BASa+Sa2−A2sinψ−BASa−BA22σRp2+3Sa+A3+sinψ23−ABSaA1+sinψ+2Sa2+BSa2=0,where(18)A=P01+λ+2P01−λcos2θ,B=27+η.According to equation (17), σRp can be obtained numerically by using appropriate algorithm such as Cardano formula. Then, σθ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 = Rp) and the radius of plastic area (Rp) are also function of θ and is not a constant. Actually, equations (16) and (17) can also explain this. Therefore, 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].Figure 2
Normalized plastic zone analysis model.The normalized radius,ρi, is used to represent the inner radius of ith annulus and can be expressed as(19)ρi=riRp.The normalized radii of the elastic-plastic boundary can be depicted byρ0, and there is(20)ρ0=RpRp=1.Therefore, the radial stress increment ofith annulus can be given by(21)Δσr=Pi−σRPn.By substituting equation (10) into equation (15), there is(22)2+M+Nσr+1+Mrdσr/dr+3Sa3σr+rdσr/dr+SaMσr+Mrdσr/dr+Nσr+Saσr+Sa27−η=0,where(23)M=1+sinψ2,N=1−sinψ2.Based on equation (15), the tangential stress σθ(i) of ith annulus can be obtained by(24)σθi=σri−σri−1ρi−ρi−1ρi+σri.To the plane strain problem, the elastic strain increment ofith annulus can be written by Hooke’s law as follows:(25)ΔεrieΔεθie=1+νE1−ν−ν−ν1−νΔσriΔσθi.The total strains in the plastic zone can be resolved into elastic and plastic strain, and there is(26)εrεθ=εreεθe+εrpεθp,while the compatibility equation must be satisfied as(27)dεθdr+εθ−εrr=0.Based on equations (26) and (27), the compatibility equation can be rewritten as(28)dεθpdρ+εθp−εrpρ=−dεθedρ−εθe−εreρ.Combining with equations (13), (25), (24), and (28), Δεθip can be obtained by the approximating differential formation, and there is(29)Δεθip=−Δεθieρ¯i−Δρiεθie−εrie−Δρiεθi−1p−εri−1pρ¯i+1+kγiΔρi,where(30)ρ¯i=ρi+ρi−12,(31)γi=γi−1+Δεθip−Δεrip.Based on equation (26), the total strain of ith annulus in the plastic zone can be expressed as(32)εriεθi=εri−1εθi−1+ΔεrieΔεθie+ΔεripΔεθip.So, the radial displacementu(i) of ith annulus can be calculated from the normalized radial displacement, and there is(33)ui=εθiρiRp.BecauseRp is not a constant, different Rp can be obtained based on different θ with equations (19)–(30). For example, when θ = 0°(marked as θ|θ = 0°), Rp can be calculated by the proposed method and marked as Rp|θ = 0°. In order to calculate Rp|θ = 0°, the plastic area can be divided into n (n ≥ 500) concentric annuli with the same radial stress (Δσr|θ=0°), and then, Rp|θ = 0° can be obtained by equations (19)–(30), as shown in Figure 3(a). Then, Rp|θ = 1°, Rp|θ = 2°, …, Rp|θ = 90° can be calculated with the same method. The plastic zones were plotted in Figures 3(b)–3(d) to make it easy to understand when θ = 20°, θ = 40°, and θ = 90°. The flowchart for the sequence of calculations is shown in Figure 4.Figure 3
Rp with different θ: (a) θ = 0°,(b) θ = 20°, (c) θ = 40°, and (d) θ = 90°.
(a)(b)(c)(d)Figure 4
Flowchart for the sequence of calculations.


## 2.3.1. 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. Therefore, the MLC not only can take the effect of intermediate principal stress into account but also is easier to be applied in engineering, which is expressed as(7)I1′3I3′=27+η,(8)I1′=σ1+Sa−Pp+σ2+Sa−Pp+σ3+Sa−Pp,I3′=σ1+Sa−Ppσ2+Sa−Ppσ3+Sa−Pp,η=4tan2φ9−7sinφ1−sinφ,Sa=ctanφ,where I1′ and I3′ are the appropriate first and third stress tensor invariants, σ1, σ2, and σ3 are the principal stresses, Sa and η are material constants, which can be determined by the cohesion, c, and internal friction angle, φ, from the M-C criterion, and Pp 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]:(9)σ2=1+sinψσ1+1−sinψσ32,where ψ is the dilation angle.Substituting equation (9) into equation (7), the MLC can be reformulated as(10)3+sinψ/2σ1+3−sinψ/2σ3+3Sa3σ1+Sa1+sinψσ1+1−sinψσ3/2+Saσ3+Sa=27+η.


## 2.3.2. Plastic Potential Function and Softening Parameter

The plastic potential function can be written as(11)Gσθ,σr,γp=σθ−kγpσr,(12)kγp=1+sinψγp1−sinψγp,where γp represents the deviatoric plastic strain and k(γ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(13)dεrp=−kγpdεθp.Considering the strength parameters as a function of deviatoric plastic strainγp, it is reasonable to assume the bilinear functions of γp can reflect the change of strength parameters as(14)ωγp=ωp−ωp−ωrγpγp∗,0<γp<γp∗,ωr,γp≥γp∗,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 γp∗ represents the critical deviatoric plastic strain.


## 2.3.3. Stress and Displacement Solution in the Plastic Zone

The 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(15)rdσrdr+σr−σθ=0.Based on equation (3), the stress components of the boundary between elastic and plastic zone (r = Rp) can be given as(16)σrb=σRp,σθb=P01+λ+2P01−λcos2θ−σRp,where Rp 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(17)−sin3ψ−BsinψσRp3+−3sin2ψ3Sa+3+sinψ2+3BAsinψ2+BA2+BSaσRp2+−3sinψ3Sa+A3+sinψ22+BASa+Sa2−A2sinψ−BASa−BA22σRp2+3Sa+A3+sinψ23−ABSaA1+sinψ+2Sa2+BSa2=0,where(18)A=P01+λ+2P01−λcos2θ,B=27+η.According to equation (17), σRp can be obtained numerically by using appropriate algorithm such as Cardano formula. Then, σθ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 = Rp) and the radius of plastic area (Rp) are also function of θ and is not a constant. Actually, equations (16) and (17) can also explain this. Therefore, 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].Figure 2
Normalized plastic zone analysis model.The normalized radius,ρi, is used to represent the inner radius of ith annulus and can be expressed as(19)ρi=riRp.The normalized radii of the elastic-plastic boundary can be depicted byρ0, and there is(20)ρ0=RpRp=1.Therefore, the radial stress increment ofith annulus can be given by(21)Δσr=Pi−σRPn.By substituting equation (10) into equation (15), there is(22)2+M+Nσr+1+Mrdσr/dr+3Sa3σr+rdσr/dr+SaMσr+Mrdσr/dr+Nσr+Saσr+Sa27−η=0,where(23)M=1+sinψ2,N=1−sinψ2.Based on equation (15), the tangential stress σθ(i) of ith annulus can be obtained by(24)σθi=σri−σri−1ρi−ρi−1ρi+σri.To the plane strain problem, the elastic strain increment ofith annulus can be written by Hooke’s law as follows:(25)ΔεrieΔεθie=1+νE1−ν−ν−ν1−νΔσriΔσθi.The total strains in the plastic zone can be resolved into elastic and plastic strain, and there is(26)εrεθ=εreεθe+εrpεθp,while the compatibility equation must be satisfied as(27)dεθdr+εθ−εrr=0.Based on equations (26) and (27), the compatibility equation can be rewritten as(28)dεθpdρ+εθp−εrpρ=−dεθedρ−εθe−εreρ.Combining with equations (13), (25), (24), and (28), Δεθip can be obtained by the approximating differential formation, and there is(29)Δεθip=−Δεθieρ¯i−Δρiεθie−εrie−Δρiεθi−1p−εri−1pρ¯i+1+kγiΔρi,where(30)ρ¯i=ρi+ρi−12,(31)γi=γi−1+Δεθip−Δεrip.Based on equation (26), the total strain of ith annulus in the plastic zone can be expressed as(32)εriεθi=εri−1εθi−1+ΔεrieΔεθie+ΔεripΔεθip.So, the radial displacementu(i) of ith annulus can be calculated from the normalized radial displacement, and there is(33)ui=εθiρiRp.BecauseRp is not a constant, different Rp can be obtained based on different θ with equations (19)–(30). For example, when θ = 0°(marked as θ|θ = 0°), Rp can be calculated by the proposed method and marked as Rp|θ = 0°. In order to calculate Rp|θ = 0°, the plastic area can be divided into n (n ≥ 500) concentric annuli with the same radial stress (Δσr|θ=0°), and then, Rp|θ = 0° can be obtained by equations (19)–(30), as shown in Figure 3(a). Then, Rp|θ = 1°, Rp|θ = 2°, …, Rp|θ = 90° can be calculated with the same method. The plastic zones were plotted in Figures 3(b)–3(d) to make it easy to understand when θ = 20°, θ = 40°, and θ = 90°. The flowchart for the sequence of calculations is shown in Figure 4.Figure 3
Rp with different θ: (a) θ = 0°,(b) θ = 20°, (c) θ = 40°, and (d) θ = 90°.
(a)(b)(c)(d)Figure 4
Flowchart for the sequence of calculations.


## 3. Verification Examples

### 3.1. 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. The dataset was taken from [9] as input data: R0 = 5.0 m, P0 = 3 MPa, Pi = 0 MPa, E = 10 GPa, v=0.2, φp = 30°, φr = 26°, cp = 0.5 MPa, and cr = 0.2 MPa.Figure5 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.Figure 5
Comparison between this study and literature results. (a) Radial and tangential stress. (b) Radial displacement.
(a)(b)


### 3.2. 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 R0 = 2.965 m, P0 = 10 MPa, Pi = 0 MPa, E = 20 GPa, v=0.3 , φp = 30°, φr = 30°, cp = 0.3 MPa, cr = 0.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.The results obtained by this study and published literatures are shown in Table1. 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.Table 1
Ratio results of the plastic zone radius and the tunnel radius (Rp/R0).
λ[33][34][35]This studyθ = 0°θ = 90°θ = 0°θ = 90°θ = 0°θ = 90°θ = 0°θ = 90°1.03.183.183.183.183.183.183.153.151.52.854.222.784.162.783.542.764.102.02.625.112.304.962.303.862.294.84


## 3.1. 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. The dataset was taken from [9] as input data: R0 = 5.0 m, P0 = 3 MPa, Pi = 0 MPa, E = 10 GPa, v=0.2, φp = 30°, φr = 26°, cp = 0.5 MPa, and cr = 0.2 MPa.Figure5 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.Figure 5
Comparison between this study and literature results. (a) Radial and tangential stress. (b) Radial displacement.
(a)(b)


## 3.2. 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 R0 = 2.965 m, P0 = 10 MPa, Pi = 0 MPa, E = 20 GPa, v=0.3 , φp = 30°, φr = 30°, cp = 0.3 MPa, cr = 0.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.The results obtained by this study and published literatures are shown in Table1. 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.Table 1
Ratio results of the plastic zone radius and the tunnel radius (Rp/R0).
λ[33][34][35]This studyθ = 0°θ = 90°θ = 0°θ = 90°θ = 0°θ = 90°θ = 0°θ = 90°1.03.183.183.183.183.183.183.153.151.52.854.222.784.162.783.542.764.102.02.625.112.304.962.303.862.294.84


## 4. Parametric Analysis and Discussion

### 4.1. Influence of Lateral Pressure Coefficient,λ

The 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. The plastic zone and relationship between λ and Rmp/R0 are shown in Figure 6, where Rmp represents the maximum plastic radius.Figure 6
The plastic zone and relationship betweenλ and Rmp/R0.From Figure6, when λ < 1, Rmp/R0 decreases gradually with the increase of λ. The 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, Rmp/R0 presents positive logarithm relevant to λ. The elliptical plastic zone appears again, while the major axis turns into the vertical direction. The plastic zone develops mainly along the top and bottom when λ = 3.0. Therefore, λ 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.


### 4.2. Influence of Critical Deviatoric Plastic Strain,γp∗

γp∗ represents the start point of residual behavior of surrounding rock, which can indicate the softening degree to a certain extent. The surrounding rock enters the residual state quickly with no softening behavior when γp∗ is very small. To study the influence of γp∗ on the radius of the plastic zone, the datasets from Section 3.2 was taken as input data with λ = 1.0, and γp∗=0.005,0.01,0.05,0.1,0.2,and0.3, respectively. Relationship between Rp/R0 and γp∗ is plotted in Figure 7.Figure 7
Relationship betweenRp/R0 and γp∗.Figure7 shows that, with the increase of γp∗, Rp/R0 decreases by means of a negative power function. Especially, Rp/R0 drops from 8.73 to 3.2 rapidly when γp∗ raises from 0.005 to 0.05. Then, with the increase of γp∗ continuously, the rate of decrease gradually reduces. In this way, it is pretty significant to determine the correct γp∗ for the calculation of the plastic zone.


### 4.3. 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 Rp/R0, radial displacement, and different ψ is shown in Figure 8.Figure 8
Relationship betweenRp/R0, radial displacement, and different ψ. (a) Rp/R0. (b) Radial displacement.
(a)(b)From Figure8(a), with the increase of ψ, Rp/R0 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°, Rp/R0 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 ψ. Rp/R0 increases more than 3 times when ψ raises from 5° to 30°. That 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.


## 4.1. Influence of Lateral Pressure Coefficient,λ

The 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. The plastic zone and relationship between λ and Rmp/R0 are shown in Figure 6, where Rmp represents the maximum plastic radius.Figure 6
The plastic zone and relationship betweenλ and Rmp/R0.From Figure6, when λ < 1, Rmp/R0 decreases gradually with the increase of λ. The 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, Rmp/R0 presents positive logarithm relevant to λ. The elliptical plastic zone appears again, while the major axis turns into the vertical direction. The plastic zone develops mainly along the top and bottom when λ = 3.0. Therefore, λ 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.


## 4.2. Influence of Critical Deviatoric Plastic Strain,γp∗

γp∗ represents the start point of residual behavior of surrounding rock, which can indicate the softening degree to a certain extent. The surrounding rock enters the residual state quickly with no softening behavior when γp∗ is very small. To study the influence of γp∗ on the radius of the plastic zone, the datasets from Section 3.2 was taken as input data with λ = 1.0, and γp∗=0.005,0.01,0.05,0.1,0.2,and0.3, respectively. Relationship between Rp/R0 and γp∗ is plotted in Figure 7.Figure 7
Relationship betweenRp/R0 and γp∗.Figure7 shows that, with the increase of γp∗, Rp/R0 decreases by means of a negative power function. Especially, Rp/R0 drops from 8.73 to 3.2 rapidly when γp∗ raises from 0.005 to 0.05. Then, with the increase of γp∗ continuously, the rate of decrease gradually reduces. In this way, it is pretty significant to determine the correct γp∗ for the calculation of the plastic zone.


## 4.3. 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 Rp/R0, radial displacement, and different ψ is shown in Figure 8.Figure 8
Relationship betweenRp/R0, radial displacement, and different ψ. (a) Rp/R0. (b) Radial displacement.
(a)(b)From Figure8(a), with the increase of ψ, Rp/R0 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°, Rp/R0 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 ψ. Rp/R0 increases more than 3 times when ψ raises from 5° to 30°. That 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.


## 5. 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. The 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)
Theλ 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. Therefore, 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γ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.



---

*Source: 4872558-2021-10-25.xml*