# Study on Energy Release of Surrounding Rock under the Multiple Unloading Disturbance during Tunnel Excavation

Authors: Yong Fan; Junwei Zheng; Xinyi Hu; Xianze Cui; Weiping He
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

During excavation of deep rock, the release of strain energy plays an important role in geologic hazards caused by excavation. However, in the previous studies, the influence of transient unloading of in situ stress caused by blasting excavation has been ignored, and the blasting excavation of a tunnel is regarded as a single blast process. In this paper, the dynamic adjustment process of strain energy and the energy storage limit of surrounding rocks caused by transient unloading of in situ stress under elastic conditions were firstly analyzed. Then, the brittle-ductile-plastic transition model based on the Hoek–Brown strength criterion was simulated in FLAC3D. Finally, the dynamic release process of strain energy of surrounding rocks caused by multiple unloading disturbances was analyzed during the excavation of 2# Underground Laboratory of Jinping II Hydropower Station employing the newly proposed index, energy release coefficient (ERC). Results show that the strain energy of surrounding rock masses firstly decreases, then increases, next reduces, and finally stabilizes under the transient unloading of in situ stress. In the process of dynamic change of strain energy, when the strain energy exceeds its storage limit, a large amount of strain energy will be released and thus will lead to damage of the surrounding rock masses. Because the cut holes and the first circle of breaking holes are far away from the final excavation boundary, the unloading disturbance to the strain energy of surrounding rock masses is small. Furthermore, the energy release of surrounding rock masses is mainly caused by the unloading of the last circle of breaking holes and perimeter holes, and the closer to the final excavation boundary, the more intense the energy release.


---

## Body

## 1. Introduction

In the western region of China, the high in situ stress induced by the large buried depth becomes one of the typical geological features of underground engineering, which is a threat to the safety and stability of surrounding rocks of a deep tunnel after excavation. The environment of high in situ stress will lead to the storage of high strain energy in rock masses. During excavation, as the initial stress constraint of excavated rock masses on the excavation surface is released, the unloading phenomenon of in situ stress occurs, which causes the initial stress and energy state of the surrounding rock masses to be broken, inducing a series of geological disasters such as damage, microearthquake, rockburst, and deformation in the surrounding rock masses [1–4]. A large number of field tests and theoretical studies showed that energy release was the essential cause of failure of rock masses [5–7]. It can be said that the rock mass is constantly exchanging energy with the outside world in the whole process of deformation, damage, and destruction. Thus, the study of energy release law of rock masses in the process of excavation plays an important role in revealing the formation mechanism of geological disasters.As early as the last century, some scholars analyzed the law of energy variation of surrounding rock masses in the excavation process of a deep tunnel from the perspective of statics. For instance, Cook et al. [8] proposed the concept and calculation formula of the energy release rate; afterwards, Walsh and Salamon [9, 10] studied the energy redistribution process of surrounding rocks during excavation according to the law of the conservation of energy. On this basis, Napier [11] employed the unified calculation method to study the energy change law in the process of intersection and rapid expansion of complex structural planes, and the research results showed that the energy release was affected not only by the excavation footage but also by the crack growth size near the excavation plane. Mitri et al. [12] applied the finite element method to calculate the energy release rate and energy storage rate of rockburst induced by excavation unloading and introduced the critical strain energy index to assess the risk of rockburst. Revuzhenkor and Klishin [13] (2009) employed the energy flow vector index and energy streamline to analyze the energy flow characteristics of surrounding rocks during the excavation of a deep-buried tunnel. During the excavation process, the in situ stress or strain energy of rock mass would adjust dramatically, and the strain energy would flow along the energy streamline from the outer boundary to the excavation surface, inducing the accumulation and release of strain energy of rock masses [14–16].However, most of the abovementioned studies regard the unloading of in situ stress on the excavation surface as a process of quasistatic unloading. In fact, this kind of approximate treatment method may be acceptable when the in situ stress is low, but when the local stress reaches 20 MPa–50 MPa or higher, the duration of unloading of in situ stress on the excavation face is only a few milliseconds, and the influences of high-speed (transient) unloading on the stress field and energy release of the surrounding rocks should not be ignored [17–21]. Furthermore, it should be mentioned that the abovementioned studies are based on the assumption that the tunnel was formed by a single blast process. In practical engineering, the cavern and tunnel blasting operation is a multiple blasting process, and millisecond delay blasting is generally employed in every cycle of footage. For instance, from the centre of the excavation surface to the periphery, it is often employed the detonation of cut holes, breaking holes, and perimeter holes in turn to complete a blasting operation. The detonation of each round of blastholes will release the stress constraint of the excavated rock mass on the original rock mass. Therefore, in the process of blasting excavation, the transient unloading of in situ stress constraint is actually a multiple unloading process. Also, the energy release law of surrounding rock masses will be different under the influence of multiple unloading disturbances. Thus, it is crucial to study the energy release law of surrounding rock masses under multiple unloading disturbances.In the present study, the dynamic evolution law of strain energy of surrounding rock masses during the transient unloading of in situ stress was firstly analyzed. After successfully describing the postpeak mechanical property of brittle-ductile-plastic transition of Jinping marble, taking the Jinping 2# Underground Laboratory as the engineering case, the release law of surrounding rock energy under multiple unloading disturbance was analyzed in FLAC3D. The analysis results can provide theoretical reference for understanding the geologic hazards and stability analysis during construction of a deep-buried tunnel.


## 2. Energy Release Condition of Surrounding Rock in Tunnel Excavation

