Calculation Method of Underground Cavern Loosening Pressure Based on Limit Analysis Using Upper-Bound Theory and Its Engineering Application

Abstract

From the point of view of the failure mechanism of the disturbed zone, this paper uses the limit analysis upper-bound theory to analyze the calculation formula of the loosening pressure, distinguish the difference between the vertical pressure and the horizontal pressure in the underground cavern, combine the loosening characteristics of the disturbed zone with the open-type disturbed zone and the annular disturbed zone, and construct the multirigidity slider translation and rotation failure mode to discuss the calculation method of surrounding rock loosening pressure of underground caverns in upper soft and hard rock stratum. The relevant calculation examples are given, and the application of the upper-bound theory of limit analysis is demonstrated in detail. Based on the actual engineering background, the calculation results of the calculation method of the loosening pressure of the cavity based on the upper-bound theory of the limit analysis are analyzed and compared for the different depths and different types of caverns. The difference, rationality, and applicability of the calculation results of this method are analyzed and discussed.

1. Introduction

When calculating the loose surrounding rock pressure of an underground cavern in a rock mass with joints, the influence of the distribution characteristics and mechanical properties of various structural surfaces in the rock mass on the shape and size of the loose zone around the underground cavern will be ignored, which will make the calculation result produce a large deviation. Therefore, for underground caverns in jointed rock masses, it has become an urgent problem to study the surrounding rock deformation and failure mechanisms that can reflect the actual characteristics of the rock mass and to develop a more reasonable method for calculating the loosening pressure of surrounding rock masses.

Using reasonable mechanical models, using analytical methods to solve the surrounding rock pressure of underground caverns has always been the focus of the researcher’s research. Up to now, the commonly used analytical calculation formulas are mostly based on the circular underground cavern model. These formulas include Fenner’s formula, modified Fenner’s formula, Casaer’s formula, and Caquot’s formula. The objects of calculation are mostly deep underground caverns. In the actual geotechnical engineering, when calculating the problems such as the surrounding rock pressure of the underground cavern, the most concern is not the specific size and distribution of the internal stress field and displacement field of the surrounding rock; usually, only the failure load or stability degree corresponding to the final plastic flow state can be calculated.

Fenner’s formula, modified Fenner’s formula, and Casaer’s formula all assume that the cavern is circular and in infinite homogeneous stratum and that rock and soil mass are regarded as an ideal elastic-plastic medium [1]. Since the above assumptions are inconsistent with the stress-strain characteristics of actual geomaterials, some scholars have deduced a computational model that can consider the strain softening of surrounding rock [2]. The Caquot formula is based on Mohr-Coulomb strength theory. For circular tunnels with uniform plasticity surrounding medium, the weight of rock and soil mass in the plastic zone is derived from the pressure of surrounding rock based on the theory of secondary stress elasticity analysis.

The limit analysis method was first created in the 1920s. In 1975, Chen [3] published the book Limiting Analysis and Soil Plasticity. The book systematically clarified the application of limit analysis theory in geotechnical engineering problems. Sloan et al. [4] combined the finite element method of numerical analysis algorithm with the limit analysis theory for the first time and effectively solved the optimization problem in the process of Solving Limit Analysis theory. Sloan et al.’s research on limit analysis upper-bound method and limit analysis lower bound method promoted the application of limit analysis method in geotechnical and rock mechanics. Fraldi and Guarracino [5–8] used the variational method to analyze the collapse loads of circular and rectangular caverns by introducing the modified Hoek-Brown (H-B) criterion and Greenberg’s optimization method pioneering. Atkinson and Potts [9], Davis et al. [10], and Subrin and Wong [11] have done a lot of work on the stability of soil tunnels. Fraldi used a variational approach, pioneering the introduction of the modified Hoek-Brown (H-B) criterion and Greenberg’s optimization method to analyze the collapse load of circular and rectangular caverns.

