Many metal mines in China will gradually enter the depth range of 1000–2000 m during the 14^th Five Year Plan Period, suggesting that deep mining will become widespread in the coming decades. The hard and brittle rocks at great depths often exhibit significantly different mechanical properties from rocks at shallow depths (Tian et al. 2020; Li et al. 2017). The frequent disasters induced by excavation activities in roadways and caverns restrict the safe and efficient operation of deep underground engineering and cause tremendous casualties and economic losses (Feng et al. 2021a; Wu et al. 2021; Han et al. 2020). In these deep hard metal mines, special failure phenomena such as slabbing, rockburst, buckling, and roof caving are frequently encountered around caverns or chambers where stresses are always concentrated, which potentially affect the stability and increase the cost of rock supporting (Ortlepp and Stacey 1994; Martin and Christiansson 2009; Gong et al. 2019; Kaiser and Cai 2012).

It has been demonstrated that the stress conditions of deep surrounding rocks are significantly affected by self-excavation unloading activities, in which the tangential stress increases and the radial stress decreases under high static stresses. In practice, the surrounding rocks of deep roadways and caverns are inevitably affected by external dynamic disturbances, such as ore collapse, ore caving, mechanical drilling, and blasting of adjacent stopes during the exploitation of solid mineral resources (Li et al. 2011; Yang et al. 2018). The two aforementioned types of disturbances are collectively denoted as “mining effects”. Therefore, the stress state of the deep engineering rock mass can be attributed to “high initial static stress + external dynamic disturbance”, and the analysis and interpretation of rock mechanical behaviors in deep hard rock mining—particularly for rockburst, slabbing, and zonal disintegration—should be focused on the coupled static and dynamic loads in the mining process (Li et al. 2017). There are two types of dynamic loads in deep mining: low-frequency disturbances generated by regional stress adjustments (Carter and Booker 1990; Su et al. 2017) and strong dynamic disturbances induced by blasting loads (Jiang et al. 2021). Considerable effort has been directed toward investigating the failure properties of rocks subjected to these two types of dynamic loads. For example, Du et al. (2016) performed true-triaxial unloading and local dynamic disturbance tests on three rock types. They found that under local dynamically disturbed loading, none of the specimens displayed obvious fracturing at low-amplitude local dynamic loading, while the degree of rock failure increased as the local dynamic loading amplitude increased. Li et al. (2008) developed a modified split Hopkinson pressure bar (SHPB) for a coupling load experiment at medium-to-high strain rate, and investigated the mechanical properties and energy dissipation law of rock under different pre-stresses and impact disturbances.

In addition to the stress conditions in deep mining activities, the presence of geological discontinuities significantly affects the stability of the engineering structure (Chen et al. 2020; Dou et al. 2021; Zhang et al. 2021). A structural plane (mainly refers to a discontinuity joint, hard structural plane, discontinuity, or other structural plane in the rock mass) is a geological discontinuity that cuts the stratum and destroys its continuity and integrity in deep mining. Researchers have investigated the mechanical behaviors and failure properties of caverns with nearby structural planes at the laboratory and engineering scales under such conditions. For example, for the former one, Lin et al. (2021) investigated the mechanical behavior of a jointed rock mass with a circular hole under compression-shear loading by means of DIC (Digital Image Correlation) and DEM (Discrete Element Method) modeling, and the crack evolution process was investigated and four types of crack coalescence among joints were summarized. Chen et al. (2017) studied the fracture evolution characteristics of sandstone containing double fissures and a single circular hole under uniaxial compression. Yang et al. (2019) conducted a series of laboratory testing to investigate the failure behavior and crack evolution mechanism of a non-persistent jointed rock mass containing a circular hole. For the latter one, Feng et al. (2019b) used Elfen software to modeling the hard rock failure induced by structural planes around deep circular tunnels. They found that that the failure intensity of rock tunnel is a function of both dip angles and frictional coefficients of structural planes. Ma et al. (2020) investigated the distance effects of the fault on the rock mass stability of the main powerhouse at the Huanggou pumped-storage power station, and concluded that when the distance is low, the damage of the surrounding rock masses subjected to excavation evolved from the spandrel toward the fault, eventually causing a V-shaped collapse if the strength reduction continues.

To better understand and prevent the aforementioned dynamic disasters, it is necessary to consider the combined effects of the external dynamic disturbances and the geological structure on the failure properties of deep underground structures. For example, Li et al. (2020) studied the failure mechanism of rock with fissure-opening flaws under coupled static-dynamic loads. However, the opening was pre-fabricated during the preparation of specimens, indicating that the excavation unloading effect prior to be impacted by dynamic loads was ignored in their research. Manouchehrian and Cai (2017; 2018) analyzed the rockburst in tunnels subjected to static and dynamic loads, with particular attention paid to the effects of nearby faults. The modeling results confirm that the presence of geological structures in the vicinity of deep excavations could be one of the major influence factors for the occurrence of rockburst. Because the dip angles of fault were not fully considered in their research, the conclusions still needs further validation and evaluation. Moreover, the dynamic stress direction was not taken into account in the numerical simulation.