The excavation and unloading of deep-buried tunnels will lead to the aggregating phenomenon of strain energy of surrounding rock. When the accumulated strain energy exceeds a certain threshold, it will cause the release of strain energy and lead to damage of surrounding rock masses [14]. Therefore, the key to study the release process of surrounding rock strain energy is to determine this limit, which is the energy storage limit of surrounding rock.Firstly, it was assumed that a circular tunnel was excavated under the condition of hydrostatic stress field. Before excavation, the constraint of excavated rock masses on the surrounding rock masses was regarded as reverse load acting on the excavation boundary which was equal to the in situ stress. Under the effect of the excavation, the load will decrease rapidly from the initial value to 0 and generate the phenomenon of transient unloading of in situ stress.To simplify the calculations, the unloading process of in situ stress on the excavated boundary is divided into an initial in situ stress state and excavation state. After the treatments mentioned above, the unloading process of in situ stress is equivalent to the superposition of the initial stress state and dynamic tensile load on the excavation boundary [22]. Based on the abovementioned analysis, the unloading problem of in situ stress can be further transformed into the problem of cavity excitation: at the beginning, there is time-dependent tensile stress acting on the cylindrical cavity in the stress-free state.In the cylindrical coordinate system, the governing equation of the elastic wave excited by the cylindrical cavity is shown in equation (1), the boundary condition is shown in equations (2), (3), and (4), in which the stress is positive in pressure and negative in tension, and the radius of the cylindrical cavity is assumed to be a0:(1)∂2u∂r2+1r∂u∂r−ur2=1Cp2∂2u∂t2,σrr,t=λ+2μ∂u∂r+λur,σθr,t=λ∂u∂r+λ+2μur,(2)ur,t=∂ur,t∂r=0,r≥a0,t<0,(3)σrr,tr=a0=0t≤0,−P0tt0,−P0t≥t0,0<t≤t0,(4)limr⟶∞ur,t=0,t>0,where Cp is the elastic longitudinal wave velocity, unit: m/s, Cp = λ+2μ/ρ, ρ is the density, unit: kg/m3; r is the distance from the centre of the column tunnel to the medium; σr and σθ are the radial stress and tangential stress caused by the unloading stress wave in the medium, respectively, unit: Pa; u(r, t) is the radial displacement, unit: m; λ and μ both are the Lame constant; P0 is the initial in situ stress, unit: Pa; and t0 is the duration of transient unloading, unit: s.After equation (1) is transformed by Laplace, the Laplace space solutions of radial stress and tangential stress are calculated according to (2), (3), and (4). Then, the inverse transformation of the space solution is obtained by employing contour integral. Finally, the dynamic stress field induced by transient unloading can be calculated by superposition of the original in situ stress field:(5)σrr,tP0=1,t<r−a0Cp,1−a022t0r2t−r−a0Cp−1t0π∫0∞Q1η,tdη,r−a0Cp≤t≤r−a0Cp+t0,1−a022r2−1t0π∫0∞Q2η,tdη,t>r−a0Cp,(6)σθr,tP0=1,t<r−a0Cp,1+a022t0r2t−r−a0Cp−1t0π∫0∞Q3η,tdη,r−a0Cp≤t≤r−a0Cp+t0,1+a022r2−1t0π∫0∞Q4η,tdη,t>r−a0Cp,where(7)Q1=Ar+Brcosr−a0/Cpη−costηE+Fη2+Cr−Drsintη−sinr−a0/CpηE+Fη2,Q2=Ar+Brcost−t0η−costηE+Fη2+Cr−Drsintη−sint−t0ηE+Fη2,Q3=Aθ+Bθcostη−cosr−a0/CpηE+Fη2+Cθ−Dθsintη−sinr−a0/CpηE+Fη2,Q4=Aθ+Bθcostη−cost−t0ηE+Fη2+Cθ−Dθsintη−sint−t0ηE+Fη2,where Ar=α1α2, Br=α3α4, Cr=α1α4, Dr=α2α3, Aθ=α5α2, Bθ=α6α4, Cθ=α5α4, Dθ=α2α6, E=α22, and F=α42.(8)G=λ+2μμ,α1=2G2rJ1x−ηCpJ0x,α2=2G2a0J1x0−ηCpJ0x0,α3=2G2rY1x−ηCpY0x,α4=2G2a0Y1x0−ηCpY0x0,α5=−2G2rJ1x−1−2G2ηCpJ0x,α6=−2G2rY1x−1−2G2ηCpY0x,where J0, J1, Y0, and Y1 are the Bessel functions; η is the integral path of the contour integral.Substituting the dynamic stress field calculated by equations (5) and (6) into equation (9), the principal stresses of surrounding rock masses under transient unloading are obtained [23].(9)σ1=σθ,σ2=νσr+σθ+1−2νP0,σ3=σr,where v is the Poisson ratio.Substituting the main stress calculated by equation (9) and the initial stress value of the original rock (σ1 = σ2 = σ3 = P0) into equation (10), the elastic strain energy density of the surrounding rock Uw(r, t) and the original strain energy density at different times U0 could be calculated, respectively.(10)U=σ12+σ22+σ32−2νσ1σ2+σ2σ3+σ1σ32E0,where E0 is the elastic modulus.Xie et al. [24] indicated that when the elastic strain energy stored in the rock masses exceeded the storage limit, it would release along the direction of one of the principle stress. At the same time, once the elastic strain energy in this direction exceeded its damage limit, the rock masses will be damaged. Considering that the surrounding rock masses of the tunnel are in hydrostatic compressive stress field in this study, according to the analysis in article [24], it can be concluded that the energy release Gi of the compressed rock masses in the direction of principle stress σi is as follows:(11)Gi=Kiσ1−σiUc,here Ki is the material constant.According to equation (11), it can be inferred that the energy release of surrounding rock masses in the direction of the minimum compressive stress is the largest. Now, the damage of the surrounding rock masses satisfies:(12)G3=K3σ1−σ3Uc,where Gc is the material parameter and the energy release peak of the surrounding rock masses which can be determined by a uniaxial compression test.Supposeσ1 = σc and σc is the uniaxial compressive strength of rock masses; σ2 = σ3 = 0, According to equation (10),(13)Uc=σc22E0.Based on the assumption of the abovementioned equations, equation (14) can be obtained by substituting equation (13) into equation (12).(14)Gc=K3σc32E0.Then, equation (14) is substituted into equation (12) to eliminate K3, and the strain energy storage limit of compressed surrounding rock masses during energy releasing can be obtained:(15)Uc=σc32E0σ1−σ3.Taking the surrounding rock masses ofr = 1.0a0 and r = 1.4a0 as examples, the ratios of the strain energy density of these two positions [Uw(r, t)] to the initial strain energy density of the surrounding rock masses (U0) under the transient unloading of in situ stress were calculated by equation (10); according to equation (15), the ratios of the energy storage limit of the surrounding rock masses [Uc(r, t)] to the initial strain energy density (U0) of these two positions under the transient unloading of in situ stress were calculated, respectively. At last, these four ratios are as shown in Figure 1.Figure 1
Adjustment history curve of the energy storage limit and strain energy density. The figure on the left is the schematic diagram of the cylindrical cavity, and the red arrow is the propagation of the unloading wave. The abscissa of the graph on the right is the time step.It can be seen from Figure1 that the unloading stress wave reaches the position of r = 1.0a0 at the time of t1 = 0, which causes the initial energy state of the surrounding rock masses to be broken and leads to the strain energy density of the surrounding rock masses experience the dynamic variation process that first decrease, then increase, and eventually, tend to stabilize. At the time of t2 = 0.4a0/Cp, the unloading stress wave reaches the position of r = 1.4a0, causing the strain energy density of the surrounding rock masses to experience the similar dynamic process as the position of r = 1.0a0. The difference between them is that the variation extent of the strain energy density caused by the transient unloading of in situ stress at different positions is different. It can be further seen that, from t1 to t3, with the propagation of unloading stress wave, the strain energy density of the surrounding rock masses at position of r = 1.0a0 under transient unloading is less than the initial strain energy density, and the strain energy density of surrounding rock masses is less than the energy storage limit of surrounding rock masses (Uw(r, t)/U0<Uc(r, t)/U0) during this time period; thus, there will be no energy release. From time t3 to t5, the strain energy density of surrounding rock masses at the position of r = 1.0a0 increases, which indicates that the surrounding rock masses is accumulating energy. When the strain energy exceeds its storage limit (t > t4, Uw(r, t)/U0>Uc(r, t)/U0), the strain energy will be released, resulting in the reduction of strain energy density of surrounding rock masses. For surrounding rock masses at the position of r = 1.4a0, from time t2 to t3, the strain energy density of surrounding rock masses is less than that of the original rock masses. In this time period, the strain energy density of surrounding rock masses is less than the energy storage limit (Uw(r, t)/U0<Uc(r, t)/U0), so it will not cause the energy release. Therefore, in the process of unloading stress wave propagating from the excavation boundary into the surrounding rocks, the strain energy of the surrounding rock masses first decreases, and then, absorbs the strain energy from the surrounding rock at a distance to promote the increase of strain energy. When the accumulated strain energy exceeds its storage limit, the energy will be released, and the strain energy of the surrounding rock masses will decrease again. The decrease and increase of the strain energy of the surrounding rock masses in the early stage are caused by the elastic unloading wave, which will not lead to the damage of the surrounding rock masses. The decrease of the strain energy of the surrounding rock masses in the second time is caused by the accumulated strain energy of the surrounding rock masses exceeds its energy storage limit, which will undoubtedly lead to the damage of the surrounding rock masses.In fact, during the process of excavation of a deep-buried tunnel, once the strain energy of surrounding rock masses exceeds its storage limit, the energy will be released, resulting in the decrease of strain energy. The variation process is shown in Figure2. During this process, there are four important eigenvalues: Uy, the initial value of the surrounding rock strain energy; Ug, the valley value of dynamic fluctuation; Uf, the peak value of dynamic fluctuation; and Us, the final stable value.Figure 2
Adjustment history curve of the energy storage limit and strain energy density.In view of the dynamic characteristics of surrounding rock energy changes under the excavation mentioned above, the energy release rate is proposed [25]:(16)ERRi=Uimax−Uis,where ERRi is the energy release rate of the i-th unit; Uimax is the maximum strain energy density of the i-th unit before energy release; and Uis is the final stable value of strain energy density after the energy release of the i-th unit.An energy release coefficient is proposed in order to more intuitively represent the energy release of surrounding rocks after transient unloading:(17)ERCi=Uimax−UisUimax.The energy release coefficient ERRi is between 0 and 1. The larger the value is, the more intense is the energy release, and the smaller the value is, the smoother is the energy release.However, the specific energy release process after the strain energy of the surrounding rock masses exceeding its limit (the black dotted line afterUw(r, t)/U0 exceeds Uc(r, t)/U0 in Figure 2) cannot be accurately obtained by theoretical analysis. Therefore, the release process of strain energy of surrounding rock masses under transient unloading is analyzed by the numerical simulation method in the following part.


