# Comparative Study on the Plastic Zone of Circular Hole Surrounding Rock in Anisotropic In Situ Stress Conditions

Authors: Pei Zhou; Peng Wu
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Anisotropic in situ ground stress is an objective stress state of rock mass that should be taken into account when calculating the plastic zone of the circular hole surrounding rock. The point stress method and approximate plastic condition method for calculating plastic zone are derived, and the finite element numerical calculation is carried out. Different analytical and finite element approaches are employed for calculation based on the representative parameters. The outcomes of the calculations are thoroughly examined. There are four shapes of plastic zone, among which the butterfly plastic zone can be obtained only by approximate plastic condition method and finite element method. The distribution of plastic zone calculated by using the point stress method and the Ruppneyt formula is identical, with only small deviations. The modified Fenner formula should not be used to calculate the plastic zone under in situ stress anisotropy since it would result in an underestimation of the maximum plastic radius. The detailed calculation and results comparison of this paper can provide a reference for more comprehensive and reasonable evaluation of the surrounding rock plastic zone.


---

## Body

## 1. Introduction

The in situ stress state of surrounding rock changes as a result of the excavation of circular hole (e.g., roadway and tunnel). When the redistributed stress meets the yield condition, the surrounding rock enters the plastic state [1, 2]. The distribution shape and size of the plastic zone serve as the foundation for evaluating the stability of the surrounding rock, as well as an essential basis for supporting design [3, 4].The in situ stress field was simplified to uniform distribution in early elastic-plastic analyses of surrounding rock, and the theory represented by the modified Fenner formula was proposed. Many scholars later optimized and improved this method in combination with the physical and mechanical properties of geotechnical materials, resulting in the development of a series of analytical calculation methods [5–8]. It is able to obtain a detailed analytical calculation theory of stress, strain, displacement, and plastic zone distribution. The plastic zone of circular hole surrounding rock is axisymmetric under isotropic stress conditions.The true in situ stress field, on the other hand, is generally anisotropic. Due to the stress redistribution effect, it is difficult to solve the exact mechanical field under the condition of anisotropic in situ stress after the surrounding rock enters the plastic state. Obviously, the distribution of the plastic zone is axi-asymmetric, and researchers have proposed some approximate solutions [9–13]. The shape and size of the plastic zone obtained by different methods differ due to different assumptions and mathematical solutions.In practice, the choice of calculation methods is currently subjective, and there is a lack of comparative analysis of different methods. Considering anisotropic in situ stress state, the point stress method and approximate plastic condition method for plastic zone calculation are derived. Based on the representative parameters, calculation is performed using four analytical methods: modified Fenner formula, Ruppneyt solution, point stress method, and approximate plastic condition method. Moreover, the nonlinear finite element software ABAQUS is used for comprehensive calculations in this paper.


## 2. Analytical Solution

The analytical solution is based on the reasonable mechanical model which is an idealized analysis model based on practical engineering problems that grasps the mechanical essence and major principles while making appropriate assumptions and simplifications. The basic assumptions for elastic-plastic analysis of circular hole surrounding rock are as follows: the cross section of hole is circular, the length of hole is much larger than its diameter, and the mechanical model is simplified to plane strain problem; geotechnical materials are homogeneous and isotropic; they are infinite boundary media; the gravity stress gradient is not considered, and in situ stress is regarded as the initial pressure acting on the boundary. The Mohr–Coulomb strength criterion is adopted as the condition for judging plastic yield. The basic mechanical model is shown in Figure1.Figure 1
Basic mechanical model.Kirsch was the first to investigate the elastic stress distribution of a tensile-stressed infinite plate, while other scholars updated and revised it, which later became well-known as the Kirsch equation [14–16]. Airy function trial methodology or the complex variable approach can be used to obtain Kirsch equation; detailed solution process is given in references [17, 18].The stress distribution of the hole surrounding rock in the elastic state is as follows based on Kirsch equation:(1)σr=12σv1+K01−Rr2−1−K01−4Rr2+3Rr4cos2θ+PiRr2,σθ=12σv1+K01+Rr2+1−K01+3Rr4cos2θ−PiRr,τrθ=12σv1−K01+2Rr2−3Rr4sin2θ.


### 2.1. Modified Fenner Formula

It is the elastic solution of isotropic stress condition when the lateral pressure coefficient in Kirsch equation is equal to 1. If the actual support force is less than the critical support force, the surrounding rock becomes plastic, and the radius of the plastic zone is [16](2)Rp=Rσv+ccotφ1−sinφPi+ccotφ1−sinφ/2sinφ.The modified Fenner formula, also known as the Kastner formula, is shown above. This approach is widely used and has a straightforward calculation. This assumption of in situ stress isotropy, on the other hand, is clearly incompatible with actual situations.


