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

Authors: Zhiwei Wang; Zhenghao Liu; Weibin Ma; Shaoshuai Shi; Lin Bu
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 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 pressureq, 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 drumIn 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 Figure1, 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 = axb, 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 H0. The line of GAm is a horizontal line connecting CA and GAm. The ∠ACAm is divided into m parts, and the angle is β. ∠AmGF 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, A1, A2, ...Am…Am+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, AA1…AiAi+1…Am+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 V0, the direction is vertical downward, the speed of △CDG is Vm+n+1, the direction is vertical downward, and the speed of the triangular block at the slip line is V1…Vm+n from top to bottom. The relative speeds of the respective sliding blocks are V0,1…Vm+n-1,m+n, Vm+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 Figure3, 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 H0. The GAm line is a horizontal straight line connecting CB and CAm. The ∠BCAm is divided into m parts, the angle is β, and the ∠AmGF 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, A1, A2, ...Am…Am+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 H0. Therefore, as long as the calculation formula of the surrounding rock pressure represented by H0 is obtained, the extremum point can be obtained by deriving H0. 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, AA1…AiAi+1…Am+n-1F on the curve are connected by a straight line, and the shape of each triangular rigid slider is determined by β, ε, and ai. The velocity of △BCA is V0. If the value of m is large enough, the BA segment is approximately horizontal, and then, the angle between V0 and the vertical direction is infinitely close to φJ. This paper assumes the horizontal of the BA segment; the velocity of △CDG is Vm+n+1, and the direction is vertical. Downward, the velocity of the triangular block at the slip line is V1…Vm+n from top to bottom, and the relative speed of each sliding block is V0,1…Vm+n-1,m+n, Vm+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 isy = axb. 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 Ai and Ai+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 γ1, the internal friction angle is φJ1, the cohesive force is CJ1, the bulk density of the lower hard formation is γ2, the internal friction angle is φJ2, and the cohesion is CJ2. 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 = axb. 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 ModeUsePr3=∑k=1iSΔAkCAk+1⋅γ1⋅Vk⋅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. The analytical formula of the FA curve is y = axb.Let|BC| = H0, |BO| = h, and tunnel span W; then,(1)OC=h−H0,AC=H02+H+Ta1/b+W22,MAm=Ta1/b+W2,Fromtanθ1=AB/H0, tanθ3=CM/MA5, tanθ4=CM/MG, and θ1, θ3, θ4 expression(2)OO1=OCtanθ5OO1=Wh−H02H−H0.(3)The coordinates ofπ2−β−φJ2>0,π−αi+1+θ2>0,π−αs−β−2φJ1>0,s=1…i,π−αt−β−2φJ2>0,s=i…m+n,αk−m−kβ+θ3+φJ1−π2>0,k=1…i,π−αv−ε−2φJ2>0,v=m+1…m+n,αl−m−lβ+θ3+φJ2−π2>0,l=i+1…m,αe−e−mε+φJ2−π2>0,e=m+1…m+n,are(4)O1Wh−H02H−H0−W2,T+H−h.The slope ofCR=H−H0,BC=V1−H+H0 and tanθ1=W/2H−H0 is tanθ3=2H−H0/L; then,(5)θ1=arctanW2H−H0,θ3=arctan2H−H0L.The analytical expression is(6)y=x+W2×tanπ2−θ1−β+T+H−H0.Similarly, the analytical formula of theCAi line can be obtained:(7)y=x+W2×tanπ2−θ1−iβ+T+H−H0.The CA1 linear equation is connected with the FA equation y = axb of the slip line, and the coordinate A1(x1,y1) of the A1 point can be obtained. Similarly, A2(x2,y2), Ai (xi,yi), Ai+1(xi+1,yi+1), and Am (xm,ym) can be obtained.The equation of the upper soft layer and the lower hard layer boundary lineOOm-i+3 is y = T + H-h, which is in line with the linear equation of the CAi+1 linear equation, and the coordinates of the Om-i+2 point are Om-i+2H0−h/tanπ/2−θ1−i+1β−W/2,T+H−h.The equation of the upper soft layer and the lower hard layer boundary lineOOm-i+3 is y = T + H-h, which is in conjunction with the FA equation y = axb of the slip line, and the coordinates of the Om-i+3 point can be obtained as Om-i+3T+H−h/a1/b,T+H−h, from O3 coordinates:(8)Om−i+2Om−i+3=T+H−ha1/b−H0−htanπ/2−θ1−i+1β+W2,available fromAH+T/a1/b+W/2,H+T, A1x1,y1(9)AA1=x1−H+Ta1/b−W22+y1−H−T2.Similarly, it is available that(10)AiAi+1=xi+1−xi2+yi+1−yi2,i=2…m+n−1,available from−W/2,T+H−H0, A1x1,y1(11)CAi=xi+W22+yi−T−H+H02.By Sine theorem, it is available that(12)Vk−1,k=sinm−iβ+θ3+φJ1+π/2sin3π/2−αk−β−3φJ1Vk−1,k=1…i,Vl−1,l=sinm−lβ+θ3+φJ2+π/2sin3π/2−αl−β−3φJ2Vl−1,l=i+1…m,Ve−1,e=sinαe−e−mε+φJ2+π/2sin3π/2−αe−ε−3φJ2Ve−1,e=m+1…m+n,αi=arcsinCAiAi−1Aisinβ,i=2…m.By Sine theorem, it is available fromGAm=T/a1/b, AmAm+1, ε that(13)αm+1=arcsinGAmAmAm+1sinε.(2)
Calculation of Velocity FieldFrom Figure3, the speed vector size of each slider is(14)Vk=sinm+1−kβ−θ3+φJ1+π/2sinπ−αk−β−2φJ1Vk−1,k=2…i,Vk−1,k=sinm−iβ+θ3+φJ1−π/2sinπ−αk−β−2φJ1Vk−1,k=2…i,Vl=sinm+1−lβ−θ3+φJ2+π/2sinπ−αl−β−2φJ2Vl−1,l=i+1…m,Vl−1,l=sinm−lβ+θ3+φJ2−π/2sinπ−αl−β−2φJ2Vl−1,l=i+1…m,Ve=sinβ+e−mε+φJ2+π/2sinπ−αe−β−2φJ2Ve−1,e=m+1…m+n,Ve−1,e=sinαe−e−mε+φJ2−π/2sinπ−αe−β−2φJ2Ve−1,e=m+1…m+n,Vm+n+1=sinπ−αm−1+θ2sinπ/2−ε−φJ2Vm+n,Vm+n,m+n+1=sinαm+n+θ3+φJ2−θ2−π/2sinπ/2−ε−φJ2Vm+n.(3)
Gravity power calculationThe area of each rigid block is(15)SΔABC=H02×H+Ta1/b+W2,SΔAkCAk+1=12AkC⋅CAk+1⋅sinβk=1…m,SΔOm−i+2Om−i+3Ai+1=Om−i+2Om−i+32×T+H−h−yi+1,SΔAiCOm−i+2Om−i+3=SΔCAiAi+1−SΔOm−i+2Om−i+3Ai+1,SΔAeGAe+1=12GAe⋅GAe+1⋅sinεe=m…m+n−2,SΔAm+n−1GF=12GAm+n−1⋅T⋅sinε,SΔCGAm=CG2⋅W2⋅sinθ2+θ3,SΔCGM=H−H02⋅W2,SΔCGAmAi+1=SΔCAmAi+1+SΔCAmG.The gravity power of each sliding block of the damage mode (i + 3) is(16)Pr1=SΔABC⋅γ1⋅V0,Pr2=SΔACA1⋅γ1⋅V1⋅cosα1−m−1β+θ3+φJ1−π2,Pr3=∑k=1iSΔAkCAk+1⋅γ1⋅Vk⋅cosαk−m−kβ+θ3+φJ1−π2,Pr4=SΔAiCOm−i+2Om−i+3⋅γ1⋅Vi+1⋅cosαi+1−m−i−1β+θ3+φJ1−π2,Pr5=SΔOm−i+3Om−i+2Ai+1⋅γ2⋅Vi+1⋅cosαi+1−m−i−1β+θ3+φJ2−π2,Pr6=∑k=i+1m−1SΔAkCAk+1⋅γ2⋅Vk⋅cosαk−m−kβ+θ3+φJ2−π2,Pr7=SΔAm−1CGAm⋅γ2⋅Vm⋅cosαm+θ3+φJ2−π2,Pr8=∑e=mm+n−1SΔAeGAe+1⋅γ2⋅Ve⋅cosαe+1−e−mε+φJ2−π2,Pr9=SΔAm+n−1GF⋅γ2⋅Vm+n⋅cosαm+n−nε+φJ2−π2.From the above, the damage mode shown by the damage mode (i + 3) is half of the gravity power.(17)Pr=∑c=114Prc.(4)
Internal energy dissipation powerHalf of the energy dissipation in the damage mode (i + 3) is(18)PW=CJ1⋅cosφJ1∑k=1iAAk⋅Vk+AC⋅V0,1+∑k=1iAkC⋅Vk,k+1+CJ2cosφJ2Om−i+2C⋅Vi+1,i+2+∑e=i+1m+n−1Ae−1Ae⋅Ve+∑k=i+1mAkC⋅Vk,k+1+Am+n−1F⋅Vm+n+∑e=mm+n−1AeG⋅Ve,e+1.(5) Supporting reaction power(19)PF=W⋅q⋅Vm+n+1−T⋅k⋅q⋅Vm+n⋅sinαm+n+φJ2.The undetermined parameterk is the ratio of the horizontal support reaction force to the vertical support reaction force, k = e/q.(6)
Calculation of surrounding rock pressureThe virtual power equation:(20)Pr+PF=PW⟺∫Ω∗W∗υ∗dΩ∗+∫A∗T∗υ∗dA∗=∫Ω∗σij∗ε˙ij∗dΩ∗+∫S∗τJ−σnJtanφJΔvt∗dS∗,kPa top plate vertical support reaction force:(21)q=Pr−PWW⋅V0+T⋅k⋅sinαm+n+φJ2⋅Vm+n.In the above formula,q is a function of the variable H0; that is,(22)q=fH0.It is known from the upper limit theorem that the vertical support reaction force of the top plate is the minimum value ofq = f (H0). When q = f (H0) takes the minimum value, the following conditions should be met:(23)dqdH0=0.Through equation (22), the value of H0 when q takes the extreme value can be obtained. It is also necessary to verify whether the obtained H0 satisfies the constraint condition corresponding to the failure mode:(24)π2−β−φJ2>0,π−αi+1+θ2>0,π−αs−β−2φJ1>0,s=1…i,π−αt−β−2φJ2>0,s=i…m+n,αk−m−kβ+θ3+φJ1−π2>0,k=1…i,π−αv−ε−2φJ2>0,v=m+1…m+n,αl−m−lβ+θ3+φJ2−π2>0,l=i+1…m,αe−e−mε+φJ2−π2>0,e=m+1…m+n.If the above constraint is satisfied, the value ofH0 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(25)e=kq.


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

