Numerical Simulation of the Dynamic Responses and Cumulative Damage of Underground Caverns under Multiple Explosions

Attacking underground caverns with earth-penetrating bombs usually involves multiple explosions in succession. To assess the dynamic responses and cumulative damage of underground caverns under multiple explosions, based on a reduced-scale physical model test, the modified Riedel–Hiermaier–+oma (RHT) model in the finite-element software LS-DYNA is used to build an underground cavern model that encounters four explosions above the vault. +e characteristics of the stress wave attenuation and the evolution laws for the cumulative damage of the surrounding rock in the process of the four explosions are presented. Also, the displacement of the vault, the strain of the cavern wall, and the damage of a rock bolt-supported cavern and an unanchored cavern are compared. +e results indicate that the peak pressure is attenuated increasingly in the latter three explosions. +e circumferential strain of the cavern wall changes from tensile to compressive from the vault to the corner. +e damage of the surrounding rock on the left and right sides of the explosion source is attenuated with increasing distance from the explosion source, and the attenuation curve has a reverse “S” shape. Moreover, the attenuation rate of the curve decreases with each explosion. Multiple explosions do not affect the size of the crushed zone, but they do increase the range of the fracture zone. With each explosion, the cumulative damage of the surrounding rock increases irreversibly, but the damage increment decreases. +e cumulative damage of the surrounding rock exhibits a highly nonlinear relationship with successive explosions, and the effect of the rock bolt reinforcement becomes more obvious with successive explosions. Accordingly, the present research results offer a reference for antiexplosion design and support the optimization of underground engineering.


Introduction
As urbanization advances, ground space no longer meets the needs of human development. For environmental and economic reasons, projects are increasingly being built underground, such as subways, bomb shelters, and storage rooms. Because of its advantages of long service life, low resource consumption, and large storage capacity, underground engineering has been used widely in civil and military fields. However, earth-penetrating weapons have been threatening the safety of underground engineering. Nowadays, a penetrating bomb usually explodes multiple times when attacking a target; for example, the GBU-28 guided bomb may explode many times during the penetrating process [1]. It has been reported that surrounding rock subjected to repeated low-level blasting is damaged more than it is by a single high-level event [2]. With the increasing of repetitive impacts, Young's modulus of the rock mass slightly increases first and then decreases, resulting in more severe damage of the rock mass [3]. erefore, ensuring that underground engineering is unaffected by such explosions has received considerable attention over the past decades.
A good understanding of the dynamic responses and cumulative damage of underground engineering under explosion is an essential prerequisite for improving its safety. Fortunately, there has been much research into how the explosion affects underground engineering. Many scholars have studied the dynamic responses of underground engineering under explosion through field tests [4,5]. Ramulu et al. [6] investigated surrounding rock under repeated blast loading during tunneling works and obtained a specific threshold vibration level for safe construction; however, unfortunately, they did not study the cumulative damage of the rock mass. As is well known, field tests are dangerous and uneconomical, so instead scholars have analyzed the dynamic responses and damage of underground caverns under explosion utilizing reduced-scale physical model tests [7,8]. Using the artificial gravity field generated by a centrifuge, a large explosion can be simulated with fewer explosives [9]. erefore, scholars have established centrifuge models to evaluate how surface explosion affects underground structures [10]. Also, as the processing capability of computers has grown, numerical simulations have entered widespread use in this field. Mussa et al. [11] conducted numerical research on the damage and antiexplosion ability of an underground box tunnel. Using the discrete-element (DE) program UDEC, many scholars [12] have studied the propagation of explosion waves in jointed rock masses. Dong et al. [13] used the extended finite-element (FE) method to investigate the energy dissipation caused by damage of roof strata during longwall mining. Meanwhile, Karrech et al. [14] used a new numerical approach to research the hydraulic damage of resource reservoirs; their method uses thermodynamic principles to describe the propagation of damage in porous media. e peak particle velocity (PPV), displacement, acceleration, and stress of the surrounding rock are usually used to study the damage and dynamic response of underground projects under explosion. e PPV has been proved to be the most representative parameter for the dynamic response of underground engineering [15]. Gao et al. [16] used numerical methods to analyze the displacement and damage characteristics of the surrounding rock with different buried depths, lateral pressure coefficients, excavation methods, and section shapes during tunnel excavation. Using the parameters of stress and acceleration, Yang et al. [17] studied how the ground explosion affected a shallow-buried subway tunnel. An appropriate material model is a reliable prerequisite for numerical calculation. Currently, many material models have been developed and applied to rock mass. e soil/ crushable foam model [18] was designed to handle foam or soil-like materials with geometric boundaries; however, full plastic flow is used to approximate the yield behavior of the materials; therefore, various softening behaviors of concrete under different loading conditions cannot be observed. e viscous-plastic model is applied to the failure analysis of brittle materials under tensile and shear stresses; however, that model ignores the third stress deviator invariant, strain rate, and postpeak softening process under extreme load, so its applicability is rather limited. e kinematic hardening cap model [19] can be formulated by simple parameters, but it is not good at predicting the softening behavior of concrete, nor can it describe the inhibition of confining pressure on the expansion of concrete. None of the aforementioned material models can depict the cumulative damage of the rock mass. erefore, scholars have looked for others that can. e Holmquist-Johnson-Cook material model can represent the compression damage and dynamic responses of brittle materials such as concrete under large strain, high strain rate, and high confining pressure [20], but it ignores the tensile damage of the materials. e Taylor-Chen-Kuszmaul model [21] considers a rock mass containing many randomly distributed cracks that are reactivated under tensile load; that model can simulate well the tensile damage of concrete but cannot describe the compression damage of concrete. e continuous surface cap model describes well the damage of concrete under low confining pressure [22], but many parameters must be determined. Gaede et al. [23] proposed an anisotropic continuum damage mechanics model based on irreversible thermodynamics, which is suitable for brittle rock damage analysis. e RHT material model improves the treatment of the hydrostatic tensile and compressive zones and reflects comprehensively the tensile and compressive damage, strain rate effect, strain hardening, softening, and failure of materials [24,25]. Some authors have obtained ideal simulation results when using the RHT model [26][27][28]. Because tensile and compressive damage must be considered simultaneously for the rock mass studied herein, the RHT model is utilized to simulate the rock mass through many attempts and comparative analysis.
In general, although the dynamic responses and damage of underground caverns under explosion have been studied thoroughly and fruitful results have been obtained, most of those previous studies dealt exclusively with a single explosion. To date, there have been few studies involving multiple explosions and even fewer involving the relationship between the explosion number and the cumulative damage of the surrounding rock. However, multiple explosions are more common in underground projects. If only analyzing the dynamic response and damage of the cavern under a single explosion, the mechanism of the surrounding rock damage under multiple explosions cannot be fully revealed. To evaluate the effects of multiple explosions around an underground cavern, based on a similarity model test, the finiteelement software LS-DYNA3D and the modified RHT model are utilized to analyze the dynamic responses and cumulative damage of an underground cavern under four explosions above its vault. e validity of the numerical model is proved by comparing the stress-time curves of the surrounding rock between the numerical and test model.