For underground engineering, in practical engineering, people are not most concerned about the specific size and distribution of internal stress field and displacement field, but about the stability of surrounding rock and the safety of supporting structure after excavation. Therefore, the maximum load acting on the supporting structure when the surrounding rock reaches the limit state is more concerned in the structural design. For such problems, the ultimate load and failure modes corresponding to the failure of reaching the critical state of the rock and soil can be obtained by the method of limit analysis in plastic mechanics [12]. On this basis, according to the law of conservation of work, the virtual work equation and the virtual power equation are established to solve the physical quantity obtained. Therefore, for the calculation of surrounding rock pressure in shallow tunnels, the use of the limit analysis method is an effective means.

2. Limit Analysis Upper-Bound Method for Solving Surrounding Rock Pressure

When applying the limit analysis upper-bound method, the following three basic assumptions are made:(1)The surface is assumed to be a horizontal plane. The surrounding rock of the tunnel is assumed to be an ideal elastoplastic body, and the structural plane obeys the surface contact slip constitutive model(2)The tunnel vault is subjected to the vertical uniform pressure q, and the tunnel sidewall is subjected to the average cloth pressure e. Due to the poor applicability of the optimized solution engineering involved in the unknown number, in order to reduce the difficulty of the solution and reduce the unknown number, suppose that q is in a proportional relationship with e, and the scale factor is k; that is, e = kq(3)We do not consider the influence of the construction process and the tunnel bottom drum

In this paper, the minimum buried depth corresponding to the circular loose zone is used as the boundary standard between deep and shallow burial of underground caverns.

2.1. Destruction Mode and Velocity Field

For the underground cavern with the loose zone as the open type, as shown in Figure 1, the tangential line is made at the top of the tunnel. The tunnel section is approximated by a rectangular DEFG. The boundary curve of the loose zone is y = ax^b, H is the tunnel depth, and T is the height of the tunnel. Point C is located on the centerline of the tunnel excavation face. The height of point NA from the point C to the ground surface is H₀. The line of GA_m is a horizontal line connecting CA and GA_m. The ∠ACA_m is divided into m parts, and the angle is β. ∠A_mGF is divided into n equal parts, the angle size is ε, and the intersections of the boundary curves of the two sides of each angle and the loose zone are A, A₁, A_2, ...A_m…A_m+n-1, respectively. Assume that the roof and sidewalls of the cave are uniformly stressed, respectively, q and e.

Figure 1 

The failure mode schematic of the ringent trumpet-shaped loosening zone.

According to the associated flow rule, the velocity vector direction on the velocity discontinuity line between the rigid sliders should be φ_J angle with the discontinuity line, and each velocity vector should satisfy the vector closure condition. The velocity field corresponding to the failure of the tunnel under the failure mode shown in Figure 2 is obtained. Because of the symmetry, the right side is taken. Within the range allowed by the error, AA₁…A_iA_i+1…A_m+n-1F on the curve are connected by a straight line, and the shape of each triangular rigid slider is determined by β, ε, and α_i. The speed of △ANC is V₀, the direction is vertical downward, the speed of △CDG is V_m+n+1, the direction is vertical downward, and the speed of the triangular block at the slip line is V₁…V_m+n from top to bottom. The relative speeds of the respective sliding blocks are V_0,1…V_m+n-1,m+n, V_m+n,m+n+1.

Figure 2 

The velocity field schematic diagram of the ringent trumpet-shaped loosening zone.