In the failure mode shown in Figure3, the BC spacing H0 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 Ai and Ai+1. Let the upper softer formation have a bulk density of γ1, an internal friction angle of φJ1, a cohesive force of CJ1, a lower hard formation density of γ2, an internal friction angle of φJ2, and a cohesive force of CJ2. 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)
Figure7 shows the solution of each geometric quantity in the loose failure mode.As shown in Figure7, 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| = V1, |NmAm| = L, and tunnel span W, then(26)CR=H−H0,BC=V1−H+H0.Fromtanθ1=W/2H−H0 and tanθ3=2H−H0/L, the following can be obtained:(27)θ1=arctanW2H−H0,θ3=arctan2H−H0L.Fromθ2 = π/2-θ1-θ3, the following can be obtained:(28)β=π/2+θ3m,ε=π2n.The coordinates of pointC are C[0, (H-H0-V1)], and the slope of the line CA is cotβ; then, the straight-line equation of CA is y = x cot β + (H-H0-V1). Similarly, the linear equations of CA1, CA2, CAi, and CAi+1 are y = x cot(2β) + (H-H0-V1), y = x cot(3β)+(H-H0-V1), y = xcot [(i + 1)β] + (H-H0-V1), y = x cot[(i + 2)β] + (H-H0-V1), and y = x cot β + (H-H0-V1).Let(V1+H0-H) =  H∗; the linear equation of the simultaneous CA and the boundary equation of the loose zone can get the coordinates of the point.(29)Ax0,y0⇔AbH∗−acotβ+ξβ2bcotβ,2bcotβbH∗+acotβ+ξβ.In the formula,ξ(β) = (a cot ß + bH∗)2 + 4b cot ß.In the same way, the coordinates ofA1 and Ai can be obtained:(30)A1x1,y1⇔A1bH∗−acotβ+ξ2β2bcot2β,2bcot2βbH∗+acot2β+ξ2β,Aixi,yi⇔AibH∗−acoti+1β+ξi+1β2bcoti+1β,2bcoti+1βbH∗+acoti+1β+ξi+1β.In the formula,(31)ξ2β=acot2β+bH∗2+4bcot2β,ξ3β=acot3β+bH∗2+4bcot3β,ξi+1β=acoti+1β+bH∗2+4bcoti+1β,ξi+2β=acoti+2β+bH∗2+4bcoti+2β.The straight-line equation of the interface between the upper softer stratum and the lower harder stratum isy = H-V1-h =  h∗, and the coordinates of the O1 point can be obtained in conjunction with the CA straight line equation, O1[(h∗ + H∗)/cosβ, h∗]. Similarly, y =  h∗ and the linear equations of CA1, CA2, and CAi can be used to obtain the coordinates of O2, O3, and Oi+1 points: O2[(h∗ + H∗)/cos 2β, h∗], O3[(h∗ + H∗)/cos3β, h∗], and Oi+1[(h∗ + H∗)/cos (i + 1)β. y =  h∗ is associated with the boundary curve equation of the loose zone and the coordinates of the Oi+2 point: Oi+2 [(1−ab)/bh∗, h∗]. The distance between the points can be obtained from the coordinates of the above points. The coordinate of point G is G(W/2, −V1), the slope of the line of GAm+1 is -tanε, and the equation of the line of GAm+1 is y = y = (−tanε) (x−W/2)−V1. Similarly, the linear equations for GAm+2, GAm+j, and GAm+n-1 are y = [−tan(2ε)](x−W/2)−V1, y = [−tan(jε)](x−W/2)-V1, and y = {−tan[(m + n−1)ε)](x−W/2)−V1.LetbW−2a =  W∗; the line equation of the simultaneous GAm+1 and the boundary curve equation of the loose zone can get the coordinates of the Am+1 point:(32)Am+1W∗tanε−2bV1+ηε4btanε,4tanε4atanε+W∗tanε−2bV1+ηε.In the formula,η(ε) = b2V12 + (a2 + b2W2/4 + abW)tan2ε−(b2V12W + 2abV1 + 4b)tan ε.Similarly, the coordinates ofAm+2, Am+j, and Am+n-1 can be obtained:(33)Am+2W∗tan2ε−2bV1+η2ε4btan2ε,4tan2ε4atan2ε+W∗tan2ε−2bV1+η2ε,Am+jW∗tanjε−2bV1+ηjε4btanjε,4tanjε4atanjε+W∗tanjε−2bV1+ηjε,Am+n−1W∗tanm+n−1ε−2bV1+ηm+n−1ε4btanm+n−1ε,4tanm+n−1ε4atanm+n−1ε+W∗tanm+n−1ε−2bV1+ηm+n−1ε.In the formula,η(2ε), η(jε), and η[(m + n−1)ε] are(34)η2ε=b2V12+a2+b2W24+abWtan22ε−b2V12W+2abV1+4btan2ε,ηjε=b2V12+a2+b2W24+abWtan2jε−b2V12W+2abV1+4btanjε,ηm+n−1ε=b2V12+a2+b2W24+abWtan2m+n−1ε−b2V12W+2abV1+4btanm+n−1ε.The distance betweenAA1 can be obtained from A(x0, y0) and A1(x1, y1):(35)AA1=x1−x02+y1−y02.The same is available:(36)BA=x02+y1−1a2,CB=x02+1a−H−H0−V12,CAi=xi2+yi−1a2,i=1…m,GAj=xj−W22+yj+V12,j=m…m+n−1,AiAi+1=xi+1−xi2+yi+1−yi2,i=1…m+n−1,Am+n−1F=xm+n−1−W22+ym+n−1+V1+T2.Using the sine theorem, we can findα, α1…αm + n:(37)α=arcsinCBABsinβ,α1=arcsinCAAA1sinβ,αi=arcsinCAi−1Ai−1Aisinβ,i=2…m,αj=arcsinGAj−1Gj−1Ajsinε,j=m+1…m+n−1,αm+n=arcsinGAm+n−1FAm+n−1sinε.(2)
Speed field calculationFrom Figure3, the sine theorem can be used to find the velocity vector size of each slider:(38)Vk=sinm+1−kβ−θ3+φJ1+π/2sin3π/2−αk−β−3φJ1Vk−1,k=1…i,Vl=sinm+1−lβ−θ3+φJ2+π/2sin3π/2−αl−β−3φJ2Vl−1,l=i+1…m,Ve=sinε+e−mε+φJ2+π/2sin3π/2−αe−ε−3φJ2Ve−1,e=m+1…m+n,Vm+n+1=sinπ−αm−1+θ2sin3π/2−ε−2φJ2Vm+n,Vk−1,k=sinm−iβ+θ3+φJ1+π/2sin3π/2−αk−β−3φJ1Vk−1,k=1…i,Vl−1,l=sinm−lβ+θ3+φJ2+π/2sin3π/2−αl−β−3φJ2Vl−1,l=i+1…m,Ve−1,e=sinαe−e−mε+φJ2+π/2sin3π/2−αe−ε−3φJ2Ve−1,e=m+1…m+n,Vm+n,m+n+1=sinαm+n+θ3+φJ2−θ2+π/2sin3π/2−ε−φJ2Vm+n.(3)
Gravity power calculationThe area of each rigid block is(39)SΔABC=12CB⋅CA⋅sinβ,SΔACA1=12AC⋅CA1⋅sinβ,SΔAkCAk+1=12AkC⋅CAk+1⋅sinβ,k=1…m,SΔCOO1=12OC⋅OO1⋅sinβ,SΔCO1O2=12OC⋅O1O2⋅sinβ,SΔCO2O3=12OC⋅O2O3⋅sinβ,SΔCO3Oi+1=12OC⋅O3Oi+1⋅sinβ,SΔOi+1AiOi+2=12yi−h∗⋅Oi+1Oi+2⋅sinβ,SΔBOO1A=SΔCBA−SΔCOO1,SΔAO1O2A1=SΔCAA1−SΔCO1O2,SΔA1O2O3A2=SΔCA1A2−SΔCO2O3,SΔA2O3Oi+1Ai=SΔCA2Ai−SΔCO3Oi+1,SΔOi+1CAi+1Oi+2=SΔCAiAi+1−SΔAiOi+1Oi+2,SΔCGAm=14CR⋅L−W,SΔCRG=14CR⋅W,SΔGAjAj+1=12GAj⋅GAj+1⋅sinεj=m…m+n−1.The weight of each sliding block of the damage mode (i + 1) is(40)Pr1=S▱BOO1A⋅γ1⋅V0cosα+β−φJ1,Pr2=SΔCOO1⋅γ2⋅V0cosα+β−φJ2,Pr3=S▱AO1O2A1⋅γ1⋅V1cosα1+2β−φJ1,Pr4=SΔCO1O2⋅γ2⋅V1cosα1+2β−φJ2,Pr5=S▱A1O2O3A2⋅γ1⋅V2cosα2+3β−φJ1,Pr6=SΔCO2O3⋅γ2⋅V2cosα2+3β−φJ2,Pr7=S▱A2O3Oi+1Ai⋅γ1⋅Vicosαi+i+1β−φJ1,Pr8=SΔCO3Oi+1⋅γ2⋅Vicosαi+i+1β−φJ2,Pr9=SΔAiOi+1Oi+2⋅γ1⋅Vi+1cosαi+1+i+2β−φJ1,Pr10=S▱Oi+1CAi+1Oi+2⋅γ2⋅Vi+1cosαi+1+i+2β−φJ2,Pr11=SΔCAi+1Am⋅γ2⋅Vmcosαm+m+1β−φJ2,Pr12=∑e=mm+n−1SΔAeGAe+1⋅γ2⋅Ve+1⋅cosαe+1+e+1−mε−φJ2,Pr13=SΔGFAm+n−1⋅γ2⋅Vm+ncosαm+n+nε−φJ2,Pr14=SΔCRG⋅γ2⋅Vm+n+1.Half of the damage mode gravity power shown by the upper damage mode (i + 1) is(41)Pr=∑c=114Prc.(4)
Internal energy dissipation powerHalf of the energy dissipation power in the failure mode (i + 1) is(42)PW=CJ1⋅cosφJ1∑k=1iAAk⋅Vk+∑t=1iOtAt−1⋅Vt−1,t+BA⋅V0+AiOi+2⋅Vi+1+CJ2cosφJ2∑t=1i+1COt⋅Vt−1,t+∑p=i+1mCAp⋅Vp−1,p+Ai+1Oi+2⋅Vi+1+∑q=i+1mAqAq+1⋅Vq+1+∑e=m+1m+n−1Ae−1Ae⋅Ve+mm+n−1F⋅Vm+n+∑e=mm+n−1AeG⋅Ve,e+1+CR⋅Vm+n+1+CG⋅Vm+n,m+n+1.(5)
Support reaction powerThe undetermined parameterk 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.(43)PF=W⋅q⋅Vm+n+1−T⋅k⋅q⋅Vm+n⋅sinαm+n+φJ2.(6)
Calculation of surrounding rock pressureVirtual power equation:(44)Pr+PF=PW⇔∫Ω∗W∗υ∗dΩ∗+∫A∗T∗υ∗dA∗=∫Ω∗σij∗ε˙ij∗dΩ∗+∫S∗τJ−σnJtanφJΔvt∗dS∗,V0,1=sinα1−4β+θ3+φ1−π/2sinπ−α1−β−2φ1V0.Roof vertical support reaction:(45)q=Pr−PWW⋅V0+T⋅k⋅sinαm+n+φJ2⋅Vm+n.In the above formula,q is a function of the variable H0, which is(46)q=fH0.It is known from the upper limit theorem that the vertical support reaction force of the top plate is the minimum value ofq = f (H0). When q = f (H0) takes the minimum value, the following conditions should be met:(47)dqdH0=0.Through equation (47), the value of H0 when q takes the extreme value can be obtained. It is also necessary to verify whether the obtained H0 satisfies the constraint condition corresponding to the failure mode:(48)π−αi+1+θ2>0,π−β−2φJ2>0,3π2−αs−β−3φJ1>0,s=1…i,3π2−αt−β−3φJ2>0,s=i+1…m,3π2−αu−ε−3φJ2>0,u=m+1…m+n,αk−m−kβ+θ3+2φJ1−π>0,k=1…i,αl−m−lβ+θ3+2φJ2−π>0,i=i+1…m,αe−e−mε+2φJ2−π>0,e=m+1…m+n−1.If the above constraint condition is satisfied, the value ofH0 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:(49)e=kq.


## 2.1. Destruction Mode and Velocity Field

For the underground cavern with the loose zone as the open type, as shown in Figure1, 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 = axb, 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 H0. The line of GAm is a horizontal line connecting CA and GAm. The ∠ACAm is divided into m parts, and the angle is β. ∠AmGF 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, A1, A2, ...Am…Am+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, AA1…AiAi+1…Am+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 V0, the direction is vertical downward, the speed of △CDG is Vm+n+1, the direction is vertical downward, and the speed of the triangular block at the slip line is V1…Vm+n from top to bottom. The relative speeds of the respective sliding blocks are V0,1…Vm+n-1,m+n, Vm+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 Figure3, 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 H0. The GAm line is a horizontal straight line connecting CB and CAm. The ∠BCAm is divided into m parts, the angle is β, and the ∠AmGF 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, A1, A2, ...Am…Am+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 H0. Therefore, as long as the calculation formula of the surrounding rock pressure represented by H0 is obtained, the extremum point can be obtained by deriving H0. 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, AA1…AiAi+1…Am+n-1F on the curve are connected by a straight line, and the shape of each triangular rigid slider is determined by β, ε, and ai. The velocity of △BCA is V0. If the value of m is large enough, the BA segment is approximately horizontal, and then, the angle between V0 and the vertical direction is infinitely close to φJ. This paper assumes the horizontal of the BA segment; the velocity of △CDG is Vm+n+1, and the direction is vertical. Downward, the velocity of the triangular block at the slip line is V1…Vm+n from top to bottom, and the relative speed of each sliding block is V0,1…Vm+n-1,m+n, Vm+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 isy = axb. 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 Ai and Ai+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 γ1, the internal friction angle is φJ1, the cohesive force is CJ1, the bulk density of the lower hard formation is γ2, the internal friction angle is φJ2, and the cohesion is CJ2. 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 = axb. 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 ModeUsePr3=∑k=1iSΔAkCAk+1⋅γ1⋅Vk⋅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. The analytical formula of the FA curve is y = axb.Let|BC| = H0, |BO| = h, and tunnel span W; then,(1)OC=h−H0,AC=H02+H+Ta1/b+W22,MAm=Ta1/b+W2,Fromtanθ1=AB/H0, tanθ3=CM/MA5, tanθ4=CM/MG, and θ1, θ3, θ4 expression(2)OO1=OCtanθ5OO1=Wh−H02H−H0.(3)The coordinates ofπ2−β−φJ2>0,π−αi+1+θ2>0,π−αs−β−2φJ1>0,s=1…i,π−αt−β−2φJ2>0,s=i…m+n,αk−m−kβ+θ3+φJ1−π2>0,k=1…i,π−αv−ε−2φJ2>0,v=m+1…m+n,αl−m−lβ+θ3+φJ2−π2>0,l=i+1…m,αe−e−mε+φJ2−π2>0,e=m+1…m+n,are(4)O1Wh−H02H−H0−W2,T+H−h.The slope ofCR=H−H0,BC=V1−H+H0 and tanθ1=W/2H−H0 is tanθ3=2H−H0/L; then,(5)θ1=arctanW2H−H0,θ3=arctan2H−H0L.The analytical expression is(6)y=x+W2×tanπ2−θ1−β+T+H−H0.Similarly, the analytical formula of theCAi line can be obtained:(7)y=x+W2×tanπ2−θ1−iβ+T+H−H0.The CA1 linear equation is connected with the FA equation y = axb of the slip line, and the coordinate A1(x1,y1) of the A1 point can be obtained. Similarly, A2(x2,y2), Ai (xi,yi), Ai+1(xi+1,yi+1), and Am (xm,ym) can be obtained.The equation of the upper soft layer and the lower hard layer boundary lineOOm-i+3 is y = T + H-h, which is in line with the linear equation of the CAi+1 linear equation, and the coordinates of the Om-i+2 point are Om-i+2H0−h/tanπ/2−θ1−i+1β−W/2,T+H−h.The equation of the upper soft layer and the lower hard layer boundary lineOOm-i+3 is y = T + H-h, which is in conjunction with the FA equation y = axb of the slip line, and the coordinates of the Om-i+3 point can be obtained as Om-i+3T+H−h/a1/b,T+H−h, from O3 coordinates:(8)Om−i+2Om−i+3=T+H−ha1/b−H0−htanπ/2−θ1−i+1β+W2,available fromAH+T/a1/b+W/2,H+T, A1x1,y1(9)AA1=x1−H+Ta1/b−W22+y1−H−T2.Similarly, it is available that(10)AiAi+1=xi+1−xi2+yi+1−yi2,i=2…m+n−1,available from−W/2,T+H−H0, A1x1,y1(11)CAi=xi+W22+yi−T−H+H02.By Sine theorem, it is available that(12)Vk−1,k=sinm−iβ+θ3+φJ1+π/2sin3π/2−αk−β−3φJ1Vk−1,k=1…i,Vl−1,l=sinm−lβ+θ3+φJ2+π/2sin3π/2−αl−β−3φJ2Vl−1,l=i+1…m,Ve−1,e=sinαe−e−mε+φJ2+π/2sin3π/2−αe−ε−3φJ2Ve−1,e=m+1…m+n,αi=arcsinCAiAi−1Aisinβ,i=2…m.By Sine theorem, it is available fromGAm=T/a1/b, AmAm+1, ε that(13)αm+1=arcsinGAmAmAm+1sinε.(2)
Calculation of Velocity FieldFrom Figure3, the speed vector size of each slider is(14)Vk=sinm+1−kβ−θ3+φJ1+π/2sinπ−αk−β−2φJ1Vk−1,k=2…i,Vk−1,k=sinm−iβ+θ3+φJ1−π/2sinπ−αk−β−2φJ1Vk−1,k=2…i,Vl=sinm+1−lβ−θ3+φJ2+π/2sinπ−αl−β−2φJ2Vl−1,l=i+1…m,Vl−1,l=sinm−lβ+θ3+φJ2−π/2sinπ−αl−β−2φJ2Vl−1,l=i+1…m,Ve=sinβ+e−mε+φJ2+π/2sinπ−αe−β−2φJ2Ve−1,e=m+1…m+n,Ve−1,e=sinαe−e−mε+φJ2−π/2sinπ−αe−β−2φJ2Ve−1,e=m+1…m+n,Vm+n+1=sinπ−αm−1+θ2sinπ/2−ε−φJ2Vm+n,Vm+n,m+n+1=sinαm+n+θ3+φJ2−θ2−π/2sinπ/2−ε−φJ2Vm+n.(3)
Gravity power calculationThe area of each rigid block is(15)SΔABC=H02×H+Ta1/b+W2,SΔAkCAk+1=12AkC⋅CAk+1⋅sinβk=1…m,SΔOm−i+2Om−i+3Ai+1=Om−i+2Om−i+32×T+H−h−yi+1,SΔAiCOm−i+2Om−i+3=SΔCAiAi+1−SΔOm−i+2Om−i+3Ai+1,SΔAeGAe+1=12GAe⋅GAe+1⋅sinεe=m…m+n−2,SΔAm+n−1GF=12GAm+n−1⋅T⋅sinε,SΔCGAm=CG2⋅W2⋅sinθ2+θ3,SΔCGM=H−H02⋅W2,SΔCGAmAi+1=SΔCAmAi+1+SΔCAmG.The gravity power of each sliding block of the damage mode (i + 3) is(16)Pr1=SΔABC⋅γ1⋅V0,Pr2=SΔACA1⋅γ1⋅V1⋅cosα1−m−1β+θ3+φJ1−π2,Pr3=∑k=1iSΔAkCAk+1⋅γ1⋅Vk⋅cosαk−m−kβ+θ3+φJ1−π2,Pr4=SΔAiCOm−i+2Om−i+3⋅γ1⋅Vi+1⋅cosαi+1−m−i−1β+θ3+φJ1−π2,Pr5=SΔOm−i+3Om−i+2Ai+1⋅γ2⋅Vi+1⋅cosαi+1−m−i−1β+θ3+φJ2−π2,Pr6=∑k=i+1m−1SΔAkCAk+1⋅γ2⋅Vk⋅cosαk−m−kβ+θ3+φJ2−π2,Pr7=SΔAm−1CGAm⋅γ2⋅Vm⋅cosαm+θ3+φJ2−π2,Pr8=∑e=mm+n−1SΔAeGAe+1⋅γ2⋅Ve⋅cosαe+1−e−mε+φJ2−π2,Pr9=SΔAm+n−1GF⋅γ2⋅Vm+n⋅cosαm+n−nε+φJ2−π2.From the above, the damage mode shown by the damage mode (i + 3) is half of the gravity power.(17)Pr=∑c=114Prc.(4)
Internal energy dissipation powerHalf of the energy dissipation in the damage mode (i + 3) is(18)PW=CJ1⋅cosφJ1∑k=1iAAk⋅Vk+AC⋅V0,1+∑k=1iAkC⋅Vk,k+1+CJ2cosφJ2Om−i+2C⋅Vi+1,i+2+∑e=i+1m+n−1Ae−1Ae⋅Ve+∑k=i+1mAkC⋅Vk,k+1+Am+n−1F⋅Vm+n+∑e=mm+n−1AeG⋅Ve,e+1.(5) Supporting reaction power(19)PF=W⋅q⋅Vm+n+1−T⋅k⋅q⋅Vm+n⋅sinαm+n+φJ2.The undetermined parameterk is the ratio of the horizontal support reaction force to the vertical support reaction force, k = e/q.(6)
Calculation of surrounding rock pressureThe virtual power equation:(20)Pr+PF=PW⟺∫Ω∗W∗υ∗dΩ∗+∫A∗T∗υ∗dA∗=∫Ω∗σij∗ε˙ij∗dΩ∗+∫S∗τJ−σnJtanφJΔvt∗dS∗,kPa top plate vertical support reaction force:(21)q=Pr−PWW⋅V0+T⋅k⋅sinαm+n+φJ2⋅Vm+n.In the above formula,q is a function of the variable H0; that is,(22)q=fH0.It is known from the upper limit theorem that the vertical support reaction force of the top plate is the minimum value ofq = f (H0). When q = f (H0) takes the minimum value, the following conditions should be met:(23)dqdH0=0.Through equation (22), the value of H0 when q takes the extreme value can be obtained. It is also necessary to verify whether the obtained H0 satisfies the constraint condition corresponding to the failure mode:(24)π2−β−φJ2>0,π−αi+1+θ2>0,π−αs−β−2φJ1>0,s=1…i,π−αt−β−2φJ2>0,s=i…m+n,αk−m−kβ+θ3+φJ1−π2>0,k=1…i,π−αv−ε−2φJ2>0,v=m+1…m+n,αl−m−lβ+θ3+φJ2−π2>0,l=i+1…m,αe−e−mε+φJ2−π2>0,e=m+1…m+n.If the above constraint is satisfied, the value ofH0 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(25)e=kq.


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

In the failure mode shown in Figure3, the BC spacing H0 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 Ai and Ai+1. Let the upper softer formation have a bulk density of γ1, an internal friction angle of φJ1, a cohesive force of CJ1, a lower hard formation density of γ2, an internal friction angle of φJ2, and a cohesive force of CJ2. 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)
Figure7 shows the solution of each geometric quantity in the loose failure mode.As shown in Figure7, 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| = V1, |NmAm| = L, and tunnel span W, then(26)CR=H−H0,BC=V1−H+H0.Fromtanθ1=W/2H−H0 and tanθ3=2H−H0/L, the following can be obtained:(27)θ1=arctanW2H−H0,θ3=arctan2H−H0L.Fromθ2 = π/2-θ1-θ3, the following can be obtained:(28)β=π/2+θ3m,ε=π2n.The coordinates of pointC are C[0, (H-H0-V1)], and the slope of the line CA is cotβ; then, the straight-line equation of CA is y = x cot β + (H-H0-V1). Similarly, the linear equations of CA1, CA2, CAi, and CAi+1 are y = x cot(2β) + (H-H0-V1), y = x cot(3β)+(H-H0-V1), y = xcot [(i + 1)β] + (H-H0-V1), y = x cot[(i + 2)β] + (H-H0-V1), and y = x cot β + (H-H0-V1).Let(V1+H0-H) =  H∗; the linear equation of the simultaneous CA and the boundary equation of the loose zone can get the coordinates of the point.(29)Ax0,y0⇔AbH∗−acotβ+ξβ2bcotβ,2bcotβbH∗+acotβ+ξβ.In the formula,ξ(β) = (a cot ß + bH∗)2 + 4b cot ß.In the same way, the coordinates ofA1 and Ai can be obtained:(30)A1x1,y1⇔A1bH∗−acotβ+ξ2β2bcot2β,2bcot2βbH∗+acot2β+ξ2β,Aixi,yi⇔AibH∗−acoti+1β+ξi+1β2bcoti+1β,2bcoti+1βbH∗+acoti+1β+ξi+1β.In the formula,(31)ξ2β=acot2β+bH∗2+4bcot2β,ξ3β=acot3β+bH∗2+4bcot3β,ξi+1β=acoti+1β+bH∗2+4bcoti+1β,ξi+2β=acoti+2β+bH∗2+4bcoti+2β.The straight-line equation of the interface between the upper softer stratum and the lower harder stratum isy = H-V1-h =  h∗, and the coordinates of the O1 point can be obtained in conjunction with the CA straight line equation, O1[(h∗ + H∗)/cosβ, h∗]. Similarly, y =  h∗ and the linear equations of CA1, CA2, and CAi can be used to obtain the coordinates of O2, O3, and Oi+1 points: O2[(h∗ + H∗)/cos 2β, h∗], O3[(h∗ + H∗)/cos3β, h∗], and Oi+1[(h∗ + H∗)/cos (i + 1)β. y =  h∗ is associated with the boundary curve equation of the loose zone and the coordinates of the Oi+2 point: Oi+2 [(1−ab)/bh∗, h∗]. The distance between the points can be obtained from the coordinates of the above points. The coordinate of point G is G(W/2, −V1), the slope of the line of GAm+1 is -tanε, and the equation of the line of GAm+1 is y = y = (−tanε) (x−W/2)−V1. Similarly, the linear equations for GAm+2, GAm+j, and GAm+n-1 are y = [−tan(2ε)](x−W/2)−V1, y = [−tan(jε)](x−W/2)-V1, and y = {−tan[(m + n−1)ε)](x−W/2)−V1.LetbW−2a =  W∗; the line equation of the simultaneous GAm+1 and the boundary curve equation of the loose zone can get the coordinates of the Am+1 point:(32)Am+1W∗tanε−2bV1+ηε4btanε,4tanε4atanε+W∗tanε−2bV1+ηε.In the formula,η(ε) = b2V12 + (a2 + b2W2/4 + abW)tan2ε−(b2V12W + 2abV1 + 4b)tan ε.Similarly, the coordinates ofAm+2, Am+j, and Am+n-1 can be obtained:(33)Am+2W∗tan2ε−2bV1+η2ε4btan2ε,4tan2ε4atan2ε+W∗tan2ε−2bV1+η2ε,Am+jW∗tanjε−2bV1+ηjε4btanjε,4tanjε4atanjε+W∗tanjε−2bV1+ηjε,Am+n−1W∗tanm+n−1ε−2bV1+ηm+n−1ε4btanm+n−1ε,4tanm+n−1ε4atanm+n−1ε+W∗tanm+n−1ε−2bV1+ηm+n−1ε.In the formula,η(2ε), η(jε), and η[(m + n−1)ε] are(34)η2ε=b2V12+a2+b2W24+abWtan22ε−b2V12W+2abV1+4btan2ε,ηjε=b2V12+a2+b2W24+abWtan2jε−b2V12W+2abV1+4btanjε,ηm+n−1ε=b2V12+a2+b2W24+abWtan2m+n−1ε−b2V12W+2abV1+4btanm+n−1ε.The distance betweenAA1 can be obtained from A(x0, y0) and A1(x1, y1):(35)AA1=x1−x02+y1−y02.The same is available:(36)BA=x02+y1−1a2,CB=x02+1a−H−H0−V12,CAi=xi2+yi−1a2,i=1…m,GAj=xj−W22+yj+V12,j=m…m+n−1,AiAi+1=xi+1−xi2+yi+1−yi2,i=1…m+n−1,Am+n−1F=xm+n−1−W22+ym+n−1+V1+T2.Using the sine theorem, we can findα, α1…αm + n:(37)α=arcsinCBABsinβ,α1=arcsinCAAA1sinβ,αi=arcsinCAi−1Ai−1Aisinβ,i=2…m,αj=arcsinGAj−1Gj−1Ajsinε,j=m+1…m+n−1,αm+n=arcsinGAm+n−1FAm+n−1sinε.(2)
Speed field calculationFrom Figure3, the sine theorem can be used to find the velocity vector size of each slider:(38)Vk=sinm+1−kβ−θ3+φJ1+π/2sin3π/2−αk−β−3φJ1Vk−1,k=1…i,Vl=sinm+1−lβ−θ3+φJ2+π/2sin3π/2−αl−β−3φJ2Vl−1,l=i+1…m,Ve=sinε+e−mε+φJ2+π/2sin3π/2−αe−ε−3φJ2Ve−1,e=m+1…m+n,Vm+n+1=sinπ−αm−1+θ2sin3π/2−ε−2φJ2Vm+n,Vk−1,k=sinm−iβ+θ3+φJ1+π/2sin3π/2−αk−β−3φJ1Vk−1,k=1…i,Vl−1,l=sinm−lβ+θ3+φJ2+π/2sin3π/2−αl−β−3φJ2Vl−1,l=i+1…m,Ve−1,e=sinαe−e−mε+φJ2+π/2sin3π/2−αe−ε−3φJ2Ve−1,e=m+1…m+n,Vm+n,m+n+1=sinαm+n+θ3+φJ2−θ2+π/2sin3π/2−ε−φJ2Vm+n.(3)
Gravity power calculationThe area of each rigid block is(39)SΔABC=12CB⋅CA⋅sinβ,SΔACA1=12AC⋅CA1⋅sinβ,SΔAkCAk+1=12AkC⋅CAk+1⋅sinβ,k=1…m,SΔCOO1=12OC⋅OO1⋅sinβ,SΔCO1O2=12OC⋅O1O2⋅sinβ,SΔCO2O3=12OC⋅O2O3⋅sinβ,SΔCO3Oi+1=12OC⋅O3Oi+1⋅sinβ,SΔOi+1AiOi+2=12yi−h∗⋅Oi+1Oi+2⋅sinβ,SΔBOO1A=SΔCBA−SΔCOO1,SΔAO1O2A1=SΔCAA1−SΔCO1O2,SΔA1O2O3A2=SΔCA1A2−SΔCO2O3,SΔA2O3Oi+1Ai=SΔCA2Ai−SΔCO3Oi+1,SΔOi+1CAi+1Oi+2=SΔCAiAi+1−SΔAiOi+1Oi+2,SΔCGAm=14CR⋅L−W,SΔCRG=14CR⋅W,SΔGAjAj+1=12GAj⋅GAj+1⋅sinεj=m…m+n−1.The weight of each sliding block of the damage mode (i + 1) is(40)Pr1=S▱BOO1A⋅γ1⋅V0cosα+β−φJ1,Pr2=SΔCOO1⋅γ2⋅V0cosα+β−φJ2,Pr3=S▱AO1O2A1⋅γ1⋅V1cosα1+2β−φJ1,Pr4=SΔCO1O2⋅γ2⋅V1cosα1+2β−φJ2,Pr5=S▱A1O2O3A2⋅γ1⋅V2cosα2+3β−φJ1,Pr6=SΔCO2O3⋅γ2⋅V2cosα2+3β−φJ2,Pr7=S▱A2O3Oi+1Ai⋅γ1⋅Vicosαi+i+1β−φJ1,Pr8=SΔCO3Oi+1⋅γ2⋅Vicosαi+i+1β−φJ2,Pr9=SΔAiOi+1Oi+2⋅γ1⋅Vi+1cosαi+1+i+2β−φJ1,Pr10=S▱Oi+1CAi+1Oi+2⋅γ2⋅Vi+1cosαi+1+i+2β−φJ2,Pr11=SΔCAi+1Am⋅γ2⋅Vmcosαm+m+1β−φJ2,Pr12=∑e=mm+n−1SΔAeGAe+1⋅γ2⋅Ve+1⋅cosαe+1+e+1−mε−φJ2,Pr13=SΔGFAm+n−1⋅γ2⋅Vm+ncosαm+n+nε−φJ2,Pr14=SΔCRG⋅γ2⋅Vm+n+1.Half of the damage mode gravity power shown by the upper damage mode (i + 1) is(41)Pr=∑c=114Prc.(4)
Internal energy dissipation powerHalf of the energy dissipation power in the failure mode (i + 1) is(42)PW=CJ1⋅cosφJ1∑k=1iAAk⋅Vk+∑t=1iOtAt−1⋅Vt−1,t+BA⋅V0+AiOi+2⋅Vi+1+CJ2cosφJ2∑t=1i+1COt⋅Vt−1,t+∑p=i+1mCAp⋅Vp−1,p+Ai+1Oi+2⋅Vi+1+∑q=i+1mAqAq+1⋅Vq+1+∑e=m+1m+n−1Ae−1Ae⋅Ve+mm+n−1F⋅Vm+n+∑e=mm+n−1AeG⋅Ve,e+1+CR⋅Vm+n+1+CG⋅Vm+n,m+n+1.(5)
Support reaction powerThe undetermined parameterk 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.(43)PF=W⋅q⋅Vm+n+1−T⋅k⋅q⋅Vm+n⋅sinαm+n+φJ2.(6)
Calculation of surrounding rock pressureVirtual power equation:(44)Pr+PF=PW⇔∫Ω∗W∗υ∗dΩ∗+∫A∗T∗υ∗dA∗=∫Ω∗σij∗ε˙ij∗dΩ∗+∫S∗τJ−σnJtanφJΔvt∗dS∗,V0,1=sinα1−4β+θ3+φ1−π/2sinπ−α1−β−2φ1V0.Roof vertical support reaction:(45)q=Pr−PWW⋅V0+T⋅k⋅sinαm+n+φJ2⋅Vm+n.In the above formula,q is a function of the variable H0, which is(46)q=fH0.It is known from the upper limit theorem that the vertical support reaction force of the top plate is the minimum value ofq = f (H0). When q = f (H0) takes the minimum value, the following conditions should be met:(47)dqdH0=0.Through equation (47), the value of H0 when q takes the extreme value can be obtained. It is also necessary to verify whether the obtained H0 satisfies the constraint condition corresponding to the failure mode:(48)π−αi+1+θ2>0,π−β−2φJ2>0,3π2−αs−β−3φJ1>0,s=1…i,3π2−αt−β−3φJ2>0,s=i+1…m,3π2−αu−ε−3φJ2>0,u=m+1…m+n,αk−m−kβ+θ3+2φJ1−π>0,k=1…i,αl−m−lβ+θ3+2φJ2−π>0,i=i+1…m,αe−e−mε+2φJ2−π>0,e=m+1…m+n−1.If the above constraint condition is satisfied, the value ofH0 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:(49)e=kq.


## 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.
Rock formationBulk weightModulus of deformationCohesionInternal friction angleHorizontal bed coefficientVertical bed coefficientPoisson’s ratio(kN/m3)(GPa)(MPa)(o)(MPa/m)(MPa/m)Strongly weathered granite18.51.30.2222301800.33Moderately weathered granite21.06.50.4315505000.38


#### 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.
Surrounding rock categoryElastic modulusES (GPa)Poisson’s ratioLateral pressure coefficientCohesion (KPa)Internal friction angle (o)Bulk weightγ (kN/m3)Strongly weathered carbonaceous slate0.80.401.0552017.5Moderately weathered carbonaceous slate3.00.311.01353023.5


### 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 Table3.Table 3
The relaxation pressures of Jiangxi Road metro station.
Calculation methodq/kPae1/kPae2/kPae0/kPae0/qRemarkWeighted Terzaghi’s formula181.8090.12175.76135.960.75q is the vertical load of the tunnel; e1 is the load at the top of the sidewall; e2 is the load at the bottom of the sidewall; e0 is the mean value of the top and bottom loads of the sidewall, e0=(e1 + e2)/2。Modified Terzaghi’s formula162.8498.41211.23130.630.80Weighted Bierbäumer’s formula198.03101.26186.66149.650.76Modified Bierbäumer’s formula168.37119.53233.81143.950.85Weighted Jiaxiao Xie’s formula174.6372.51157.68123.570.71Modified Jiaxiao Xie’s formula153.1976.32189.40114.760.75Field measured data135.6578.12192.01106.890.79Upper-bound theorem of limit analysisk = 0.40269.41107.760.40The upper-bound theorem of limit analysis assumes that the sidewall load is a uniform loade1 = e2 = e0, where k is the ratio of the horizontal support reaction force to the vertical support reaction force, k = e0/q.k = 0.60217.82130.690.60k = 0.80179.16143.330.80k = 1.00159.67159.671.00k = 1.20146.31175.571.20k = 1.40138.20193.481.40It is known from Table3 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.
Calculation methodq/kPae1/kPae2/kPae0/kPae0/qRemarkPlatts pressure arch theory240.2370.19101.6985.940.36q is the vertical load of the tunnel; e1 is the load at the top of the sidewall; e2 is the load at the bottom of the sidewall; e0 is the mean value of the top and bottom loads of the side wall, e0 = (e1 + e2)/2。Modified Platts pressure arch theory225.4682.34113.5897.960.43Railway specification method301.9190.57150.96120.760.40Fraldi limit analysis method310.25115.28160.32137.800.44Field measured data200.1066.3895.2080.790.40Upper-bound theorem of limit analysisk = 0.40368.97147.590.40The upper-bound theorem of limit analysis assumes that the sidewall load is a uniform loade1 = e2 = e0, where k is the ratio of the horizontal support reaction force to the vertical support reaction force, k = e0/q。k = 0.60323.64194.180.60k = 0.80284.36227.490.80k = 1.00248.79248.791.00k = 1.20220.56264.671.20k = 1.40204.68286.551.40Table 5
The relaxation pressures of Zhifang tunnel when the depth of the tunnel is 75 m.
Calculation methodq/kPae1/kPae2/kPae0/kPae0/qRemarkPlatts pressure arch theory240.2370.19101.6985.940.36q is the vertical load of the tunnel; e1 is the load at the top of the sidewall; e2 is the load at the bottom of the sidewall; e0 is the mean value of the top and bottom loads of the sidewall, e0 = (e1 + e2)/2。Modified Platts pressure arch theory225.4682.34113.5897.960.43Railway specification method301.9190.57150.96120.760.40Fraldi limit analysis method310.25115.28160.32137.800.44Field measured data210.5668.3898.2083.290.40Upper-bound theorem of limit analysisk = 0.40379.56151.820.40The upper-bound theorem of limit analysis assumes that the sidewall load is a uniform loade1 = e2 = e0, where k is the ratio of the horizontal support reaction force to the vertical support reaction force, k = e0/q。k = 0.60333.01199.810.60k = 0.80294.76235.810.80k = 1.00259.87259.871.00k = 1.20231.51277.811.20k = 1.40215.34301.481.40


## 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.
Rock formationBulk weightModulus of deformationCohesionInternal friction angleHorizontal bed coefficientVertical bed coefficientPoisson’s ratio(kN/m3)(GPa)(MPa)(o)(MPa/m)(MPa/m)Strongly weathered granite18.51.30.2222301800.33Moderately weathered granite21.06.50.4315505000.38


### 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.
Surrounding rock categoryElastic modulusES (GPa)Poisson’s ratioLateral pressure coefficientCohesion (KPa)Internal friction angle (o)Bulk weightγ (kN/m3)Strongly weathered carbonaceous slate0.80.401.0552017.5Moderately weathered carbonaceous slate3.00.311.01353023.5


## 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.
Rock formationBulk weightModulus of deformationCohesionInternal friction angleHorizontal bed coefficientVertical bed coefficientPoisson’s ratio(kN/m3)(GPa)(MPa)(o)(MPa/m)(MPa/m)Strongly weathered granite18.51.30.2222301800.33Moderately weathered granite21.06.50.4315505000.38


## 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.
Surrounding rock categoryElastic modulusES (GPa)Poisson’s ratioLateral pressure coefficientCohesion (KPa)Internal friction angle (o)Bulk weightγ (kN/m3)Strongly weathered carbonaceous slate0.80.401.0552017.5Moderately weathered carbonaceous slate3.00.311.01353023.5


## 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 Table3.Table 3
The relaxation pressures of Jiangxi Road metro station.
Calculation methodq/kPae1/kPae2/kPae0/kPae0/qRemarkWeighted Terzaghi’s formula181.8090.12175.76135.960.75q is the vertical load of the tunnel; e1 is the load at the top of the sidewall; e2 is the load at the bottom of the sidewall; e0 is the mean value of the top and bottom loads of the sidewall, e0=(e1 + e2)/2。Modified Terzaghi’s formula162.8498.41211.23130.630.80Weighted Bierbäumer’s formula198.03101.26186.66149.650.76Modified Bierbäumer’s formula168.37119.53233.81143.950.85Weighted Jiaxiao Xie’s formula174.6372.51157.68123.570.71Modified Jiaxiao Xie’s formula153.1976.32189.40114.760.75Field measured data135.6578.12192.01106.890.79Upper-bound theorem of limit analysisk = 0.40269.41107.760.40The upper-bound theorem of limit analysis assumes that the sidewall load is a uniform loade1 = e2 = e0, where k is the ratio of the horizontal support reaction force to the vertical support reaction force, k = e0/q.k = 0.60217.82130.690.60k = 0.80179.16143.330.80k = 1.00159.67159.671.00k = 1.20146.31175.571.20k = 1.40138.20193.481.40It is known from Table3 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.
Calculation methodq/kPae1/kPae2/kPae0/kPae0/qRemarkPlatts pressure arch theory240.2370.19101.6985.940.36q is the vertical load of the tunnel; e1 is the load at the top of the sidewall; e2 is the load at the bottom of the sidewall; e0 is the mean value of the top and bottom loads of the side wall, e0 = (e1 + e2)/2。Modified Platts pressure arch theory225.4682.34113.5897.960.43Railway specification method301.9190.57150.96120.760.40Fraldi limit analysis method310.25115.28160.32137.800.44Field measured data200.1066.3895.2080.790.40Upper-bound theorem of limit analysisk = 0.40368.97147.590.40The upper-bound theorem of limit analysis assumes that the sidewall load is a uniform loade1 = e2 = e0, where k is the ratio of the horizontal support reaction force to the vertical support reaction force, k = e0/q。k = 0.60323.64194.180.60k = 0.80284.36227.490.80k = 1.00248.79248.791.00k = 1.20220.56264.671.20k = 1.40204.68286.551.40Table 5
The relaxation pressures of Zhifang tunnel when the depth of the tunnel is 75 m.
Calculation methodq/kPae1/kPae2/kPae0/kPae0/qRemarkPlatts pressure arch theory240.2370.19101.6985.940.36q is the vertical load of the tunnel; e1 is the load at the top of the sidewall; e2 is the load at the bottom of the sidewall; e0 is the mean value of the top and bottom loads of the sidewall, e0 = (e1 + e2)/2。Modified Platts pressure arch theory225.4682.34113.5897.960.43Railway specification method301.9190.57150.96120.760.40Fraldi limit analysis method310.25115.28160.32137.800.44Field measured data210.5668.3898.2083.290.40Upper-bound theorem of limit analysisk = 0.40379.56151.820.40The upper-bound theorem of limit analysis assumes that the sidewall load is a uniform loade1 = e2 = e0, where k is the ratio of the horizontal support reaction force to the vertical support reaction force, k = e0/q。k = 0.60333.01199.810.60k = 0.80294.76235.810.80k = 1.00259.87259.871.00k = 1.20231.51277.811.20k = 1.40215.34301.481.40


## 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 ofk, 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



---

*Source: 5096931-2020-12-30.xml*