Overview of Test Model
e test model simulates a straight-wall arched-top cavern with a buried depth of 20 m and a span of 4 m. e test devices have been described previously by Wang et al. [29,30]; therefore, only an overview is given here. e surrounding rock of the cavern is fixed by four lateral limit devices that can move back and forth; being aluminum plates with high porosity, they eliminate the reflection of the stress wave and accurately restore the infinite boundary conditions of the surrounding rock. Because the test is a similarity model test, the size and material of the model must be converted according to the Froude similarity criterion. e similarity coefficients of the density, length, and stress  [29]. e surrounding rock consists of sand, cement, water, and accelerating agent at a mixed proportion of 15 : 1 : 1.6 : 0.0166. As shown in Figure 1(a), the cavern is supported by single-row full-length bonded rock bolts with a spacing and row spacing of 4 cm and a length of 18 cm; these are simulated as aluminum rods. e test assumes that the underground cavern is attacked by an earth-penetrating bomb, with four explosions occurring above the vault. Figure 1(b) presents the locations of the four explosions.
To ensure that the buried depth of the explosives in the model does not cause the throwing phenomenon of the overlying rock mass [31], the scaled depth h/W 1/3 (where h is the buried depth of the explosive and W is the charge mass) of each explosive is set as 1.71 m· kg −1/3 . e buried depths and charge masses of the four explosives are given in Table 1.

Finite-Element Model.
Considering the symmetries of the model, to save computation time, a model with a thickness of 4 cm (the row spacing of the rock bolts) is established after approximately two-dimensional processing. Referring to the locations and charge masses of the explosives in the test model, four explosives are arranged at different depths in the calculation model. To analyze fully the dynamic responses and damage of the cavern after each explosion, the detonation interval of the explosion source was set as 4 ms. To compare and analyze the antiexplosion abilities of anchored and unanchored caverns, a cavern without rock bolt support is also simulated. e dimensions of the calculation model and the locations of the explosives are shown in Figure 2(a).