Numerical simulation is an effective method for analyzing the failure behaviors of rock materials subjected to complex stress conditions, stress paths, and geological structures at the engineering scale (Zhou et al. 2022). In the present study, the combined finite–discrete element method (FDEM) was adopted to simulate the failure process of a deep underground cavern. This method can reveal the failure properties of deep hard rock, which are affected by key factors such as the dynamic stress-wave amplitudes, disturbance direction, and dip angles of the structural plane. Typical hard marble extracted from Hongling Lead–Zinc mine in Inner Mongolia, China was selected for the numerical simulation, and the stress-field distribution, displacement, velocity of failed rock, and failure zone around the circular cavern were analyzed. This study provides insight into failure mechanism of hard surrounding rock, which is affected by the external dynamic disturbances (mainly refers to blasting loads), excavation unloading effects, and geological conditions in deep mining engineering.

The FDEM is a numerical method that combines continuum mechanics principles with dis-continuum algorithms to simulate multiple interacting deformable solids (Li et al. 2018; Feng et al. 2019c). In the current study, Elfen software, developed by Rockfield Software, was adopted to reproduce entire failure process of an underground cavern. This method can solve two-dimensional (2D) and 3D models problems implicitly and explicitly, and is able to simulate the rock failure process under static and dynamic loads, which shows great advantage in the numerical simulation. Elfen has a large number of constitutive models for selection. Among these models, the Mohr–Coulomb constitutive model with rotating crack was commonly selected because it can simulate rock fracturing due to both tension and compression. A detailed description of the fundamental principles of the FDEM approach and applicability of constitutive model on the basis of Elfen can be found in previous literature (Li et al. 2018; Feng et al. 2019a; Cai 2013).

2.1 Numerical simulation parameter validation

A uniaxial compressive test based on the FDEM was conducted to validate the simulation with laboratory testing. The objectives were to calibrate the selected parameters. For 2D simulations, plane strain analysis assumes a unit thickness of 1 mm in current laboratory testing. The dimensions of the cylindrical specimen used in the uniaxial test were 100 mm × Φ50 mm. The basic physical and mechanical and discrete contact parameters used in this study are presented in Table 1. For details regarding the selection of the discrete contact parameters, the reader is referred to the works of Feng et al. (2019b) and Chen et al. (2020).

There are five basic material parameters in the Mohr–Coulomb constitutive model: cohesion C, frictional angle φ, dilation angle ψ, tensile strength σ_t, and fracture energy G_f. The residual cohesion and frictional angle are calculated in accordance with CWSF (cohesion weakening and frictional strengthening) model proposed by Hajiabdolmajid et al. (2002). The relationship between fracture toughness in mode I and fracture energy can be defined as Mitelman and Elmo (2014):

where K_IC (MPa/m^−1/2) denotes the fracture toughness in mode I, and E (GPa) is the Young’s modulus. In the current study, the fracture toughness K_IC can be estimated using an empirical relationship with the rock tensile strength σ_t  (Mitelman and Elmo 2014; Zhang 2002):

Based on above equations, the fracture energy is calculated to be about 5 N/m.

Figure 1 presents the stress–strain curves and final failure modes of marble under the uniaxial compressive test based on the numerical simulation and laboratory testing. Macro-shear bands coupled with longitudinal splitting cracks appeared in the specimen, which is a typical hard/brittle property of marble. The final failure mode was tension with shear. Observation and comparison of the crack distribution within the marble specimen revealed good agreement between the two cases.

Typical stress–strain curves and failure modes of marble specimen under uniaxial compressive test based on numerical simulation and laboratory testing

The stress–strain curves differed slightly with regard to the total strain and pre-peak stage. This is attributed to the compaction and closure of primitive cracks within the rock specimens, which led to larger deformation and strain values in the early stage compared with the result of the numerical simulation. This phenomenon could not be reproduced in the current numerical simulation analysis, because the micro-pore structure was not fully considered when specifying material parameters in the software. Nevertheless, the typical hard-brittle property in which the yield deformation stage representative of a nonlinear segment (plastic strain) was extremely limited can be observed in Fig. 1. Moreover, the peak strengths for the two cases were almost identical (approximately 70 MPa), indicating that the selected material parameters and discrete contact parameters can accurately characterize the failure property of marble specimens in numerical modeling. The results confirm the applicability and accuracy of the FDEM.

2.2 Model setup and boundary conditions

The established 2D models are shown in Fig. 2. The rock mass is represented by a single face, and the tunnel is represented by a circular cavern. The external boundary geometry length of the model was 200 m, and the cavern diameter was 20 m.

Geometry size and established model in the current study

In the present study, models without a structural plane (Fig. 2b) and with a structural plane (Fig. 2c) were considered. For the case of a nearby structural plane with a constant length of 20 m, four dip angles were modeled: 0º, 30º, 60º, and 90º (angle between the horizontal direction and the structural plane). The total number of unstructured triangular elements established by the model was 31,574. Dense meshes with a length of 1 m were arranged around the central cavern, and the meshes were sparser (length of 2 m) for the remaining regions. The objective was to reduce the computational burden to the greatest extent possible (Feng et al. 2021b). The bottom boundary was fixed in the vertical direction, and the top boundary and both sides were unconstrained.

