# Elastic Stress and Plastic Zone Distributions around a Deeply Buried Tunnel under the Nonhydrostatic Pressure

Authors: Shuhong Liu; Jinlu Fan; Chenhao Wu; Yongquan Zhu
Journal: Geofluids
 (2022)

---

## Abstract

The elastic stress and the plastic zone are the important mechanical parameters to determine the tunnel support design. Based on Muskhelishvili’s complex variable function, the analytical solution for the elastic stress around a deeply buried noncircular tunnel under the nonhydrostatic pressure is firstly derived. The shape and size of the plastic zone of the surrounding rock mass are then determined by substituting the elastic stresses into the Drucker-Prager yield criterion. Finally, taking a horseshoe-shaped tunnel as an example, the analytical solutions of elastic stress distribution and plastic zone shape around the tunnel under different lateral pressure coefficients are in good agreement with ANSYS numerical solutions. The calculation results show that if the vertical in situ stress exceeds the critical value, with increasing lateral pressure coefficient, the hoop stress at the roof and floor of the tunnel increases significantly and the shape and size of the plastic zone change obviously.


---

## Body

## 1. Introduction

With the rapid development of modern roads, railways, and mining, the role of tunnels is becoming more and more important. The safety of tunnel construction has attracted attention in both academia and engineering circles. For the plane problem of a hole or a crack in an infinite body, the most effective method is using the complex variable theory. Closed-form solutions for simple shapes, such as an elliptical hole, a crack, and a square hole, have been discussed in detail. The exact solutions for cracks emanating from an elliptical hole were also studied [1–3]. This method has been successfully used to solve the tunnel problem. A deeply buried tunnel can be regarded as a hole in an infinite elastic body. For unlined tunnels, Exadaktylos et al. [4, 5] gave analytical solutions of stress and displacement distributions around the semicircular or notched circular tunnels and compared with numerical results of FLAC2D. Guan et al. [6] and Liu et al. [7, 8] studied horseshoe-shaped tunnels and compared analytical solutions with FLAC3D and ANSYS numerical results, respectively. For a lined straight wall arch tunnel, Kargar et al. [9, 10] and Lu et al. [11] provided analytical solutions for stress according to the different contact condition between lining and surrounding rock; Lu et al. [12] gave analytic solutions for stress and displacement considering support delay. Li and Chen [13] and Liu et al. [14] obtained the analytical solutions for stress and displacement of lined horseshoe-shaped tunnels in isotropic and orthotropic surrounding rock, respectively.The plastic zone around a tunnel is induced simultaneously when tunnels are excavated in the surrounding rock mass. The extension of the plastic zone may increase pressure of the lining and thus inevitably reduce the safety of the tunnel. Analytical solutions for the elastic–plastic stresses and the plastic zone around a circular tunnel under an axial in situ stress were investigated [15–17]. While the surrounding rock is often under unequal biaxial in situ stresses in actual geological conditions, analytical solutions for the plastic zone in the vicinity of a circular tunnel were also obtained. Behnam et al. [18], Deng et al. [19], and Guo et al. [20] gave the shape of the plastic zone for different lateral pressure coefficient. Guo et al. [21] and Dong et al. [22] verified the correctness of the plastic zone predicted by theoretical calculation through numerical simulation.Considering the pressure relief effect, Xu et al. [23] proposed a method of grooving slots in the roof and floor of a circular tunnel under the nonuniform stress; then, a circular tunnel actually becomes a circular tunnel with two cracks. Based on the different yield criterion, the shape and size of the crack tip plastic zone were analyzed [24–26]. Li et al. [27] obtained the numerical distribution of the plastic zone in the surrounding rock with a large-span rectangular roadway using FLAC3D software. Using the Griffith yield criterion, Shi and Bai [28] studied the plastic zone between a horseshoe-shaped tunnel and a concealed cave to ensure rock safety thickness. Guo et al. [29] obtained the explicit form of analytical solution caused by a shallow circular tunnel and analyzed the plastic zone around the tunnel with pile load combined with Mindlin’s solution. Zou et al. [30] considered the effect of the intermediate principal stress and interaction between the surrounding rock and support structure on the plastic zone of the surrounding rock around the circular tunnel. Shi et al. [31] and Zhou and Wu [32] calculated the plastic zones of circular tunnel surrounding rock by different strength criteria. Ma et al. [33, 34] determined the plastic zones around two circular tunnels. Hu et al. [35] investigated the influence of left and right tunnels on the plastic zone shape of the surrounding rock around the center circular tunnel.The above-described literature shows that many studies focus on analytical solutions for elastic stresses of the surrounding rock with deep-buried tunnels. Additionally, more results of the plastic zone around circular tunnels are presented. The analytical distribution of the plastic zone around noncircular tunnels in engineering is very few due to the complexity of mathematics.In this paper, based on Muskhelishvili’s complex variable function theory, the analytical solution for the elastic stress of a deeply buried noncircular tunnel under the nonhydrostatic pressure is derived. Subsequently, the elastic stresses are taken into the Drucker-Prager yield criterion; then, the plastic zone can be estimated in the surrounding rock mass. Taking a horseshoe-shaped tunnel as an example, the analytical solution of elastic stress distributions and plastic zone shapes for different lateral pressure coefficients is given. The finite element software ANSYS is used to verify the accuracy of the proposed analytical solutions.