Material Model and Parameters.
e left, right, and lower boundaries of the model are set as nonreflection conditions to meet the infinite boundaries in the actual situation.
e three-dimensional Solid164 element is adopted for the rock mass and explosive, and the Beam161 element is used for the rock bolts. e rock mass is divided into a 1 cm grid, with that around the explosion source being divided more finely, as shown in Figure 2(b).

Surrounding Rock.
e RHTconcrete model is utilized to simulate the rock mass, which takes into account the failure surface, residual surface, and elastic limit surface related to pressure [32], as indicated in Figure 3(a). e equivalent force strength of the failure surface σ * eq can be expressed by the product of three fractional functions related to pressure P, Lode angle θ, and equivalent strain rate _ ε. e relationship is as follows: where Y * TXC (P) is the compressive meridian strength, R 3 (θ) is the Lode angle factor, and F rate (_ ε) is the strain rate dynamic enhancement factor. ey are given by the following equations: where P * spall is the normalized spallation strength; r t is the deviator stress at the tensile meridian; r c is the deviator stress at the compression meridian; B Q is the pressure influence parameter; J 2 and J 3 are the second and third invariants, respectively; and A, N, α, δ, and Q 2,0 are material constants. e equivalent stress of the elastic limit surface is derived from the failure surface, as defined in the following equation: F el is the elastic scaling factor, which can be expressed as follows:    where R t and R c are material parameters, f t,el is the elastic limit of uniaxial tension, and f c,el is the elastic limit of uniaxial compression. F cap can be expressed by an elliptical cap function: where P 0 is the initial pore compress pressure and f c is the uniaxial compressive strength. e residual strength is related to the hydrostatic pressure, and its function expression is in exponential form: where P * is the normalized hydrostatic pressure and B and M are material constants. e relationship between the pressure P and porosity α in the RHT model needs to be calculated by the equation of state (EOS). Figure 3(b) depicts the EOS of P-α in the RHT model. e model is elastic when the pressure is lower than the pore crushing pressure. Once the pressure exceeds the pore crushing pressure, the stiffness of the material decreases, and the plastic strain begins to appear in the model. e cumulative damage does not occur either in the elastic stage or in the plastic strain stage. Only when the equivalent stress exceeds the failure stress does the material begin to accumulate damage. e damage index D represents the damage degree of the material, which ranges from zero to one. e closer the damage index D is to one, the higher the damage degree is. e damage index D is given as follows: where Δε P is the equivalent plastic strain increment and ε failure P is the final equivalent plastic strain, and it is defined as follows: among them, D 1 and D 2 are material parameters and E f, min is the minimum equivalent plastic strain when the material fails. A large number of parameters need to be determined in the RHT model, and most of them are determined through experiments. Some parameters are more difficult to obtain, and the original parameters provided by the model or modified parameters are used in the research process [33,34]. In this paper, some parameters of the RHT model are obtained by experiments, and the rest parameters are referred to the article [32] and modified appropriately, as illustrated in Table 2.

TNT.
e TNT is modeled as a high-energy explosive material ( * MAT_HIGH_EXPLOSIVE_BURN) [35]. e Jones-Wilkins-Lee EOS [35] is used to define the detonation pressure P, which is expressed as follows: where A, B, R 1 , R 2 , and ω are material constants, V is the current volume divided by the initial volume, and E 0 is the initial specific internal energy. e TNT parameters are selected from Meyer [36] and are listed in Table 3.

Rock Bolts.
e aluminum rods that simulate the rock bolts in the test are elastoplastic. erefore, the isotropic elastoplastic model * MAT_PLASTIC_KINEMATIC [35] is adopted to simulate the rock bolts in the numerical model; this describes the isotropic hardening and the follow-up hardening plastic of the rock bolts. Table 3 lists the rock bolts parameters. Shock and Vibration