2.3 Stress path

In current numerical modeling, face loading was employed to apply geostatic stresses at the model boundaries. A vertical stress was applied to the top domain, and horizontal stresses were applied on both sides of the model. To avoid unexpected fracturing within the model due to instantaneous loading, the loading was conducted linearly, with a total duration of 1 s. The in-situ stress remained constant until the end of the numerical simulation. To reproduce the excavation activities, a linear relaxation path over 1 s was formulated for the central cavern. When the cavern was deactivated, the relaxation was automatically applied to release the stresses according to the specified curve, as shown in Fig. 3a. The monitored in situ stress fields of the Hongling Lead–Zinc mine were calculated as follows (Feng et al. 2019b):

where σ_v represents the vertical stress, and σ_h max and σ_h min represent the maximum horizontal principal stress and minimum horizontal principal stress, respectively. In this study, a circular cavern at a burial depth of approximately 1100 m was simulated to reflect the higher geostatic stress field. The h in Eq. (3) was regarded as 1000 m because there was a depth of 100 m between the top domain and the cavern roof. Therefore, the vertical and horizontal stress components were σ_v = 27.9 MPa and σ_h max = 35.8 MPa (σ_h min is not considered in 2D simulations), respectively.

Stress path applied in the model. a Loading and unloading path; b Triangular load function. Note: σ_s indicates the initial static stress, and σ_d means the dynamic stress

External dynamic disturbances can be produced by various sources, such as explosions, vibrations, earthquakes, ore drawing, and adjacent mining activities (Manouchehrian and Cai 2017). In this study, the dynamic disturbances were assumed to be generated by adjacent blasting loads, and the dynamic loads were applied to the model boundary via face loading (see Fig. 4), where the stress path was designed on the basis of the literature (Yang et al. 2020). The dynamic disturbances imposed on the top domain and right boundary were modeled separately. A triangular blast load function, which included a rising stage (1 ms) and descending stage (2 ms), was selected to reflect the dynamic stress-wave induced by the blasting operation. The load function of the dynamic stress-wave is shown in Fig. 3b. An important parameter during the analysis of dynamics is the damping value. A numerical viscous damping (i.e., proportional to velocity) is introduced in the FDEM simulation. The value of damping prescribed is defined relative to the automatically estimated lowest frequency of vibration for the application and the damping may be prescribed differently for each freedom of the node. The factor ranges from 0 (i.e. no damping) to 1 (i.e. critical damping). Actually, a damping value of 0.1 or less should be used for dynamic relaxation simulation of linear problems (Rockfield 2013), and a value of 0.015 is recommended in the current numerical simulation. The influencing factors considered in this model include: (1) stress wave amplitude: the indicating values are 20, 40, 60, 80 and 120 MPa. (2) disturbance direction. (3) dip angle of structural plane which may be 0º, 30º, 60ºor 90º. (4) Relative location between structural plane and dynamic disturbance source: (1) same side: dynamic direction from vertical direction; (2) opposite side: dynamic direction from hozirontal direction; (3) Same and opposite side: dynamic direction from both vertical and horizontal direction.

Combined static and dynamic loads in the model. a Model without structural plane; b Model with structural plane. Note: σ_vs indicates the vertical static stress, and σ_hs indicates the horizontal static stress

The failure around circular cavern caused by the excavation process was examined. The objectives were to validate the accuracy of numerical results and to perform a comparison with the subsequent simulation results obtained with consideration of either the external dynamic disturbance or the structural-plane effect. This is important, as in several previous studies, the excavation unloading process was ignored during the analysis of the cracking behavior and mechanism of deep caverns influenced by various factors, i.e., the stress conditions, geological discontinuity, and stress paths.

Figure 5a plots the results of the circular cavern induced by excavation unloading effect after it was entirely unloaded according to the linear relaxation function. As shown in Fig. 5a, the stress redistribution (maximum principal stress contours) around the cavern did not reach the model boundaries, confirming the rationality of previously selected model size.

Comparison of results of the circular cavern induced by excavation unloading effect. a Maximum principal stress contours derived from numerical simulation; b Crack distribution derived from numerical simulation; c Rockburst simulation experiment on circular hole sidewalls (Gong et al. 2019); d Final shape of v-shaped breakout notches in the crown and invert of the tunnel (Read 2004)

As shown in Fig. 5b, several cracks appeared around the circular cavern. The current results exhibit discrepancies compared with those of a previous experimental study (Li et al. 2020; Zeng et al. 2018), because the hole or cavern was prefabricated manually prior to testing. In the previous study, unloading-induced cracks were not generated, because the rock was fractured by persistent loading, which was attributed to the limited experimental techniques used for reproducing the unloading process. Generally, the failure behavior around deep tunnels and chambers is induced by the excavation unloading effect rather than monotonically compressive loads (Feng et al. 2019a). This may lead to different crack-propagation patterns and locations of the failure zone between two stress paths.

