# Elastoplastic Analysis of Circular Tunnel Based on Drucker–Prager Criterion

Authors: Xiongfei Yang; Hong Yuan; Jiayu Wu; Shanqing Li
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

Based on the Drucker–Prager yield criterion, the theoretical solution of stratigraphic deformation in tunnel excavation process is deduced by the cavity expansion theory. In view of soil loosening around the tunnel caused by the tunnel excavation process, the internal friction angle of the surrounding soil is not a constant but a function of normal stress. The piecewise linearization of the nonlinear yield function is used to analyze the elastoplastic solution of the cylindrical hole shrinkage. A comparison is conducted with a plastic zone in which the internal friction angle of the soil remains unchanged. It can be concluded that the radial stress, the tangential stress, the radial strain, and the tangential strain around the inner wall calculated from the former are smaller.


---

## Body

## 1. Introduction

After subway tunnel excavation, the stress of surrounding soil will be redistributed. When the stress on the surrounding soil exceeds the elastic limit of rock mass and enters into the plastic state, elastoplastic analysis of the surrounding soil must be conducted to determine the stability of the subway tunnel and to provide a basis for the quantitative design of roadway support [1,2]. For deep tunnels, to analyze the surrounding soil stress and displacement generated by tunnel excavation and support, it is not necessary to consider the surface boundary conditions and the surrounding soil can be regarded as an infinite domain and can be solved [3].There are many examples of prediction of ground displacement caused by tunnel excavation based on the theory of cavity expansion. Mair and Taylor [4] used the Tresca criterion to predict the deformation of soil around the tunnel in clay and compared them with those displacements measured in Green Park and Regent Park in London. Yu and Rowe [5] proposed theoretical and semitheoretical solutions to the expansion of the borehole and the expansion of the sphere in the unloading process during drained and undrained conditions. Pinto and Whittle [6] predicted and expound the deformation of shallowly buried tunnel in the soft soil layer, based on the cavity expansion theory. Mo et al. [7] deduced the theoretical solution of the concentric regions of two different rock soils and also carried out the finite element simulation. For a long time, the elastoplastic analysis of the tunnel has been conducted by adopting the Mohr–Coulomb criterion [4, 8–10] or Hoek–Brown criterion [11–13], based on the fact that the internal friction angle of the surrounding soil is a constant. However, the internal friction angle is a function of normal stress rather than a constant. Kennedy and Lindberg [14] derived small strain solutions based on the segmented Mohr–Coulomb yield surface. Florence and Schwer [15] extended the linear Mohr–Coulomb criterion solution. Recently, Zou et al. [16, 17] focused on employing the elastoplastic theories of circular tunnel to solve geotechnical engineering problems and achieved outstanding results.In this paper, the Drucker–Prager criterion is used as the yield condition, and the plastic analysis of deep-buried circular subway tunnels is performed by using the piecewise linearization of nonlinear yield function. A comparison is conducted with a plastic zone in which the internal friction angle of the soil remains unchanged.


## 2. Drucker–Prager Yield Criterion