Comparison between Test and Numerical Results.
To validate the accuracy of the numerical analysis, the stresstime curves of measuring points at the same locations in the simulation and test model are compared. Measuring points P1, P2, and P3 are arranged at 20, 35, and 61 cm, respectively, directly below the first explosion source, as shown in Figure 4. Figure 5 compares the stress-time curves at P1, P2, and P3 after the first explosion in the numerical calculation and test; the stress-time curves of the test are those recorded by Wang et al. [29]. Figure 5 shows that the stress-time curves at the same location in the test and numerical calculation have similar trends: they rise rapidly to a peak initially and then fall gradually to a relatively stable value. e stress peaks at the same measuring point are relatively close and have the same order of magnitude. However, the stress curves in the simulation rise earlier than those in the test, which indicates that the stress in the simulation reaches the measuring point earlier than it does in the test. Besides, the residual stress of the stress-time curve in the numerical calculation is less than that in the test. e main reason for this phenomenon is that the surrounding rock in the numerical simulation is regarded as an isotropic medium, while the model in the test is made of cement mortar, which is formed by layered compaction. Additionally, the rock mass layer and joint surface were reserved manually during the construction in the test, which can slow down the propagation of the stress wave and increase its attenuation. erefore, the stress wave in the simulation arrives at each measuring point before that in the test. Due to the existence of layer and joint surface in the test surrounding rock, the unrecoverable residual deformation occurs after the explosion, and the residual deformation causes the residual stress. Because of the isotropy, the residual stress of the surrounding rock in the numerical TNT U1 Figure 4: Arrangement of measuring points (unit: cm).
calculation is relatively small. Meanwhile, the residual stress at the measuring points P1, P2, and P3 in the test becomes smaller and smaller, mainly because with increasing distance from the explosion source, the explosive power on the surrounding rock decreases. erefore, the residual stress of the surrounding rock at P1, P2, and P3 decreases gradually. e stress wave attenuation formula given in [37] has been used widely to investigate the propagation characteristics of stress waves. e attenuation formula involves a power-law relationship between the peak pressure and the scaled distance R/W 1/3 (where R is the distance from the explosion source to the measuring point and W is the charge mass), namely, where P max is the peak pressure, c is a typical value that depends on the charge material and soil properties, and n is an attenuation factor that depends mainly on the material properties. e attenuation formulas of the peak pressure and scaled distance of the measuring points are plotted in Figure 6. e shape of the peak pressure attenuation curve in the test is similar to that in the simulation. e curve is steep at first and then becomes slower. Otherwise, the attenuation factor of the calculation model is 1.207, which is less than the attenuation factor of 1.345 in the test, indicating that the attenuation of the stress wave in the test is faster than that in the numerical calculation. is is consistent with the fact that the rock mass layer and joint surface in the test model can increase the attenuation of the stress wave. Considering what has been discussed, it can be seen that the results of the numerical simulation have a high degree of credibility.

Attenuation Characteristics of Stress Wave.
e detonation pressure is a major factor in the damage of an underground cavern, from which we can study the law governing the stress wave propagation. To study the attenuation characteristics of the stress wave, three measuring points are arranged at a distance R directly below each explosion source. e peak pressure and scaled distance R/ W 1/3 are given in Table 4. Figure 7 shows the peak pressure versus the scaled distance of the measuring points, as well as the fitted attenuation formulas of the peak pressure and scaled distance after each explosion. According to the attenuation formulas, except for the fact that the attenuation factor of the second explosion is smaller than that of the first explosion, the other two explosions are larger than the first explosion. e reason is that the primary functions of the first explosion are to (i) densify the surrounding rock, (ii) close the original structure surface of the rock mass, and (iii) narrow the cracks. erefore, the attenuation of the peak pressure in the second explosion slows down, and the attenuation factor becomes smaller. e second, third, and fourth explosions increase the damage degree of the surrounding rock, and the number and widths of the cracks in the rock mass increase gradually.
is accelerates the attenuation of the peak pressure in the latter three explosions. As can be seen, the attenuation factors of the peak pressure in the last two explosions do not differ by much, but they are quite different from that of the second explosion. e reason is that there is only one measuring point in the anchored zone in the second explosion but more in the latter two explosions. As we know, the reinforcement effect of the rock bolts improves the strength of the rock mass; therefore, the attenuation factor of the stress wave in the last two explosions does not change much, increasing only slightly. Shock and Vibration

