Study on Rock Damage Mechanism for Lateral Blasting under High In Situ Stresses

A 3D numerical model was presented to investigate the blast-induced damage characteristics of highly stressed rock mass. The RHT (Riedel, Hiermaier, and Thoma) model in LS-DYNA was used to simulate the blast-induced damage and its parameters were calibrated by a physical model test. Based on the calibrated numerical model, the influences of confining pressure and free surface span on the blast-induced damage characteristics were investigated. The results show that under uniaxial loading, the crater volume increases with confining pressure increases. The uniaxial static load can change the optimal burden and the critical embedding depth of charge. In stressed rock, the variation law of the crater shape affected by radial tensile fractures is opposite to that affected by reflected tensile fractures. Under the biaxial static load, the crater volume of the borehole placed on the side of the max static load is greater than the other side. The explosion crater can be improved by increasing the free surface span on the same side. Finally, it is suggested that the blasting efficiency can be improved by preferentially detonating the charge on the side of the max static load, and then the charge on the other side can be detonated with a wider free surface span.


Introduction
With the increase of excavation depth, the in situ stress increases gradually and plays an increasingly important role in the rock breaking by blasting. Due to the existence of in situ stress, blast-induced damage characteristics are different from those in surface and shallow sub-surface rock blasting, especially considering the role of lateral free surface in production blasting. In order to determine the parameters of borehole layout and blasting parameters in highly stressed rock mass as well as for safe and efficient production, it is necessary to investigate the blast-induced damage characteristics in the static-dynamic stress field.
A lot of studies have been done in the area of rock breaking due to blasting considering the effect of in situ stress. Kutter et al. [1] analytically and experimentally investigated the influence of in situ stress on the blast-induced rock fracture. The results showed that the cracks induced by blasting stress wave and gas pressure grow preferably in the direction of maximum principal stress of superimposed stress fields. Zhang and Peng [2][3][4][5] studied the crater blasting under different confining pressures via theoretical analysis and physical model tests. The outcomes showed that the crater shape becomes oval with the long axis aligned on the loading direction, and the open angle in this direction and the crater volume is greater with the increase of confining pressure under uniaxial static load. Based on the fracture mechanics and the rock damage failure criterion, Xiao et al. [6] calculated rock fragmentation induced by blasting under high stress. It is concluded that the release of strain energy in the highly stressed rock mass is helpful to improve the breaking effect.
Yang et al. [7] conducted caustics experiments to investigate the propagation characteristics of blast-induced cracks in the dynamic-static stress field. The results indicated that the in situ stress has an important effect on crack propagation induced by blasting and the crack propagation is restrained when the crack propagation direction is perpendicular to the direction of in situ stress. Hu and Lu [8,9] studied the formation and propagation of crack induced by presplitting blasting in highly stressed rock mass via a mathematical model and concluded that the in situ stress can restrain the development of cracks between the presplitting holes when the in situ stress is perpendicular to the crack face. Yang and He [10,11] experimentally investigated the influence of confining pressure and ratios of horizontal-to-vertical pressure on the blast-induced rock fracture. The results showed that the direction of crack growth was largely controlled by the hoop tensile stress and biaxial pre-pressure ratio.
As a research tool, the numerical modeling method has been widely used to investigate the blast-induced damage characteristics of rock. Donzé et al. [12] used the discrete element method (DEM) to study the blast-induced radial fractures under confining pressure and found that the radial fractures induced by blasting tend to grow in the direction of maximum principal stress. Yilmaz et al. [13] investigated the blast-induced damage characteristics under different in situ stresses via a 3D FLAC (Fast Lagrangian Analysis of Continua) analysis. The results indicated that the development of fractures around the borehole is governed by the maximum principal stress and it is more obvious with the increase of the difference between the two principal stresses. Xie [14] used LS-DYNA to study the damage characteristics in cutting blasting under different in situ stresses. The results showed that with the increase of in situ stress, the damage zone becomes smaller. With the increase of the lateral pressure coefficient, the extending direction of the tensile damage zone becomes more obvious, which causes a great challenge to the cutting blasting excavation in deep rock masses. Yi, Jayasinghe, Ma and Li [15][16][17][18] used LS-DYNA to investigate the influence of in situ stress on the blast-induced cracks. Their results showed that the crack propagation trends towards the direction of maximum compressive pressure. Han, Wei and Deng [19] used a numerical model to study the contour control blasting under different in situ stresses. The result indicated that the in situ stress could affect the crack evolution and direction, and the quality of the contour surface is hard to control in highly stressed rock masses.
The studies mentioned above mainly focus on the plane problems of blasting under static load. However, the three-dimension propagation of stress wave induced by explosives, the charge length, and the detonation velocity of explosive cannot be considered in these 2D plane strain models. The above factors can be involved in a 3D model analysis to obtain more realistic results. Additionally, the key factor of the free surface is rarely considered, especially the lateral free surface, which plays an important role in the rock breaking by blasting. In this study, a 3D blasting model of coupling static and dynamic loads is developed in LS-DYNA and the model parameters are calibrated by the physical model test. Subsequently, the calibrated numerical model is used to simulate the blastinduced damage considering the roles of in situ stress and lateral free surface. Based on the damage distribution, the blast-induced damage characteristics and the explosion craters under different static loads and free surface spans are analyzed. To verify the material model and apply it to the subsequence simulation of lateral blasting under static load, a physical single-hole crater blasting model test was conducted firstly to calibrate the material parameter, as shown in Figure 1a. The whole model is 400 × 400 × 200 mm. A borehole with a diameter of 8.0 mm and a length (L) of 40 mm is drilled in the center and the explosive with a diameter (d) of 8.0 mm and a length (l) of 12 mm is charged in the borehole. The cemented sand, which is composed of ordinary Portland cement (PC32.5), uniform-grained sand and water in the mass ratio of 3:3:1, is used as the model material to study the blast-induced damage of rock. The material mechanical parameters are determined by averaging the measured data from six mortar cubic blocks. The density ρ 0 is 2456 kg/m 3 ; the compressive strength f c is 48.3 MPa and the elastic modulus E is 32.36 GPa; Poisson's ratio µ is 0.24; P-wave velocity v p is 3828 m/s. According to the physical model, the single-hole crater blasting numerical model was developed for comparison with the test results, as shown in Figure 1b. The model consists of rock, explosive, and stemming. The size of the numerical model is same as the physical model and the total number of the meshed elements is 0.56 million, where the numerical convergence tests has been carried out and the calculation results of the model are convergent and accurate.