## 3. Release Process of Strain Energy of Surrounding Rock Masses under the Multiple Excavation Disturbances

### 3.1. Project Background

Jinping II hydropower station is located on the bend of the Yalong River in Sichuan province, China. Taking the advantage of high natural drop to generate power, the total installed capacity is 4800 MW. There are several parallel tunnels across the Jinping Mountains to straighten the river bend. The average length of each tunnel is about 16.6 km, and the maximum buried depth is up to 2525 m. The first unit started operation in August 2013, and by the end of 2014, all units were put into operation to generate power.In order to explore the various rock mechanic problems in the long-term operation of deep-buried tunnels, the underground laboratories were planned to be arranged in the tunnels. The location of the underground laboratory is shown in Figure3. The underground laboratories utilize the auxiliary construction branch tunnel to excavate 4 sets of experimental caverns which are 130 m in length, 14 m in width, and 14 m in height. Two laboratories which are 65 m in length, 14 m in width, 14 m in height, and in the shape of city gate are set in each experimental cavern, forming a layout of 4 caverns and 8 laboratories. The location, in situ stress, and specific arrangement of the laboratories are shown in Figure 3. The laboratories were excavated in three layers, and the heights from the top layer to the bottom layer were 7.0 m, 5.5 m, and 1.5 m. The top layer adopted the excavation procedure that the middle drift excavation was advanced, and the two sides were expanded and followed. The arrangement of blast holes is shown in Figure 4, and the blasting parameters are shown in Table 1. The numbers ms1–ms13 in the Figure 4 are the nonelectric millisecond detonator sections.Figure 3
Adjustment history curve of the energy storage limit and strain energy density. (a) Layout of the Jinping Underground Laboratory. (b) In situ stress field of the Jinping Underground Laboratory.
(a)(b)Figure 4
Blastholes layout of middle drift excavation in the Jinping 2# Laboratory.Table 1
Blasting design parameters of 2# underground laboratory.
Hole typeCharge weightQuantityHole countDelayCut hole0.2128ms1Breaking hole0.21920ms30.21717ms50.23519ms70.22923ms9Periphery hole0.22350ms11