The Mohr–Coulomb failure criterion can be adopted to describe the yield and failure mechanisms of geomaterials. However, singularity will occur when the directional derivative is on the yield surface at the apex or ridgeline of the hexagonal cone [18, 19]. In order to overcome the shortcoming of the Mohr–Coulomb failure criterion, Drucker and Prager originally proposed the lower limit of the inscribed Mohr–Coulomb yield surface in 1952 [20], which can be expressed as(1)f=J2−αI1−Y=0,I1=σ1+σ2+σ3,J2=16σ1−σ22+σ2−σ32+σ3−σ12,where I1 is the first invariant of the Cauchy stress and J2 is the second invariant of the deviator stress tensor. In the formula, the material constants α and Y can be written as(2)α=sinφ33+sinφ,Y=3ccosφ3+sinφ,where φ and c are the internal friction angle and the cohesion of soil, respectively.Assuming the compressive stress as positive and the tensile stress as negative in the calculation, we can obtain the following relation on the basis of the cylindrical cavity expansion theory:(3)σ1=σθ,σ2=σr,σ3=0.Thus, the first invariant of the Cauchy stress and the second invariant of the deviator stress tensor can be simplified as(4)I1=σθ+σr,J2=σθ−σr26.By substituting Equation (4) into Equation (1), the Drucker–Prager yield criterion can be given by(5)16−ασθ−16+ασr=Y.By letting(6)λ=6+6α6−6α,γ=66−6α,and substituting Equation (6) into Equation (5), σθ can be transformed as a function of σr, which can be written as(7)σθ=γY+λσr.Based on the positional relationship between the Drucker–Prager criterion and the Mohr–Coulomb criterion on theπ plane (Figure 1), a variety of expressions for the parameters α and Y can be obtained. The parameters α and Y for various yield criteria are listed in Table 1.Figure 1
Yield surface projected onπ plane.Table 1
The parameters of the common Drucker–Prager yield criterion [21, 22].
No.
Criterion types
α
Y
DP1
The Mohr–Coulomb criterion based with external corner circumscribed, a circle yield criterion
2
sin
φ
/
3
3
−
sin
φ
6
c
cos
φ
/
3
3
−
sin
φ
DP2
The Mohr–Coulomb criterion based with inner corner circumscribed, a circle yield criterion
2
sin
φ
/
3
3
+
sin
φ
6
c
cos
φ
/
3
3
+
sin
φ
DP3
The Mohr–Coulomb criterion based with equivalent area, a circle yield criterion
2
3
sin
φ
/
2
3
π
9
−
sin
2
φ
6
3
c
cos
φ
/
2
3
π
9
−
sin
2
φ
DP4
The Mohr–Coulomb criterion based with matching circles, for plain strain problems with associated flow rules
sin
φ
/
3
3
+
sin
φ
k
=
3
c
cos
φ
/
3
+
sin
φ


## 3. Governing Equation

Before analyzing the elastic and plastic zones of the tunnel, the governing equations should be derived firstly. According to the theory of elasticity, the equilibrium equation of the infinitesimal body around the tunnel can be expressed as(8)rdσrdr+σr−σθ=0,with boundary conditions(9)σr=patr=a,σr=p0atr=∞.In addition, strain-displacement relations of the soil can be given by(10)εr=dudr,εθ=ur,(11)dεθdr+εθ−εrr=0.After the strain components and displacement are obtained, the stress-strain relation involved in plane strain problem can be given by(12)εr−εrp=εre=1E′σr−ν′σθ,(13)εθ−εθp=εθe=1E′σθ−ν′σr,where(14)E′=E1−ν2,ν′=ν1−ν.


## 4. Multisection Drucker–Prager Criterion Solution

In this paper, the elastoplastic closed-form solution of cylindrical hole contraction in the unrestricted Drucker–Prager medium is analyzed by piecewise linearization of nonlinear yield function. The pore radius is set asa, and the hydrostatic pressure p0 is uniformly applied to isotropic soil. Tunnel wall pressure is reduced very slowly, so the power effect can be ignored. As the internal pressure decreases, the inner wall of the hole will fail first. As the pressure of the inner wall of the hole decreases further, the plastic zone will expand outward until it reaches a certain radius c1. Eventually, the material in the annular area a≤r≤c1 is in the plastic state, and the other area c1≤r≤∞ is purely elastic.


### 4.1. Stress Distribution in the Elastic Zone

As the tunnel pressurep decreases from p0, the deformation of the soil around the tunnel is purely elastic at the beginning. According to the theory of elasticity, the elastic solutions of stress components can be expressed as(15)σr=p0−p0−par2,(16)σθ=p0+p0−par2.According to Equations (5) and (6), the yield equation can be obtained as follows(17)σθ−λσr=γY.With the continuous reduction of the pore pressure, the inner wall will trigger an initial yielding when pore pressure reaches the following condition:(18)py=2p0−γY1+λ.


### 4.2. Stress Distribution in the Plastic Zone

