Numerical Simulation Analysis Method of the Surrounding Rock and Support Bearing Capacity in Underground Cavern

Abstract

After the excavation of an underground cavern, how the surrounding rock and the support work together to bear the excavation load is an important prerequisite to correctly analyze the joint force characteristics; effectively play the role of support; and ensure the safety, efficiency, and economy of underground cavern construction. Starting from the elastic-plastic load release characteristics of surrounding rock, this paper proposes a calculation method of the elastic load coefficient of surrounding rock and a graded release method of plastic load, which ensures the actual effect of the synergistic action of the first support and surrounding rock. Based on the elastic-plastic damage evolution characteristics of surrounding rock, a weighted iterative calculation method of elastic-plastic damage is proposed, and an evaluation method of load release ultimate bearing capacity of surrounding rock is determined. By monitoring the change law of rock acoustic wave velocity with surrounding rock damage, the relationship between the wave velocity and the damage coefficient of the surrounding rock in the excavation process is deduced, and it is proposed to determine the latest support time for first support by using the measured rock damage wave velocity. Through the numerical simulation analysis of a diversion tunnel excavation and support, the damage evolution law of the surrounding rock with the release of the excavation load is studied. The ultimate bearing capacity of various surrounding rocks and supporting opportunity is determined. The results demonstrate the validity and practicality of the analysis and calculation methods in this paper, which provide a new idea and analysis method for quantifying the bearing capacity of surrounding rock and determining the support timing in the excavation and support design of underground caverns.

1. Introduction

A reasonable analysis of the surrounding rock and support collaborative bearing capacity of underground caverns is the basic premise to ensure the safety, efficiency, and economy of underground engineering design and construction. Since the underground cavern is an underground space excavated in a complex geological environment, the collaborative bearing capacity of the surrounding rock and support is affected by many environmental factors [1,2]. How to analyze the collaborative bearing capacity of surrounding rock and support after excavation disturbance is the focus and difficulty of academic and engineering communities [3].

The current analysis method of the surrounding rock and support bearing capacity by scholars can be divided into two categories: One is the quantitative analysis of the limited shape of the underground cavern (such as a circular tunnel) or limited boundary conditions (such as only considering the action in the plane range) [4,5]. It is mainly based on the load structure method and is established by the calculation method of resistance coefficient K of the surrounding rock of the tunnel through various theoretical derivations. Then, K is modified according to the bearing capacity of the rock, and the time-varying bearing coefficient is proposed to evaluate the bearing capacity of the surrounding rock when it is self-stabilized. The other category is qualitative analysis based on numerical methods and model tests [6,7,8,9]. Qualitative analysis is mainly based on basic theories of the joint bearing of surrounding rock and support. According to numerical methods and model tests, the damage characteristics of surrounding rock and the stress of support are discussed, and the bearing capacity of the surrounding rock is qualitatively evaluated based on the stress state, surrounding rock deformation control, and the failure range of rock. For example, Zhang et al. [10] analyzed the load transfer mechanism of the tunnel pressure arch and the progressive damage law of the surrounding rock based on the composite structure of the surrounding rock. Yu et al. [11] revealed the time-varying evolution law of the bearing capacity of the surrounding rock. In terms of mechanical characteristics of supporting structures, Liu et al. [12] put forward a finite element stress analysis method for the fully grouted rock bolts, obtained the stress distribution along the length of the bolt, and could reasonably simulate the bond and slip state of the bolt. Shi et al. [13] used PFC2D to explore the meso interaction mechanism between the surrounding rock and bolt, simulated the crack expansion process of the bolt in a pullout test, and revealed the failure mechanism of the bolt. Li et al. [14] proposed a constitutive model based on the Mohr–Coulomb criterion to describe the mechanical behavior of the anchor cable pullout, which reasonably described the failure mechanism between the anchor cable and grouting ring.

To sum up, the research results of domestic and foreign scholars have effectively promoted the development of the research on the stability of surrounding rock in underground engineering. However, the current quantitative evaluation method is limited to some specific conditions, and its application is greatly restricted. The qualitative analysis is mainly based on numerical analysis and test methods, and comprehensive analysis and evaluation of the failure zone, displacement, stress, and strain are carried out. It is limited by theoretical methods and experimental conditions, which exist in large human evaluation factors, so its application is greatly influenced by the level of human cognition. In this paper, the collaborative bearing capacity of the surrounding rock and support is researched according to the law of load release of the surrounding rock and the evolution law of rock damage in underground cavern excavation. The numerical analysis and evaluation method of the ultimate bearing capacity of the surrounding rock and the idea of determining the latest time of application of first support are proposed based on the principle of the rock acoustic test method.