For the underground cavern with the loose zone as the ring type, a similar treatment method is adopted. As shown in Figure 3, the tangential line is made at the top of the tunnel hole. The tunnel section is approximated by a rectangular DEFG, and the boundary curve function of the loose zone is in the form of y = 1/(a + bx). H is the tunnel depth and T is the tunnel height. Point C is located on the centerline of the tunnel excavation face. The height from point C to the ground surface is H₀. The GA_m line is a horizontal straight line connecting CB and CA_m. The ∠BCA_m is divided into m parts, the angle is β, and the ∠A_mGF is equally divided into n equal parts, the angle size is ε, and the intersections of the boundary curves of the two sides of each corner and the loose zone are A, A₁, A_2, ...A_m…A_m+n-1, respectively. It is assumed that the vertical and horizontal forces received by the roof and the sidewall of the cavern are evenly distributed, respectively, q and e. Due to the aliquot angle, the boundary equation y = 1/(a + bx) of the loose region can be obtained. It can be known that the unknown quantity is only H₀. Therefore, as long as the calculation formula of the surrounding rock pressure represented by H₀ is obtained, the extremum point can be obtained by deriving H₀. In turn, the minimum value of the surrounding rock pressure calculated by the upper limit method can be obtained, and finally, it is necessary to verify whether the speed vector field satisfies the closing condition. Similarly, according to the associated flow rule, the velocity vector direction on the velocity discontinuity line between the rigid sliders should be φ_J angle with the discontinuity line, and each velocity vector should satisfy the vector closure condition. The velocity field corresponding to the failure of the tunnel under the failure mode shown in Figure 3 is obtained (Figure 4). Due to the symmetry, the right side is taken. Within the range allowed by the error, AA₁…A_iA_i+1…A_m+n-1F on the curve are connected by a straight line, and the shape of each triangular rigid slider is determined by β, ε, and a_i. The velocity of △BCA is V₀. If the value of m is large enough, the BA segment is approximately horizontal, and then, the angle between V₀ and the vertical direction is infinitely close to φ_J. This paper assumes the horizontal of the BA segment; the velocity of △CDG is V_m+n+1, and the direction is vertical. Downward, the velocity of the triangular block at the slip line is V₁…V_m+n from top to bottom, and the relative speed of each sliding block is V_0,1…V_m+n-1,m+n, V_m+n,m+n+1.

Figure 3 

The failure mode schematic of the arch-shaped loosening zone.

Figure 4 

The velocity field schematic diagram of the arch-shaped loosening zone.

2.2. Derivation Process of Loosening Pressure of Underground Cavern in Open-Type Loose Area

When the stratum above the vault is heterogeneous, the geometrical resolution of the boundary of the rupture zone is y = ax^b. The surrounding rock pressure is calculated by the limit theoretical upper limit method.

Taking the failure mode (i+1) as an example (the failure mode shown in Figure 5), when the intersection of the upper soft layer and the lower hard layer boundary and the fracture boundary curve is between A_i and A_i+1, the upper limit method in the limit analysis theory is used to calculate this. Surrounding rock pressure under normal conditions. In the failure mode (i+1), the bulk density of the upper softer formation is γ₁, the internal friction angle is φ_J1, the cohesive force is C_J1, the bulk density of the lower hard formation is γ₂, the internal friction angle is φ_J2, and the cohesion is C_J2. F is the origin of the coordinate, and a Cartesian coordinate system is established. The horizontal right is positive and the vertical upward is positive. The analytical formula of the FA curve is y = ax^b. The corresponding velocity field at this time is shown in Figure 6.

Figure 5 

The (i + 1)-th failure.

Figure 6 

The velocity field schematic diagram of the (i + 1)-th failure.

(1)Solving the Geometric Quantities in the Failure Mode Use as the origin of coordinates. Establish a Cartesian coordinate system with a positive horizontal to positive and positive vertical. The analytical formula of the FA curve is y = ax^b. Let | BC | = H₀, | BO | = h, and tunnel span W; then, From , , , and , , expression are The slope of and is ; then, The analytical expression is Similarly, the analytical formula of the CA_i line can be obtained: The CA₁ linear equation is connected with the FA equation y = ax^b of the slip line, and the coordinate A₁(x₁,y₁) of the A₁ point can be obtained. Similarly, A₂(x₂,y₂), A_i (x_i,y_i), A_i+1(x_i+1,y_i+1), and A_m (x_m,y_m) can be obtained. The equation of the upper soft layer and the lower hard layer boundary line OO_m-i+3 is y = T + H-h, which is in line with the linear equation of the CA_i+1 linear equation, and the coordinates of the O_m-i+2 point are O_m-i+2. The equation of the upper soft layer and the lower hard layer boundary line OO_m-i+3 is y = T + H-h, which is in conjunction with the FA equation y = ax^b of the slip line, and the coordinates of the O_m-i+3 point can be obtained as O_m-i+3, from O₃ coordinates: available from ,  Similarly, it is available that available from ,  By Sine theorem, it is available that By Sine theorem, it is available from , , that(2)Calculation of Velocity Field From Figure 3, the speed vector size of each slider is(3)Gravity power calculation The area of each rigid block is The gravity power of each sliding block of the damage mode (i + 3) is From the above, the damage mode shown by the damage mode ( + 3) is half of the gravity power.(4)Internal energy dissipation power Half of the energy dissipation in the damage mode ( + 3) is (5) Supporting reaction power The undetermined parameter k is the ratio of the horizontal support reaction force to the vertical support reaction force, k = e/q.(6)Calculation of surrounding rock pressure The virtual power equation:  top plate vertical support reaction force: In the above formula, q is a function of the variable H₀; that is,