When the pore wall reaches the initial yielding, the plastic area will be formed near the inner wall in the regiona≤r≤c1 with decrease of the pressure p. As mentioned above, the plastic area is divided into m annular zone, as shown in Figure 2. For different plastic zones, there are different yield conditions represented by different stress functions; hence, each zone can be analyzed separately.Figure 2
Plastic area is divided into several zones.For the plastic zone 1, the yield criterion is given by(19)f=f1=σθ−λ1σr−γ1Y1=0.Substituting this yield criterion Equation (19) into the equilibrium equation (Equation (8)), we can obtain the governing equation for the radial stress, which can be expressed as(20)rdσrdr+1−λ1σr=γ1Y1.By solving the Euler equation, the general solution for stress components in the plastic zone 1 can be obtained:(21)σr=γ1Y11−λ1+Arλ1−1,(22)σθ=γ1Y11−λ1+Aλ1rλ1−1.Additionally, the stress components in the outer elastic zone can be expressed as(23)σr=p0−Br−2,σθ=p0+Br−2.The integral constantsA and B can be solved by using the boundary and continuous conditions. Finally, we can obtain(24)A=21−λ1p0−2γ1Y11−λ12c11−λ1,(25)B=γ1Y1−1−λ1p01+λ1c12.By substituting Equation (24) into Equation (21), the radial stress at the elastoplastic interface can be obtained(26)p1=2p0−γ1Y11+λ1.Similarly, the equilibrium equation for theith plastic zone can be written as(27)rdσrdr+1−6+6αi6−6αiσr=66αi+6Yi,where(28)λi=6+6αi6−6αi,γi=66−6αi.The general solution can be obtained as follows by solving Equation (27).(29)σr=γiYi1−λi+Arλi−1.According to the continuous condition(30)σr|r=ci=pi.The stress expressions for other plastic zones can be derived. For the zonei, stress components can be obtained as(31)σr=γiYi1−λi+γiYiλi−1+pirciλi−1,σθ=γiYi1−λi+λiγiYiλi−1+pirciλi−1.The boundary radiusc2,c3,…,cm−1,cm of each region can be determined by the radial stress continuous conditions at the junction of the zone i−1 and the zone i, which yield(32)cm=aγmYm+pλm−1γmYm+pmλm−11/1−λm,cm−1=cmγm−1Ym−1+pmλm−1−1γm−1Ym−1+pm−1λm−1−11/1−λm−1,cm−2=cm−1γm−2Ym−2+pm−1λm−2−1γm−2Ym−2+pm−2λm−2−11/1−λm−2,c1=c2γ1Y1+p2λ1−1γ1Y1+p1λ1−11/1−λ1.


### 4.3. Displacement Analysis

This paper mainly discusses the strain and displacement in the elastic and plastic zones. It is convenient to work out the strain expression in the elastic zone by solving the simultaneous Equations (12), (13), (15), and (16). However, the law of plastic flow is required for the calculation strain and displacement in the plastic zone, which will be analyzed by following the associated flow rule.Firstly, we calculate the strain and displacement distributions of the plastic zone 1 in the regionc2≤r≤c1, which is close to the outer elastic area. The relationship between two plastic strain components can be obtained by referring to the associated flow rule:(33)εrp+λ1εθp=0.Equation (33) can be alternatively written as(34)εr=εre−λ1εθp,(35)εθ=εθe+εθp.where the tangential plastic strain εθp is undetermined. When Equations (34) and (35) are substituted into Equation (11), the governing equation for εθp can be obtained as follows(36)rdεθpdr+1+λ1εθp=g1r,where(37)g1r=−λ12−1E′γ1Y1λ1−1+p1rc1γ1−1.The general solution can be derived as follows by solving Equation (36):(38)εθp=−λ12−12λ1E′γ1Y1λ1−1+p1rc1λ1−1+Br−λ1−1,where the integral constant B is determined by the elastoplastic boundary condition:(39)εθp|r=c1=0.Thus, the integral constantB can be solved:(40)B=λ12−12λ1E′γ1Y1λ1−1+p1c1λ1+1.Thus, the strain components and displacement in zone 1 are expressed as(41)εr=1E′1−ν′1−λ1γ1Y1+1−ν′λ1+λ12−12γ1Y1λ1−1+p1rc1λ1−1−λ12−12E′γ1Y1λ1−1+p1c1rλ1+1,εθ=1E′1−ν′1−λ1γ1Y1−λ12−12λ1−λ1+ν′γ1Y1λ1−1+p1rc1λ1−1+λ12−12λ1E′γ1Y1λ1−1+p1c1rλ1+1,u=rεθ.Similarly, strain and displacement expressions in other plastic zones can be obtained. The solution of the integral constants can be determined by displacement continuous conditions.