### 3.2. Numerical Model and Parameters

According to the analysis mentioned above, when the deep rock masses are affected by the effect of transient unloading of in situ stress caused by blasting excavation in the hydrostatic stress field, the strain energy of surrounding rock masses will generate a dynamic fluctuation phenomenon. However, in the case of a nonhydrostatic stress field or irregular excavation boundary, it is hard to obtain the dynamic release process of strain energy by theoretical calculation. Also, under the transient unloading condition, the storage limit and residual strain energy density are more difficult to be obtained than quasistatic unloading. Moreover, Cai [26] reported that the transient unloading of in situ stress accompanied by blasting excavation will produce a large unbalanced force on the excavation boundary and induce a significant dynamic effect. Since the Fast Lagrangian Analysis of Continua (FLAC) is a numerical method based on the explicit finite difference method, it can accurately simulate the dynamic mechanical process induced by transient unloading. Using FLAC, Cai [26] successfully simulated the stress adjustment process on the excavation boundary induced by transient unloading in tunnel excavation (Figure 5). As can be seen from the Figure 5, the minimum principal stress of the surrounding rock unit on the excavation boundary fluctuates dynamically at the instant of the formation for the excavation surface, which indicates that FLAC can capture the dynamic effect induced by transient unloading. Thus, FLAC is employed in this paper to simulate the dynamic energy release process of surrounding rock masses induced by transient unloading during tunnel excavation.Figure 5
Stress change path of surrounding rock on the excavation boundary induced by transient unloading [26].Before numerical calculation, the constitutive model matching the postpeak mechanical property of rock masses should be established. The lithology of the rock mass excavated in the Jinping Underground Laboratory is mainly marble. A large number of studies and tests have shown that the marble exhibits the mechanical properties of brittle-ductile-plastic transition under different confining pressures. For instance, Chu obtained marble samples from the east end of the Jinping auxiliary tunnel with more than 2000 m depth for the triaxial compression test [27]. The test result was performed as shown in Figure 6. It can be seen from Figure 6 that when the confining pressure is 2 MPa, the Jinping marble shows brittleness characteristic after peak. With the confining pressure increasing (2–8 MPa), the Jinping marble shows ductile characteristic after peak, and the residual strength is still high. While the confining pressure reaches 40 MPa, the residual strength after peak does not decrease, showing plastic characteristic. Therefore, the Jinping marble in different parts of the surrounding rocks will be in different stress environments after excavation, resulting in different mechanical properties. Thus, the numerical description method of the postpeak mechanical property of Jinping marble should be established at first.Figure 6
Triaxial test results of Jinping marble [27].In the equation of the Hoek–Brown strength criterion, strength parameters such asmb, s, and a can change with the accumulation of plastic strain ε3p after rock masses yielding. Thus, the Hoek–Brown strength criterion in FLAC3D can be applied to describe the strengthening and softening behaviour of the Jinping marble after yielding. The yielding equation is(18)σ1=σ3+σcmbσ3σc+sa,where mb, s, and a are the parameters associated with the rock mass quality GSI and the rock mass material parameter mi.In order to accurately describe the changes ofmb, s, a, and other parameters caused by the accumulation plastic strain ε3p under different confining pressures, Cundall et al. [28] proposed a factor μ1 related to the minimum principal stress σ3 in the description of softening and hardening of materials after yielding. The factor μ1 is employed to describe the change of strength parameters such as mb, s, and a with plastic strain ε3p under different confining pressure levels. Essentially, μ1 is a scaling factor. By properly setting the value of μ1 under different confining pressures, the brittle-ductile-plastic transition characteristics can be described. When determining μ1, there were still 8 parameters to describe the brittle-ductile-plastic transition characteristic of Jinping marble: 4 residual strength parameters and 4 peak strength parameters [29].The Hoek–Brown mechanical parameters and in situ stress field measured by the triaxial compression test and field tests are shown in Table2. However, the determination of the residual strength and scaling factor is a troublesome process. First, several groups of parameters were assumed for numerical calculation and, then, compared with the excavation damage extents detected by field tests. When the numerical calculation results were consistent with the field test results, the calculation parameters adopted could be considered as the ideal parameters. Through multiple numerical inversion tests, the ideal residual strength parameters and scaling factors were finally determined, as shown in Tables 3 and 4, respectively.Table 2
Hoek–Brown mechanical parameters and in situ stress field.
Compressive strengthElasticity modulusGSImimb140 MPa45 GPa7093.0827saStrength of extensionPoisson ratioBulk modulus0.03570.51.6 MPa0.2428.8 GPaShear elasticityThe horizontal direction:σxxThe axial direction:σyyThe vertical direction:σzz18.1 GPa39.7 MPa34.9 MPa58.3 MPaTable 3
Strength parameters of marble changing with plastic strain.
Signalmentε3p (10−3)σci (MPa)mbsaPeak strength0.001403.08270.03570.50Ductility section1.20903.08270.03570.50Residual strength2.30403.800000.70Table 4
Scaling factor related to the minimum principal stress.
σ3 (MPa)02481824.925μ11.000.850.750.650.450.350The numerical simulation curves of the axial stress changing with strain of marble under different confining pressures are shown in Figure7. According to the inversion stress-strain curve, when the confining pressure is low, the rock mass mainly shows the brittle feature. With the increase of the confining pressure, the ductile feature of the rock mass gradually occurs. Finally, when the confining pressure increases to 25 MPa, the rock mass shows the plastic feature. The correctness of the inversion results will be verified in the following paper by comparing the damage extents of surrounding rock masses simulated by employing the parameters in Tables 2–4 with the damage extents of surrounding rock masses detected in the field tests.Figure 7
Axial stress-strain curve of marble by numerical simulation.


### 3.3. Numerical Simulation of Energy Release