Displacement of Vault.
Being closest to the explosion source, the vault is the most vulnerable to damage. erefore, to study the safety of the cavern, it is essential to analyze the deformation of the vault, and to do so, the displacement measuring point U1 is selected on the vault (as shown in Figure 4). Figure 8 shows the displacement-time curves at U1 in the anchored and unanchored cavern. It can be seen that the vault displacement increases in a stepwise manner; that is to say, the vault displacement increases by a certain amount after each explosion. In the first two explosions, the peak displacement of the anchored cavern vault is 0.797 and 1.168 mm, respectively. Excluding the fact that the power of the second explosion is greater than that of the first explosion, the difference in the peak displacement caused by each explosion is very small. e main reason is that the first explosion acts mainly to make the surrounding rock mass more compact, and its strength is increased to a certain extent. erefore, the vault displacement increases little in the second explosion. With decreasing distance from the explosion source and increasing charge mass, the latter two explosions are more powerful. e large explosive power and cumulative damage effect mean that the peak displacement of the vault in the latter two explosions increases markedly, reaching 3.620 and 7.062 mm, respectively. e displacement-time curves of the first two explosions rise to a maximum initially, then fall slowly, and finally stabilize, but this behavior is not obvious in the last two explosions. e reason for this phenomenon is that when the deformation of the rock mass is small, the elastic deformation that can be recovered accounts for a large proportion of the whole deformation; with increasing deformation, however, the plastic deformation increases and the proportion of the elastic deformation decreases. Meanwhile, it can be observed that the peak displacement of the anchored cavern is smaller than that of the unanchored cavern, and the reduction ratios of the four explosions are 7.6%, 15.5%, 23.8%, and 26.1%, respectively, indicating that the reinforcement effect of the rock bolts is more obvious when the deformation of the rock mass is larger. Figure 9 shows the peak displacement and scaled distance of measuring point U1 on the vault after each explosion. e scaled distances of the four explosions are 2.84, 2.19, 1.6, and 1 m·kg −1/3 , respectively, and the vault displacement increases with decreasing scaled distance. e slope of the straight line from the second to the fourth explosion is much larger than that from the first to the second explosion, indicating that damage to the vault begins when the scaled distance is 2.19 m·kg −1/3 (the second explosion). Additionally, the slopes of the straight lines with the anchored cavern are smaller, revealing that the deformation of the anchored cavern is slower.

Strain of Cavern Wall.
e stress wave is scattered and diffracted when it propagates to the cavern wall, which may cause stress concentration around the cavern wall. e circumferential strain reflects the magnitude of the   Shock and Vibration circumferential stress. erefore, 10 strain measuring points (ε1-ε10) are arranged around the cavern wall, as shown in Figure 4. Figure 10 presents the peak circumferential strains at the measuring points around the anchored and unanchored cavern wall after each explosion. e values at the measuring points represent the strain magnitude, where a positive value indicates a tensile strain and a negative value indicates a compressive strain. e strain at the measuring points from the vault to the corner is positive first and then negative; that is to say, the circumferential stress from the vault to the corner is first tensile and then compressive. e tensile part is concentrated mainly between the vault and the side arch, and the tensile range increases with successive explosions. Additionally, the circumferential tensile strain decreases gradually from the vault to both sides. Analysis suggests that is because the strength of the tensile wave decreases gradually from the vault to both sides with increasing distance from the explosion source. Meanwhile, the largest circumferential compressive strain is that near the spandrel, this being because the spandrel is located at the junction where stress concentration is prone to occur. With successive explosions, the strain increases more dramatically, not only because of the increasing explosive power but also because of the cumulative damage effect. In several explosions, the circumferential strain of the surrounding rock in the middle of the sidewall is small, indicating that these places are relatively safe when subjected to explosions above the vault, and important equipment can be installed here. Comparing the strain at the same location in the anchored and unanchored caverns shows that the strain in the anchored cavern is significantly smaller than that in the unanchored cavern. Furthermore, the average strain reductions of the 10 measuring points in the four explosions are 11.7%, 14.5%, 22.8%, and 26.6%, respectively. e results indicate that the rock bolts are effective at reducing the strain of the cavern wall, and the greater the strain, the more obvious their reinforcement effect.

Damage Evolution Process of Underground Cavern.
To analyze the damage process of the underground cavern during multiple explosions, the damage of the cavern at different times is shown in Figure 11. Figure 11(a) indicates the damage evolution process of the first explosion. e high temperature and pressure of the explosive gas mean that compression damage occurs near the explosion source at 0.3 ms. At 0.6 ms, the surrounding rock near the ground surface begins to be damaged, mainly because the stress wave propagates to the ground surface and generates tensile stress, and the rock mass has weak resistance to tensile stress, resulting in tensile damage. e closer the damage zone near the ground surface is to the explosion source, the smaller the distribution range. When the stress wave propagates to the free surface of the cavern, tensile damage also occurs.
Subsequently, damage appears on the contact surface between the rock bolts and the surrounding rock, and it is distributed in a strip along the direction of the rock bolts. Conversely, the surrounding rock between the rock bolts is less damaged. Figure 11(b) shows the damage evolution process of the second explosion. After the second explosion, the damage zone around the explosion source expands to a certain extent, and a divergent strip damage zone forms. e damage zone near the ground surface then extends to the left and right sides of the model as well as to the explosion source, and it connects with the damage zone near the explosion source. At 8 ms, the damage degree of the surrounding rock near the vault increases significantly. In terms of the damage range, besides the increased width of the damaged strip between the three rock bolts above the vault, there are also arch-shaped areas of slight damage outside the unanchored zone.
With the occurrence of the third explosion, the damage zone around the explosion source increases remarkably at 8.3 ms, presenting an inverted triangular distribution, and the damage zone narrows from the first explosion source to the third explosion source. e reason for this phenomenon is that the first two explosions bring new free surfaces in the explosion cavity, and when the compressive stress generated by the third explosion propagates to these free surfaces, the tensile wave formed by reflection generates tensile damage. erefore, the damaged area around the first two explosion sources is wider. Because a large number of cracks are produced in the first two explosions, the stress wave generated by the third explosion has lost most of its energy when it propagates to the ground surface, and it is difficult to damage the rock mass near the ground surface. erefore, the damage range near the ground surface does not change greatly; meanwhile, the damage degree of the surrounding rock increases slightly, and the closer to the ground surface, the less the damage increases. At 12 ms, part of the rock mass between the rock bolts has been seriously damaged, which is different from the previous two explosions. e damage degree of the arch-shaped damage areas outside the anchored zone increases surely. Additionally, the spandrel is also seriously damaged.