### 2.2. Ruppneyt Solution

Due to the stress redistribution effect, it is difficult to compute the exact mechanical field when in situ stress anisotropy exists. There are some approximate answers available at present, among which the Ruppneyt solution is a representative calculating approach [19]:(3)Rp=Rσv1+K0+2ccotφ1−sinφ2Pi+2ccotφ1−sinφ/2sinφ×1+σv1−K01−sinφcos2θ3σv1+K0+2ccotφsinφ.Ruppneyt equation considering anisotropic in situ stress conditions has an explicit expression and has been widely applied in engineering [20].


### 2.3. Point Stress Method

A plastic radius equation was proposed by Cai and Cai [21], and the results of the original literature are modified in this study.Assume that the circular hole surrounding rock is an axisymmetric plane strain solution after entering the plastic state, and the circumferential stress equation is as follows:(4)σθp=APi+ccotφrRA−1−ccotφ.The Kirsch equation is satisfied by the stress solution in the elastic region. The radial stress at the elastic-plastic interface is utilized to replace the support force, and the boundary radius of the plastic zone is used to replace the hole radius. The following equation can be obtained:(5)σθ=12σv1+K01+Rpr2+1−K01+3Rpr4cos2θ−σv1−sinφ−ccosφRpr2.The radius of the plastic zone is calculated using simultaneous (4) and (5):(6)Rp=R1−sinφ2σvcos2θ1−K0+σvK0+sinφ+ccosφ+cotφ1+sinφPi+ccotφ1−sinφ/2sinφ.This approach is called point stress method since it is based on the equal circumferential stress of individual points.


### 2.4. Approximate Plasticity Condition Method

Kastner proposed to ignore the stress redistribution caused by elastic-plastic deformation, i.e., the stress field of the hole surrounding rock is still the same as Kirsch equation during plastic deformation, and the boundary line equation of approximate plastic zone is obtained based on the geometric relationship of stress Mohr circle [22]. This method has recently been used in the stability analysis of roadway surrounding rock [23–28]. Taking the supporting force into account, the boundary line equation is rederived in this study.The following equation can be obtained after a series of operations such as the combination of similar terms and the transformation of trigonometric functions (the specific calculation and simplification process are omitted):(7)M8Rr8+M6Rr6+M4Rr4+M2Rr2+M0=0,in which(8)M8=9K0−12,M6=−62K0−12+cos2θK02−1+12Piσvcos2θK0−1,M4=K0−12cos4θ6−2A−12A+12+4−2A−12A+12+4cos2θK02−1+K0+12−4PiσvK0+1+2cos2θK0−1+4Piσv2,M2=−4cos4θK0−12−2cos2θK02−11−2A−12A+12+4cos2θK0−1Piσv+2BA−1A+121σv,M0=K0−12−K0+12A−12A+12−4BK0+1A−1A+121σv−4B2A+121σv2,A=1+sinφ1−sinφ,B=2ccosφ1−sinφ.The issue with this method is that the stress distribution in the plastic stage is identical to that in the elastic stage, making it a rough estimating method. Therefore, it is called approximate plastic condition method.


### 2.5. Calculation Results of Analytical Methods

Typical calculation parameters are adopted [29]: R = 2 m, γ = 25 kN/m3, H = 600 m, = 30°, and c = 1 MPa. The plastic zone boundary lines are depicted in Figures 2–4 using the Ruppneyt solution, point stress technique, and approximated plastic condition method, respectively.Figure 2
Boundary lines of plastic zone based on Ruppneyt solution.Figure 3
Boundary lines of plastic zone based on the point stress method.Figure 4
Boundary lines of plastic zone based on the approximate plastic condition method.


## 2.1. Modified Fenner Formula

It is the elastic solution of isotropic stress condition when the lateral pressure coefficient in Kirsch equation is equal to 1. If the actual support force is less than the critical support force, the surrounding rock becomes plastic, and the radius of the plastic zone is [16](2)Rp=Rσv+ccotφ1−sinφPi+ccotφ1−sinφ/2sinφ.The modified Fenner formula, also known as the Kastner formula, is shown above. This approach is widely used and has a straightforward calculation. This assumption of in situ stress isotropy, on the other hand, is clearly incompatible with actual situations.


## 2.2. Ruppneyt Solution

Due to the stress redistribution effect, it is difficult to compute the exact mechanical field when in situ stress anisotropy exists. There are some approximate answers available at present, among which the Ruppneyt solution is a representative calculating approach [19]:(3)Rp=Rσv1+K0+2ccotφ1−sinφ2Pi+2ccotφ1−sinφ/2sinφ×1+σv1−K01−sinφcos2θ3σv1+K0+2ccotφsinφ.Ruppneyt equation considering anisotropic in situ stress conditions has an explicit expression and has been widely applied in engineering [20].