## 2. Problem and Basic Equations

The schematic diagram of the noncircular tunnel is presented in Figure1. Assuming that the buried depth is great enough, the effect of gravity can be neglected. The surrounding rock mass is subjected to the nonhydrostatic pressure P and λP acting along Ox- and Oy-axes, respectively. λ is the lateral pressure coefficient.Figure 1
Schematic diagram of the noncircular tunnel.Muskhelishvili’s method has obvious advantages in solving plane problems involving holes with complex geometrical shapes. The mapping functionz=ωζ transforms the outer region of the tunnel in the physical z-plane with z=x+iy into the outer region of a unit circle in the mathematical ζ-plane with ζ=ρeiθ, where ρ and θ are polar coordinates in the ζ plane. The general form of the conformal mapping function can be expressed as
(1)z=ωζ=∑k=1nCkζ−k−2,where the coefficients Ck and m are known constants for the given tunnel shape [36]. In general, only a few terms are sufficient to satisfy the accuracy of mapping.For plane strain problems, stresses and displacements in the surrounding rock mass are determined by the following stress functionsφζ and ψζ:
(2)φζ=Bωζ+φ0ζ,(3)ψζ=B′+iC′ωζ+ψ0ζ,(4)φ0ζ=12πi∫σωσω′σ¯φ0′σ¯φ−ζdσ−12πi∫σf0σdσφ−ζ,(5)ψ0ζ=12πi∫σωσ¯ω′σφ0′σσ−ζdσ−12πi∫σf0σ¯dσφ−ζ,(6)f0σ=−2Bωσ−B′−iC′ωσ¯,where ζ=σ=eiθ denotes an arbitrary point of the unit circle. The constants B, B′+iC′, and B′−iC′ are relevant with far-field stresses P and λP:
(7)B=P41+λ,B′+iC′=B′−iC′=−P21−λ.Substituting equations (1) and (7) into equation (6), f0σ can be rewritten as
(8)f0σ=−P21+λ∑k=3nCkσ−k−2+P21−λC1σ−1−P21+λC1σ+C2+P21−λ∑k=2nCkσk−2.The stress componentsσx, σy, and τxy in Cartesian coordinates can be represented as
(9)σx+σy=4ReΦζ,σy−σx+2iτxy=2ω′ζωζ¯Φ′ζ+ω′ζΨζ,where Φζ=φ′ζ/ω′ζ, Ψζ=ψ′ζ/ω′ζ.The stress componentsσθ and σρ in curvilinear coordinates can be represented as
(10)σθ+σρ=4ReΦζ.When there is no surface force on the boundary of the tunnel,σρ=0, then, equation (10) is simplified as
(11)σθ=4ReΦζ.


## 3. Solution for Stress Functions