In the current study, the typical slabbing failure was observed around the circular cavern (Fig. 5b). Due to the unloading effect, the rock near the free surface fails in tension. During the development of V-shaped notch, the surrounding rock near the free surface is mainly damaged by tension failure. By comparison, the rock failure mode in the current simulation is found to be consistent with those of V-shaped notch formed on the surrounding rocks in the field and rockburst simulation test (Gong et al. 2019; Read 2004) (Fig. 5c, d).

4.1 Stress-wave propagation in model with circular cavern

To analyze the stress-wave propagation law for the established model, the maximum principal stress contours for cases of an amplitude of 80 MPa were examined, as shown in Fig. 6. The current initial time corresponds to the beginning of dynamic disturbance (t = 2.5 s). As indicated by Fig. 6, the stress-wave exhibited significant attenuation characteristics as the stress propagated from the top domain to the bottom domain. The superimposed maximum principal stress value was attenuated from 101 MPa (corresponding to t = 2.1 ms) to approximately 45 MPa (corresponding to t = 39.57 ms). When the stress-wave was transmitted to the cavern boundaries, large number of visible cracks were generated around the roof (after t = 25.22 ms). These cracks were mainly induced by the wave reflection when it encountered the free face. Therefore, the roof was fractured by tensile failure instead of compression failure. The failure intensities and damage ranges on both sides were not as severe as those around the roof—particularly for the case of the floor (almost no change occurred compared with the case of t = 2.1 ms). This is primarily attributed to the persistent wave attenuation, which reduced the dynamic disturbance, although the surrounding rocks in these areas were still subjected to the tension effect. The numerical results indicated that the mechanical response and failure characteristics of the cavern exhibited typical temporal and spatial changes when subjected to the external dynamic disturbance. For verification, the radial stress distribution at different locations was also examined.

Dynamic stress-wave propagation process in the numerical model. Note: The dynamic stress-wave was imposed on top domain, and the amplitude is 80 MPa. The initial time starts from 2.5 s

Several monitoring points were selected to study the radial stress evolution laws in the roof, floor, and right side of the circular cavern, as shown in Fig. 7. As shown in Fig. 7a (cavern roof), the radial stresses at points 1 and 2 did not change over time, remaining at approximately 0. This indicates that the surrounding rock mass was severely damaged by the excavation unloading prior to the external dynamic disturbance. The radial stress evolution laws differed among points 3, 4, and 5. The radial stresses at these points prior to t = 2.5 s remained at approximately 20 MPa. After t = 2.52 s, the stress values suddenly increased to approximately 50 MPa. Subsequently, the radial stress values rapidly decreased, particularly for points 3 and 4 (they reached almost zero). This indicates that the rock mass was fractured by the reflected tensile wave, which was expected, because the tensile strength of the rock was significantly lower than the compressive strength. Additionally, the peak stress values for monitoring points 3, 4, and 5 were 60, 46, 44 MPa, respectively. These results reflect the attenuation of the wave along the propagation route.

Radial stress-time evolution curves for corresponding monitoring points in different locations of cavern. Note: the negative value indicates that the rock is subjected to compressive stress condition

Figure 7b shows the radial stress variation trend for the cavern floor. The radial stress at point 6 was always zero, similar to the case of points 1 and 2, indicating that these rock masses were damaged prior to the dynamic disturbance. For the remaining four points (7, 8, 9, and 10), although the radial stress fluctuated owing to the dynamic stress-wave after t = 2.5 s, the amplitude was smaller than those for points 3, 4, and 5. This is attributed to the wave attenuation, which weakened the dynamic response of the rock mass in these areas. Moreover, the radial stress did not decrease after the fluctuation and remained close to its initial value, indicating that these areas were not fractured by the external dynamic disturbance. This also reveals why the failure intensity was lower for the floor than for the roof.

The radial stress-time evolution law for the right sidewall shows similar variation trend compared with the case of roof and floor (Fig. 7c). For rock masses farther from the cavern boundaries, the original radial stress value was higher, indicating that the rock masses were not damaged before the external dynamic disturbances. In contrast, the radial stress for the rock masses near the roadway boundaries (point 11) was always zero. The radial stress increased monotonically with an increase in the distance from the cavern boundaries prior to the dynamic disturbance. Similarly, the stress amplification increased monotonically with an increase in the distance from the cavern boundaries. The stress amplifications for points 12, 13, and 14 were 3.5, 6, and 9 MPa, respectively. This phenomenon was not determined by the wave attenuation (same propagation distance); rather, it was determined mainly by the wave reflection and scattering. When the rock mass was close to the cavern boundaries (free surface), the reflection and scattering inhibited the propagation of the stress-wave. Therefore, the stress amplifications for points 14 and 13 exceeded that for point 12. However, this does not imply that a point with larger stress amplification was more prone to fracture, as points 14 and 13 remained in compressive stress states. This is why damage due to the external dynamic disturbance occurred at point 12 (the radial stress at point 12 decreased to 0, while it remained at 10 and 15 MPa for points 13 and 14, respectively).