## 2.3. Point Stress Method

A plastic radius equation was proposed by Cai and Cai [21], and the results of the original literature are modified in this study.Assume that the circular hole surrounding rock is an axisymmetric plane strain solution after entering the plastic state, and the circumferential stress equation is as follows:(4)σθp=APi+ccotφrRA−1−ccotφ.The Kirsch equation is satisfied by the stress solution in the elastic region. The radial stress at the elastic-plastic interface is utilized to replace the support force, and the boundary radius of the plastic zone is used to replace the hole radius. The following equation can be obtained:(5)σθ=12σv1+K01+Rpr2+1−K01+3Rpr4cos2θ−σv1−sinφ−ccosφRpr2.The radius of the plastic zone is calculated using simultaneous (4) and (5):(6)Rp=R1−sinφ2σvcos2θ1−K0+σvK0+sinφ+ccosφ+cotφ1+sinφPi+ccotφ1−sinφ/2sinφ.This approach is called point stress method since it is based on the equal circumferential stress of individual points.


## 2.4. Approximate Plasticity Condition Method

Kastner proposed to ignore the stress redistribution caused by elastic-plastic deformation, i.e., the stress field of the hole surrounding rock is still the same as Kirsch equation during plastic deformation, and the boundary line equation of approximate plastic zone is obtained based on the geometric relationship of stress Mohr circle [22]. This method has recently been used in the stability analysis of roadway surrounding rock [23–28]. Taking the supporting force into account, the boundary line equation is rederived in this study.The following equation can be obtained after a series of operations such as the combination of similar terms and the transformation of trigonometric functions (the specific calculation and simplification process are omitted):(7)M8Rr8+M6Rr6+M4Rr4+M2Rr2+M0=0,in which(8)M8=9K0−12,M6=−62K0−12+cos2θK02−1+12Piσvcos2θK0−1,M4=K0−12cos4θ6−2A−12A+12+4−2A−12A+12+4cos2θK02−1+K0+12−4PiσvK0+1+2cos2θK0−1+4Piσv2,M2=−4cos4θK0−12−2cos2θK02−11−2A−12A+12+4cos2θK0−1Piσv+2BA−1A+121σv,M0=K0−12−K0+12A−12A+12−4BK0+1A−1A+121σv−4B2A+121σv2,A=1+sinφ1−sinφ,B=2ccosφ1−sinφ.The issue with this method is that the stress distribution in the plastic stage is identical to that in the elastic stage, making it a rough estimating method. Therefore, it is called approximate plastic condition method.


## 2.5. Calculation Results of Analytical Methods

Typical calculation parameters are adopted [29]: R = 2 m, γ = 25 kN/m3, H = 600 m, = 30°, and c = 1 MPa. The plastic zone boundary lines are depicted in Figures 2–4 using the Ruppneyt solution, point stress technique, and approximated plastic condition method, respectively.Figure 2
Boundary lines of plastic zone based on Ruppneyt solution.Figure 3
Boundary lines of plastic zone based on the point stress method.Figure 4
Boundary lines of plastic zone based on the approximate plastic condition method.


## 3. Finite Element Method

Finite element software ABAQUS is used for numerical calculation in order to evaluate the accuracy of the analytical solutions derived by various ideas and methodologies. There are two stages of the numerical simulation: in situ stress balance and hole excavation [30].The calculation parameters are identical to those used in the previous analytical calculations. The model’s width and height are both 40 m, and the hole is located in the center. The horizontal displacement at the left and right sides of the model is zero, and the vertical displacement at the bottom is also zero. The gravity of the rock mass is taken into account in the numerical simulation. As a result, a pressure of 14.5 MPa is applied on the top boundary, and the pressure on the left and right sides is determined according to lateral pressure coefficient. The element type of surrounding rock in the model is CPE4 solid element. To achieve excellent precision, the model has a total of 201,550 elements. Eight groups of numerical simulation are carried out in this research, with the calculated plastic zone distribution illustrated in Figure5.Figure 5
Distribution of plastic zone based on the finite element method: (a)K0 = 0.3; (b) K0 = 0.4; (c) K0 = 0.5; (d) K0 = 0.6; (e) K0 = 0.7; (f) K0 = 0.8; (g) K0 = 0.9; (h) K0 = 1.
(a)(b)(c)(d)(e)(f)(g)(h)


## 4. Comparison and Analysis