After successfully describing the postpeak mechanical properties of brittle-ductile-plastic transformation of Jinping marble, the parameters and constitutive model mentioned above were employed to simulate the energy variation process of surrounding rock masses induced by transient unloading of in situ stress during the blasting excavation of central pilot tunnel (7 m × 7 m, gate shape) in the Jinping Underground Laboratory. The numerical results are shown in Figure8. In Figure 8, Figure 8(a) is the calculation model, 22 m in width and 23 m in height. Also, the nonreflecting boundary condition (energy absorbing boundary condition) was adopted in the model. The axial excavation footage is 3 m and 26760 units in total.Figure 8
Calculation model and energy change of surrounding rock masses under multiple unloading (J/m3). (a) Calculation model. (b) Initial energy state. (c) ms1. (d) ms3. (e) ms5. (f) ms7. (g) ms9. (h) ms11.
(a)(b)(c)(d)(e)(f)(g)(h)From Figure8(b), it can be seen that the original rock masses will store a certain amount of strain energy before excavation, which is caused by the elastic deformation of the original rock under the effect of the initial in situ stress. From Figures 8(c)–8(h), it can be found that the energy state of the surrounding rock masses changes continuously with the sequential initiation of the cut holes, breaking holes, and the perimeter holes. In fact, with the transient unloading of the in situ stress on the excavation surface, the strain energy of the surrounding rock masses decreases rapidly and absorbed from the nearby rock masses. When the absorbed strain energy exceeds the energy storage limit, it will lead to the rapid release of strain energy and eventually cause damage to surrounding rock masses. After the release of strain energy of surrounding rock masses, it is stable in the state, as shown in Figures 8(c)–8(h).In order to monitor the dynamic release process of strain energy of surrounding rock masses in the process of multiple unloading disturbance, the variation processes of strain energy of units at the top, bottom, upper right corner, lower right corner, and right position of the central pilot tunnel in Figure8(a) were also recorded in the numerical calculation, taking the right side as an example, as shown in Figure 9.Figure 9
Energy release process of surrounding rock on the right side of the excavation section. (a) Rock mass unit at 0.2 m away from the excavation boundary. (b) Rock mass unit at 1.2 m away from the excavation boundary. (c) Rock mass unit at 2.2 m away from the excavation boundary. (d) Rock mass unit at 3.2 m away from the excavation boundary.
(a)(b)(c)(d)It can be seen from Figure9 that because the cut holes and the first round of breaking holes are far away from the final excavation boundary, the disturbance of excavation unloading of ms1 and ms3 to the surrounding rock energy is small, and the strain energy of surrounding rock masses almost increases linearly. However, the disturbance of transient unloading of in situ stress in ms5, ms7, ms9, and ms11 to the strain energy of surrounding rock masses increases gradually, and the closer to the final excavation boundary, the greater the disturbance. Under the influence of transient unloading, the energy of rock mass elements will experience the process of firstly reducing, then increasing, and finally, reducing. The first reduction is caused by the unloading wave, which will not cause damage and failure of surrounding rock. The second reduction is caused by the large release of strain energy which is induced by the accumulated strain energy in surrounding rock masses which exceeds its storage limit and leads to the damage or failure of surrounding rock masses.In addition, we can also find from Figure9 that the closer the surrounding rock masses are to the final excavation boundary, the greater the energy release coefficient and energy release rate. At 0.2 m away from the final excavation boundary, the energy release is mainly caused by the excavation unloading of ms5, ms7, and ms9, among which the energy release of surrounding rock masses caused by the excavation unloading of ms7 is the largest and fastest. At 1.2 m away from the final excavation boundary, the energy release is mainly caused by the excavation unloading of ms7, ms9, and ms11, among which the energy release of surrounding rock masses caused by the unloading of ms9 is the largest and the fastest. The energy release of the rock mass unit at 2.2 m away from the final excavation boundary is mainly caused by the unloading of ms9 and ms11, among which the unloading of ms11 causes the largest and the fastest energy release of the surrounding rock masses. However, there is almost no energy release of rock mass unit at 3.2 m away from the final excavation boundary.


### 3.4. Damage to Surrounding Rock Masses Caused by Energy Release

In the numerical calculation, the damage distributions of surrounding rock masses caused by the release of strain energy were recorded as shown in Figure10.Figure 10
Damage distributions of surrounding rock caused by transient unloading of in situ stress. Both shear-n and shear-p in the figure indicate that the rock mass is damaged. (a) ms1. (b) ms3. (c) ms5. (d) ms7. (e) ms9. (f) ms11.
(a)(b)(c)(d)(e)(f)It can be seen from Figures8 and 10 that there is a closely corresponding relationship between the variation law of strain energy of surrounding rock masses and the damage distribution of surrounding rock masses. After the unloading of ms1, ms3, and ms5, the strain energy of surrounding rock masses accumulates. Because the accumulated strain energy only exceeds the storage limit of surrounding rock masses in a small scale, there is no obvious damage in surrounding rock masses (the damage here only considers the surrounding rock masses beyond the final excavation boundary). However, after the unloading of ms7, ms9, and ms11, the strain energy gathered in the surrounding rock masses is further increasing, and the value exceeds the energy storage limit of the surrounding rock masses, resulting in a large amount of release of strain energy and damage to the surrounding rock masses, the main form of damage is shear failure.In order to verify the correctness of the model and the calculation methods in this study, two groups of acoustic holes were set at the right wall of the Jinping 2# Laboratory for acoustic detection. There were 3 holes in each group, which were arranged in a position of a triangle, 0.5 m∼1.5 m higher than the bottom of the cavern. The holes in each group were numbered A, B, and C, with a diameter of 76 mm, a depth of 6 m, and a spacing of about 1.0 m. The positions and arrangement of the acoustic holes in the cavern are shown in Figures10(f) and 11, respectively, and the testing results of two groups of acoustic holes after excavation are shown in Figure 12.Figure 11
Layout of acoustic holes on the right side of the excavation section in the Jinping 2# Laboratory.Figure 12
Testing results of acoustic holes. (a) Acoustic detection results of group 1. (b) Acoustic detection results of group 2.
(a)(b)Due to the influence of surrounding rock fracture on the transmission of acoustic wave, the wave velocity will be reduced. When the velocity of acoustic wave reduces rapidly, the rock mass can be regarded as being damaged. According to the Construction Specification on Underground Excavating Engineering of Hydraulic Structures (DL/T 5099-1999), when the acoustic velocity of the surrounding rock masses of the current depth is reduced more than 10% of the average velocity of the undisturbed rock masses, these rock masses can be regarded as being damaged. The greater the velocity reduction, the more the intensity of damage. Finally, from the acoustic data in Figure12, the damage depths of surrounding rock masses corresponding to each detection hole were obtained. The comparison between the testing damage ranges of the surrounding rocks and the numerical simulation of the damage range of the surrounding rocks is shown in Table 5. From the comparison, we can find that the average value of the measured result is consistent with the result of numerical simulation, which fully illustrates accuracy of the constitutive model selected in this paper.Table 5
Comparison of the numerical calculation result and acoustic detecting result of the damage zone of the Jinping 2# Laboratory.
Group numberDetection hole numberSurrounding rock damage depth (m)1A2.5B2.3C1.92A2.3B2.3C2.1Average depth of the testing result—2.2Numerical simulation result—2.3