From equation (1), the following formula may be derived:
(12)ωσω′σ¯=∑k=1nCkσ−k−2−k−2∑k=1nCkσk−1=∑k=1n−2Lkσ−k+∑k=0∞Lk′σk,where the coefficients Lk1≤k≤n−2 are expressed by the known coefficients Ck1≤k≤n−2 of equation (1).
(13)Ln−k−1=Cn−k+1k−2Ln−2Ck+k−3Ln−3Ck−1+⋯+Ln−k+1C3C1,1≤k≤n−2.The analytic functionφ0ζ can be expressed as the following series form:
(14)φ0ζ=∑k=1∞akζ−k,where ak are unknown real coefficients obtained from the following boundary conditions.Combining equations (12) and (14) yields
(15)ωσω′σ¯φ0′σ¯=−∑k=1n−2Lkσ−k+∑k=0∞Lk′σk∑k=1∞kakσk+1=∑k=1n−4Skσ−k+∑k=0∞Sk′σk,where
(16)Sk=−Ln−2a1,k=n−4,−2Ln−2a2−Ln−3a1,k=n−5,−n−3−kLn−2an−3−k−⋯−Lk+2a1,2≤k≤n−6,−n−4Ln−2an−4−⋯−2L4a2−L3a1,k=1.The real coefficientsSk1≤k≤n−4 are expressed by the known coefficients Ck1≤k≤n−2 of equation (1) and the still unknown coefficients ak1≤k≤n−4 of equation (14). Lk′ and Sk′ are not given because they are not used in the following derivation.Substituting equations (8), (14), and (15) into (4) and calculating the Cauchy integral on the right-hand side of this equation yield
(17)∑k=1∞akζ−k=−∑k=1n−4Skζ−k−P21+λ∑k=3nCkζ−k−2+P21−λC1ζ−1.In equation (17), it can be seen that ak=0k≥n−1. Because the coefficients of the same negative exponents of ζ on both sides of equation (17) are equal, the following equations are obtained:
(18)ak=−P21+λCn,k=n−2,−P21+λCn−1,k=n−3,−Sk−P21+λCk+2,2≤k≤n−4,−S1−P21+λC3+P21−λC1,k=1.The coefficientsak1≤k≤n−2 solved by the equation (18) give the analytic function φ0ζ:
(19)φ0ζ=∑k=1n−2akζ−k.Substituting equations (1), (7), and (19) into (2), the first stress function φζ is given as follows:
(20)φζ=P41+λ∑k=1nCkζ−k−2+∑k=1n−2akζ−k.Next, another analytic functionψ0ζ will be derived. The conjugate of complex numbers in equations (8) and (15) may be written:
(21)f0σ¯=−P21+λ∑k=2nCkσk−2+P21−λC1σ+C2−P21+λC1σ−1+P21−λ∑k=3nCkσ−k−2,ωσ¯ω′σφ0′σ=∑k=1n−4Skσk+∑k=0∞Sk′σ−k.Similarly, substituting equation (21) into (5) and calculating the Cauchy integral on the right-hand side of this equation, ψ0ζ can be obtained:
(22)ψ0ζ=ω¯1/ζω′ζ∑k=1n−2kakζ−k−1+∑k=1n−4Skζk−P21+λC1ζ−1+P21−λ∑k=3nCkζ−k−2.Substituting equations (1), (7), and (22) into (3), then, the second stress function ψζ is given:
(23)ψζ=−P1−λ2∑k=1nCkζ−k−2+ω¯1/ζω′ζ∑k=1n−2kakζ−k−1+∑k=1n−4Skζk−P21+λC1ζ−1+P21−λ∑k=3nCkζ−k−2.


## 4. Plastic Zone

For the noncircular tunnel under different in situ stresses, the exact closed-form solution for the plastic zone in the surrounding rock mass cannot be obtained. A suitable yield criterion for the surrounding rock mass is introduced to estimate the plastic zone. In this section, the Drucker-Prager (D-P) yield criterion is adopted:(24)αI1+J2−k=0,where I1 is the first stress invariant and J2 is the second invariant of the deviator stress tensor. α and k are internal friction angle ϕ and adhesion C of the surrounding rock mass, respectively. The following formulas may be written:
(25)I1=σx+σy+σz,J2=16σx−σy2+σy−σz2+σz−σx2+6τxy2+τyz2+τxz2,α=2sinϕ33sinϕ,k=6Ccosϕ33−sinϕ.For plane strain problem,σz=μσx+σy,τyz=τxz=0, and μ represents the Poisson’s ratio. Substituting equations (9) and (25) into (24), the relationship between ρ and θ can be obtained; then, the shape and size of the plastic zone are determined from equation (1).


## 5. Results and Discussion