The numerical results clearly show that the plastic zone has four shapes: butterfly, curved rectangle with concave horizontal direction and convex vertical direction, approximate ellipse, and circle. Specifically, the plastic zone is butterfly whenK0 = 0.3∼0.5; when K0 = 0.6, it is a curved rectangle with concave horizontal direction and convex vertical direction; when K0 = 0.7∼0.9, the plastic zone is approximate ellipse; and when K0 = 1, it is round.Among various analytic approaches, only the approximate plastic condition method can calculate the butterfly-shaped plastic zone. AtK0 = 0.3∼0.4, there is an obvious butterfly plastic zone which is not noticeable at K0 = 0.5. Ruppneyt solution and point stress method, although considering the anisotropy of in situ stress, cannot reflect the distribution characteristics of butterfly plastic zone.WhenK0 = 0.7∼0.9, the plastic zone obtained from the two analytical methods (Ruppneyt solution and point stress method) and numerical simulation are approximate elliptical. Specifically, when K0 = 0.7, the plastic radius calculated by Ruppneyt solution is 3.45∼4.72 m and the result of the point stress method is 2.99∼4.77 m; the numerical simulation solution is 3.24∼4.73 m. When K0 = 0.8, the plastic radii calculated by these three methods (Ruppneyt solution, point stress method, and finite element method) are 3.78∼4.61 m, 3.52 ∼ 4.65 m, and 3.67∼4.66 m, respectively. These three methods are 4.09∼4.50 m, 3.98∼4.52 m, and 4.03∼4.59 m, respectively, at K0 = 0.9. Obviously, the results of these three methods are very close.Although the approximate plastic condition method can also obtain the elliptical plastic zone whenK0 = 0.8 ∼ 0.9, the plastic radius is significantly smaller than the results of the previous three calculation methods. When K0 = 0.8, the plastic radius is 2.59 ∼2.73 m and it is 2.64∼2.70 m when K0 = 0.9.The distribution of the plastic zone calculated using the point stress method and the Ruppneyt solution is nearly identical, with only minor differences. Specifically, whenK0 = 0.3, the vertical direction of the inner boundary of circular hole has an elastic range (85°∼ 95°and 25°∼275°) according to Ruppneyt solution.WhenK0 = 0.3∼0.5, there are elastic zones in this position calculated by the point stress method, and the range is larger than results of Ruppneyt solution. When K0 = 0.3, the range is 59°∼121°and 239°∼301°. The range is 64°∼ 116°and 244°∼296°when K0 = 0.4; the range is 73°∼ 107°and 253°∼ 287°when K0 = 0.5.When the lateral pressure coefficient is equal to 1, the boundary line of the plastic zone is circular. The plastic radii calculated using the modified Fenner formula, Ruppneyt solution, point stress method, and finite element method are 4.39 m, 4.40 m, 2.68 m, and 4.32∼4.51 m, respectively. The plastic radius computed by other methods is quite close, with the exception of the approximate plastic condition method.When the anisotropy of in situ stress is strong, the maximum plastic radius calculated by numerical simulation is much larger than the result of modified Fenner formula which is based on the assumption of in situ stress isotropic.


## 5. Conclusions

In this research, the plastic zone of circular hole surrounding rock is calculated and compared using four analytical approaches (modified Fenner formula, Ruppneyt solution, point stress method, and approximate plasticity condition method) and finite element method. Some key conclusions are reached through in-depth comparison and analysis with various methods.(1)
Four shapes of plastic zone (butterfly, curved rectangle with concave horizontal direction and convex vertical direction, approximate ellipse, and circle) can be obtained by using the finite element method and approximate plastic condition method. Other analytic approaches cannot get butterfly shape. Therefore, the approximate plastic condition method has advantages in describing the shape of plastic zone.(2)
In the distribution and evolution of the plastic zone, the degree of in situ stress anisotropy is critical. The butterfly distribution appears when the horizontal stress differs significantly from the vertical stress (whenK0 < 0.5 or K0 > 2), which is the key shape of the plastic radius with a significant abrupt change. The modified Fenner formula will severely underestimate the plastic zone range at this time.(3)
When lateral pressure coefficient is close to 1, Ruppneyt solution, point stress method, and finite element method calculation results are similar. The shape of the plastic zone calculated by the approximate plastic condition method is similar to that calculated by other methods, but the plastic radius is significantly underestimated.(4)
In the supporting design, the anisotropic in situ stress conditions need to be carefully considered. While using an anchor, the length of anchor should be greater than the maximum plastic radius. When in situ stress anisotropy is considerable, the appropriate plastic zone calculation method and support scheme must be carefully chosen.



---

*Source: 2208528-2022-06-06.xml*