## 3.1. Project Background

Jinping II hydropower station is located on the bend of the Yalong River in Sichuan province, China. Taking the advantage of high natural drop to generate power, the total installed capacity is 4800 MW. There are several parallel tunnels across the Jinping Mountains to straighten the river bend. The average length of each tunnel is about 16.6 km, and the maximum buried depth is up to 2525 m. The first unit started operation in August 2013, and by the end of 2014, all units were put into operation to generate power.In order to explore the various rock mechanic problems in the long-term operation of deep-buried tunnels, the underground laboratories were planned to be arranged in the tunnels. The location of the underground laboratory is shown in Figure3. The underground laboratories utilize the auxiliary construction branch tunnel to excavate 4 sets of experimental caverns which are 130 m in length, 14 m in width, and 14 m in height. Two laboratories which are 65 m in length, 14 m in width, 14 m in height, and in the shape of city gate are set in each experimental cavern, forming a layout of 4 caverns and 8 laboratories. The location, in situ stress, and specific arrangement of the laboratories are shown in Figure 3. The laboratories were excavated in three layers, and the heights from the top layer to the bottom layer were 7.0 m, 5.5 m, and 1.5 m. The top layer adopted the excavation procedure that the middle drift excavation was advanced, and the two sides were expanded and followed. The arrangement of blast holes is shown in Figure 4, and the blasting parameters are shown in Table 1. The numbers ms1–ms13 in the Figure 4 are the nonelectric millisecond detonator sections.Figure 3
Adjustment history curve of the energy storage limit and strain energy density. (a) Layout of the Jinping Underground Laboratory. (b) In situ stress field of the Jinping Underground Laboratory.
(a)(b)Figure 4
Blastholes layout of middle drift excavation in the Jinping 2# Laboratory.Table 1
Blasting design parameters of 2# underground laboratory.
Hole typeCharge weightQuantityHole countDelayCut hole0.2128ms1Breaking hole0.21920ms30.21717ms50.23519ms70.22923ms9Periphery hole0.22350ms11


## 3.2. Numerical Model and Parameters

According to the analysis mentioned above, when the deep rock masses are affected by the effect of transient unloading of in situ stress caused by blasting excavation in the hydrostatic stress field, the strain energy of surrounding rock masses will generate a dynamic fluctuation phenomenon. However, in the case of a nonhydrostatic stress field or irregular excavation boundary, it is hard to obtain the dynamic release process of strain energy by theoretical calculation. Also, under the transient unloading condition, the storage limit and residual strain energy density are more difficult to be obtained than quasistatic unloading. Moreover, Cai [26] reported that the transient unloading of in situ stress accompanied by blasting excavation will produce a large unbalanced force on the excavation boundary and induce a significant dynamic effect. Since the Fast Lagrangian Analysis of Continua (FLAC) is a numerical method based on the explicit finite difference method, it can accurately simulate the dynamic mechanical process induced by transient unloading. Using FLAC, Cai [26] successfully simulated the stress adjustment process on the excavation boundary induced by transient unloading in tunnel excavation (Figure 5). As can be seen from the Figure 5, the minimum principal stress of the surrounding rock unit on the excavation boundary fluctuates dynamically at the instant of the formation for the excavation surface, which indicates that FLAC can capture the dynamic effect induced by transient unloading. Thus, FLAC is employed in this paper to simulate the dynamic energy release process of surrounding rock masses induced by transient unloading during tunnel excavation.Figure 5
Stress change path of surrounding rock on the excavation boundary induced by transient unloading [26].Before numerical calculation, the constitutive model matching the postpeak mechanical property of rock masses should be established. The lithology of the rock mass excavated in the Jinping Underground Laboratory is mainly marble. A large number of studies and tests have shown that the marble exhibits the mechanical properties of brittle-ductile-plastic transition under different confining pressures. For instance, Chu obtained marble samples from the east end of the Jinping auxiliary tunnel with more than 2000 m depth for the triaxial compression test [27]. The test result was performed as shown in Figure 6. It can be seen from Figure 6 that when the confining pressure is 2 MPa, the Jinping marble shows brittleness characteristic after peak. With the confining pressure increasing (2–8 MPa), the Jinping marble shows ductile characteristic after peak, and the residual strength is still high. While the confining pressure reaches 40 MPa, the residual strength after peak does not decrease, showing plastic characteristic. Therefore, the Jinping marble in different parts of the surrounding rocks will be in different stress environments after excavation, resulting in different mechanical properties. Thus, the numerical description method of the postpeak mechanical property of Jinping marble should be established at first.Figure 6
Triaxial test results of Jinping marble [27].In the equation of the Hoek–Brown strength criterion, strength parameters such asmb, s, and a can change with the accumulation of plastic strain ε3p after rock masses yielding. Thus, the Hoek–Brown strength criterion in FLAC3D can be applied to describe the strengthening and softening behaviour of the Jinping marble after yielding. The yielding equation is(18)σ1=σ3+σcmbσ3σc+sa,where mb, s, and a are the parameters associated with the rock mass quality GSI and the rock mass material parameter mi.In order to accurately describe the changes ofmb, s, a, and other parameters caused by the accumulation plastic strain ε3p under different confining pressures, Cundall et al. [28] proposed a factor μ1 related to the minimum principal stress σ3 in the description of softening and hardening of materials after yielding. The factor μ1 is employed to describe the change of strength parameters such as mb, s, and a with plastic strain ε3p under different confining pressure levels. Essentially, μ1 is a scaling factor. By properly setting the value of μ1 under different confining pressures, the brittle-ductile-plastic transition characteristics can be described. When determining μ1, there were still 8 parameters to describe the brittle-ductile-plastic transition characteristic of Jinping marble: 4 residual strength parameters and 4 peak strength parameters [29].The Hoek–Brown mechanical parameters and in situ stress field measured by the triaxial compression test and field tests are shown in Table2. However, the determination of the residual strength and scaling factor is a troublesome process. First, several groups of parameters were assumed for numerical calculation and, then, compared with the excavation damage extents detected by field tests. When the numerical calculation results were consistent with the field test results, the calculation parameters adopted could be considered as the ideal parameters. Through multiple numerical inversion tests, the ideal residual strength parameters and scaling factors were finally determined, as shown in Tables 3 and 4, respectively.Table 2
Hoek–Brown mechanical parameters and in situ stress field.
Compressive strengthElasticity modulusGSImimb140 MPa45 GPa7093.0827saStrength of extensionPoisson ratioBulk modulus0.03570.51.6 MPa0.2428.8 GPaShear elasticityThe horizontal direction:σxxThe axial direction:σyyThe vertical direction:σzz18.1 GPa39.7 MPa34.9 MPa58.3 MPaTable 3
Strength parameters of marble changing with plastic strain.
Signalmentε3p (10−3)σci (MPa)mbsaPeak strength0.001403.08270.03570.50Ductility section1.20903.08270.03570.50Residual strength2.30403.800000.70Table 4
Scaling factor related to the minimum principal stress.
σ3 (MPa)02481824.925μ11.000.850.750.650.450.350The numerical simulation curves of the axial stress changing with strain of marble under different confining pressures are shown in Figure7. According to the inversion stress-strain curve, when the confining pressure is low, the rock mass mainly shows the brittle feature. With the increase of the confining pressure, the ductile feature of the rock mass gradually occurs. Finally, when the confining pressure increases to 25 MPa, the rock mass shows the plastic feature. The correctness of the inversion results will be verified in the following paper by comparing the damage extents of surrounding rock masses simulated by employing the parameters in Tables 2–4 with the damage extents of surrounding rock masses detected in the field tests.Figure 7
Axial stress-strain curve of marble by numerical simulation.