Compared with the traditional method, the method in this paper divides the excavation release load into the elastic load and plastic load and releases the plastic load in a graded manner, which can reasonably reflect the damage evolution process of rock during tunnel excavation. The quantifiable acoustic test value is used to measure the damage degree of surrounding rock and determine the latest application time of first support, which can objectively evaluate the bearing capacity of surrounding rock, facilitate engineering monitoring, and determine the support opportunity. It provides an effective analysis and judgment method for the collaborative bearing of the first support and surrounding rock, determining the latest support time after the excavation of the underground cavern.

2. Determination of Load Release of Surrounding Rock in Underground Cavern

After the excavation of the underground cavern, the total load provided to the surrounding rock by the initial geo-stress and support is composed of three parts: the released engineering force due to the deformation of the surrounding rock, the surrounding rock resistance force borne by its strength, and the externally applied anchor support force [15], giving,

$$PT=PD+PR+PS$$

(1)

where PT is total load, PD is engineering force, PR is resistance force, and PS is support force.

Underground engineering construction generally adopts the principle of New Austrian Tunnelling Method (NATM) [16,17], and the first support should be applied promptly after excavation. As the rock is an elastic-plastic medium affected by excavation disturbance, a part of the elastic deformation load, which is difficult to monitor, is quickly released before the PS has been applied. The plastic deformation load is gradually released with the different characteristics of the surrounding rock. Thus, how much of the surrounding rock release load can actually be borne by the first support is crucial to a reasonable support design.

Generally speaking, the joint bearing effect of surrounding rock and support only works in the process of plastic load release. Therefore, if the elastic load in the excavation is determined, the released engineering force PD due to the lagging support during the deformation process of the surrounding rock can be determined. After the PD is determined, the load jointly borne by the anchorage support and the surrounding rock can be reasonably calculated.

2.1. Elastic Load Coefficient of Surrounding Rock

During the excavation process, the total released load R of the surrounding rock can be decomposed into two parts: the elastic load R_e and the plastic load R_p (Figure 1). Elastic load refers to the load bearing when the surrounding rock just reaches the yield surface and will generally be released before the support is applied, while plastic load refers to the load borne after the surrounding rock exceeds the yield surface, which is jointly borne by the surrounding rock and the support. The elastic load coefficient of the surrounding rock refers to the load distribution ratio coefficient that shows that the stress state of the surrounding rock just reaches the yield surface in the entire excavation release load.

It is assumed that the surrounding rock stress increment generated by the total release load R is $\left\{\mathbf{\Delta }\mathit{\sigma }\right\}$. If the rock element transitions from the elastic stress state (point A in Figure 2) to the plastic stress state (point B in Figure 2) during loading, there exists a critical stress state (point C in Figure 2) corresponding to the elastic release load R_e, which makes the yield function F equal to zero, giving [18],

$$F(\{{\mathit{\sigma }}_c\})=F(\{{\mathit{\sigma }}_0\}+S\{\mathbf{\Delta }\mathit{\sigma }\})=0$$

(2)

where S is the elastic load coefficient of the rock element.

As the yield function $F(\{{\mathit{\sigma }}_c\})$ is generally an expression of a complex stress state, it is usually impossible to solve the S according to Equation (2). In the finite element numerical calculation, the following equation can be used to approach point C step by step from the initial stress state (point A) and the fully loaded stress state (point B) for iterative solution:

$$S_{i+1}=S_i-(S_i-S_j)\times F(\{{\mathit{\sigma }}_i\})/[F(\left\{{\mathit{\sigma }}_i\right\})-F(\left\{{\mathit{\sigma }}_j\right\})]$$

(3)

When i = 1, take $S_i$ = 1, corresponding to the yield function taken as $F(\{{\mathit{\sigma }}_1\})=F(\{{\mathit{\sigma }}_0\}+\{\mathbf{\Delta }\mathit{\sigma }\})$; $S_j$ = 0, corresponding to the yield function taken as $F(\{{\mathit{\sigma }}_j\})=F(\{{\mathit{\sigma }}_0\})$.

When i > 1, in the subsequent iterations, according to the $S_{i+1}$ and the corresponding yield function $F_{i+1}(\left\{{\mathit{\sigma }}_0\right\}+S_{i+1}\left\{\mathbf{\Delta }\mathit{\sigma }\right\})$ calculated by Equation (3), the value opposite to F_i+1 in the previous step of F_i and F_j is selected as the next iteration F_j and corresponding $S_j$, and then substituted into Equation (3) for iterative calculation. The coefficient value $S_{i+1}$ calculated by several iterations can quickly satisfy Equation (2), and its value is taken as the elastic load coefficient of the element.