The vibration velocity–time evolution curves for the different positions of cavern were plotted, as shown in Fig. 8. The maximum velocities for points A and C (approximately 1.5 m/s) were significantly higher than that for point B (0.75 m/s) during the excavation unloading process. After reaching its maximum value, the velocity rapidly decreased to 0 until the occurrence of the dynamic disturbance. When the stress-wave arrived at the cavern boundaries, the vibration velocities at three monitoring points (roof, floor, and right side of cavern) all increased—particularly at points A and B (v_A > v_B > v_C). This variation trend differed from that of the excavation unloading stage (v_C > v_A > v_B). When the stress-wave reached point C, the stress amplification was significantly reduced after the wave had passed through a series of free surfaces and long propagation route, which led to a smaller increase of vibration velocity for the floor due to the external dynamic disturbance.

Vibration velocity–time evolution curves for three monitoring points in different locations of cavern (dynamic disturbances in vertical direction)

4.2 Amplitude of dynamic load

The effects of the amplitude of the dynamic load on the failure characteristics of an excavated circular cavern were investigated. The external dynamic stress-wave was imposed on the top domain. Stress amplitudes of 0, 20, 40, 60, 80, and 120 MPa were considered in the analysis. Figure 9 presents the final failure modes of the circular cavern subjected to various stress-wave amplitudes. As shown, the failure intensity and damage range increased monotonically with an increase in the stress-wave amplitude σ_d. When the stress-wave amplitude was 20 MPa, less newly generated cracks were observed around the circular cavern compared with the case where there was no dynamic disturbance (σ_d = 0). Because of the wave attenuation and the low stress-wave amplitude, the dynamic response was weak when the wave reached the cavern boundaries. When the stress-wave amplitude was increased to 40 MPa, several visible cracks were generated around the surrounding rock of the circular cavern—particularly near the roof. Under this condition, although the rock mass was fractured by the reflected tensile wave, the extent of this influence remained at a limited level. As the stress-wave amplitude continued to increase, the failure intensity and damage range of the surrounding rock gradually increased, and visible cracks emerged and propagated around the cavern roof. When σ_d was greater than 60 MPa, newly generated cracks were observed in the floor and both sides of the circular cavern (Fig. 9e, f). Under this condition, the large stress-wave amplitude led to severe fracture of the rock mass in these areas, even though wave attenuation occurred. Additionally, the severe fracture and persistent propagation of cracks near the cavern roof caused stress redistribution and transference when the dynamic stress-wave amplitude was large. This leads to further increase in the tangential stresses around both sides and the floor. Therefore, the rock masses in these areas were fractured.

Final failure modes of circular cavern subjected to various stress-wave amplitudes. Note: σ_dy represents that the external dynamic disturbances is generated from top domain

Figure 10 presents the maximum kinetic energy per unit area (kJ/m²) and fracture damage zone distribution with respect to the dynamic stress-wave amplitude. Both these indicators increased with the dynamic stress-wave amplitude. When σ_d reached 40 MPa, the maximum kinetic energy per unit and the size of the fracture damage zone were significantly reduced. For example, the maximum kinetic energy per unit increased from 1.63 at σ_d = 40 MPa to 3.63 at σ_d = 60 MPa, and the size of the fracture damage zone increased from 228 m² at σ_d = 40 MPa to approximately 304 m² at σ_d = 60 MPa. The overall variation trends were consistent with those for the final failure mode and crack-propagation pattern, verifying the accuracy of the numerical simulation results.

Maximum kinetic energy per unit area and fracture damage zone distribution as functions of dynamic stress-wave amplitude (dynamic disturbances in vertical direction)

4.3 Direction of dynamic disturbance

The effect of the direction of the dynamic disturbance on the cavern failure was analyzed. For comparison with the numerical results presented in Sects. 4.1 and 4.2 (stress imposed from the vertical direction), the dynamic stress-wave was generated from the right domain (imposed from the horizontal direction) of the model. Figure 11 presents the maximum principal stress contours for cases of an amplitude of 80 MPa. Similar to the case of the vertical direction, the wave was attenuated, as indicated by an increase in the spread distance. The superimposed maximum principal stress value was attenuated from 111 MPa at t = 2.1 ms to approximately 52 MPa at t = 39.57 ms. However, when the wave was transmitted near the cavern, visible cracks and obvious damage occurred on the right side, whereas almost no obvious damage occurred at the other positions of the cavern (particularly near the left side). Tao et al. (2018) investigated the failure characteristics and stress distribution of a pre-stressed rock specimen with a circular cavity subjected to dynamic loading. In their study, the stress-wave propagation through the tunnel and the deformation and failure characteristics of the tunnel were reproduced using ANSYS software. They found that the plastic deformation was primarily distributed at the sides—particularly the side where the dynamic stress-wave was incident, as well as the roof and floor, where several points with slight deformation were observed. Their conclusions are consistent with our results, validating our findings.

Dynamic stress-wave propagation process in the model. Note: The dynamic stress-wave was imposed on right domain, and the amplitude is 80 MPa. The initial time starts from 2.5 s