## 4.1. Stress Distribution in the Elastic Zone

As the tunnel pressurep decreases from p0, the deformation of the soil around the tunnel is purely elastic at the beginning. According to the theory of elasticity, the elastic solutions of stress components can be expressed as(15)σr=p0−p0−par2,(16)σθ=p0+p0−par2.According to Equations (5) and (6), the yield equation can be obtained as follows(17)σθ−λσr=γY.With the continuous reduction of the pore pressure, the inner wall will trigger an initial yielding when pore pressure reaches the following condition:(18)py=2p0−γY1+λ.


## 4.2. Stress Distribution in the Plastic Zone

When the pore wall reaches the initial yielding, the plastic area will be formed near the inner wall in the regiona≤r≤c1 with decrease of the pressure p. As mentioned above, the plastic area is divided into m annular zone, as shown in Figure 2. For different plastic zones, there are different yield conditions represented by different stress functions; hence, each zone can be analyzed separately.Figure 2
Plastic area is divided into several zones.For the plastic zone 1, the yield criterion is given by(19)f=f1=σθ−λ1σr−γ1Y1=0.Substituting this yield criterion Equation (19) into the equilibrium equation (Equation (8)), we can obtain the governing equation for the radial stress, which can be expressed as(20)rdσrdr+1−λ1σr=γ1Y1.By solving the Euler equation, the general solution for stress components in the plastic zone 1 can be obtained:(21)σr=γ1Y11−λ1+Arλ1−1,(22)σθ=γ1Y11−λ1+Aλ1rλ1−1.Additionally, the stress components in the outer elastic zone can be expressed as(23)σr=p0−Br−2,σθ=p0+Br−2.The integral constantsA and B can be solved by using the boundary and continuous conditions. Finally, we can obtain(24)A=21−λ1p0−2γ1Y11−λ12c11−λ1,(25)B=γ1Y1−1−λ1p01+λ1c12.By substituting Equation (24) into Equation (21), the radial stress at the elastoplastic interface can be obtained(26)p1=2p0−γ1Y11+λ1.Similarly, the equilibrium equation for theith plastic zone can be written as(27)rdσrdr+1−6+6αi6−6αiσr=66αi+6Yi,where(28)λi=6+6αi6−6αi,γi=66−6αi.The general solution can be obtained as follows by solving Equation (27).(29)σr=γiYi1−λi+Arλi−1.According to the continuous condition(30)σr|r=ci=pi.The stress expressions for other plastic zones can be derived. For the zonei, stress components can be obtained as(31)σr=γiYi1−λi+γiYiλi−1+pirciλi−1,σθ=γiYi1−λi+λiγiYiλi−1+pirciλi−1.The boundary radiusc2,c3,…,cm−1,cm of each region can be determined by the radial stress continuous conditions at the junction of the zone i−1 and the zone i, which yield(32)cm=aγmYm+pλm−1γmYm+pmλm−11/1−λm,cm−1=cmγm−1Ym−1+pmλm−1−1γm−1Ym−1+pm−1λm−1−11/1−λm−1,cm−2=cm−1γm−2Ym−2+pm−1λm−2−1γm−2Ym−2+pm−2λm−2−11/1−λm−2,c1=c2γ1Y1+p2λ1−1γ1Y1+p1λ1−11/1−λ1.