2.2. Release Load Coefficient Jointly Borne by the Support and Surrounding Rock

Since the elastic property of the rock is very small, the elastic load will be quickly released after the cavern is excavated. If the support can be applied in time, the released engineering force PD can be assumed to be an elastic load. Therefore, the elastic release load R_e and the plastic release load R_p jointly borne by the support and the surrounding rock can be calculated by the total release load R according to Equation (4).

$$\{\begin{array}{l}{R}_e=SR\\ R_p=(1-S)R=\beta R\end{array}$$

(4)

where $\beta =1-S$ is the release load coefficient jointly borne by the support and surrounding rock.

Due to the spatial structure, storage environment, excavation method, and other factors of the underground cavern, the released load coefficients of each part are different, so the released load borne by the support is also different. In order to reasonably simulate the release load borne by the anchorage support, the following method can be used to obtain ${\beta }_i$ of each node after the unit excavation in 3D finite element analysis.

Firstly, according to Equation (2), the elastic coefficient $S_j$ of each unit j is calculated, and then the elastic coefficient $S_i$ of node i is calculated according to the average value method around the node element. If the node i has n units associated with it, then $S_i$ can be expressed as

$$S_i={\displaystyle \sum _{k=1}^n\frac{S_{jk}}{n}}$$

(5)

Then, according to the application of support on surrounding rock, the release load coefficient of jointly bearing at node i can be calculated as follow:

$${\beta }_i=1-S_i$$

(6)

Substitute Equation (6) into Equation (4) to calculate the joint load of anchorage support and surrounding rock at each node. This method can reasonably respond to the reinforcing effect of anchorage support application timing on surrounding rock. At the same time, it effectively reflects the spatial effect that supports in different positions bearing different release loads according to the actual deformation of surrounding rock.

3. Determination of Ultimate Bearing Capacity of Surrounding Rock in Underground Cavern

The resistance force PR provided by the strength of the surrounding rock is the ultimate bearing capacity without considering the external support force. Examining the ultimate bearing capacity of the surrounding rock is a prerequisite for the reasonable design of anchorage support parameters and determination of the latest time to apply support. This paper proposes an iterative analysis method of damage evolution for underground cavern excavation according to the damage evolution characteristics of the surrounding rock and proposes the criteria for discriminating the local damage of the surrounding rock according to its local damage characteristics.

3.1. Damage Evolution Calculation Method of Surrounding Rock Load Release

The elastic load can be released at one time in numerical calculation, and the plastic load can be applied in a weighted and graded manner after the anchorage support is applied. That is, when there are few surrounding rocks entering plastic failure in the early stage, the applied plastic release load is larger. Then, the graded load decreases gradually when the plastic failure zone gradually increases. The iterative release process is shown in Figure 1. The nonlinear damage constitutive relation of rock can be calculated according to the damage coefficient D, and can be expressed as [19]:

$${({\mathit{\sigma }}_{ij})}_D=(1-D){\mathit{\sigma }}_{ij}+\frac{D}{3}{\mathit{\sigma }}_{kk}{\mathit{\delta }}_{ij}$$

(7)

where ${\mathit{\delta }}_{ij}$ is the Kronecker delta function, ${\mathit{\sigma }}_{ij}$ is stress conforming to the elastic-plastic constitutive relation, and D is the damaged internal variable describing the propagation of rock microfractures. It can be calculated by Equation (8) according to the plastic strain partial tensor after the rock enters plastic failure [20].

$$D={\mathrm{D}}_{\mathrm{n}}[1-\exp (-\mathrm{K}{\mathit{\xi }}^{\mathrm{a}})]$$

(8)

where D_n, K and a are damage material constants, which can be determined according to material tests. $\mathit{\xi }=\sqrt{{\mathit{e}}_{\mathit{i}\mathit{j}}·{\mathit{e}}_{\mathit{i}\mathit{j}}}$, where ${\mathit{e}}_{\mathit{i}\mathit{j}}$ is the plastic strain partial tensor.

The plastic damage stress matrix [H_D] can be calculated according to the rock damage coefficient in Equation (8) in each level of the plastic load iteration, giving,

$$[{\mathit{H}}_D]=(1-D+D{\mathit{\delta }}_{ij}/3)[{\mathit{H}}_{\mathit{p}}]+(D-D{\mathit{\delta }}_{ij}/3)[{\mathit{H}}_{\mathit{e}}]$$

(9)

where $[{\mathit{H}}_p]$ and $[{\mathit{H}}_{\mathit{e}}]$ are plastic and elastic stress matrix of the rock, respectively.

The damage stiffness matrix [K_D] of the plastic damage element is modified by the following equation,

$$[{\mathit{K}}_D]=\int {[\mathit{B}]}^{\mathrm{T}}[{\mathit{H}}_D][\mathit{B}]dv$$

(10)

The plastic stiffness can be kept constant during the iterative process in each level of plastic load increment $\{\mathbf{\Delta }{R}_{P_i}\}$, and the displacement increment generated by the plastic load at level i can be iterated according to Equation (11).

$$[{\mathit{K}}_{\mathit{e}}]{\left\{\mathbf{\Delta }{\mathit{\delta }}_i\right\}}_j=\left\{\mathbf{\Delta }{R}_i\right\}+[{\mathit{K}}_{\mathit{D}}]{\left\{\mathbf{\Delta }{\mathit{\delta }}_i\right\}}_{j-1}$$

(11)

The damage iteration by means of graded weighting can reasonably simulate the magnitude of the release load borne by the anchoring support when it is applied and has advantages of stable iteration and fast calculation speed. At the same time, it can reasonably reflect the development process of failure evolution with the joint action of the anchoring support and the surrounding rock.

3.2. Evaluation Method for Ultimate Bearing Capacity of Surrounding Rock Based on Excavation Load Release

In actual engineering construction, due to the influence of construction conditions, equipment, organization, and other factors in the excavation of underground caverns, the anchorage support cannot be applied in time after excavation. Therefore, in addition to the elastic load, before the support is applied, the release load also includes a part of the plastic load. At this time, the release load coefficient independently borne by the surrounding rock can be expressed as

$${\alpha }_i=S_i+b{\beta }_i=S_i+b(1-S_i)$$

(12)

where b is a ratio coefficient of plastic load release before support is applied.

By continuously adjusting the b, the release load before support applied $R_b=(S_i+b{\beta }_i)R$ is applied to the surrounding rock in stages, and the damage evolution and failure process of the surrounding rock are analyzed. In the elastic-plastic numerical calculation, after the surrounding rock enters into the plastic failure zone, it will produce irreversible plastic deformation, but still has a certain bearing capacity. When the ratio of tensile strain to ultimate tensile strain $[{\mathit{\epsilon }}_a]$ of the surrounding rock is greater than 1, it is considered that the surrounding rock has reached the failure standard, and the corresponding load is the ultimate load of the surrounding rock.

$$\frac{{\left\{{\mathit{\epsilon }}_c\right\}}_i}{[{\mathit{\epsilon }}_a]}>1$$

(13)

The elastic-plastic damage iterative calculation method can reasonably analyze the damage evolution process of the surrounding rock. Continuously adjust the release load ratio coefficient of the surrounding rock by Equation (12), and when the release load makes the surrounding rock appear in the damage criteria in Equation (13), it can be considered that the release coefficient at this time is the release coefficient corresponding to the ultimate bearing capacity of the surrounding rock. Using this load release coefficient search method, the ultimate bearing capacity of the surrounding rock can be evaluated under various conditions.

4. Determination of the Latest Application Time of First Support

According to the theory of NATM, after the underground cavern is excavated, the first support should be applied in time to effectively play the joint effect of surrounding rock and support. However, in the actual project by the construction conditions or some objective reasons, the first support cannot be applied in time. Therefore, how one determines the latest time to apply the first support is the key to ensure the safety of the project and make a reasonable construction organization plan. In this paper, by studying the methods of determining the ultimate bearing capacity and load release coefficient of the surrounding rock, the correlation between the damage coefficient and the acoustic wave velocity of the rock is used, and the method of determining the latest time of applying the first support is proposed.

4.1. Relationship between Rock Acoustic Wave Velocity and Damage Coefficient

According to the propagation principle of elastic medium in infinite three-dimensional space, the velocity of longitudinal wave V_p and transversal wave V_s of rock can be expressed as [21]:

$$\{\begin{array}{l}{V}_p=\sqrt{\frac{E(1-\mu )}{\rho (1+\mu )(1-2\mu )}}\\ V_s=\sqrt{\frac{E}{2\rho (1+\mu )}}\end{array}$$

(14)

where E is the elastic modulus of the rock, $\mu$ is the Poisson’s ratio, and $\rho$ is the density.

In the excavation process with the release of load, the surrounding rock will suffer plastic damage. According to the elastic-plastic damage iterative calculation method, when the surrounding rock enters the plastic state, its damage degree can be described by damage coefficient D. The elastic modulus after damage is calculated by Equation (15).

$$\tilde{E}=E(1-D)$$

(15)

Ignoring the density change during the damage process of surrounding rock, and substituting Equation (15) into Equation (14), after the damage of the surrounding rock and excavation load is released, the acoustic wave velocity of the rock can be expressed as

$$\{\begin{array}{l}\widetilde{V_p}=\sqrt{\frac{E(1-\mu )(1-D)}{\rho (1+\mu )(1-2\mu )}}\\ \widetilde{V_s}=\sqrt{\frac{E(1-D)}{2\rho (1+\mu )}}\end{array}$$

(16)

Using the above elastic-plastic damage iterative evolution calculation method, the ultimate bearing capacity and failure state with the release of the surrounding rock load during the excavation are determined by Equation (13). The damage coefficient D, which corresponds to ultimate bearing capacity of the surrounding rock, is calculated according to Equation (8). Substitute D into Equation (16) to calculate the damage wave velocity, which is the rock acoustic wave velocity corresponding to the ultimate bearing capacity of the surrounding rock. Comparing the acoustic wave velocity before and after excavation, the relative reduction coefficient of acoustic wave of damaged rock after excavation can be obtained, giving

$$\tilde{\gamma }=\widetilde{V_p}/\widetilde{V_p}=\sqrt{1-D}$$

(17)

4.2. Determination of Latest Application Time of First Support According to Rock Acoustic Wave Velocity

The acoustic wave velocity of surrounding rock refers to the propagation velocity of elastic waves in the surrounding rock. Because of the damage of the rock elastic modulus during the excavation process, its wave velocity also decreases continuously [22]. The wave velocity after the surrounding rock damage can be measured by acoustic testing technology. It mainly includes the single-hole method and cross-hole method [23]. The cross-hole method is more complex, so in engineering, the earth’s surface excitation and in-hole reception method (single-hole method) is often used to measure the surrounding rock acoustic waves [24]. This principle is shown in Figure 3. The seismic source on the earth’s surface generates elastic waves, which are received by the geophone in the hole, and the acoustic wave velocity of the rock can be calculated by using the time difference between the elastic waves reaching different depths. After the cavern is excavated, the longitudinal and transverse waves of the rock can be calculated by Equation (18).

$$\{\begin{array}{l}\widetilde{V_{ip}}=\frac{h_{i+1}-h_i}{t_{i+1}^p-t_i^p}\\ \widetilde{V_{is}}=\frac{\sqrt{h_{i+1}{}^2+l^2}-\sqrt{h_i{}^2+l^2}}{t_{i+1}-t_i}\end{array}$$

(18)

where $h_i$ and $h_{i+1}$ are the vertical distance from the measuring point i and i + 1 to the orifice, respectively; l is the horizontal distance from the seismic source to the orifice; $t_i$ and $t_{i+1}$ are the time of the transversal wave reaching the geophone; and $t_i^p$ and $t_{i+1}^p$ are the time of longitudinal wave reaching the geophone.

The acoustic wave velocity of rock, which is continuously decreasing in the excavation process, can be measured. When the ratio of the damage acoustic wave velocity obtained from Equation (18) to the original rock wave velocity is close to the value calculated in Equation (17), this is the latest time that the first support must be applied.

5. Engineering Application and Analysis

5.1. Calculation Model and Parameters

The diversion tunnel of a project is buried at 220 m below the surface. It is excavated by drilling and blasting method. The excavation section is circular, and the tunnel diameter is 7.92 m. The surrounding rock type in the project area is mainly III~IV, and the part of the tunnel crossing the fw8-1 fault is selected for analysis and calculation. The fault width is 15 m, and the physical and mechanical parameters of surrounding rock are considered as type V. The three-dimensional finite element model of this section tunnel and the analysis model of excavation and support of each typical section are shown in Figure 4. As seen in the figure, it is completely discretized by an eight-node hexahedron, consisting of 82,491 nodes and 77,100 elements. The x-axis of the model is perpendicular to the tunnel axis along horizontal direction. The y-axis coincides with the tunnel axis, with positive direction going along the water flow. The z-axis coincides with the geodetic coordinate. The dimensions of the model along the three directions are 154.4 m, 200.0 m, and 254.3 m, respectively.

For type III surrounding rock, anchoring support is carried out only within 120° of the top arch of the cavern. The support parameters are φ22, L = 2.5 m, and @1.2 × 1.2 m, with full-section sprayed C30 concrete, and the sprayed layer thickness is 0.1 m. The support parameters of type IV surrounding rock are φ22, L = 2.5 m, and @1.0 × 1.0 m, and the sprayed layer thickness is 0.15 m. The support parameters of the type V surrounding rock are φ25, L = 3.0 m, and @0.8 × 0.8 m, and the sprayed layer thickness is 0.15 m.

The three-dimensional initial geo-stress can be obtained by stress inversion. The initial maximum principal stress of surrounding rock is between −8.9 and −11.0 MPa, and the minimum principal stress is between −4.6 and −6.2 MPa. The physical and mechanical parameters of the materials are listed in

Table 1. Physical and mechanical parameters of materials.

Materials		Elastic Modulus/GPa	Poisson’s Ratio	Cohesion Force/MPa	Friction Angle/°	Tensile Strength/MPa	Compressive Strength/MPa	Density/g·cm−3
Rock	III	6.50	0.29	0.90	41.99	3.00	50.00	2.50
	IV	2.00	0.30	0.40	26.57	0.10	15.00	2.20
	V	0.80	0.35	0.10	16.70	0.08	12.00	2.00
Sprayed layer C30		30.0	0.167	2.0	46.0	1.43	14.3	2.50
Bolt		210	0.30	400

.

5.2. Calculation Results and Analysis

In order to study the stability characteristics and bearing capacity of surrounding rock in different parts of the section during the tunnel excavation and support, characteristic sections are set up in each type of rock. Each characteristic section is arranged with eight monitoring points, and the layout of monitoring points around the tunnel is shown in Figure 5.

5.2.1. Distribution of Elastic Load Release Coefficient

The load release calculation method can reasonably simulate the anchorage support effect. According to this method, the elastic load release coefficients of each monitoring point of the tunnel characteristic section in type III, IV, and V surrounding rocks can be iterated by Equation (3), as shown in

Table 2. Elastic load release coefficient of monitoring points around the tunnel (%).

Rock	1	2	3	4	5	6	7	8
III	93	94	96	94	93	94	96	94
IV	46	57	70	58	46	58	70	57
V	12	25	44	26	12	26	44	25

.

From the distribution law of elasticity coefficient around the tunnel, the haunch of the tunnel is large, and the top arch first enters the plastic, then the bottom arch enters, and the failure zone gradually develops to the middle parts. From the distribution of the elastic coefficient of each type of surrounding rock, the elastic load coefficient of type III is above 93%, and the load tends to be completely released with the completion of the excavation. The elastic load release coefficient of type IV surrounding rock is distributed in 46~70%, among which the top arch is about 50% and the haunch is about 70%. The elastic load release of type V surrounding rock is relatively small, where the haunch is about 40% and the top arch is only 12%. It can be seen that the section with better surrounding rock conditions has a larger proportion of elastic load release, because the better rock required a larger load to reach the plastic state.

After the elastic load is released, the joint bearing action of anchorage support and surrounding rock is applied. The stress distribution of the bolts for each type of surrounding rock are as follows: type III surrounding rock is 40~100 MPa, type IV is 150~250 MPa, and type V is 200~400 MPa (most of the bolts are in yield). The worse the surrounding rock type is, the greater the release load borne by its bolts is. The calculation result shows that the self-stabilization capacity of type III surrounding rock is better after excavation. The type IV surrounding rock can basically ensure the stability of the surrounding rock under anchorage support. The self-stabilization capacity of type V surrounding rock is poor, and the load of anchorage support is relatively large. In this situation, it is difficult for the first support to meet the requirements of surrounding rock stability, and the use of secondary support is necessary.

5.2.2. Analysis of the Surrounding Rock Ultimate Bearing Capacity

According to the calculation results of elastic load release coefficient, the load of type III surrounding rock is basically released completely, and the residual load can be carried by the surrounding rock itself, which has a good self-stabilization capacity. Therefore, this section only studies the ultimate bearing capacity of type V surrounding rock because its bearing capacity is the poorest.

The damage iterative calculation is carried out by constantly adjusting the value of the plastic load release ratio coefficient b by Equation (12). The graded loading is carried out on the sections of tunnels with surrounding rock type V to calculate and analyze the distribution of the tensile strain ratio of the surrounding rock under all levels of load. Studying the change law of failure zone, stress, and displacement of the surrounding rock with the load release coefficient is used to analyze the ultimate bearing capacity of the surrounding rock.

Ultimate Bearing Capacity Analysis of Surrounding Rock

The distribution of the tensile strain ratio of type Ⅴ surrounding rock under plastic loading at all levels is shown in Figure 6. As can be seen from the figure, with the gradual release of plastic loading, the value of the tensile strain ratio gradually increases and develops to deep rock. From its distribution, it is not uniform in the same section, and the tensile strain at the haunch is relatively small, but the depth and distribution range are larger. The tensile strain at the top and bottom of the tunnel has a larger value and a faster growth rate, which are the control parts of the bearing capacity of surrounding rock. When the plastic load is released to 30%, the first area with tensile strain ratio greater than 1 appears at the bottom of the tunnel, which means that the tensile stress here has exceeded the ultimate tensile strength of the surrounding rock. The surrounding rock starts to produce a disturbed zone, and the support should be applied in time.

• Development law of surrounding rock failure zone

In the process of plastic load release, the surrounding rock is gradually damaged, and the volume and depth of the failure zone are increasing. The development law of the surrounding rock failure zone under plastic load at all levels is shown in Figure 7. In general, the surrounding rock failure is mainly caused plastic and tensile failure, and the volume of all failure types are accumulated continuously. In terms of development law, several failure types are also basically the same, only differing in values. When the plastic load release coefficient is less than 0.3, the surrounding rock is mainly in the plastic failure zone and the growth rate is small. After the plastic load is released to 30%, the tensile failure zone starts to appear in the surrounding rock, and the volume of the failure zone also starts to grow rapidly. It can be considered that the released load has reached the maximum allowable bearing capacity of the surrounding rock at this time. Support should be applied in time to ensure the stability of the surrounding rock, which is relatively consistent with the analysis result of the tensile strain ratio in the previous section.

2.

Development law of surrounding rock stress

According to the development law of surrounding rock stress (Figure 8), with the release of the plastic load around the tunnel, the compressive stress at the top and bottom of the tunnel increases continuously, and the stress at the haunch decreases first and then increases slightly. This is because, with the excavation of the tunnel, the stress is adjusted and redistributed, and the load at the haunch is transferred to the top and bottom of the tunnel, resulting in the surrounding rock at the top arch and bottom reaching the bearing capacity first.

3.

Development law of displacement

With the release of plastic load around the tunnel, the displacement of the surrounding rock around the tunnel increases continuously (Figure 9), indicating that the plastic displacement is accumulated gradually in the process of load release. From the distribution of values, there is little difference between the top and bottom of the tunnel, and the displacement of the left and right haunch is basically the same, showing a symmetrical distribution. In the same characteristic section, the displacement at the haunch is larger than that at the top arch and bottom, and the maximum displacement is 98.23 mm after all the plastic load is released. It can also be seen from the contour map of tensile strain ratio that the distribution depth at the haunch is obviously larger. Therefore, the surrounding rock load can be released in the way that the deformation is acting.

In summary, through the analysis of the development law of the tensile strain ratio, the failure zone, stress, and displacement of the surrounding rock with the change of plastic load release coefficient can be determined. It can be seen that the tensile strain ratio of the surrounding rock starts to be greater than 1, and tensile failure occurs when the plastic load is released to 30%, the volume of the failure zone increases rapidly, and the plastic displacement accumulates continuously. Therefore, it can be considered that the maximum allowable bearing capacity of the surrounding rock is all of the elastic load and 30% of the plastic load (total release coefficient ${\alpha }_{\alpha }$ = 12% + 0.3 × 88% = 0.384), which means that the ultimate capacity of type V surrounding rock is only 38.4% of the total released load. It can be seen that type III surrounding rock has complete self-stabilization ability and type IV surrounding rock can bear most of excavation release load, while type V surrounding rock can only bear a small part of excavation release load and its self-stabilization ability is poor, so more external force support needs to be provided.

5.2.3. Analysis of Surrounding Rock Supporting Time

The damage coefficient corresponding to the ultimate load of the surrounding rock can effectively determine the latest time to apply the support. For type V surrounding rock, the distribution law of damage coefficient with plastic load release coefficient around the tunnel is shown in Figure 10. In the process of tunnel excavation unloading, the damage degree and damage coefficient of the surrounding rock increase continuously. The growth rate of damage degree at different parts is inconsistent. The damage degree at the haunch increases faster than that at the top and bottom, and finally tends to be stable. When the plastic load release coefficient b = 0.3, it can be seen that D at the top and haunch of the tunnel are 0.77 and 0.54, respectively.

The distribution law of acoustic wave velocity reduction coefficient of damaged rock under different plastic load release coefficients is calculated according to the corresponding damage coefficients, as shown in Figure 11.

After the plastic load is released, the surrounding rock starts to show plastic damage, which also means that $\tilde{\gamma }$ decreases, especially at the top and bottom of the tunnel, but the haunch is relatively small, while the depth of the reduction coefficient is smaller, and the overall range of $\tilde{\gamma }$ change is “X”. With the further release of plastic load, the descent range develops deeper into the surrounding rock, and the descending range continues to expand. The acoustic wave velocity ratio of the surrounding rock at the haunch decreases faster. After the load release is completed, the acoustic wave velocity ratio of the surrounding rock at the haunch is already close to the top and bottom of the tunnel, which is about 0.2, indicating that the wave velocity of the surrounding rock at this time has already decreased to 20% before excavation. When b = 0.3, the acoustic wave velocity reduction coefficients at the top and haunch of the tunnel are 0.479 and 0.68, respectively. It indicates that the anchorage support must be applied when the acoustic wave velocity reduction coefficient of damaged rock at the top and haunch of the tunnel reaches 47.9% and 68%, respectively. So, as long as the acoustic wave velocity is monitored around the tunnel during construction, the acoustic wave reduction coefficient is calculated, and the latest time of anchor support application at this point can be determined according to the distribution law in Figure 11.

6. Discussion

Theoretical analysis and various numerical simulation methods have been widely used in the study of underground cavern surrounding rock stability [5,25]. Since the theoretical analysis is mainly aimed at circular tunnels to obtain accurate analytical solutions, and the non-homogeneity and discontinuity of the rock cannot be considered, this paper adopts a numerical simulation method to study the bearing capacity of surrounding rock and support and determine the support timing. In the excavation process of cavern construction, the load release coefficient method has been commonly used to control the excavation process of caverns [26], but a uniform load release coefficient is used in most numerical simulations. The value of this coefficient relies mostly on experience, and the calculation method is still unclear, without considering the space-time effect of load release at the excavation surface. The research of some scholars shows that the load release process at different locations on the same section varies considerably [27]. Therefore, this paper proposes a calculation method of the elastic load coefficient of surrounding rock and a graded release method of plastic load based on the excavation load release characteristics, which reasonably reflect the space-time effect in the process of surrounding rock load release. In terms of the bearing capacity of the surrounding rock, this paper obtains the proportion of the excavation load to be borne when the bearing capacity of the surrounding rock is exerted to the maximum according to the damage characteristics and failure criterion of the rock. It provides a new analytical calculation idea for the load distribution problem in the joint bearing of the surrounding rock and support. Combined with the rock acoustic wave velocity testing technology, the relationship between the wave velocity and damage coefficient of the surrounding rock during construction excavation is derived, and the method of using the measured damaged rock acoustic wave velocity to determine the latest support timing applied for the first support is proposed. In general, the method in this paper can reasonably reflect the actual release characteristics of loads in cavern excavation, quantify the proportion of loads borne by the surrounding rock and support, and obtain the latest application time of the first support intuitively and conveniently, which has certain guiding significance for the design and construction of underground caverns.

However, some problems are still found in the research worthy of further discussion: (1) In this paper, the ultimate tensile strain criterion is used to judge the bearing capacity of the surrounding rock, but it is often difficult to obtain the ultimate tensile strain, and the applicability of this method in other projects needs to be further explored. (2) The calculation of rock acoustic wave velocity only considers the influence of elastic modulus and density and ignores the influence of Poisson’s ratio. In fact, Poisson’s ratio also changes during the load release process, which in turn affects the surrounding rock acoustic wave velocity. The above two points are also the direction and focus of future research.

7. Conclusions

According to the elastic-plastic damage evolution characteristics of surrounding rock excavation unloading, combined with 3D numerical calculation and acoustic monitoring methods, based on the analysis theory of the bearing capacity of the surrounding rock and support of the underground cavern, and through the simulation analysis of the underground cavern excavation, the following innovative results were obtained.

(1)

Based on the characteristics of surrounding rock load release in the excavation, the method of determining the elastic load according to the time when the surrounding rock enters the critical stress state is proposed. The method can reasonably calculate the released load jointly borne by the first support and surrounding rock; correctly reflect the applied effect of anchorage support and the temporal-spatial stress distribution characteristics after excavation of the cavern; and provide a theoretical basis for safe, efficient, and economic design of anchorage support parameters.

(2)

Based on the elastic-plastic damage theory of rock, the damage evolution analysis method of underground cavern excavation with variable rock element damage stiffness is proposed. It has fast iterative convergence speed and can effectively reflect the gradual damage evolution process of the surrounding rock in the underground cavern excavation process, which provides an efficient analysis method for reasonably determining the ultimate bearing capacity of different surrounding rocks and judging the stability state of them in various conditions.

(3)

Based on the damage characteristics of the surrounding rock and the propagation law of acoustic wave in the surrounding rock, the relationship between the damage coefficient and the acoustic wave velocity reduction coefficient of damaged rock is studied, and a method for determining the latest application time of anchorage support is proposed. According to the test results of the damage surrounding rock acoustic wave velocity during the construction and the attenuation distribution law of acoustic wave velocity analyzed by the numerical method, this method can determine the latest application time of anchorage support for each part of the surrounding rock, which has good operability and practicability for actual engineering construction.