Shock and Vibration
As shown in Figure 11(d), the distribution of the inverted triangular damage zone is increased further after the fourth explosion, but it extends only to the free surface of the cavern and both sides of the explosion source. e maximum damage between the explosion source and the anchored zone has been connected. Similar to the third explosion, the range of the damage zone near the ground surface shows no obvious change, and the damage degree is increased only slightly. What is different from the previous explosions is that the reinforced area has gone tremendous damage. Although the strength of the rock mass in the reinforced area near the vault increases, its tensile strength is still weak. Furthermore, with decreasing distance from the explosion source and increasing charge mass, the tensile stress intensity of the reinforced area near the vault increases, so the reinforced area has gone tremendous damage. From the perspective of depth, the width of the damaged area from the first explosion source to the reinforced area becomes smaller with increasing depth, but the damage at the junction of the anchored and unanchored zones is serious, and the damage range from the junction down to the vault is getting larger.
In summary, with successive explosions, the damage range and degree around the explosion source increase. e damage degree near the ground surface increases with successive explosions, but the range of the damage zone first increases and then remains unchanged. e damage in the anchored zone is distributed mainly along the direction of the rock bolts in the first two explosions, and the rock mass damage between the rock bolts is very small. Subsequently, with successive explosions, the rock mass between the rock bolts also breaks down. In the whole process, the sidewall and floor are damaged only slightly.

Influence of Multiple Explosions on Explosion Damage
Zone.
ere has been much research into the explosion damage zone of an infinite rock mass after an explosion. Whittaker et al. [38] divided the explosion damage zone into the crushed zone, fracture zone, and elastic vibration zone and concluded that the radius of the fracture zone is 10-15 times that of the charging radius. Moreover, Donze et al. [39] found that the radius of the crushed zone is roughly five times that of the charging radius. In the RHT model, when the damage index D of the rock surrounding the explosion source reaches one, the surrounding rock is in the crushed zone; when D is between zero and one, the surrounding rock is in the fracture zone; and when D is equal to zero, the surrounding rock is in the elastic vibration zone. e analysis of the explosion damage zone requires the rock mass to have an infinite boundary. However, both the ground and the cavern surface are free surfaces in this model. To reduce the influence of the free surfaces, only the surrounding rock on the left and right sides of the explosion Damage index 1.000e + 00 9.000e -01 8.000e -01 7.000e -01 6.000e -01 5.000e -01 4.000e -01 3.000e -01 2.000e -01 1.000e -01 0.000e + 00 t = 12.0 ms t = 9.0 ms t = 8.6ms t = 8.3ms source is analyzed. From the explosion source, a measuring point is set every 1 cm at the left and right sides of the explosion source horizontally, and the average value of the damage index of the two measuring points with the same distance from the explosion source is taken as D at distance R (from the explosion source to the measuring point). Figure 12 shows the relationship between D and R after each explosion. According to the attenuation trend of the curves, the Boltzmann function of the "S" curve is used for fitting. e functional expression relating D and R is as follows: where D 1 is the initial damage, D 2 is the final damage, R 0 is the corresponding value of R when D(D 1 + D 2 )/2, and K is a constant that determines the inclination degree of the function. e smaller the value of K, the steeper the curve. Figure 12 shows that the fitting curve has an obvious reverse "S" shape, and the fitting degree of the four curves reaches 0.99, indicating good fits. e curvature of the fitting curve shows that with increasing distance from the explosion source, the damage of the rock mass decreases slowly initially and then rapidly and finally stabilizes. Comparing the fitting formulas of the four curves shows that K increases with successive explosions; that is to say, the damage curve becomes increasingly gentle, indicating that the attenuation rate of the curves decreases with successive explosions.
In the first explosion, Figure 12 shows that the damage index D can be taken as being one when R is less than 4 cm and zero when R is greater than 25 cm. After the second explosion, D can be taken as one when R is less than 6 cm and zero when R exceeds 35 cm. For the third and fourth explosions, D can be taken as 1 when R is less than 9 and 13 cm, respectively, but it does not stabilize around zero with increasing distance from the explosion source. e charge radii of the four explosion sources are 2.21, 2.77, 3.56, and 4.81 cm, respectively, and the damage index D of the surrounding rock can be taken as being one when the distance from the explosion source is less than two to three times the charge radius, namely, the surrounding rock is in the crushed zone. erefore, multiple explosions do not affect the distribution of the crushed zone. In the first two explosions, the damage index D of the rock mass tends to be stable near zero when the distance from the explosion source exceeds 13 times the charge radius, namely, the rock mass is in the elastic vibration zone. However, in the third and fourth explosions, the surrounding rock is still in the fracture zone when the distance from the explosion source exceeds 21 and 16 times the charge radius, respectively, which indicates that multiple explosions increase the range of the fracture zone.