## 4.3. Displacement Analysis

This paper mainly discusses the strain and displacement in the elastic and plastic zones. It is convenient to work out the strain expression in the elastic zone by solving the simultaneous Equations (12), (13), (15), and (16). However, the law of plastic flow is required for the calculation strain and displacement in the plastic zone, which will be analyzed by following the associated flow rule.Firstly, we calculate the strain and displacement distributions of the plastic zone 1 in the regionc2≤r≤c1, which is close to the outer elastic area. The relationship between two plastic strain components can be obtained by referring to the associated flow rule:(33)εrp+λ1εθp=0.Equation (33) can be alternatively written as(34)εr=εre−λ1εθp,(35)εθ=εθe+εθp.where the tangential plastic strain εθp is undetermined. When Equations (34) and (35) are substituted into Equation (11), the governing equation for εθp can be obtained as follows(36)rdεθpdr+1+λ1εθp=g1r,where(37)g1r=−λ12−1E′γ1Y1λ1−1+p1rc1γ1−1.The general solution can be derived as follows by solving Equation (36):(38)εθp=−λ12−12λ1E′γ1Y1λ1−1+p1rc1λ1−1+Br−λ1−1,where the integral constant B is determined by the elastoplastic boundary condition:(39)εθp|r=c1=0.Thus, the integral constantB can be solved:(40)B=λ12−12λ1E′γ1Y1λ1−1+p1c1λ1+1.Thus, the strain components and displacement in zone 1 are expressed as(41)εr=1E′1−ν′1−λ1γ1Y1+1−ν′λ1+λ12−12γ1Y1λ1−1+p1rc1λ1−1−λ12−12E′γ1Y1λ1−1+p1c1rλ1+1,εθ=1E′1−ν′1−λ1γ1Y1−λ12−12λ1−λ1+ν′γ1Y1λ1−1+p1rc1λ1−1+λ12−12λ1E′γ1Y1λ1−1+p1c1rλ1+1,u=rεθ.Similarly, strain and displacement expressions in other plastic zones can be obtained. The solution of the integral constants can be determined by displacement continuous conditions.


## 5. Case Study