It is known from the upper limit theorem that the vertical support reaction force of the top plate is the minimum value of q = f (H₀). When q = f (H₀) takes the minimum value, the following conditions should be met:

Through equation (22), the value of H₀ when q takes the extreme value can be obtained. It is also necessary to verify whether the obtained H₀ satisfies the constraint condition corresponding to the failure mode:

If the above constraint is satisfied, the value of H₀ when the extreme value is obtained is substituted into equation (21) to obtain the maximum value of the upper limit solution of the vertical support reaction force of the top plate.

Horizontal support reaction is

2.3. Detonation Pressure Detonation Process of Underground Cavern in Annular Loose Zone

In the failure mode shown in Figure 3, the BC spacing H₀ is unknown. Considering the layer stratification, the values of the external force power and the internal energy loss power are different depending on the thickness of the softer upper layer, due to the calculation of the upper limit method. Both are related to the value of the formation bulk density, the value of the internal friction angle, and the value of the cohesion force, so the final calculation formula of the surrounding rock pressure must be different.

Taking the failure mode (i + 1) as an example (Figure 7), the intersection of the upper softer layer and the lower harder layer boundary line and the fracture boundary curve is located between A_i and A_i+1. Let the upper softer formation have a bulk density of γ₁, an internal friction angle of φ_J1, a cohesive force of C_J1, a lower hard formation density of γ₂, an internal friction angle of φ_J2, and a cohesive force of C_J2. F is the origin of the coordinate, and a Cartesian coordinate system is established. The horizontal right is positive and the vertical upward is positive. The FA curve analytical expression is y = 1/(a + bx). The corresponding velocity field at this time is shown in Figure 8.

Figure 7 

The (i + 1)-th failure mode.

Figure 8 

The velocity field schematic diagram of the (i + 1)-th failure mode.

