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

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 diﬀerence 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 diﬀerent depths and diﬀerent types of caverns. The diﬀerence, rationality, and applicability of the calculation results of this method are analyzed and discussed.


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. erefore, 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. ese formulas include Fenner's formula, modified Fenner's formula, Casaer's formula, and Caquot's formula. e 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]. e 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. e limit analysis method was first created in the 1920s. In 1975, Chen [3] published the book Limiting Analysis and Soil Plasticity. e 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][6][7][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. erefore, 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. erefore, for the calculation of surrounding rock pressure in shallow tunnels, the use of the limit analysis method is an effective means.

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) e surface is assumed to be a horizontal plane. e 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) e 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.

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. e tunnel section is approximated by a rectangular DEFG. e 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. e height of point NA from the point C to the ground surface is H 0 . e line of GA m is a horizontal line connecting CA and GA m . e ∠ACA m is divided into m parts, and the angle is β. ∠A m GF 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 1 , 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.
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. e 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 1 . . .A i A i+1 . . .A m+n-1 F on the curve are connected by a straight line, and the shape of each triangular rigid slider is determined by β, ε, and α i . e speed of △ANC is V 0 , 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 1 . . .V m+n from top to bottom. e relative speeds of the respective sliding blocks are V 0,1 . . .V m+n-1,m+n , V m+n,m+n+1 .
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.
e 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. e height from point C to the ground surface is H 0 . e GA m line is a horizontal straight line connecting CB and CA m . e ∠BCA m is divided into m parts, the angle is β, and the ∠A m GF 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 1 , 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 0 . erefore, as long as the calculation formula of the surrounding rock pressure represented by H 0 is obtained, the extremum point can be obtained by deriving H 0 . 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. e 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 1 . . .A i A i+1 . . .A m+n-1 F on the curve are connected by a straight line, and the shape of each triangular rigid slider is determined by β, ε, and a i . e velocity of △BCA is V 0 . If the value of m is large enough, the BA segment is approximately horizontal, and then, the angle between V 0 and the vertical direction is infinitely close to φ J . is 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 1 . . .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 .

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 . e surrounding rock pressure is calculated by the limit theoretical upper limit method.

Mathematical Problems in Engineering
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 c 1 , the internal friction angle is φ J1 , the cohesive force is C J1 , the bulk density of the lower hard formation is c 2 , 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. e horizontal right is positive and the vertical upward is positive. e analytical formula of the FA curve is y � ax b . e corresponding velocity field at this time is shown in Figure 6.
(1) Solving the Geometric Quantities in the Failure Mode Use P r3 � i k�1 S Δ A k CA k+1 · c 1 · V k · cos[α k − (m − k)β + θ 3 + φ J1 − (π/2)]} as the origin of coordinates. Establish a Cartesian coordinate system with a positive horizontal to positive and positive vertical. e analytical formula of the FA curve is The coordinates of Mathematical Problems in Engineering e analytical expression is Similarly, the analytical formula of the CA i line can be obtained: e CA 1 linear equation is connected with the FA equation y � ax b of the slip line, and the coordinate A 1 (x 1 ,y 1 ) of the A 1 point can be obtained. Similarly, e equation of the upper soft layer and the lower hard Similarly, it is available that By Sine theorem, it is available that (2) Calculation of Velocity Field From Figure 3, the speed vector size of each slider is (3) Gravity power calculation e area of each rigid block is e gravity power of each sliding block of the damage mode (i + 3) is Figure 6: e velocity field schematic diagram of the (i + 1)-th failure.

Mathematical Problems in Engineering
From the above, the damage mode shown by the damage mode (i + 3) is half of the gravity power.
(4) Internal energy dissipation power Half of the energy dissipation in the damage mode (i + 3) is e 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 e virtual power equation: kPa top plate vertical support reaction force: In the above formula, q is a function of the variable H 0 ; 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 0 ). When q � f (H 0 ) takes the minimum value, the following conditions should be met: rough equation (22), the value of H 0 when q takes the extreme value can be obtained. It is also necessary to verify whether the obtained H 0 satisfies the constraint condition corresponding to the failure mode: If the above constraint is satisfied, the value of H 0 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

Detonation Pressure Detonation Process of Underground
Cavern in Annular Loose Zone. In the failure mode shown in Figure 3, the BC spacing H 0 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 c 1 , an internal friction angle of φ J1 , a cohesive force of C J1 , a lower hard formation density of c 2 , 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. e horizontal right is positive and the vertical upward is positive. e FA curve analytical expression is y � 1/(a + bx). e corresponding velocity field at this time is shown in Figure 8.
(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. e 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.
In the same way, the coordinates of A 1 and A i can be obtained: In the formula, η(ε) 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 Mathematical Problems in Engineering e distance between AA 1 can be obtained from A(x 0 , y 0 ) and A 1 (x 1 , y 1 ): e same is available:

|CB| �
��������������������� � 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 e area of each rigid block is (39) e 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 Mathematical Problems in Engineering 13 (5) Support reaction power e undetermined parameter k is the ratio of the horizontal support reaction force to the vertical support reaction force, k � e/q. e sum of the powers of the vertical support reaction force q and the horizontal support reaction force e of the tunnel roof is PF.
(43) (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 0 ). When q � f (H 0 ) takes the minimum value, the following conditions should be met: rough equation (47), the value of H 0 when q takes the extreme value can be obtained. It is also necessary to verify whether the obtained H 0 satisfies the constraint condition corresponding to the failure mode: If the above constraint condition is satisfied, the value of H 0 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:

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. e 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. e method calculates the vertical pressure and horizontal pressure of the tunnel.

Qingdao Jiangxi Road Metro
Station. e geology of Qingdao city is a typical upper soft and hard rock formation.
is 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. e ground level of the metro station is about 11.71 m. e overlying strongly weathered 1d granite has a calculated depth of 5.5 m, and the harder stratum is the late Yanshanian granite. e mid-weathered rock face has a large undulation, and the calculated depth is 6.  Table 1.

Zhifang Tunnel of Lanzhou-Chongqing Railway.
Zhifang tunnel of Lanzhou-Chongqing Railway [16] is a double-line tunnel in a carbonaceous slate. e 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. e tunnel section height is 12.  Table 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.
It is known from Table 3 that when k � 0.8-1.0, for the vertical pressure value, the calculation result of the upperbound 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 upperbound 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. e 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 deepburied 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][6][7][8]

Conclusion
is 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. e 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. e analysis process is cumbersome, and it is necessary to use MATLAB or other types of commercial mathematical software to solve the equation. e 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
e authors declare that they have no conflicts of interest.