## 3.3. Numerical Simulation of Energy Release

After successfully describing the postpeak mechanical properties of brittle-ductile-plastic transformation of Jinping marble, the parameters and constitutive model mentioned above were employed to simulate the energy variation process of surrounding rock masses induced by transient unloading of in situ stress during the blasting excavation of central pilot tunnel (7 m × 7 m, gate shape) in the Jinping Underground Laboratory. The numerical results are shown in Figure8. In Figure 8, Figure 8(a) is the calculation model, 22 m in width and 23 m in height. Also, the nonreflecting boundary condition (energy absorbing boundary condition) was adopted in the model. The axial excavation footage is 3 m and 26760 units in total.Figure 8
Calculation model and energy change of surrounding rock masses under multiple unloading (J/m3). (a) Calculation model. (b) Initial energy state. (c) ms1. (d) ms3. (e) ms5. (f) ms7. (g) ms9. (h) ms11.
(a)(b)(c)(d)(e)(f)(g)(h)From Figure8(b), it can be seen that the original rock masses will store a certain amount of strain energy before excavation, which is caused by the elastic deformation of the original rock under the effect of the initial in situ stress. From Figures 8(c)–8(h), it can be found that the energy state of the surrounding rock masses changes continuously with the sequential initiation of the cut holes, breaking holes, and the perimeter holes. In fact, with the transient unloading of the in situ stress on the excavation surface, the strain energy of the surrounding rock masses decreases rapidly and absorbed from the nearby rock masses. When the absorbed strain energy exceeds the energy storage limit, it will lead to the rapid release of strain energy and eventually cause damage to surrounding rock masses. After the release of strain energy of surrounding rock masses, it is stable in the state, as shown in Figures 8(c)–8(h).In order to monitor the dynamic release process of strain energy of surrounding rock masses in the process of multiple unloading disturbance, the variation processes of strain energy of units at the top, bottom, upper right corner, lower right corner, and right position of the central pilot tunnel in Figure8(a) were also recorded in the numerical calculation, taking the right side as an example, as shown in Figure 9.Figure 9
Energy release process of surrounding rock on the right side of the excavation section. (a) Rock mass unit at 0.2 m away from the excavation boundary. (b) Rock mass unit at 1.2 m away from the excavation boundary. (c) Rock mass unit at 2.2 m away from the excavation boundary. (d) Rock mass unit at 3.2 m away from the excavation boundary.
(a)(b)(c)(d)It can be seen from Figure9 that because the cut holes and the first round of breaking holes are far away from the final excavation boundary, the disturbance of excavation unloading of ms1 and ms3 to the surrounding rock energy is small, and the strain energy of surrounding rock masses almost increases linearly. However, the disturbance of transient unloading of in situ stress in ms5, ms7, ms9, and ms11 to the strain energy of surrounding rock masses increases gradually, and the closer to the final excavation boundary, the greater the disturbance. Under the influence of transient unloading, the energy of rock mass elements will experience the process of firstly reducing, then increasing, and finally, reducing. The first reduction is caused by the unloading wave, which will not cause damage and failure of surrounding rock. The second reduction is caused by the large release of strain energy which is induced by the accumulated strain energy in surrounding rock masses which exceeds its storage limit and leads to the damage or failure of surrounding rock masses.In addition, we can also find from Figure9 that the closer the surrounding rock masses are to the final excavation boundary, the greater the energy release coefficient and energy release rate. At 0.2 m away from the final excavation boundary, the energy release is mainly caused by the excavation unloading of ms5, ms7, and ms9, among which the energy release of surrounding rock masses caused by the excavation unloading of ms7 is the largest and fastest. At 1.2 m away from the final excavation boundary, the energy release is mainly caused by the excavation unloading of ms7, ms9, and ms11, among which the energy release of surrounding rock masses caused by the unloading of ms9 is the largest and the fastest. The energy release of the rock mass unit at 2.2 m away from the final excavation boundary is mainly caused by the unloading of ms9 and ms11, among which the unloading of ms11 causes the largest and the fastest energy release of the surrounding rock masses. However, there is almost no energy release of rock mass unit at 3.2 m away from the final excavation boundary.