Numerical
Appl. Sci. 2021, 11, x FOR PEER REVIEW 3 of 22 in the center and the explosive with a diameter (d) of 8.0 mm and a length (l) of 12 mm is charged in the borehole. The cemented sand, which is composed of ordinary Portland cement (PC32.5), uniform-grained sand and water in the mass ratio of 3

RHT Material Model for Rock
In this study, the dynamic response and damage process of rock mass were simulated by Riedel, Hiermaier, and Thoma (RHT) in LS-DYNA, which were widely used in the

RHT Material Model for Rock
In this study, the dynamic response and damage process of rock mass were simulated by Riedel, Hiermaier, and Thoma (RHT) in LS-DYNA, which were widely used in the numerical simulation of concrete and rock [3,[20][21][22][23][24][25]. The damage degree of the RHT material model is given by D = ∑ ∆ε p ε f , in which ∆ε p is the accumulated plastic strain and ε f is the failure strain. The detailed introduction of the RHT constitutive model can be found in [20]. Based on the tested mechanical parameter, the following sections obtain the other material parameters through empirical formulas or related literatures.

Strain rate parameters
The effect of strain rate on the rock strength is clear. The strain rate strength factor F r . ε p is expressed as [21]: where . ε p is the strain rate; . ε c 0 and . ε t 0 are the compressive and tensile reference strain rates, which are 3 × 10 −5 s −1 and 3 × 10 −6 s −1 , respectively; P is the hydrostatic pressure; f c and f t are the uniaxial strengths in compression and tension. The strain rates in compression (β c ) and in tension (β t ) are constant for the material and can be calculated by where the unit is MPa, and then β c and β t are determined as 0.024 and 0.029.

Failure surface parameters
A and N are the constants in the failure surface and can be obtained by: where σ * f (P * , F r ) is the normalized strength and can be calculated by σ * f = σ f f c ; P * is the normalized hydrostatic pressure and can be calculated by P * = P f c ; F r is the strain rate strength factor.
When the rock is in a quasi-static state, . ε p = 3.0 × 10 −5 s −1 , and then F r = 1 can be calculated by Equation (1). The rock strength under different confining pressures can be calculated by empirical equation of Hoek-Brown [26] and the fitting equations for rock material are as follows [27]: The axial stress at failure (σ 1 ) under different confining pressures (σ 2 = σ 3 ) can be calculated by Equation (4). The results are shown in Table 1. Subsequently, the pressure P = The P-α compaction of RHT is given by where B 0 and B 1 are the material constants; α and ρ are the initial porosity and density; e is the specific internal energy; µ 0 is the volumetric strain; A 1 , A 2 , A 3 are the polynomial coefficients. A 1 , A 2 , and A 3 can be calculated by formulas in [24], where c is the wave speed; k is the material constant. A 1 , A 2 and A 3 can be calculated as 36.0, 40.4 and 4.8 GPa, respectively. The minimum damaged residual strain (ε m p ) can be determined by the calibration of the physical model test. The remaining model parameters in this study, which are not sensitive to the numerical results, are referred to the cemented sand parameters in the literatures [28,29]. The determined RHT parameters are listed in Table 2. Compaction pressure P co 6 MPa Break compressive strain rate . ε c 3.0 × 10 25 Porosity exponent N p 3 Break tensile strain rate . ε t 3.0 × 10 25 Initial porosity α 0 1 1 ONEMPA is the unit conversion factor defining 1 MPa in the pressure units used.

Material Identifications of Charge, Air and Stemming
The charge is modeled by MAT_HIGH_EXPLOSIVE_BURN in LS-DYNA [30]. The JWL(Jones-Wilkins-Lee) EOS are given by: where P is the pressure, A, B, R 1 , R 2 , and ω are constants, V is the specific volume, and E is the internal energy with an initial value of E 0 . In this study, the explosive is a mixture of RDX (Hexogen), PETN, DDNP, et al. The estimation of JWL parameters of explosive is complex and costly [31], so the parameters refer to similar explosive parameters [32]: A = 524 GPa, B = 7.68 GPa, R 1 = 4.2, R 2 = 1.1, ω = 0.34, E 0 = 8.5 GPa. The charge density is 1.6 × 10 3 kg/m 3 and the detonation velocity is 6950 m/s. The air is modeled by MAT_NULL in LS-DYNA, and the corresponding EOS is given by [30]: where C 0 , C 1 , C 2 , C 3 , C 4 , C 5 and C 6 are polynomial coefficients; µ = ρ ρ 0 − 1 is specific volume; e is the internal energy per volume and has the unit of pressure, Pa. In this study, the air is modeled as an ideal gas by setting C 0 = C 1 = C 2 = C 3 = C 4 = 0 and C 5 = C 6 = 0.4, and the initial internal energy per volume is set to 0.25 J/cm 3 [25].
The stemming is modeled by MAT_SOIL_AND_FOAM in LS-DYNA and its parameters are shown in Table 3.  Figure 1c shows the comparison of explosion craters between the physical test results and the simulated results. In the simulated results, the critical damage is set to 0.6 or 0.7 according to the previous studies [15,22,23,33]. In this study, a critical value D of 0.7 is reasonable for the consistency between the physical test results and the numerical results. It can be found that the crater boundary in the simulated results is similar to that in the test results. Therefore, the calibrated numerical model is able and feasible to study the blast-induced damage characteristics.

Numerical Model for Lateral Blasting under Static Load
A numerical model with dimensions of 400 × 400 × 200 mm was built to simulate the dynamic response and damage evolution of lateral blasting, as shown in Figure 2. It is commonly seen that rectangular cavern is usually generated by production boreholes with an excavation method of lateral caving with large-diameter long-hole blasting or by tunneling excavation for its advantages of simple procedures, high excavation efficiency and convenient support measures [22,34]. Therefore, the cavern prototype was set as a rectangle in this study, which also has important enlightening significances for other shapes. In the center of the model, a rectangular cavern with a size of span X × span Y × 200 mm is placed to form the lateral free surface and a borehole with a diameter of 8.0 mm is placed near the cavity with a distance of W. It should be noted that the span X and span Y are no more than 1/3 of the model size of 400 mm to decrease the influence of boundary on the stress distribution. The explosive with a diameter of 6.3 mm and a length of 20 mm is charged in the hole centrally. The ends of the borehole are filled with stemming. Static stresses, P1 and P2, are applied to the four external boundaries of the model in X and Y directions respectively using a dynamic relaxation scheme, and the four sides inside the model are specified as free surfaces. After the stress initialization, the charge is loaded in the model and detonated. The numerical model is meshed by hexahedral elements, with a size of 4 mm, which is small enough to avoid any wave distortion [22]. The total number of meshed elements is 0.5 million. In this study, in order to monitor the damage distribution in the rock mass, cut Y1 and cut Z1 are selected, as shown in Figure 2. The evolution process of blast-induced damage is completed before 100 µs, which is set to the calculation termination time.
model are specified as free surfaces. After the stress initialization, the charge is loaded in the model and detonated. The numerical model is meshed by hexahedral elements, with a size of 4 mm, which is small enough to avoid any wave distortion [22]. The total number of meshed elements is 0.5 million. In this study, in order to monitor the damage distribution in the rock mass, cut Y1 and cut Z1 are selected, as shown in Figure 2. The evolution process of blast-induced damage is completed before 100 μs, which is set to the calculation termination time.

Influence of Uniaxial Static Load on the Damage Distribution under Different Burdens
In this section, a span X of 100 mm and a span Y of 100 mm are applied to the numerical model, and then the simulations of uniaxial loading of P1, uniaxial loading of P2, and biaxial loading were conducted. The damage contours of cut Z1 at 100 μs were extracted from the numerical results. The variation characteristics of damage distribution, the volume V and the shape of the crater were studied in detail in this section.
According to the elastic mechanics, on the right side of the rectangle cavern, the Xdirection stress is small, and the stress field is governed by the Y-direction stress, especially near the free surface. Therefore, the blast-induced damage is mainly affected by the original Y-direction stress field. Figure 3 shows the Y-direction elastic stress fields ( y σ ) under different static loads (W = 4 cm). For the stress y σ , its sign is positive in tension and negative in compression. As seen from the Y-direction stress contours, the rectangle cavern induces stress concentration near the free surface. For P1 = 5 MPa, as shown in Figure 3a, there is a large tensile stress zone on the left side of the borehole, especially near the free surface, and a small tensile stress zone on the right side of the borehole. The maximum tensile stress is 6.6 MPa near the free surface and y σ decreases to around 2 MPa

Influence of Uniaxial Static Load on the Damage Distribution under Different Burdens
In this section, a span X of 100 mm and a span Y of 100 mm are applied to the numerical model, and then the simulations of uniaxial loading of P1, uniaxial loading of P2, and biaxial loading were conducted. The damage contours of cut Z1 at 100 µs were extracted from the numerical results. The variation characteristics of damage distribution, the volume V and the shape of the crater were studied in detail in this section.
According to the elastic mechanics, on the right side of the rectangle cavern, the Xdirection stress is small, and the stress field is governed by the Y-direction stress, especially near the free surface. Therefore, the blast-induced damage is mainly affected by the original Y-direction stress field. Figure 3 shows the Y-direction elastic stress fields (σ y ) under different static loads (W = 4 cm). For the stress σ y , its sign is positive in tension and negative in compression. As seen from the Y-direction stress contours, the rectangle cavern induces stress concentration near the free surface. For P1 = 5 MPa, as shown in Figure 3a, there is a large tensile stress zone on the left side of the borehole, especially near the free surface, and a small tensile stress zone on the right side of the borehole. The maximum tensile stress is 6.6 MPa near the free surface and σ y decreases to around 2 MPa on the left side of the borehole. For P2 = 5 MPa, as shown in Figure 3b, the excavation zone is in a compressive stress field. The σ y is maximum around the free surface (around 12 MPa), and it decreases to around 7 MPa on the left side of the borehole. For P1 = P2 = 5 MPa, as shown in Figure 3c, the excavation zone is also in a compressive stress field, but the compressive stress field is weakened, and the distribution changes a lot. The maximum σ y transfers from the free surface, where the σ y decreases to around 4 MPa, to the four corners of the cavern. The above static stress field analysis is beneficial to understanding the coupling mechanism of static load and blasting stress wave load on the rock damage characteristics in the subsequent dynamic analysis.
To evaluate the influence of P1 on the damage distribution due to blasting, five cases of uniaxial static loads, P1 = 0, 2, 5, 8, and 10 MPa, were first conducted in this section. Figure 4 shows the damage contours for different P1 with different burdens at cut Z1. As mentioned in Section 2.2, the elements with a damage level above 0.7 are regarded as severe damage zone and form the explosion crater. In the case of P1 = 0 MPa, blast-induced severe damage zones (D ≥ 0.7) are widely distributed and can form craters from the charge center to the free surface when the burden W is no more than 4 cm. When W is more than 5 cm, the blast-induced severe rock damage (D ≥ 0.7) mainly distributes around the explosive and little severe damage zone covers the free surface, but the two zones are not connected. Thus, in the cases of W = 5 cm and W = 6 cm, only blasting cavities are formed around the charge but no crater is formed by the blasting. Therefore, the burden should be no more than 4 cm to form an explosion crater in the case of P1 = 0 MPa. In the case of P1 = 2 MPa, the damage zones are enlarged for each burden, but the severe damage zone (D ≥ 0.7) around the charge and the damage zone near the free surface are still separated when W = 5 cm and W = 6 cm, which indicates no crater is formed. When P1 increases to 5 MPa, the damage zones are further enlarged, and the two zones begin to connect for W = 5 cm but not for W = 6 cm. In the case of P1 = 8 MPa and 10 MPa, with the increase of static load, the damage zones become larger. Especially for W = 5 cm, the severe damage zone is clearly enlarged near the free surface, thus a crater is formed. It should be noted that there is still a large low-level damage zone (D < 0.7) between the blasting cavity and the free surface for W = 6 cm, as a result, the explosion crater cannot be formed. To evaluate the influence of P1 on the damage distribution due to blasting, five cases of uniaxial static loads, P1 = 0, 2, 5, 8, and 10 MPa, were first conducted in this section. Figure 4 shows the damage contours for different P1 with different burdens at cut Z1. As mentioned in Section 2.2, the elements with a damage level above 0.7 are regarded as severe damage zone and form the explosion crater. In the case of P1 = 0 MPa, blast-induced severe damage zones (D ≥ 0.7) are widely distributed and can form craters from the charge center to the free surface when the burden W is no more than 4 cm. When W is more than 5 cm, the blast-induced severe rock damage (D ≥ 0.7) mainly distributes around the explosive and little severe damage zone covers the free surface, but the two zones are not connected. Thus, in the cases of W = 5 cm and W = 6 cm, only blasting cavities are formed around the charge but no crater is formed by the blasting. Therefore, the burden should be no more than 4 cm to form an explosion crater in the case of P1 = 0 MPa. In the case of P1 = 2 MPa, the damage zones are enlarged for each burden, but the severe damage zone (D ≥ 0.7) around the charge and the damage zone near the free surface are still separated when W = 5 cm and W = 6 cm, which indicates no crater is formed. When P1 increases to 5 MPa, the damage zones are further enlarged, and the two zones begin to connect for W = 5 cm but not for W = 6 cm. In the case of P1 = 8 MPa and 10 MPa, with the increase of In order to evaluate the explosion crater clearly, the crater volume V is measured by counting the high-level damage elements (D ≥ 0.7) and summing their volumes. The crater volumes for each burden under different static loads are shown in Figure 5. It can be found that with P1 increases, the volume of explosion crater tends to increase. When P1 is less than 8 MPa, the V for W = 2 cm is the smallest (no crater for W = 5 and 6 cm, as shown in Figure 4) due to excessive dissipation of the explosion energy into the atmosphere and the V for W = 3 cm is the largest, which indicates that the optimal burden is 3 cm. However, when P1 is more than 8 MPa, the burden of 4 cm is optimal because its corresponding crater volume is the largest. It can be clearly seen that the V for W = 2 cm is not sensitive to the static load, but the others vary greatly with changing P1, especially for W = 4 cm, W = 5 cm and W = 6 cm when P1 ≥ 5 MPa. The above results are due to the rapid expansion of the damage zone near the free surface for W = 5 cm and W = 6 cm (as shown in Figure 4). It should be noted that the V for W = 6 cm is the volume sum of the blasting cavity and the damage zone near the free surface, but not the crater volume (as shown in Figure 4). It can be concluded that the P1 can change the optimal burden of charge and increase the critical embedding depth of the charge.  In order to evaluate the explosion crater clearly, the crater volume V is measured by counting the high-level damage elements (D ≥ 0.7) and summing their volumes. The crater volumes for each burden under different static loads are shown in Figure 5. It can be found that with P1 increases, the volume of explosion crater tends to increase. When P1 is less than 8 MPa, the V for W = 2 cm is the smallest (no crater for W = 5 and 6 cm, as shown in Figure 4) due to excessive dissipation of the explosion energy into the atmosphere and the V for W = 3 cm is the largest, which indicates that the optimal burden is 3 cm. However, when P1 is more than 8 MPa, the burden of 4 cm is optimal because its corresponding crater volume is the largest. It can be clearly seen that the V for W = 2 cm is not sensitive to the static load, but the others vary greatly with changing P1, especially for W = 4 cm, W = 5 cm and W = 6 cm when P1 ≥ 5 MPa. The above results are due to the rapid expansion of the damage zone near the free surface for W = 5 cm and W = 6 cm (as shown in Figure  4). It should be noted that the V for W = 6 cm is the volume sum of the blasting cavity and the damage zone near the free surface, but not the crater volume (as shown in Figure 4). It can be concluded that the P1 can change the optimal burden of charge and increase the critical embedding depth of the charge.  In order to investigate the variation in the shape of the explosion crater, the craters (formed by the elements with D ≥ 0.7) for W = 4 cm under different P1 are plotted in Figure  6. In the XY plane, the shape of the explosion craters is similar to a triangle and expands with the increase of P1. The diameter of the crater in the Y direction also becomes larger with the increase of P1 at different depths (X direction), especially in the top of the crater, where a new damage zone is generated. The above results are induced by the combined effect of the Y-direction tensile component of the stress wave and the Y-direction tensile stress field (as shown in Figure 3a). In the XZ plane, there is a clear trend that with the increase of P1, the diameter in the Z direction becomes larger, especially when P1 = 8 MPa. This is because with the increase of Y-direction stress field ( 1  In order to investigate the variation in the shape of the explosion crater, the craters (formed by the elements with D ≥ 0.7) for W = 4 cm under different P1 are plotted in Figure 6. In the XY plane, the shape of the explosion craters is similar to a triangle and expands with the increase of P1. The diameter of the crater in the Y direction also becomes larger with the increase of P1 at different depths (X direction), especially in the top of the crater, where a new damage zone is generated. The above results are induced by the combined effect of the Y-direction tensile component of the stress wave and the Y-direction tensile stress field (as shown in Figure 3a). In the XZ plane, there is a clear trend that with the increase of P1, the diameter in the Z direction becomes larger, especially when P1 = 8 MPa. This is because with the increase of Y-direction stress field (σ y1 ) induced by P1, the combined effect of the Y-direction tensile component of incident stress wave (σ yI ) and Y-direction stress field (σ y1 ) is intensified and promotes the initiation and propagation of radial tensile fracture at point A, as shown in Figure 7a. Besides, some damage zones appear in the left side of the borehole due to the stress concentration, but they expand little in the Y direction, as shown in the YZ plane. In the YZ plane, when P1 = 5 MPa, the long axis of the bottom circle of the crater is in the Y direction. This is because the superposition of the Y-direction component of reflected tensile stress (σ yR ) and Y-direction stress field (σ y1 ) at point B (as shown in Figure 7b) increases the Y-direction dimension of reflected tensile damage zone around the free surface, which is also shown in the XY plane. However, with the increase of P1, the long axis of the bottom circle of the crater transfers from Y direction to Z direction due to the faster growth of the diameter in the Z direction, which is consistent with the results in the XZ plane and can be illustrated by Figure 7a. The result is also consistent with the law that the long axis of blast-induced damage is parallel to the max principal compressive stress (Z direction). It can be concluded that the crater shape is governed by the reflected tensile fractures when P1 ≤ 5 MPa, but governed by the radial tensile fractures when P1 ≥ 8 MPa.
Another five cases of uniaxial static loads, P2 = 0 MPa, P2 = 2 MPa, P2 = 5 MPa, P2 = 8 MPa and P2 = 10 MPa, were simulated to investigate the effect of load direction on damage distribution. Figure 8 shows the damage contours at cut Z1 and the crater volumes for different burdens under different P2. In the case of P2 = 2 MPa, the crater volumes are enlarged for W = 2 cm, 3 cm, and 4 cm, but shrunk for W = 5 and 6 cm, as shown in Figure 8b. There is still no crater formed by the blasting for W = 5 and 6 cm. as shown in Figure 8a. The increase of the crater volumes for W = 2, 3, and 4 cm is mainly induced by the increase of reflected tensile fractures around the free surface, where the combined effect of the X-direction component (σ xR ) of reflected tensile wave and the compressive stress field (σ y2 ) induced by P2 promotes the damage development at point C and D, and the increase of radial tensile damage zone, where the combined effect of the Z-direction tensile component (σ zL ) of incident wave and σ y2 promotes the damage development at point B, as shown in Figure 9. However, the reductions of the crater volumes for W = 5 cm and 6 cm are induced by the reduction of radial damage zones distributed around the charge, where the volumes of blasting cavities are mainly restrained by σ y2 , especially at point A, as shown in Figure 9. In the case of P2 = 5 MPa, the crater volumes are increased when W = 2 cm, 4 cm, 5 cm, and 6 cm but reduced when W = 3 cm. For W = 3 cm, this may be because the increase of reflected tensile fractures around the free surface is smaller than the reduction of the radial tensile damage zone around the charge. There is a clear increase of the damage zones for W = 4, 5, and 6 cm, which is induced by the great increase of fractures around the free surface. Especially for W = 5 cm, the fractures around the free surface are clearly enlarged and begin to connect with the blasting cavity formed by the damage zone around the charge. For W = 6 cm, the combined effect of σ xR and σ y2 is enhanced due to the intensification of the latter, and some damage zones extend from the free surface to the borehole. However, the damage zones only distribute along the Y direction but expand little in the Z direction, so the crater is hard to form. In the case of P2 = 8 MPa and 10 MPa, with the increase of static load, the damage zones become larger, except for W = 3 cm, where the damage zone distribution along the Z direction is reduced, as shown in Figure 8a. It should be noted that the crater is still not formed for W = 6 cm. For each W, the fractures tend to extend along the Y direction with P2 increases, which is consistent with the law that the long axis of the blast-induced damage zone is parallel to the max principal stress (Y direction). Appl. Sci. 2021, 11, x FOR PEER REVIEW 11 of 22    rection but expand little in the Z direction, so the crater is hard to form. In the case of P2 = 8 MPa and 10 MPa, with the increase of static load, the damage zones become larger, except for W = 3 cm, where the damage zone distribution along the Z direction is reduced, as shown in Figure 8a. It should be noted that the crater is still not formed for W = 6 cm. For each W, the fractures tend to extend along the Y direction with P2 increases, which is consistent with the law that the long axis of the blast-induced damage zone is parallel to the max principal stress (Y direction).  With the increase of P1 or P2, the crater volume can be enlarged. However, the increase of V with P1 is more than that with P2. For example, when W = 4 cm, the V for P1 increases by 17.1, 25.8, 26.8 and 60.3 cm 3 compared with that for P2 when the stress level is 2, 5, 8 and 10 MPa, respectively. In other words, when the static load values are the same, the V is increased by 15-31% for P1 compared with P2. The results show that the effect of P1 on the increase of V is greater than that of P2. With the increase of P1 or P2, the crater volume can be enlarged. However, the increase of V with P1 is more than that with P2. For example, when W = 4 cm, the V for P1 increases by 17.1, 25.8, 26.8 and 60.3 cm 3 compared with that for P2 when the stress level is 2, 5, 8 and 10 MPa, respectively. In other words, when the static load values are the The shapes of the craters for W = 4 cm under different P2 are plotted in Figure 10a. In the XY plane, the shape of the explosion crater gradually changes from a triangle to a trapezoid and has a significant expansion in the Y direction at different depths due to the directional effect of the Y-direction compressive stress field. However, the depth of the explosion crater is reduced a little due to the volume shrink at the top of the crater, especially in the case of P2 = 10 MPa, where the crater is mainly distributed on the left side of the borehole. In the XZ plane, near the free surface, the Z-direction diameter gradually increases under the combined effect of reflected tensile wave and σ y2 . In the YZ plane, the shape of the crater becomes an oval and its long axis is in the Y direction for P2 = 2 MPa, which obeys the law that the long axial of blast-induced damage zone is parallel to the max principal stress. However, when P2 ≥ 5 MPa, the shape tends to expand in the Z direction. This is because the reflected tensile fracture zone becomes the dominant factor affecting the crater shape. As shown in Figure 9b, σ y2 is perpendicular to σ xR and Z-direction tensile stress component (σ zR ) of reflected wave at point C and point D, and it is conductive to the growth of reflected tensile fractures induced by σ xR . However, it is opposite to the Y-direction tensile stress component (σ yR ) at point D and it will restrain the formation of reflected tensile fractures induced by σ yR . As a result, the Z-direction reflected tensile fractures are easier to propagate. To study the effect of P2 on the radial tensile fracture zone distribution, as shown in Figure 10a, section A-A at 1 cm to the left side of the borehole is selected, which is away from the free surface, and its damage zone is mainly governed by radial tensile fractures. is conductive to the growth of reflected tensile fractures induced by xR σ . However, it is opposite to the Y-direction tensile stress component ( yR σ ) at point D and it will restrain the formation of reflected tensile fractures induced by yR σ . As a result, the Z-direction reflected tensile fractures are easier to propagate. To study the effect of P2 on the radial tensile fracture zone distribution, as shown in Figure 10a, section A-A at 1 cm to the left side of the borehole is selected, which is away from the free surface, and its damage zone is mainly governed by radial tensile fractures.   In section A-A, it can be clearly found with P2 increases, the Z-direction dimension of the crater reduces (P2 ≥ 5 MPa) and the Y-direction dimension of the crater increases gradually, as shown in Figure 10b. The above results are caused by the coupling mechanism of σ y2 and the incident compressive stress wave, which is illustrated by Figure 9a. At point A, the original static compressive stress field (σ y2 ) is opposite to the Y-direction tensile stress component (σ yI ) and it will prevent the formation of radial tensile fractures induced by σ yI . However, at point B, σ y is perpendicular to σ zI , and it is conductive to the growth of reflected tensile fractures induced by σ zI . Therefore, the radial tensile fracture zone is an ellipse with a long axis in the Y direction.

Influence of Biaxial Static Load on the Damage Distribution
To investigate the characteristic of damage distribution under biaxial loading, W = 4 cm and P1 = 5 MPa were kept, 11 cases of P2, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10 MPa, were simulated to investigate the effect of the pressure coefficient λ = P2 P1 on damage distribution. The borehole is placed on the right side of the cavern, named as case I, or the top side of the cavern, named as case II. Figure 11a shows the shapes of craters for W = 4 cm under biaxial loads with different λ for case I. With the λ increases, the shapes in the three planes are shrinking in the overall trend. In the XY plane, the shape transforms from a trapezoid to a triangle, and the X-direction dimension reduces much to cause the decrease of the depth of the crater. In the XZ plane, away from the free surface, the Z-direction and X-direction dimensions are both reduced due to the increasing Y-direction compressive stress field. In the YZ plane, two damage zones distribute the leftmost side and the rightmost side along the Z direction for λ = 0, but they disappear for λ = 0.4. For λ = 0, the two damage zones are formed by the combined effect of σ yR and σ y1 induced by P1. However, for λ = 0.4, the addition of the new Y-direction compressive stress field (σ y2 ) induced by P2 will neutralize part of the tensile stress and prevent the formation of the two damage zones. When λ ≥ 0.4, the shape is getting flatter due to the increase of the additional compressive stress field. For case I, the crater volume V reduces monotonically with the increase of λ, as shown in Figure 11c.
It should be noted that the law only applies to the case of small λ. When λ is large enough, the Y-direction compressive stress field induced by lateral pressure P2 will govern the damage distribution around the surface and away from the charge, and the crater volume may be increased with the increase of λ, as described in the uniaxial loading case of P2.  Figure 11b shows the shapes of craters for W = 4 cm under biaxial loads with different λ for case II. The shape of the explosion crater is the smallest in each plane when λ = 0.4. In the XY plane, the dimension in the X direction reduces first when = 0.4. This is because the X-direction compressive stress field, which is beneficial for the damage development along the X direction, is weakened by the X-direction tensile stress field induced by P2. When λ = 1.0, the dimension in the X direction increases instead, especially the damage zones on the upside and downside of the borehole. This is because that the damage mechanism has changed, and the damage zones are mainly formed by the combined effect of the X-direction tensile stress, which is caused by the rock rebound, and the Y-direction compressive stress field induced by P2. Compared with the case of λ = 0.4, the compressive stress field reduces in the X direction but enhances in the Y direction, which is beneficial for the evolution of the damage zones on the upper and lower sides of the borehole. Therefore, the X-direction crater dimension increases at the top of the crater. When λ = 1.6 and 2.0, the crater shape is enlarged further due to the increase of P2. In the YZ plane, the variation characteristics of the crater shape are similar to those in the XY plane, and the shape tends to be a triangle, which is consistent with the characteristics in the uniaxial loading cases of P1. In the XZ plane, when λ = 0, the damage zone near the free surface is induced by the coupling effect of the reflected tensile stress wave and the X-direction compressive stress field. When λ = 0.4, the combined effect is weakened by the addition of the X-direction tensile stress field induced by P2. Afterwards, with the increase of λ, the P2 becomes the dominant factor affecting the crater shape and the damage mechanism begins to change. When λ = 1.6 and 2.0, two damage zones appear on the left side and the right side along the Z direction, which are formed by the coupling effect of the X-direction tensile component of the reflected wave and the X-direction tensile stress field. It can be seen that when λ ≥ 0.4, the variation characteristics of crater shape at each plane are similar to these in the uniaxial loading cases of P1. For case II, as shown in Figure 11c, the crater volume V reduces first before λ increases to 0.4 and then increases with λ increases. The turning point of V is mainly caused by the change of the dominant damage mechanism, as described in the above analysis. It can be found that λ = 1 is a demarcation point. The crater volume V is greater for case I (the borehole is placed on the right side of the rectangular cavity) when λ < 1, but greater for case II (the borehole is placed on the top side of the rectangular cavity) when λ > 1. There is a common feature that the V is greater when the borehole is placed on the side of the max static load. Taking the demarcation point of λ = 1 as the reference point, the increase of V with the increase of the static load on the side of the borehole is greater than that with the reduction of the static load on the other side. The result indicates that the crater volume is more sensitive to the variation of static load on the same side than the other side, which is consistent with the uniaxial load numerical result in Section 3.1 that the effect of P1 on the increase of V is greater than that of P2.

Influence of Span Ratio on the Damage Distribution
In this section, the span X of 100 mm, W = 4 cm, and P1 = 5 MPa were kept, and the span Y = 100, 110, 120 and 130 mm were considered to investigate the characteristic of the damage distribution with different span ratios k = span Y span X . Considering that the static stress field on the top side is less influenced by the variation of span Y, the borehole layout placed on the top side is not considered in this section. Figure 12 shows the Y-direction elastic stress fields with different k under P1 = P2 = 5 MPa. It can be found that with span Y increases, the σ y on the right side of the rectangle cavern reduces from around 6 MPa to around 2 MPa, but varies little on the top side. The results indicate that the Y-direction compressive stress field on the right side is weakened with the increase of k. Figure 13a shows the explosion crater with different k under P1 = P2 = 5 MPa. In the XY plane, the crater Y-direction dimensions near the free surface are enlarged with the increase of k, which is similar to the characteristics of case II when λ ≥ 0.4 in Section 3.2. The expanded damage zones are also induced by the combined effect of the Y-direction rock rebound and the X-direction compressive stress field. When k increases, the Y-direction compressive stress field reduces, and the X-direction compressive stress field increases, which can intensify the combined effect and improve the X-direction damage development. In the XZ plane, the crater shape is enlarged with the increase of k, and the variation characteristic is also consistent with the results of case II when λ ≥ 0.4 in Section 3.2. In the YZ plane, the crater shape changes from an ellipse to a circle with the increase of k, which is induced by the weakening of the Y-direction compressive field. Figure 13b shows the explosion crater with different k under P1 = 5 MPa and P2 = 10 MPa. The results show that the shape variation characteristics are similar to those under P1 = P2 = 5 MPa. The crater volumes under different k are listed in Table 4. It can be found that the crater volume increases with k increases. The results indicate that in the stressed rock mass, the explosion crater can be improved by increasing the free surface span on the side of the borehole. Especially for the case of unequal biaxial loading, the rock on the side of the maximum principle stress should be excavated first, where the rock breaking efficiency is higher than the other side with the same span, and then the span on the other side can be increased, which is beneficial to improving the explosion crater on this side. In the view of strain energy density, Yang [35] pointed out that the rock mass with a poor strain energy density should be excavated first to release the high strain energy of the adjacent rock, and then the release intensity of strain energy can be effectively controlled and the vibration induced by the instantaneous unloading can be reduced. In our study, the borehole on the side of P2, where the strain energy density is poor (as shown in Figure 14), should also be detonated firstly to improve the crater volume and the critical embedding depth of the charge. Besides, as shown in Figure 14, with the excavation of rock mass on the side of P2, the strain energy density on the other side (P1) will be reduced. The result is beneficial to the control of vibration induced by the instantaneous unloading and the increase of the rock breaking efficiency in the high strain energy zone. results indicate that the Y-direction compressive stress field on the right side is weakened with the increase of k. Figure 13a shows the explosion crater with different k under P1 = P2 = 5 MPa. In the XY plane, the crater Y-direction dimensions near the free surface are enlarged with the increase of k, which is similar to the characteristics of case II when λ ≥ 0.4 in Section 3.2. The expanded damage zones are also induced by the combined effect of the Y-direction rock rebound and the X-direction compressive stress field. When k increases, the Y-direction compressive stress field reduces, and the X-direction compressive stress field increases, which can intensify the combined effect and improve the X-direction damage development. In the XZ plane, the crater shape is enlarged with the increase of k, and the variation characteristic is also consistent with the results of case II when λ ≥ 0.4 in Section 3.2. In the YZ plane, the crater shape changes from an ellipse to a circle with the increase of k, which is induced by the weakening of the Y-direction compressive field. Figure 13b shows the explosion crater with different k under P1 = 5 MPa and P2 = 10 MPa.
The results show that the shape variation characteristics are similar to those under P1 = P2 = 5 MPa. The crater volumes under different k are listed in Table 4. It can be found that the crater volume increases with k increases. The results indicate that in the stressed rock mass, the explosion crater can be improved by increasing the free surface span on the side of the borehole. Especially for the case of unequal biaxial loading, the rock on the side of the maximum principle stress should be excavated first, where the rock breaking efficiency is higher than the other side with the same span, and then the span on the other side can be increased, which is beneficial to improving the explosion crater on this side. In the view of strain energy density, Yang [35] pointed out that the rock mass with a poor strain energy density should be excavated first to release the high strain energy of the adjacent rock, and then the release intensity of strain energy can be effectively controlled and the vibration induced by the instantaneous unloading can be reduced. In our study, the borehole on the side of P2, where the strain energy density is poor (as shown in Figure  14), should also be detonated firstly to improve the crater volume and the critical embedding depth of the charge. Besides, as shown in Figure 14, with the excavation of rock mass on the side of P2, the strain energy density on the other side (P1) will be reduced. The result is beneficial to the control of vibration induced by the instantaneous unloading and the increase of the rock breaking efficiency in the high strain energy zone.

Discussion and Conclusions
This study mainly investigates the blast-induced damage characteristics considering the lateral free surface in highly stressed rock mass by using a 3D numerical model. Firstly, the numerical model is calibrated by comparing the results of single-hole crater blasting model test and the numerical simulation. The results of the test and simulation are in good agreement in the upper part of the crater boundary, but the test result gives a slightly larger contour in the lower part. This may be due to the existence of an uneven weakening part in the lower part during the pouring process of cemented sand, which is more likely to be damaged and not considered in the numerical simulation. On the whole, the numerical model is reasonably accurate to study the blast-induced damage characteristics. And then, the influence of uniaxial static load under different burdens, biaxial static load and span ratio on the damage distribution were evaluated with the calibrated numerical model. When the free surface exists, the damage mechanisms are clearly different from these of the plane problem without considering the free surface and the 3D stress state. The development of the blast-induced damage away from the charge is governed by the static load, especially near the free surface.
For uniaxial loading, the damage zone and the crater volume V increase with the increase of uniaxial static load P1 or P2. The V is greater when the borehole is placed on the side of the static load. For example, when W = 4 cm, compared with P2, the V is increased by 15-31% for P1 with the same static load value. The variations of crater volume also show that the uniaxial static load can change the optimal burden of charge and increase the critical embedding depth of the charge, especially for P1.
Static load significantly affects the blast-induced damage distribution, especially for the radial tensile fractures zone and the reflected tensile fractures zone. For the radial tensile fractures zone, the long axis of the bottom circle of the crater turns parallel to the max principal compressive stress, which has been proposed by many researchers [2][3][4]8,[13][14][15]. However, the law only applies to the radial tensile fractures but not to the reflected tensile fractures. For the latter, the opposite law that the long axis of the reflected tensile fracture zone is perpendicular to the max principal compressive stress will be obtained.
The variation law of the crater volume is different for the case of the biaxial static load. When P1 = 5 MPa is kept and the borehole is placed on the side of P1, the V reduces monotonously with the increase of P2 before λ increases to 2.0. When the borehole is placed on the other side, the V first reduces until λ increases to 0.4 and then increases with the increase of P2 before λ increases to 2.0. The turning point at λ = 0.4 is mainly caused by the transformation of the dominant factor affecting the crater shape and the damage mechanism from P1 to P2. The crater volume is greater when the borehole is placed on the side of the max static load, where the strain energy density is lower. Meanwhile, the V is also more sensitive to the variation of static load on the same side than the other side.
The crater volume increases with k increases, which indicates that in the stressed rock mass, the explosion crater can be improved by increasing the free surface span on the side of the borehole. Especially for the case of unequal biaxial loading with the same span (λ = 1 and k = 1), the rock on the side of max principle stress should be excavated first, where the rock breaking efficiency is higher than the other side, and then the span on the other side can be increased (k increases), which is beneficial to improving the blasting efficiency on this side. Meanwhile, when the charge on the side of the max static load (lower strain energy density) is detonated first and the charge on the other side (higher strain energy density) is detonated later, the transient unloading induced vibration can be reduced [35].

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.