There is a tunnel with the radiusr=a, the cohesion of the surrounding soil C=5MPa, the internal friction angle of the surrounding soil φ=30∘, and the initial hydrostatic pressure p0=28.2MPa acting uniformly on the soil. After tunnel excavation, the pressure of the inner wall of the tunnel will decrease slowly. When the pressure p on the inner wall of the tunnel decreases to 18MPa, 12.63MPa, and 3.17MPa, respectively, the distributions of radial stress and tangential stress (dimensionless) around the tunnel are shown in Figure 3.Figure 3
Variation of radial stress and tangential stress of soil around the tunnel.It can be seen from Figure3 that the tunnel is in an elastic state when the pore pressure is greater than 12.63 MPa, and the radial stress increases with the increase of r. In contrast, the tangential stress has opposite tendency. Once the inner wall pressure p decreases to 12.63 MPa, plasticity commences around the inner wall of the tunnel, so part of surrounding soil enters the plastic state. The tangential stress in the plastic area increases with the increase of r, whereas the tangential stress in the elastic zone decreases with the increase of r. In addition, both the radial and tangential stresses far from the tunnel are almost equal to the hydrostatic pressure.In the plastic area, by substituting Equation (24) into Equation (21), the relationship between the pore pressure p and the plastic area radius c can be given by(42)p=γY1−λ+21−λp0−2γY1−λ2ca1−λ.When the tunnel wall pressurep decreases to 3.17 MPa, the radius of the plastic area varies with the internal friction angle φ and cohesion C of the soil, as shown in Figures 4 and 5, respectively. It can be noted from Figures 4 and 5 that the radius of the plastic area decreases with the increase of the inner friction angle φ and cohesion C of the surrounding soil. It indicates that it is essential to use varying inner friction angle φ and cohesion C of the surrounding soil to characterize the mechanical behavior of tunnel excavation.Figure 4
Relationship between the radius of plastic zone and internal friction angleφ.Figure 5
Relationship between the radius of plastic zone and cohesionC.Notably, it is an efficient way to calculate the elastoplastic solution of tunnel excavation by using piecewise linear yield function to represent the nonlinear yield criterion. Based on the influence of the normal stress on the internal friction angle of the soil, we divided the plastic area into seven zones. The yield criterion for each zone with various internal friction angles can be expressed as(43)φ1=30∘,f1=0when0.9p1≤σr<p1φ2=28.5∘,f2=0when0.8p1≤σr<0.9p1φ3=27∘,f3=0when0.7p1≤σr<0.8p1φ4=25.5∘,f4=0when0.6p1≤σr<0.7p1φ5=24∘,f5=0when0.5p1≤σr<0.6p1φ6=22.5∘,f6=0when0.4p1≤σr<0.5p1φ7=21∘,f7=0whenp≤σr<0.4p1.Comparison of stress components between one plastic area with internal friction angleφ=30∘ and seven separated plastic zones with various internal friction angles is shown in Figure 6, and the corresponding strain components’ comparison is shown in Figure 7.Figure 6
Stress distributions divided into a plastic zone and multiple plastic zones.Figure 7
Strain distributions divided into a plastic zone and multiple plastic zones.In practical engineering, the closer the distance from the tunnel is, the looser the soil becomes during tunnel excavation and the smaller the internal friction angle becomes. Figures6 and 7 depict that the radius of the plastic area increases when the plastic area is divided into seven zones. The stress and strain components around the inner wall of the tunnel decrease. However, the values of stress and strain components in the elastic area are very close.To have a better understanding of the stress and strain distributions around the tunnel, finite element software ANSYS is used for numerical simulation [23]. Two-dimensional element PLANE42 is used to simulate surrounding soil. Under symmetrical conditions, a-quarter mode with 1,600 elements and 1,683 nodes is constructed, as depicted in Figure 8. The radius of the surrounding soil is 100 m, and the radius a of the tunnel is equal to 3 m. Comparison between finite element results and analytical results is illustrated in Figures 9 and 10, respectively.Figure 8
Finite element model of a-quarter structure.Figure 9
Comparison of strain distribution between analytical solution and finite element solution.Figure 10
Comparison of stress distribution between analytical solution and finite element solution.According to Figures9 and 10, we consider tunnel excavation as the two-dimensional plane strain problem in the theoretical derivation, without considering the influence of strain εz along the tunnels. This leads to a larger analytical solution of tangential stress σθ in plastic area and a smaller radius in plastic area. Notably, by using the associated flow rule to analyze the plastic strain, the analytical solution of εr near the inner wall of the tunnel becomes larger. The good agreement in elastic area between finite element results and theoretical results validates the model developed in the current study. Furthermore, our model results are in good agreement with the numerical results proposed by Zhang et al. [24].


## 6. Conclusions

Taking into account loosening of the soil around the tunnel caused by the tunnel excavation, the internal friction angle of rock and soil is not a constant but a function of normal stress. The piecewise linearization of nonlinear yield function was used to analyze the elastoplastic solution of shrinkage of cylindrical pores based on the Drucker–Prager yield criterion. A comparison is conducted with a plastic zone in which the internal friction angle of the soil remains unchanged. Based on this study, the following conclusions can be drawn:(1)
The tunnel excavation process leads to decrease of the pressure on the inner wall of the tunnel and increase of the radius of the plastic area. The radius of the plastic area is related to the internal friction angle and cohesion of the soil in the tunnel.(2)
The radius of the plastic area decreases with increasing of the inner friction angle and cohesion of the surrounding soil.(3)
Compared with only one plastic zone (internal friction angleφ=30∘), when the plastic area is divided into seven zones, the radius of the plastic area increases and the stress and strain components around the inner wall of tunnel decrease, but the distributions of stress and strain in the elastic area are very close.



---

*Source: 5149789-2018-09-18.xml*