## 3.4. Damage to Surrounding Rock Masses Caused by Energy Release

In the numerical calculation, the damage distributions of surrounding rock masses caused by the release of strain energy were recorded as shown in Figure10.Figure 10
Damage distributions of surrounding rock caused by transient unloading of in situ stress. Both shear-n and shear-p in the figure indicate that the rock mass is damaged. (a) ms1. (b) ms3. (c) ms5. (d) ms7. (e) ms9. (f) ms11.
(a)(b)(c)(d)(e)(f)It can be seen from Figures8 and 10 that there is a closely corresponding relationship between the variation law of strain energy of surrounding rock masses and the damage distribution of surrounding rock masses. After the unloading of ms1, ms3, and ms5, the strain energy of surrounding rock masses accumulates. Because the accumulated strain energy only exceeds the storage limit of surrounding rock masses in a small scale, there is no obvious damage in surrounding rock masses (the damage here only considers the surrounding rock masses beyond the final excavation boundary). However, after the unloading of ms7, ms9, and ms11, the strain energy gathered in the surrounding rock masses is further increasing, and the value exceeds the energy storage limit of the surrounding rock masses, resulting in a large amount of release of strain energy and damage to the surrounding rock masses, the main form of damage is shear failure.In order to verify the correctness of the model and the calculation methods in this study, two groups of acoustic holes were set at the right wall of the Jinping 2# Laboratory for acoustic detection. There were 3 holes in each group, which were arranged in a position of a triangle, 0.5 m∼1.5 m higher than the bottom of the cavern. The holes in each group were numbered A, B, and C, with a diameter of 76 mm, a depth of 6 m, and a spacing of about 1.0 m. The positions and arrangement of the acoustic holes in the cavern are shown in Figures10(f) and 11, respectively, and the testing results of two groups of acoustic holes after excavation are shown in Figure 12.Figure 11
Layout of acoustic holes on the right side of the excavation section in the Jinping 2# Laboratory.Figure 12
Testing results of acoustic holes. (a) Acoustic detection results of group 1. (b) Acoustic detection results of group 2.
(a)(b)Due to the influence of surrounding rock fracture on the transmission of acoustic wave, the wave velocity will be reduced. When the velocity of acoustic wave reduces rapidly, the rock mass can be regarded as being damaged. According to the Construction Specification on Underground Excavating Engineering of Hydraulic Structures (DL/T 5099-1999), when the acoustic velocity of the surrounding rock masses of the current depth is reduced more than 10% of the average velocity of the undisturbed rock masses, these rock masses can be regarded as being damaged. The greater the velocity reduction, the more the intensity of damage. Finally, from the acoustic data in Figure12, the damage depths of surrounding rock masses corresponding to each detection hole were obtained. The comparison between the testing damage ranges of the surrounding rocks and the numerical simulation of the damage range of the surrounding rocks is shown in Table 5. From the comparison, we can find that the average value of the measured result is consistent with the result of numerical simulation, which fully illustrates accuracy of the constitutive model selected in this paper.Table 5
Comparison of the numerical calculation result and acoustic detecting result of the damage zone of the Jinping 2# Laboratory.
Group numberDetection hole numberSurrounding rock damage depth (m)1A2.5B2.3C1.92A2.3B2.3C2.1Average depth of the testing result—2.2Numerical simulation result—2.3


## 4. Results and Discussion

### 4.1. Conclusions

The conclusions are as follows:(1)
With the transient unloading of the in situ stress on the excavation surface, the strain energy of the surrounding rock masses near the excavation surface reduces rapidly and, then, absorbs the strain energy from the surrounding rock masses. When the accumulated strain energy exceeds the energy storage limit, it will cause a large amount of strain energy release and damage to the rock masses. In fact, with the sequential initiation of cut holes, breaking holes, and perimeter holes, the transient unloading effect causes the accumulation of strain energy in the surrounding rock masses near the excavation surface. When the accumulated strain energy exceeds the energy storage limit of the surrounding rock masses, a large amount of strain energy is released, resulting in the damage of the surrounding rock masses, and the main form of the damage is shear failure.(2)
The closer the rock masses are to the excavation surface, the greater the energy release rate and the release coefficient are. If the rock mass units are close to the final excavation boundary, the energy release is mainly caused by the unloading of the breaking holes and the perimeter holes. Here, the energy release of surrounding rock masses caused by the unloading of the breaking holes in the second round is the largest and fastest. However, the strain energy of rock mass units far away from the final excavation boundary is almost not released.


## 4.1. Conclusions

The conclusions are as follows:(1)
With the transient unloading of the in situ stress on the excavation surface, the strain energy of the surrounding rock masses near the excavation surface reduces rapidly and, then, absorbs the strain energy from the surrounding rock masses. When the accumulated strain energy exceeds the energy storage limit, it will cause a large amount of strain energy release and damage to the rock masses. In fact, with the sequential initiation of cut holes, breaking holes, and perimeter holes, the transient unloading effect causes the accumulation of strain energy in the surrounding rock masses near the excavation surface. When the accumulated strain energy exceeds the energy storage limit of the surrounding rock masses, a large amount of strain energy is released, resulting in the damage of the surrounding rock masses, and the main form of the damage is shear failure.(2)
The closer the rock masses are to the excavation surface, the greater the energy release rate and the release coefficient are. If the rock mass units are close to the final excavation boundary, the energy release is mainly caused by the unloading of the breaking holes and the perimeter holes. Here, the energy release of surrounding rock masses caused by the unloading of the breaking holes in the second round is the largest and fastest. However, the strain energy of rock mass units far away from the final excavation boundary is almost not released.



---

*Source: 6486815-2020-08-26.xml*