The vibration velocity–time evolution curves for three monitoring points at different positions of the cavern (dynamic disturbance in the horizontal direction) are presented in Fig. 12. The variation trend after the dynamic disturbance was analyzed. The maximum velocity amplifications for the points decreased in the following order: v_B > v_C > v_D. The main reason for this is consistent with that in the situation of stress imposition in the vertical direction (stress-wave attenuation and reflection) and will not be discussed here. Additionally, the maximum vibration velocity in the horizontal direction was almost identical to that in the vertical direction (2.6 and 2.9 m/s, respectively). This indicates that the direction of the dynamic disturbance hardly affected the failure intensity of the deep cavern.

Vibration velocity–time evolution curves for three monitoring points in different locations of cavern (dynamic disturbances in horizontal direction)

Figure 13 shows a plot of the final failure modes of the circular cavern subjected to different stress-wave amplitudes (in the horizontal direction). The failure intensity and damage range increased monotonically with the increasing stress-wave amplitude. This trend is consistent with the results for stress imposition in the vertical direction (Fig. 9). However, the main damage range changed from the roof to the right side of the circular cavern when the disturbance source was from the top domain to the right domain. We also found that the differences in the overall failure intensity around the cavern were insignificant for identical stress-wave amplitudes.

Final failure modes of circular cavern subjected to various stress-wave amplitudes. Note: σ_dx represents that the external dynamic disturbances is generated from right domain

To validate the aforementioned conclusions, the maximum kinetic energy per unit and the fracture damage zone for the two disturbance directions were plotted, as shown in Fig. 14. The maximum kinetic energy was always identical for same stress-wave amplitude and different disturbance directions, similar to the fracture damage zone. This was particularly true when the dynamic stress-wave amplitude was less than 60 MPa. Du et al. investigated the failure behaviors of different rock types (granite, red sandstone, and cement mortar) using a novel testing system coupled to true-triaxial static loads and local dynamic disturbances. The experimental results indicated that the fracture process did not depend on the direction of local dynamic loading, and the degree of rock failure increased as the local dynamic loading amplitude increased. This indicates that the direction of the dynamic load had little effect on the development of rockburst. The experimental results of Du et al. (2016) are consistent with our numerical results. Additionally, we found that the fracture intensity was higher for stress imposition in the vertical direction than for stress imposition in the horizontal direction. This is attributed to the initial geostatic stress states. Because more visible cracks were generated in the roof and floor than in both sides of the cavern during the excavation unloading stage, the final failure intensity became more severe in these areas after the dynamic disturbance.

Maximum kinetic energy per unit and fracture damage zone distribution as functions of dynamic stress-wave amplitudes. Note: y direction indicates vertical direction, and x direction indicates horizontal direction)

Figure 15 shows the displacement distribution for different monitoring points around the cavern with respect to the stress-wave amplitude (in both the vertical and horizontal directions). The maximum displacement value for point A under a vertical disturbance was 0.1729 m, and the maximum displacement value for point B under a horizontal disturbance was 0.1644 m. The maximum displacement value for point C under a vertical disturbance was 0.0259 m, and the maximum displacement value for point D under a horizontal disturbance was 0.0255 m. These results indicate that for the same stress-wave amplitude, the relative location between the disturbance direction and the maximum tangential stress had little effect on the damage degree of the circular cavern in the specific static stress environment. These conclusions are based on the small difference between the vertical and horizontal static stresses (27.9 and 35.8 MPa, respectively). Therefore, the applicability of the foregoing conclusions to cases with large differences between the maximum and the minimum principal stresses must be further analyzed.

Displacement distribution for different monitoring points around cavern as functions of stress-wave amplitudes in vertical and horizontal direction

Previous findings have indicated that geological discontinuities and weakened planes significantly affect the mechanical properties of deep underground rock masses (Feng et al. 2019b; Manouchehrian and Cai 2017, 2018; Zhou et al. 2015; Yang and Zhang 2020). The effects of the structural plane and external dynamic disturbance on the dynamic response of the deep circular cavern are investigated in this section.

5.1 Dip angle of structural plane

Four different dip angles of the structural plane were considered in the analysis: 0º, 30º, 60º, and 90º. The friction coefficient of the structural plane was set as 0.1. The dynamic stress-wave was imposed on the top domain, with an amplitude of 60 MPa. The failure characteristics of the circular cavern with a nearby structural plane before and after the external dynamic disturbance (vertical direction) were analyzed, as shown in Fig. 16. Before the cavern was affected by the dynamic stress-wave, the surrounding rock was fractured, and the failure behavior was determined by the excavation unloading process and the structural-plane effect. As shown in Fig. 16a, the failure intensity and damage degree first increased and then decreased. When the dip angle θ was 0º, visible cracks were generated in the roof. Careful observation revealed that these cracks were merely distributed between the cavern boundaries and structural plane, and the failure mode resembled a pan shape. This result is consistent with Zhou’s description (2015). When θ was 30º, numerous cracks were formed near the structural plane, and new cracks appeared at the tip of the structural plane. Analogous slabbing failure was observed around the cavern boundary, which is similar to the previous findings (Feng et al. 2019c; Sinha and Walton 2019; Cho et al. 2010). When θ increased to 60º, the failure intensity and crack propagation became more severe. Under this condition, the rock masses on both sides and the floor of the cavern were fractured violently. As reported by Feng et al. (2019b), gradual slabbing failure may contribute to the stress concentration and dynamic response around the structural plane; the activated structural plane begins to shear and slip and may in turn cause slabbing failure near excavation boundaries. Therefore, the extensive and violent coupled reaction may lead to a serious failure and dynamic response, even rockburst. When θ was 90º, the failure intensity was reduced instead—particularly near the cavern roof.