(1)Figure 7 shows the solution of each geometric quantity in the loose failure mode. As shown in Figure 7, the vertex of the boundary curve of the loose zone is taken as the coordinate origin, and a Cartesian coordinate system is established, with the horizontal rightward being positive and the vertical upward being positive. The analytical formula of the BF curve is y = 1/(a + bx), the points B and R are located on the central axis of the tunnel, the loose failure mode is the axis of symmetry of the BR, and the right part of the BR is analyzed. From | BO | = h, | CR | = V₁, | N_mA_m| = L, and tunnel span W, then From and , the following can be obtained: From θ₂ = π/2-θ₁-θ₃, the following can be obtained: The coordinates of point C are C[0, (H-H₀-V₁)], and the slope of the line CA is cotβ; then, the straight-line equation of CA is y = x cot β + (H-H₀-V₁). Similarly, the linear equations of CA₁, CA₂, CA_i, and CA_i+1 are y = x cot(2β) + (H-H₀-V₁), y = x cot(3β)+(H-H₀-V₁), y = xcot [(i + 1)β] + (H-H₀-V₁), y = x cot[(i + 2)β] + (H-H₀-V₁), and y = x cot β + (H-H₀-V₁). Let(V₁+H₀-H) =  ; the linear equation of the simultaneous CA and the boundary equation of the loose zone can get the coordinates of the point. In the formula, ξ(β) = (a cot ß + b)² + 4b cot ß. In the same way, the coordinates of A₁ and A_i can be obtained: In the formula, The straight-line equation of the interface between the upper softer stratum and the lower harder stratum is y = H-V₁-h =  , and the coordinates of the O₁ point can be obtained in conjunction with the CA straight line equation, O₁[( + )/cosβ, ]. Similarly, y =  and the linear equations of CA₁, CA₂, and CA_i can be used to obtain the coordinates of O₂, O₃, and O_i+1 points: O₂[( + )/cos 2β, ], O₃[( + )/cos3β, ], and O_i+1[( + )/cos (i + 1)β. y =  is associated with the boundary curve equation of the loose zone and the coordinates of the O_i+2 point: O_i+2 [(1−ab)/b, ]. The distance between the points can be obtained from the coordinates of the above points. The coordinate of point G is G(W/2, −V₁), the slope of the line of GA_m+1 is -tanε, and the equation of the line of GA_m+1 is y = y = (−tanε) (x−W/2)−V₁. Similarly, the linear equations for GA_m+2, GA_m+j, and GA_m+n-1 are y = [−tan(2ε)](x−W/2)−V₁, y = [−tan(jε)](x−W/2)-V₁, and y = {−tan[(m + n−1)ε)](x−W/2)−V_1. Let bW−2a =  ; the line equation of the simultaneous GA_m+1 and the boundary curve equation of the loose zone can get the coordinates of the A_m+1 point: In the formula, η(ε) = b²V₁² + (a² + b²W²/4 + abW)tan²ε−(b²V₁² W + 2abV₁ + 4b)tan ε. Similarly, the coordinates of A_m+2, A_m+j, and A_m+n-1 can be obtained: In the formula, η(2ε), η(jε), and η[(m + n−1)ε] are The distance between AA₁ can be obtained from A(x₀, y₀) and A₁(x₁, y₁): The same is available: Using the sine theorem, we can find α, α1…αm + n:(2)Speed field calculation From Figure 3, the sine theorem can be used to find the velocity vector size of each slider:(3)Gravity power calculation The area of each rigid block is The weight of each sliding block of the damage mode (i + 1) is Half of the damage mode gravity power shown by the upper damage mode (i + 1) is(4)Internal energy dissipation power Half of the energy dissipation power in the failure mode (i + 1) is(5)Support reaction power The undetermined parameter k is the ratio of the horizontal support reaction force to the vertical support reaction force, k = e/q. The sum of the powers of the vertical support reaction force q and the horizontal support reaction force e of the tunnel roof is PF.(6)Calculation of surrounding rock pressure Virtual power equation: Roof vertical support reaction: In the above formula, q is a function of the variable H0, which is It is known from the upper limit theorem that the vertical support reaction force of the top plate is the minimum value of q = f (H₀). When q = f (H₀) takes the minimum value, the following conditions should be met: Through equation (47), the value of H₀ when q takes the extreme value can be obtained. It is also necessary to verify whether the obtained H₀ satisfies the constraint condition corresponding to the failure mode: If the above constraint condition is satisfied, the value of H₀ when the extreme value is obtained is substituted into equation (45) to obtain the maximum value of the upper limit solution of the vertical support reaction force of the top plate. Horizontal support reaction:

3. Case Analysis and Engineering Application

Based on two specific engineering cases, this paper uses 6 classic solutions before and after correction and the upper limit method of limit analysis based on the boundary of the loose zone for the shallow tunnel [13, 14] combined with the field measured data. A total of 7 methods are used to calculate the vertical pressure and horizontal pressure of the tunnel. The method is used to calculate the vertical pressure and horizontal pressure of the tunnel. For the deep-buried tunnel combined with the field measured data, the Platts pressure arch theory method, the standard method, the limit analysis upper limit method based on the loose zone boundary, and the Fraldi limit analysis method are used, respectively. The method calculates the vertical pressure and horizontal pressure of the tunnel.

3.1. Engineering Application Background

3.1.1. Qingdao Jiangxi Road Metro Station

The geology of Qingdao city is a typical upper soft and hard rock formation. This section uses the Jiangxi Road metro station as an example to calculate the surrounding rock pressure and safety factor of a single-arched large-span metro station excavated in such a stratum.

Jiangxi Road metro station [15] is located in Qingdao city, and the station is constructed by shallow tunneling method, with a width of 20.6 m and a height of 15.5 m. The ground level of the metro station is about 11.71 m. The overlying strongly weathered 1d granite has a calculated depth of 5.5 m, and the harder stratum is the late Yanshanian granite. The mid-weathered rock face has a large undulation, and the calculated depth is 6.5 m. The tunnel depth is 12 m. The values of physical and mechanical parameters of the softer stratum in the upper part are taken according to the parameters of the surrounding rock of Grade IV. The values of the physical and mechanical parameters of the lower hard stratum are taken according to the parameters of the surrounding rock of the V-class, and the pressure coefficients of the stratum are 0.56 and 0.58, respectively. Using C45 concrete to form the second lining, the thickness of the top arch is 0.65 m, the thickness of the inverted arch is 0.70 m, and the thickness of the sidewall is 0.70 m. The concrete has an elastic modulus of 32.5 GPa, a Poisson’s ratio of 0.2, and a bulk density of 25 kN/m³. The physical and mechanical parameters of the rock mass are shown in Table 1.

Table 1 

Physical and mechanical parameters of the rock mass of Jiangxi Road metro station.

3.1.2. Zhifang Tunnel of Lanzhou-Chongqing Railway

Zhifang tunnel of Lanzhou-Chongqing Railway [16] is a double-line tunnel in a carbonaceous slate. The tunnel has a total length of 5135 m, a minimum buried depth of 45 m, and a maximum buried depth of 248 m. In this paper, two sections of DK202 + 390 ∼ DK202 + 540 are taken, the buried depth is 65 m and 70 m, respectively, and the thickness of the softer upper layer is 25 m. The tunnel section height is 12.4 6 m, and the maximum span is 14.06 m. When calculating, the horizontal joint spacing is 0.6 m, the inclined transfixion joint spacing is 0.6 m, and the inclination angle is 60°. The values of the physical and mechanical parameters of the softer stratum in the upper part are taken according to the parameters of the surrounding rock of Grade IV, and the values of the physical and mechanical parameters of the lower hard stratum are taken according to the parameters of the surrounding rock of the V-class. The second lining adopts C35 waterproof concrete, the thickness of the top arch is 0.70 m, the thickness of the inverted arch is 0.80 m, the thickness of the sidewall is 0.70 m, the elastic modulus of the concrete is 30 GPa, Poisson’s ratio is 0.2, and the bulk density is 25 kN/m³. The physical and mechanical parameters of the rock mass are shown in Table 2.

Table 2 

Physical and mechanical parameters of the rock mass of Zhifang tunnel.

3.2. Loose Pressure Calculation Results and Difference Analysis

According to the standard of deep and shallow burial, the Qingdao Metro Station is a shallow tunnel. According to the structural dimensions of the cavern and the physical and mechanical parameters of the stratum, the calculation results are shown in Table 3.

Table 3 

The relaxation pressures of Jiangxi Road metro station.

It is known from Table 3 that when k = 0.8–1.0, for the vertical pressure value, the calculation result of the upper-bound theorem of limit analysis is similar to the weighted Terzaghi formula, the weighted Bierbäumer formula, and the weighted Jiaxiao Xie formula. When k = 1.2–1.4, for the vertical pressure value, the calculation result of the upper-bound theorem of limit analysis is similar to the field measured data. Since k takes different values, the calculation result of the upper-bound theorem of limit analysis can be consistent with the calculation results of the classical calculation method and can be consistent with the field measured data, so the limit analysis method is available and effective.

Compared with the vertical pressure of the tunnel, based on various factors including the joint production, the calculated geometrical quantities of each feature such as the sliding angle of the loose area and the length of the loose rock of the dome are corrected. The calculation result of the classical method is larger than the calculation result of the classical method before the correction.

According to the standard of deep and shallow burial, the two-buried-depth (65 m, 70 m) Zhifang tunnel of Lanzhou-Chongqing Railway belongs to the category of the deep-buried tunnel. According to the structural size of the cavern and the physical and mechanical parameters of the stratum, the Platts pressure arch theory algorithm before and after the correction and the limit analysis upper-bound method and the Fraldi limit analysis method based on the loose zone boundary [5–8] are used to calculate the loose pressure of the top and sidewalls of the tunnel. The calculation results are shown in Tables 4 and 5. Except for the new method proposed in this paper, the other four methods cannot consider the influence of the buried depth of the cavern, so the calculation results in the corresponding two tables are the same. It is known from Tables 4 and 5 that when k = 1.2 to 1.4, for the vertical pressure value, the calculation result of the limit analysis upper limit method is similar to the calculation result of the stress transfer method based on the loose zone boundary and the field measured data.

Table 4 

The relaxation pressures of Zhifang tunnel when the depth of the tunnel is 65 m.

Table 5 

The relaxation pressures of Zhifang tunnel when the depth of the tunnel is 75 m.

4. Conclusion

This paper discusses the research status, basic principles, and application of the plastic limit analysis principle in geotechnical engineering. It is proposed to calculate the loose rock pressure of underground caverns in upper soft and hard rock strata by using the principle of plastic limit analysis upper-bound method. We get the following conclusions:(1)One of the advantages of using the limit analysis method is that the determined loose rock failure mode can minimize the damage process under unsupported conditions after tunnel excavation(2)In the derivation process of the calculation formula of the loose surrounding rock pressure, the parameters in the shear-sliding constitutive model of the structural plane are used, and the associated flow law is used to give the objective function for calculating the loose surrounding rock pressure. In the process of analysis and derivation, in order to avoid the complexity of solving the objective function, the vertical line where the central axis of the underground cavern is located is the axis of symmetry, with a point on the axis of symmetry as the center of rotation, and the ray from the center of rotation to the boundary of the loose area, and the method of making the angle between each rigid block boundary and the center of rotation equal. With this method, the calculation method can determine the maximum value by simply calculating the maximum value in mathematics and finally obtain the surrounding rock pressure, avoiding the use of various optimization methods to determine the upper limit solution, which greatly simplifies the calculation process. Finally, the theoretical formula for calculating the surrounding rock pressure is derived. The numerical example is given to calculate the loose pressure using the limit analysis upper-bound method. It is also proved that the limit analysis method is also suitable for calculating the surrounding rock loosening pressure of underground caverns(3)Although the calculation process does not involve the optimization solution of the optimization method to solve the objective function, the relative derivation process is still complicated. In the layered formation, the corresponding failure mode changes with the change of the softer upper layer thickness. The analysis process is cumbersome, and it is necessary to use MATLAB or other types of commercial mathematical software to solve the equation. The calculation formula of the surrounding rock pressure is only applicable to the case where the height of each layer of the formation is constant, and the versatility is poor(4)As we all know, the supporting force of the sidewall is very important for the stability of the underground cavern. In the process of derivation in this paper, in order to simplify the number of unknowns, it is assumed that the horizontal force and the vertical force have a proportional force relationship. On the one hand, this assumption can reflect the relationship between vertical force and horizontal force of the lining and highlight the effect of horizontal pressure on the lining; on the other hand, in practical applications, when taking different values of k, the corresponding values of q and e are obtained, which can be compared with the calculation results of other methods, and the experience of k value is further obtained

Data Availability

Some or all data, models, or codes used during the study were provided by a third party. Direct requests for these materials may be made to the provider.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was financially supported by the Technology Research and Development Program of the China Railway Corporation (K2018G013), the International Collaborative Project of the China Academy of Railway Sciences (2019YJ027), the Key Research and Development Project of Shandong Province (Grant no. 2019GSF111018), and Open Fund of State Key Laboratory of Water Resource Protection and Utilization in Coal Mining (GJNY-18-73.3).