In order to verify the derived solution, analytical results of the elastic stress distributions and the plastic zone shapes around a horseshoe-shaped tunnel under different lateral pressure coefficients are compared with numerical results predicted by ANSYS finite element software.The tunnel is a two-lane highway tunnel located in grade V surrounding rock. The horseshoe-shaped cross-section consists of six successive arcs, as shown in Figure2. The center position and the radius of each arc are O1=0.00+0.00i, O2=1.73+1.73i, O3=−0.18+5.25i, and O4=12.80+0.00i and R1=8.33m, R2=5.88m, R3=1.88m, and R4=15.88m, respectively.Figure 2
Cross-section of the tunnel.The coefficients in equation (1) which was obtained by searching points on the boundary are as follows [36]: C1=6.7771, C2=2.3095, C3=−0.9785, C4=0.3119, C5=−0.1391, C6=0.0364, C7=0.0371, C8=−0.0374, C9=0.0181, and C10=−0.0017.The tunnel cross-section after mapping is adopted for analytical and numerical analyses.The two-dimensional finite element model of size400m×400m is much larger than the tunnel. Therefore, the numerical model can be regarded as the infinite surrounding rock mass with a horseshoe-shaped tunnel under far-field loads P and λP, just like the analytic model. The following parameters of the surrounding rock are assumed elastic modulus E=1.9GPa, Poisson’s ratio μ=0.4, adhesion C=0.1MPa, and internal friction angle ϕ=25°.


### 5.1. Elastic Stress Distribution

In this section, far-field stressP=−3MPa and the lateral pressure coefficient λ=2/3. Figure 3 displays the contours of elastic stresses σx, σy, and τxy obtained by the ANSYS finite element software. Figures 4–6 illustrate the distributions of elastic stresses σx, σy, and τxy along the tunnel boundary, the Ox-axis and the Oy-axis predicted by the proposed analytical solution and ANSYS finite element software. A positive angle α turns counterclockwise from the positive x-axis to the positive y-axis (Figure 1). In these figures, it can be seen that there is obvious stress concentration around the tunnel. Far away from the tunnel, the stresses σx, σy, and τxy tend to the constant values of −3 MPa, −2 MPa, and 0 MPa, which are the values of far-field in situ stresses. The excellent agreement of both predictions may be seen in Figures 4–6, because the same tunnel shape is adopted. These results indicate that the above-derived stress functions are accurate.Figure 3
Contours of stresses (in Pa) around the tunnel by ANSYS.
(a)
σx(b)
σy(c)
τxyFigure 4
Comparison of stresses along the tunnel boundary.Figure 5
Comparison of stresses along theOx-axis.Figure 6
Comparison of stresses along theOy-axis.Figure7 illustrates the distribution of hoop stress σθ along the tunnel boundary predicted by the proposed analytical solution for λ=0.2, 0.4, 0.6, 0.8, 1.0, and 1.2 for the case of P=−0.6MPa. An interesting phenomenon is found that there are two identical inflection points of the six curves. With increasing lateral pressure coefficient, the hoop stress at the points θ=0° and 180° changes significantly.Figure 7
Distribution of stressσθ along the tunnel boundary.


### 5.2. Plastic Zone

Figures8 and 9 illustrate the shape and size of the plastic zone around the tunnel predicted by ANSYS finite element software and the above-derived analytical solution for λ=0.2, 0.4, 0.6, 0.8, 1.0, and 1.2 for the case of P=−0.6MPa. The good agreement of both predictions may be seen from these figures. The plastic zone distribution is symmetrical to the vertical axis, and the boundary of the plastic zone is smooth.Figure 8
The plastic zone for differentλ predicted by ANSYS.
(a)
λ=0.2(b)
λ=0.4(c)
λ=0.6(d)
λ=0.8(e)
λ=1.0(f)
λ=1.2Figure 9
The plastic zone for differentλ predicted by the analytical solution.
(a)
λ=0.2(b)
λ=0.4(c)
λ=0.6(d)
λ=0.8(e)
λ=1.0(f)
λ=1.2In order to view the variation of the shape and size of the plastic zone change with the lateral pressure coefficient, the six conditions predicted by the analytical solution are plotted in Figure10. It can be seen that with the lateral pressure coefficient increases, the shape and size of the plastic zone change obviously and the location extends from the tunnel sides to the roof and floor in turn. However, the plastic zone around the tunnel sidewalls shows little change. These results indicate that the plastic zone estimation in Section 4 is accurate. In the study, we also find that if the vertical in situ stress is less than a certain critical value, such as P=−0.22MPa, the location of the plastic zone varies with the lateral pressure coefficient.Figure 10
The plastic zone predicted by the analytical solution forλ=0.2, λ=0.4, λ=0.6, λ=0.8, λ=1.0, and λ=1.2.