Final failure mode of circular cavern with a nearby structural plane before and after being subjected to external dynamic disturbances in vertical direction

After experiencing the external dynamic disturbance, the circular cavern was obviously aggravated, with more visible cracks generated in the surrounding rock, regardless of the dip angle. The overall variation trends of the failure intensity and damage range with respect to the dip angle were similar to the previous findings for the excavation unloading process (the most serious situation corresponded to θ = 60º). In contrast, at θ = 0º, no visible cracks were generated on the sides of the cavern prior to the external dynamic disturbance. This may have been due to the structural plane, which inhibited the stress propagation around the cavern wall and finally led to a relatively stable state in these areas. However, when the cavern was subjected to a dynamic stress-wave, cracks were generated in both sidewalls. This appeared to be due to the consistent fracturing at the cavern roof induced by the stress-wave, which led to stress redistribution around the sidewalls. Moreover, the cracks in the roof continued to coalesce and propagate toward both sides of the cavern. The analysis indicated that the failure of the cavern was affected by the excavation unloading, structural plane, and external dynamic disturbance (Table 2).

Figure 17 shows the fracture damage zone distribution with respect to the dip angle of the structural plane under an external dynamic disturbance in the vertical direction. The red line indicates the case without a structural plane. As indicated by Fig. 17, the damage areas for dip angles of 0º, 30º, and 60º were larger than that in the case without a structural plane, whereas for a dip angle of 90º, the point was below the red line. This indicates that the nearby structural plane had both favorable and unfavorable effects on the failure behavior of the cavern. When the dip angle was 90º, the structural plane inhibited the fracture of the surrounding rock, for the following two reasons: (1) The shear and slip failure along the structural plane were not activated by the excavation unloading or external dynamic disturbance, because the tangential stress was normal to the structural plane; thus, the mutual interaction of slabbing and shearing did not cause severe dynamic fracture. (2) The presence of the structural plane along the radial direction may have reduced the integrity of the thin rock plates around the roof.

Fracture damage zone distribution as functions of dip angle of structural plane subjected to vertical dynamic disturbances

The displacement–time evolution laws for four monitoring points around the circular cavern were analyzed, as shown in Fig. 18. Regardless of the location of the monitoring point, the displacement was always maximized when the dip angle of the structural plane was 60°, whereas at different dip angles, the displacement distribution was disorganized. This indicates that the failure intensity and damage range at θ = 60° were significantly more severe than those at the other dip angles. Moreover, all the displacement values for monitoring point C decreased (to varying degrees) at approximately t = 2 s and subsequently increased. Because point C was located on the cavern floor, the movement direction of the rock mass was opposite to the direction of the displacement generated by the initial geostatic stress during the excavation and unloading.

Displacement–time evolution curves for four monitoring points at different dip angles of structural plane. Note: Point A is located in roof, Point B is located in right side, Point C is located in floor, and Point D is located in left side

5.2 Relative location between structural plane and dynamic-disturbance source

The effect of the relative location between the structural plane and the dynamic-disturbance source (disturbance direction of the stress-wave) on the failure behavior of the circular cavern with a nearby structural plane was investigated. The dip angle of the structural plane was set as 30°, with a constant friction coefficient of 0.1. Figure 19 shows the failure process of the circular cavern after the dynamic disturbance.

Final failure mode of circular cavern with a nearby structural plane at a dip angle of 30° subjected to external dynamic disturbances with various directions. Note: The initial time in this figure starts from t = 2.5 s

For the vertical disturbance direction (structural plane and dynamic-disturbance source were on the same side), the cracks were initiated near the cavern boundaries and structural-plane tip at t = 25.89 ms and gradually propagated toward the cavern roof. At t = 45.33 ms, a few newly generated cracks appeared on the sides and the floor. The failure intensity was higher for the roof than for other locations.

For the horizontal disturbance direction (structural plane and dynamic-disturbance source were on opposite sides), the cracks were initiated on the right side of the cavern and at the structural-plane tip (t = 25.88 ms). At t = 35.35 ms, more cracks appeared around the structural plane and gradually coalesced with the cracks on the right side of the cavern. Moreover, cracks were generated around the roof. Because the structural plane was near the cavern roof, it was further activated by the dynamic stress-wave, which led to the instability of the rock mass in these areas (t = 45.33 ms).

When dynamic disturbances were imposed on both the top and right domains (structural plane and dynamic-disturbance source were on opposite sides), the failure intensity and damage range were more severe and violent than those in the two aforementioned situations. At t = 25.88 ms, the cracks were generated almost simultaneously in the roof and the right side of the cavern. Over time, the stress-waves propagated further and coalesced. The right side of the cavern, the roof, and the floor continued to deteriorate (t = 50.15 ms)—particularly near the structural plane. Under this condition, the rock mass near the structural plane was affected by both the external dynamic disturbances (tensile failure) from the two directions and the structural plane (shear failure) simultaneously. The failure intensity of the left side of the cavern was lower than those at the other positions.