Law Governing Cumulative Damage of Surrounding
Rock under Multiple Explosions. To reveal the law governing the cumulative damage of the surrounding rock under multiple explosions, only the rock mass outside the crushed zone can be analyzed. e fourth explosion source has the largest charge radius, and the radius of its crushed area is 13 cm. erefore, the damage measuring point D 1 is set at 15 cm below the fourth explosion source, and three measuring points are arranged under D 1 with a spacing of 5 cm, as shown in Figure 4. Figure 13 shows the cumulative damage curves of the four measuring points. In general, the damage index D of each measuring point increases with successive explosions, which suggests that (i) each subsequent explosion further damages the already damaged rock mass and (ii) the damage of the rock mass is irreversible.
However, not every explosion damages the rock mass further; for example, the damage index D of measuring point D 3 does not increase after the third explosion. is verifies that there is a certain energy threshold for damaging the rock mass, and only when the explosion energy exceeds that threshold will the rock mass be damaged further. e four curves show that the damage index D of each measuring point after each explosion does not decrease with increasing distance from the explosion source, mainly because the tensile stress increases with decreasing distance from the vault to measuring points D 3 and D 4 . Meanwhile, the relatively low tensile strength of the rock mass means that the damage of points D 3 and D 4 increases.
Because the surrounding rock is damaged severely when subjected to the tensile stress, more attention should be paid to the tensile zones near the ground surface and the vault. Starting from the ground and vault, a damage measuring point is located every 1 cm on a vertical line extending to the explosion source. e relationship between the damage index D of the measuring point and the distance d from the ground surface and vault after each explosion is shown in Figure 14. Figure 14(a) shows that the damage of the ground surface neither increases nor decreases regularly with decreasing distance from the explosion source. In the four explosions, the values of damage index of the measuring points at 8, 14, 15, and 18 cm from the ground surface are one, while those at the other measuring points are relatively small, indicating that the damage of the rock mass near the ground surface is stratified. In Figure 14(b), the damage of the vault also exhibits obvious stratification. e damage at 3, 10, 19, and 27 cm from the vault is significantly smaller than that of nearby measuring points. An interesting phenomenon in Figures 14(a) and 14(b) is that the complete damage zone with D � 1 is relatively thin during the first explosion but thickens gradually with subsequent explosions. Consequently, the complete damage zone of the vault is connected and the vault is destroyed completely after the fourth explosion. is indicates that the cavern arch is vulnerable to damage, so that damage must be analyzed further. erefore, measuring points D 5 -D 7 (as shown in Figure 4) are arranged at the vault, side arch, and spandrel, respectively, to measure the damage of the arch. In Figure 15, the cumulative damage curves of the four measuring points show that the vault and spandrel are damaged more than is the side arch. e reason for this is that the vault is closer to the explosion source and thus sustains greater cumulative damage. Although the spandrel is relatively far from the explosion source, stress concentration is generated there because it is located at the junction of the arch and the straight wall, and so the spandrel is also damaged more than is the side arch.
To quantify the mechanism for the cumulative damage of the rock mass under multiple explosions, further exploration is conducted with four 100 g TNTexplosions at the location of the first explosion source. When using ANSYS/LS-DYNA, it is difficult to achieve four explosions in the same place, and instead, the load is applied repeatedly on the cavity wall 10 cm from the explosion center (as plotted in Figure 16(a)). e load-time curve shown in Figure 16(b) was obtained from the test [40], and the load was applied at 0, 4, 8, and 12 ms.
Measuring points Da-Df (as shown in Figure 16(a)) are selected to record the damage of the surrounding rock at 15,20,25,30,35, and 40 cm below the explosion source to analyze the law governing the cumulative damage of the rock mass around the explosion source. Figure 17 shows fits of the functional relationship between the damage index D of the surrounding rock and the distance R from the explosion source. Comparing the four fitting functions shows that the reduction coefficient decreases from 1.225 to 1.047, indicating that the damage of the surrounding rock decreases with increasing distance from the explosion source, and the rate of decrease drops with successive explosions. Furthermore, the damage of the rock mass increases with successive explosions, whereas the damage increment decreases.
To reveal the relationship between the damage of the surrounding rock around the explosion source and the explosion number, the relationship between the cumulative damage and the explosion number at each measuring point is analyzed by regression in Figure 18(a), and the variances are all greater than 0.99. e damage of the rock mass increases with successive explosions, but rather than being a simple superposition of the damage of single explosions, the damage exhibits a highly nonlinear relationship with the number of explosions.
Because the vault damage is layered, there is no regularity between the damage of the surrounding rock near the vault and the distance from the explosion source, and so it is not analyzed herein. To explore the relationship between the cumulative damage of the vault and the explosion number, three damage measuring points are arranged at 5, 10, and 15 cm above the vault (as shown in Figure 16(a)). Figure 18(b) shows the relationship between the cumulative damage of the surrounding rock and the explosion number. It can be seen that there is a highly nonlinear relationship  Figure 19 shows the final damage of anchored and unanchored caverns after four explosions, from which the damage zone distribution range and shape are different. e most obvious difference is that the anchored cavern has the damage zone along the direction of the rock bolts, while the unanchored cavern has no such damage. e damage zone of the anchored cavern thickens from the spandrel to the vault, whereas that of the unanchored cavern remains almost uniformly thick. Additionally, compared to the anchored cavern, the unanchored cavern is damaged far greater from the side arch to the spandrel and also to a certain extent on the sidewall.