## 5.1. Elastic Stress Distribution

In this section, far-field stressP=−3MPa and the lateral pressure coefficient λ=2/3. Figure 3 displays the contours of elastic stresses σx, σy, and τxy obtained by the ANSYS finite element software. Figures 4–6 illustrate the distributions of elastic stresses σx, σy, and τxy along the tunnel boundary, the Ox-axis and the Oy-axis predicted by the proposed analytical solution and ANSYS finite element software. A positive angle α turns counterclockwise from the positive x-axis to the positive y-axis (Figure 1). In these figures, it can be seen that there is obvious stress concentration around the tunnel. Far away from the tunnel, the stresses σx, σy, and τxy tend to the constant values of −3 MPa, −2 MPa, and 0 MPa, which are the values of far-field in situ stresses. The excellent agreement of both predictions may be seen in Figures 4–6, because the same tunnel shape is adopted. These results indicate that the above-derived stress functions are accurate.Figure 3
Contours of stresses (in Pa) around the tunnel by ANSYS.
(a)
σx(b)
σy(c)
τxyFigure 4
Comparison of stresses along the tunnel boundary.Figure 5
Comparison of stresses along theOx-axis.Figure 6
Comparison of stresses along theOy-axis.Figure7 illustrates the distribution of hoop stress σθ along the tunnel boundary predicted by the proposed analytical solution for λ=0.2, 0.4, 0.6, 0.8, 1.0, and 1.2 for the case of P=−0.6MPa. An interesting phenomenon is found that there are two identical inflection points of the six curves. With increasing lateral pressure coefficient, the hoop stress at the points θ=0° and 180° changes significantly.Figure 7
Distribution of stressσθ along the tunnel boundary.


## 5.2. Plastic Zone

Figures8 and 9 illustrate the shape and size of the plastic zone around the tunnel predicted by ANSYS finite element software and the above-derived analytical solution for λ=0.2, 0.4, 0.6, 0.8, 1.0, and 1.2 for the case of P=−0.6MPa. The good agreement of both predictions may be seen from these figures. The plastic zone distribution is symmetrical to the vertical axis, and the boundary of the plastic zone is smooth.Figure 8
The plastic zone for differentλ predicted by ANSYS.
(a)
λ=0.2(b)
λ=0.4(c)
λ=0.6(d)
λ=0.8(e)
λ=1.0(f)
λ=1.2Figure 9
The plastic zone for differentλ predicted by the analytical solution.
(a)
λ=0.2(b)
λ=0.4(c)
λ=0.6(d)
λ=0.8(e)
λ=1.0(f)
λ=1.2In order to view the variation of the shape and size of the plastic zone change with the lateral pressure coefficient, the six conditions predicted by the analytical solution are plotted in Figure10. It can be seen that with the lateral pressure coefficient increases, the shape and size of the plastic zone change obviously and the location extends from the tunnel sides to the roof and floor in turn. However, the plastic zone around the tunnel sidewalls shows little change. These results indicate that the plastic zone estimation in Section 4 is accurate. In the study, we also find that if the vertical in situ stress is less than a certain critical value, such as P=−0.22MPa, the location of the plastic zone varies with the lateral pressure coefficient.Figure 10
The plastic zone predicted by the analytical solution forλ=0.2, λ=0.4, λ=0.6, λ=0.8, λ=1.0, and λ=1.2.


## 6. Conclusions

For the plane strain problem of the surrounding rock mass with a noncircular tunnel under the nonhydrostatic pressure, the analytical solution for the elastic stresses was obtained. Hence, the shape and size of the plastic zone can be estimated.(1)
The elastic stress distributions and the plastic zone shapes around the tunnel for different lateral pressure coefficients predicted by the proposed analytical model were in good agreement with those by the ANSYS numerical model(2)
There is obvious stress concentration around the tunnel, and the stress values tend to the far-field in situ stresses far away from the tunnel. With increasing lateral pressure coefficient, the hoop stress at the roof and floor of the tunnel changes significantly. For different lateral pressure coefficients, there are two identical inflection points on the distribution curves of hoop stress along the tunnel boundary(3)
If the vertical in situ stress exceeds the critical value, with increasing lateral pressure coefficient, the shape and size of the plastic zone change obviously and the location extends from the tunnel sides to the roof and floor in turn. However, the plastic zone around the tunnel sidewalls shows little change



---

*Source: 2363291-2022-07-04.xml*