The effects of the relative location between the structural plane and the dynamic-disturbance source on the displacement and vibration velocity distribution at selected monitoring points around the circular cavern were also investigated, as shown in Fig. 20. The displacement and vibration velocity were always the highest in the case of the double disturbance, regardless of the locations of points A and B. For point A, the displacement and velocity induced by the dynamic stress-wave from the vertical direction exceeded those induced by the dynamic stress-wave from the horizontal direction. This is explained as follows: the longer distance between the disturbance source and point A led to greater stress-wave attenuation for the case of the horizontal direction. For point B, the displacement and velocity induced by the dynamic stress-wave from the horizontal direction exceeded those induced by the dynamic stress-wave from the vertical direction. Additionally, another factor led to the smaller values for the case of the vertical direction: the structural plane restrained the propagation of the stress-wave toward point B, resulting in the reflection of the stress-wave. This indicates that the relative location between the structural plane and the dynamic-disturbance source significantly affected the failure intensity and damage range of the deep circular cavern.

Displacement and vibration velocity–time evolution curves for two monitoring points around circular cavern with a nearby structural plane when subjected to external dynamic disturbances. Note: y direction indicates the vertical direction, x direction indicates the horizontal direction, direction indicates both vertical and horizontal direction. Point A is located in roof, Point B is located in right side

The displacement–time history curves for different monitoring points on the structural plane under different disturbance directions were analyzed, as shown in Fig. 21. The displacements of the three monitoring points always decreased in the order of Dis_E > Dis_F > Dis_G, regardless of the disturbance direction. During the excavation unloading process (before the dynamic stress-wave was induced), the displacement was larger closer to the cavern roof. At this time, the structural plane near the cavern boundary was prone to activation and instability, and the failure behavior was determined by the excavation unloading and the structural-plane effect. After the external dynamic disturbance, because the structural plane near the roadway boundary was further activated and dislocated by the stress-wave, the failure intensities in these areas were significantly increased. Additionally, the displacements of the monitoring points were always larger for the double disturbance than for the other two cases.

Displacement velocity–time evolution curves for three monitoring points in structural plane when subjected to external dynamic disturbances. Note: y direction indicates the vertical direction, x direction indicates the horizontal direction, xy direction indicates both vertical and horizontal direction

For point E, the displacement in the case of the vertical direction exceeded that for the horizontal direction (Fig. 21). This is attributed to the longer distance between point E and the right domain, which led to greater stress-wave attenuation, for the case of the horizontal direction. Furthermore, the displacements of points F and G exhibited no obvious distinction in these two situations. Because these two monitoring points were farther from the cavern boundaries than point E, the rock masses in these areas were hardly affected by the structural-plane activation, excavation unloading, and external dynamic disturbance effect.

The effects of external dynamic disturbances and the structural plane on rock fracturing around an excavated deep underground cavern were comprehensively investigated by using a combined finite-element method. The dynamic stress-wave amplitudes, disturbance direction, and dip angle of the structural plane were considered as key factors. A typical hard rock (marble) extracted from Hongling Lead–Zinc mine in China was selected for the numerical simulation, and the stress-field distribution, displacement, velocity of failed rock, and failure zone around the circular cavern were comprehensively analyzed. The results confirmed that the dynamic response and failure properties of the deep cavern were affected by the excavation unloading, external dynamic disturbance, and structural plane. On the basis of a parametric analysis, the following conclusions are drawn.

1. (1)

   The mechanical response and failure characteristics of the cavern exhibited typical temporal and spatial differences under an external dynamic disturbance, which is primarily attributed to the stress-wave reflection that occurred when the wave was confronted with free surfaces and underwent persistent attenuation along the extended propagation route. The cracks formed by the external dynamic disturbance were tensile cracks.

2. (2)

   The failure intensity and damage range around the deep cavern increased monotonically with the increasing stress-wave amplitude. When the stress-wave amplitude was ≥ 60 MPa, the failure intensity was significantly increased. Under this condition, the sufficient stress-wave amplitude led to severe fracture of the rock mass in these areas, although wave attenuation occurred.

3. (3)

   The direction of the dynamic disturbance affected the failure position of the cavern (mainly in the roof and right side when the direction was vertical and horizontal, respectively) but hardly affected the failure intensity. The fracture intensity for the case of the vertical direction was higher than that for the horizontal direction, which is attributed to the specific initial geostatic stress states and excavation unloading disturbance prior to the external dynamic disturbance.

4. (4)

   With the increasing dip angle of the structural plane, the failure intensity and damage degree first increased and then decreased, both before (excavation unloading stage) and after the dynamic disturbance. This indicates that the nearby structural plane had both favorable and unfavorable effects on the failure behavior of the cavern. Moreover, the relative location between the structural plane and the dynamic-disturbance source significantly affected the failure intensity and damage range of the deep circular cavern.