Comparison of Cumulative Damage between Anchored and Unanchored Caverns.
To quantitatively analyze how the rock bolt reinforcement influences the damage degree, the values of the damage index D for D 5 -D 7 (as shown in Figure 4) are listed in Table 5. e anchored cavern is damaged less in the four explosions than is the unanchored cavern: the damage of the three measuring points of the anchored cavern in the four explosions is 15.9%, 23.2%, 31.4%, and 37.8% less than that of the unanchored cavern on average. Also, the damage reduction of the surrounding rock in the anchored zone increases with successive explosions; this is mainly because the rock bolt reinforcement is not fully exerted when the damage is small, but the rock bolts exert their anchoring effect fully with increasing damage. In general, the rock bolts are effective at improving the antiexplosion ability of the surrounding rock and reducing the damage degree of the anchored zone, and the greater the cumulative damage of the surrounding rock, the more obvious the anchoring effect of the rock bolts.

Conclusions
Based on a similarity model test, a numerical method is utilized to analyze the dynamic responses and cumulative damage of an underground cavern under four explosions above its vault. e accuracy of the simulation results was verified by comparison with the test, and the following useful conclusions are obtained: (1) Of the four explosions, the first acts mainly to compact the surrounding rock, and it is the second that begins the destruction of the vault. Additionally, the cavern arch is destroyed completely after the fourth explosion. (2) From the vault to the corner, the circumferential strain of the cavern wall changes from tensile to compressive, and the maximum circumferential tensile and compressive strains are located at the vault and spandrel, respectively. Moreover, the strain in the middle of the sidewall is relatively small, and important equipment can be installed here. (3) With successive explosions, the distribution range of the crushed zone does not change, but that of the fracture zone increases. (4) e damage of the surrounding rock near the explosion source decreases with increasing distance from the explosion source, and the attenuation rate decreases with successive explosions. By contrast, the damage of the surrounding rock at the vault changes irregularly with increasing distance from the explosion source. e cumulative damage of the surrounding rock near both the explosion source and the vault increases irreversibly with successive explosions, and the relationship between the cumulative damage and the explosion number is nonlinear. (5) e rock bolts are effective at reducing the displacement, strain, and cumulative damage of the surrounding rock in the anchored zone, and their reinforcement effect increases with increasing displacement, strain, and cumulative damage.
Future work should focus on how the supporting methods influence the antiexplosion performance of underground caverns. Also, the damage threshold of the surrounding rock should be explored through experiments.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.