# Numerical Study of Fracture Characteristics of Deep Granite Induced by Blast Stress Wave

Authors: Zheming Zhu; Weiting Gao; Duanying Wan; Meng Wang; Yun Shu
Journal: Shock and Vibration
 (2021)

---

## Abstract

To study the characteristics of rock fracture in deep underground under blast loads, some numerical models were established in AUTODYN code. Weibull distribution was used to characterize the inhomogeneity of rock, and a linear equation of state was applied to describe the relation of pressure and volume of granite elements. A new stress initialization method based on explicit dynamic calculation was developed to get an accurate stress distribution near the borehole. Two types of in situ stress conditions were considered. The effect of heterogeneous characteristics of material on blast-induced granite fracture was investigated. The difference between 2D models and 3D models was discussed. Based on the numerical results, it can be concluded that the increase of the magnitude of initial pressure can change the mechanism of shear failure near the borehole and suppress radial cracks propagation. When initial lateral pressure is invariable, with initial vertical pressure rising, radial cracks along the acting direction of vertical pressure will be promoted, and radial cracks in other directions will be prevented. Heterogeneous characteristics of material have an obvious influence on the shear failure zones around the borehole.


---

## Body

## 1. Introduction

Chinese Sichuan-Tibet railway, which will cross Sichuan Basin, Yunnan-Guizhou Plateau, and Qinghai-Tibet Plateau, is one of the most challenging railway lines to build in the world [1, 2], along with more than 90% of its length running on bridges or in tunnels [3]. From Lhasa to Chengdu, with an elevation drop more than 3,000 meters, as shown in Figure 1, the railway will meander through the mountains. More than anything, the accumulated height of the railway reaches more than 14,000 meters, and it will cross many high in situ stress zones, i.e., about 80% of the tunnels are buried in the range of 500–2,000 m, and the in situ stresses are in the range of 10–40 MPa. Also, such geographical conditions will bring great difficulties for drilling and blasting engineering. Therefore, it is very significant to investigate rock fracture characteristics under in situ stresses [4], which can optimize blasting design.Figure 1
A route map of the Sichuan-Tibet railway.Because the response of rocks under dynamic load is very complicated, it is challenging to conduct relative researches by experimental study directly [5, 6], and numerical experiment has gradually become the mainstream. Lu et al. [7] investigated dynamic rock response around the tunnel induced by the release of in situ stress by finite element method in ABAQUS code; similarly, Yang et al. [8–11] and Yan et al. [12] studied the rock damage near the tunnel under blast loads, and dynamic stress redistribution by plenty of numerical simulations, and the results show that the redistribution of dynamic contributes a lot to the rock damage. Xiao et al. [13] conducted numerical research for dynamic rock stress unloading near the tunnel by dynamic implicit and explicit methods in LS-DYNA code. His results show that high in situ stresses have a significant influence on rock breakage induced by blasting. Qiu et al. [14] established the numerical model of a deep-buried U-shaped tunnel in a particle flow code (PFC2D) and investigated the influence of buried depth on the tunnel dynamic stability. Hana et al. [15, 16] combined finite and discrete element methods (FDEM) to study the damage evolution of a deep tunnel by combining finite element and discrete element methods.To better understand the mechanisms of deep rock breakage by blast load, some scholars investigated rock fracture characteristics near the borehole under high in situ stresses [17]. Xie et al. [18, 19] conducted many numerical simulations to investigate the effect of a free surface, in situ stress, and lateral pressure coefficient on blast-induced rock damage zones based on ANSYS\LS-DYNA code, and the results show that initial in situ stresses will cause the propagation anisotropy of damage. Yi et al. [20] and Tao et al. [21] analyzed the transmit of blast stress wave in the prestressed rock masses based on the dynamic elastic framework and investigated the effect of in situ stresses on blast-induced damage and microcracks distributions by some numerical methods. Yilmaz and Unlu [22] studied the effect of anisotropic high in situ stresses on blast-induced rock damage by dynamic geotechnical software (FLAC 3D).So far, some scholars have researched the propagation of a single crack under high in situ stresses [23, 24]. Yang and Ding [25, 26] investigated main crack propagation under initial in situ stresses and blast loads by dynamic caustics systems and found that initial in situ stresses can change the crack propagating mode from mode-I to mode-II crack. Yang et al. [27] analyzed stress distributions near borehole and explored crack dynamic propagation behaviour under high initial static stresses, and the results show that initial stress can suppress crack propagation in the direction perpendicular to its direction.However, it is insufficient for studying rock fracture characteristics induced by blast loads in the deep underground. In one aspect, inhomogeneity of rock is rarely considered in many numerical simulations, and it is well known that lots of microcracks are irregularly distributed in the interior of rock masses, which can cause the local difference of mechanic parameters, such as density, elasticity modulus, Poisson ratio, and rock strength. It is generally considered that the microcracks conform toWeibull distribution in rock masses [28]. Therefore, it will be used to correct the numerical models in this study.In another aspect, the Damage indexD, usually used in many pieces of research, cannot explicitly express the mechanism of rock fracture, which is often used to show the damage degree. Many investigations have found two main failure types, shear failure and tensile failure when rock is subjected to dynamic loads. So, a modified max-principal failure criterion was used in the numerical models to evaluate the rock failure state [29], which can assist in studying the fracture mechanism of rock under blast loads.In investigating the deep underground blast, it is significant to realize the stress distribution as the actual conditions in numerical models. At present, a widely used method is to get the initial stress distribution by implicit dynamic solutions or static solutions [30], and then, it will be imported into the explicit dynamic model as initial stress condition, which has some shortcomings, such as difficult operating processes, high requirements for material model, and time consuming. Therefore, a simple and efficient stress initialization method was developed in this study, which is entirely based on explicit dynamic calculation and has an accurate initial stress distribution.In this study, two in situ stress conditions have been considered. The first is to explore the influence of the magnitude of initial pressure on blast-induced rock fractures when initial lateral pressure is equal to initial vertical pressure. The second is to keep initial lateral pressure stable and investigate the effect of pressure ratio on radial cracks distributions by changing initial vertical pressure. In addition, the effect of heterogeneous characteristics of material on blast-induced granite fracture, and the difference between 2D models and 3D models was discussed. Based on the numerical results in this study, the fracture characteristics of granite under different in situ stresses were analyzed, which will have some guiding significance for practical engineering.


## 2. Numerical Study

### 2.1. Dynamic Finite Volume Method

The numerical simulations conducted in this study is based on the finite volume method in AUTODYN code. In the model, each material can be designed as one subgrid, as shown in Figure2.Figure 2
Mesh of a subgrid.For nodeA in Figure 2, its accelerations u¨ in x-direction and v¨ in y-direction can be calculated as follows:(1)u¨=Fxmv¨=Fym,where Fx and Fy are nodal forces in x- and y-direction, respectively, and m is the mass of the area BDEF.At timestep (n + 1/2), the velocities of node A can be calculated by(2)u˙n+1/2=u˙n−1/2+u¨nΔtv˙n+1/2=v˙n−1/2+u¨nΔt,where n is timesteps and Δt is the time increment corresponding to every timestep. Furthermore, the displacement of node A at timestep (n+1) can be expressed as follows:(3)un+1=un+u˙n+1/2Δtvn+1=vn+v˙n+1/2Δt.Based on the velocities of the nodesA, B, C, and D, the strain rate of the element ABCD can be obtained by(4)ε˙x=∂u˙∂x=12AEn+1/2u˙A−u˙cyB−yD−u˙B−u˙DyA−yCε˙y=∂v˙∂y=12AEn+1/2v˙A−v˙cxB−xD−v˙B−v˙DxA−xCγ˙xy=∂v˙∂x+∂u˙∂y=12AEn+1/2v˙A−v˙cyB−yD−v˙B−v˙DyA−yC+u˙A−u˙cxB−xD−u˙B−u˙DxA−xC.According to the linear elastic constitutive, the components of the deviatoric stress tensor for the two-dimensional case can be calculated by(5)Sxn+1=Sxn+2Gε˙x−13e˙n+1/2ΔtSyn+1=Syn+2Gε˙y−13e˙n+1/2ΔtSxyn+1=Sxyn+2Gε˙xyn+1/2Δt,where G is the shear modulus, ε˙x, ε˙y, and ε˙xy are the components of Cauchy strain tensor, and e˙ is the volume strain rate. Therefore, the components of Cauchy stress tensor can be obtained as follows:(6)σx=P+Sxσy=P+Syτxy=Sxy,where P is the hydrostatic pressure, and it is usually calculated by the equation of state.In order to acquire the node force, integrate the wave equation in the two-dimensional case as shown in the following equation:(7)∬S∂σx∂x+dτxydydxdy=∬Sρ∂2u∂t2dxdy∬S∂σy∂y+dτxydxdxdy=∬Sρ∂2v∂t2dxdy.Applying Green’s formula in the above integrating equations, the following equation can be obtained:(8)∫Lσxdy−τxydx=∬Sρ∂2u∂t2dxdy=2m∂2u∂t2=2Fx∫Lτxydy−σydx=∬Sρ∂2v∂t2dxdy=2m∂2v∂t2=2Fy,where S is the area of the quadrilateral element BDEF, and L is the integral path BDEF, and m is the mass of the element. Based on the above equations, the node force Fx and Fy can be expressed as follows:(9)Fx=12σx1yF−yE−τxy1xF−xE+σx2yB−yF−τxy2xB−xF+σx3yD−yB−τxy3xD−xB+σx4yE−yD−τxy4xE−xDFy=12τxy1yF−yE−σy1xF−xE+τxy2yB−yF−σy2xB−xF+τxy3yD−yB−σy3xD−xB+τxy4yE−yD−σy4xE−xD.The above process is a calculating loop in AUTODYN code. By constantly calculating, the physical parameters of the numerical model can be acquired at different times.


### 2.2. Numerical Model

In deep underground blasting engineering, the displacement of rocks along the direction of a borehole is usually strongly constrained. Therefore, a numerical model based on the plane strain condition was established in this study, as shown in Figure3, which consists of four parts: granite, copper, polyethene, and PETN. The granite model is 200 mm in width and height, and the radius of the borehole is 3.5 mm. The charge method is similar to research of Banadaki and Mohanty [31]. The military explosive, PETN, is used to offer blast loads, which is a high-energy explosive and can set off the pressure of about 10 GPa on the borehole wall when detonated. The coupling material is a circle of polyethene with a thickness of 1.5 mm. A copper ring with a thickness of 1.5 mm is set on its periphery to prevent the explosion gas from going into the rock model, as shown in Figure 3. To ensure the accuracy of numerical results, a very refining meshing was applied in four parts, and all the models were discretized by two-dimensional quadrilateral elements, and elements numbers are 156400, 320, 256, and 400, respectively, for granite, copper, polyethene, and PETN.Figure 3
Numerical models proposed in this article.


#### 2.2.1. Granite Numerical Model

A suitable material model is critical for numerical simulations. Many scholars have applied the RHT model [32] and JHR model [33] in simulating rock dynamic mechanic behaviour, the excellent applicability of which has been verified by many cases. However, too many material parameters need to be considered, and their measurements are very complicated. A simple and practical way is to apply a linear elasticity model in rock material [34], which is based on the fact that rock is of high brittleness so that little plastic strain will happen when it breaks. Therefore, in this study, a linear elasticity model was used to model the relation of stress and strain, which can bring much convenience.In the granite interior, lots of microcracks are irregularly widely distributed, which will influence the homogeneity of rock mechanics properties. Because the microcracks conform toWeibull distribution in the rock mass [28], the rock mechanics properties are also considered to conform to this distribution, which can be expressed by(10)fx=mx0⋅xx0m−1⋅e−x/x0m,where fx is the probability density function of Weibull distribution, x0 is the statistical average of rock mechanical characteristics, such as the measured density and the measured elastic modulus, and m is the homogeneity coefficient, which represents inhomogeneity degree of microcracks distributions in the rock masses or dispersion degree of material properties. According to equation (10), the probability function of Weibull distribution can be obtained as follows:(11)Px=1−e−x/x0m,where Px is the probability value when the certain material property is less than the value of x. As shown in Figure 4, there are curves of the probability function with m of 20, 50, and 100.0. It can be found that the distribution of material mechanical property tends to close to its statistical average x0 with the increase of m.Figure 4
The curves of the probability of Weibull distribution.The average dynamic mechanical parameters used in the granite model has been listed in Table1 [31], and its density, bulk modulus, and dynamic shear modulus are 2.66 g/cm3, 25.59 GPa, and 21.83 GPa, respectively. Tensile strength and shear strength are set as 30 MPa and 75 MPa. In establishing models, it is considered that the distributions of the above parameters conform to Weibull distribution, and the vertical axis in Figure 4 is equally divided into 100 intervals. Then, the corresponding element numbers to one interval can be obtained by products of total elements number in the model times the probability value of the interval, and the average mechanical value at this interval was assigned to the corresponding elements. For example, in the n-th interval, the two boundary values of probability are Pn+1 and Pn, and the number of the corresponding element Nn can be calculated by(12)Nn=M⋅Pn+1−Pn,where M is the total elements number of the granite model. After determining the number of elements corresponding to every interval, they will be randomly assigned to the granite numerical model. As shown in Figure 3, elements in different colours belong to different intervals, with different dynamic mechanical parameters. In this way, the inhomogeneity of the granite model has been corrected by the Weibull distribution.Table 1
Average dynamic mechanical parameters of granite.
Densityρ (g/cm3)Dynamic shear modulusGd (GPa)Bulk modulusK (GPa)Tensile strengthσT (MPa)Shear strengthτs (MPa)2.6621.8325.5930.075.0As mentioned above, there is little plastic strain when fracture happens for granite due to brittleness characteristics. Therefore, a linear equation of state was used to express the relation of its deformation and pressure, which can be written as(13)P=k⋅ρρ0−1,where P was the hydrostatic pressure, k is the bulk modulus, and ρ/ρ0 is the ratio of the density of the deformed state and initial state. A modified principal stress failure criterion was used to evaluate the failure state of granite elements [29], which can be expressed as(14)σ1≥σTorτmax=σ1−σ32≥τS,where σT and τS are tensile strength and shear strength of granite, respectively, and σ1 and σ3 are the maximum and minimum principal stress. The modified principal stress failure criterion means that when maximum shear stress τmax exceeds the shear strength τS or maximum principal stress σ1 exceeds the tensile strength σT, the element will fail, and it will not stand any shear or tensile anymore.


#### 2.2.2. JWL Equation for PETN

The military explosive, PETN, was used to offer blast loads in this study, which can set off about 10 GPa pressure on the borehole wall and cause rock breakage. The explosion is a very complex chemical reaction, with extremely high energy release during several microseconds, and the energy will make explosion products expand rapidly. In general, the JWL equation of state is used to model the relation of blast pressure and the volume of explosion products, which can be expressed by(15)P=A1−ωR1Ve−R1V+B1−ωR2Ve−R2V+ωEV,where P is the explosive pressure, V is the special volume of explosion products, E is the initial total energy, and A, B, R1, R2, and ω are the constants of explosions. For PETN[35], A = 348.62 GPa, B = 11.29 GPa, R1 = 7.0, R2 = 2.0, and ω = 0.247.


#### 2.2.3. Shock EOS

As the copper and polyethene are directly subjected to strong blast shock waves, Mie-Gruneisenshock EOS is used to simulate their states in which the relationship of hydrostatic pressure and material density can be expressed by(16)p=ρ0c02μ1+μ1−s−1μ21−γ0μ2μ+1+γ0E,where c0 and s are the material constants, μ=ρ/ρ0−1, ρ is the current density, ρ0 is the initial density, γ0 is the Gruneisen constant γ0≈2s1−1, and E is the initial internal energy per unit volume. The material parameters of copper and polyethene are listed in Table 2.Table 2
Parameters of shock EOS of copper and polyethene.
Materialc0 (m/s)sCopper39401.489Polyethene29011.481


#### 2.2.4. Verification of Numerical Model

All the models in the article were solved in AUTODTN code, which is a highly nonlinear explicit dynamic solver and does well in modelling rock fracture under blast loads [36–39]. To verify the accuracy of the model proposed in the study, a numerical result without in situ stresses was compared with the research results [19, 40], as shown in Figures 5(a) and 5(b). According to their results, there is a very serious damage zone (called the crushed zone) caused by blasts, and some radial cracks emanating from the crushed zone can be observed. The same phenomena were acquired in our model, as shown in Figure 5(c), so the granite numerical model proposed in the article can well simulate rock fracture characteristics under blast load.Figure 5
Comparison of blast-induced fractures. (a) Numerical result by Zhu et al. [40]. (b) Numerical result by Xie et al. [19]. (c) Numerical result in this study. (d) Damage details near the borehole.
(a)(b)(c)(d)As shown in Figure5(d), it is the details near the borehole. It was divided into three zones according to the damage degree and failure types of rocks. Zone I is the pure shear failure zone caused by the severe compressive function of blast shock waves. With blast shock waves declining into stress waves of high amplitude, it just leads to some element shear failure. Due to the hoop tensile stress of blast rarefaction waves, many radial cracks appear around the shear failure elements. Therefore, zone II called the severe damage zone, includes two failure types. In zone III, blast stress waves cannot cause elements shear failure due to lower compressive stresses. Some radial crack tips in zone II will keep propagating under hoop tensile stress, so zone III is also called radial cracks zone.


### 2.3. Application of In Situ Stresses

When considering in situ stress conditions, the practical engineering problem will be simplified into the model, as shown in Figure6. The lateral pressure Px and the vertical pressure Py just need to be considered. According to Saint Venant’s Principle, there is strong stress concentration around the borehole under in situ stresses, affecting blast results. The solutions of stress distributions near the borehole can be deduced based on elasticity theory, as expressed in equations (17)–(20):(17)σr=−Py+Px21−a2r2+Py−Px21+3a4r4−4a2r2cos2θ,(18)σθ=−Py+Px21+a2r2−Py−Px21+3a4r4cos2θ,(19)σz=vσr+σθ,where a is the radius of the borehole, r is the distance from the center of the borehole, σr is initial radial stress, σθ is initial hoop stress, and σz is principal stress along the z-direction.Figure 6
Sketch of applying in situ stresses.It is significant to realize the same initial stress distributions in the numerical model as the actual state. The traditional method is to conduct stresses initialization by importing the results of the implicit dynamic calculation or static calculation to an explicit dynamic model, which is a complex operation. In this study, a new method entirely based on explicit dynamic calculation was developed to acquire actual stress distributions near the borehole. First, the four edges of the model were set as the transmit boundary to simulate infinite rock media. Then, the corresponding constant stress boundaries, such asPx = 20 MPa and Py = 20 MPa, were applied on upper-lower boundaries and right-left boundaries. Finally, the explicit dynamic solving will be carried out for the model until it reaches an equilibrium state.Figure 7
σy distributions of the granite model as a function of time during stress initialization. (a) t = 8.27 μs. (b) t = 16.57 μs. (c) t = 24.87 μs. (d) t = 33.18 μs. (e) t = 41.48 μs. (f) t = 49.78 μs.In Figure7, there are the contour plots of σy at different times during stress initialization. Figures 7(a) and 7(b) show that stress waves set off by the constant stress boundaries are propagating in the model, and in Figures 7(c) and 7(d), it can be found that stress waves have passed the borehole and caused stress concentrations near it. With stress waves going forward, they will not be reflected due to the transmit boundaries’ existence when stress waves encounter the edges of the model, as shown in Figure 7(e). Therefore, the numerical model will easily reach the statics equilibrium state, as Figure 7(f) shows. To study the change of stresses during the stress initialization, four stress gauges were set in the locations of r=a,2a,3aand4a distance from the borehole center to record histories of hoop stress, as shown in Figure 8(a). It can be found that the stress value at r=a is much higher than others, which means that initial pressures caused strongly compressive stress concentration. With time increasing, all the curves of hoop stress gradually tend to be stable, and although there still exits some vibration, it is too small to need to be considered. According to the curves in Figure 8(a), it was thought that the numerical model had reached the equilibrium state at 90 μs, so it was set as the detonation time of the explosive.Figure 8
Stress details near the borehole during applying in situ stresses. (a) History ofσθ at gauge points. (b) Distributions of σθ and σr near the borehole.
(a)(b)To verify the accuracy of the method of stress initialization developed in the article, hoop stresses and radial stresses in the range of 25 mm distance from the borehole center were acquired from numerical results corresponding tot = 90 μs, which were compared with the theoretical solutions according to equations (18) and (19), as shown in Figure 8(b). It can be observed that there is a pretty good consistency between numerical solutions and theoretical solutions, which evidences that the method in this study can do well in simulating initial stress distributions near the borehole.


## 2.1. Dynamic Finite Volume Method

The numerical simulations conducted in this study is based on the finite volume method in AUTODYN code. In the model, each material can be designed as one subgrid, as shown in Figure2.Figure 2
Mesh of a subgrid.For nodeA in Figure 2, its accelerations u¨ in x-direction and v¨ in y-direction can be calculated as follows:(1)u¨=Fxmv¨=Fym,where Fx and Fy are nodal forces in x- and y-direction, respectively, and m is the mass of the area BDEF.At timestep (n + 1/2), the velocities of node A can be calculated by(2)u˙n+1/2=u˙n−1/2+u¨nΔtv˙n+1/2=v˙n−1/2+u¨nΔt,where n is timesteps and Δt is the time increment corresponding to every timestep. Furthermore, the displacement of node A at timestep (n+1) can be expressed as follows:(3)un+1=un+u˙n+1/2Δtvn+1=vn+v˙n+1/2Δt.Based on the velocities of the nodesA, B, C, and D, the strain rate of the element ABCD can be obtained by(4)ε˙x=∂u˙∂x=12AEn+1/2u˙A−u˙cyB−yD−u˙B−u˙DyA−yCε˙y=∂v˙∂y=12AEn+1/2v˙A−v˙cxB−xD−v˙B−v˙DxA−xCγ˙xy=∂v˙∂x+∂u˙∂y=12AEn+1/2v˙A−v˙cyB−yD−v˙B−v˙DyA−yC+u˙A−u˙cxB−xD−u˙B−u˙DxA−xC.According to the linear elastic constitutive, the components of the deviatoric stress tensor for the two-dimensional case can be calculated by(5)Sxn+1=Sxn+2Gε˙x−13e˙n+1/2ΔtSyn+1=Syn+2Gε˙y−13e˙n+1/2ΔtSxyn+1=Sxyn+2Gε˙xyn+1/2Δt,where G is the shear modulus, ε˙x, ε˙y, and ε˙xy are the components of Cauchy strain tensor, and e˙ is the volume strain rate. Therefore, the components of Cauchy stress tensor can be obtained as follows:(6)σx=P+Sxσy=P+Syτxy=Sxy,where P is the hydrostatic pressure, and it is usually calculated by the equation of state.In order to acquire the node force, integrate the wave equation in the two-dimensional case as shown in the following equation:(7)∬S∂σx∂x+dτxydydxdy=∬Sρ∂2u∂t2dxdy∬S∂σy∂y+dτxydxdxdy=∬Sρ∂2v∂t2dxdy.Applying Green’s formula in the above integrating equations, the following equation can be obtained:(8)∫Lσxdy−τxydx=∬Sρ∂2u∂t2dxdy=2m∂2u∂t2=2Fx∫Lτxydy−σydx=∬Sρ∂2v∂t2dxdy=2m∂2v∂t2=2Fy,where S is the area of the quadrilateral element BDEF, and L is the integral path BDEF, and m is the mass of the element. Based on the above equations, the node force Fx and Fy can be expressed as follows:(9)Fx=12σx1yF−yE−τxy1xF−xE+σx2yB−yF−τxy2xB−xF+σx3yD−yB−τxy3xD−xB+σx4yE−yD−τxy4xE−xDFy=12τxy1yF−yE−σy1xF−xE+τxy2yB−yF−σy2xB−xF+τxy3yD−yB−σy3xD−xB+τxy4yE−yD−σy4xE−xD.The above process is a calculating loop in AUTODYN code. By constantly calculating, the physical parameters of the numerical model can be acquired at different times.


## 2.2. Numerical Model

In deep underground blasting engineering, the displacement of rocks along the direction of a borehole is usually strongly constrained. Therefore, a numerical model based on the plane strain condition was established in this study, as shown in Figure3, which consists of four parts: granite, copper, polyethene, and PETN. The granite model is 200 mm in width and height, and the radius of the borehole is 3.5 mm. The charge method is similar to research of Banadaki and Mohanty [31]. The military explosive, PETN, is used to offer blast loads, which is a high-energy explosive and can set off the pressure of about 10 GPa on the borehole wall when detonated. The coupling material is a circle of polyethene with a thickness of 1.5 mm. A copper ring with a thickness of 1.5 mm is set on its periphery to prevent the explosion gas from going into the rock model, as shown in Figure 3. To ensure the accuracy of numerical results, a very refining meshing was applied in four parts, and all the models were discretized by two-dimensional quadrilateral elements, and elements numbers are 156400, 320, 256, and 400, respectively, for granite, copper, polyethene, and PETN.Figure 3
Numerical models proposed in this article.


### 2.2.1. Granite Numerical Model

A suitable material model is critical for numerical simulations. Many scholars have applied the RHT model [32] and JHR model [33] in simulating rock dynamic mechanic behaviour, the excellent applicability of which has been verified by many cases. However, too many material parameters need to be considered, and their measurements are very complicated. A simple and practical way is to apply a linear elasticity model in rock material [34], which is based on the fact that rock is of high brittleness so that little plastic strain will happen when it breaks. Therefore, in this study, a linear elasticity model was used to model the relation of stress and strain, which can bring much convenience.In the granite interior, lots of microcracks are irregularly widely distributed, which will influence the homogeneity of rock mechanics properties. Because the microcracks conform toWeibull distribution in the rock mass [28], the rock mechanics properties are also considered to conform to this distribution, which can be expressed by(10)fx=mx0⋅xx0m−1⋅e−x/x0m,where fx is the probability density function of Weibull distribution, x0 is the statistical average of rock mechanical characteristics, such as the measured density and the measured elastic modulus, and m is the homogeneity coefficient, which represents inhomogeneity degree of microcracks distributions in the rock masses or dispersion degree of material properties. According to equation (10), the probability function of Weibull distribution can be obtained as follows:(11)Px=1−e−x/x0m,where Px is the probability value when the certain material property is less than the value of x. As shown in Figure 4, there are curves of the probability function with m of 20, 50, and 100.0. It can be found that the distribution of material mechanical property tends to close to its statistical average x0 with the increase of m.Figure 4
The curves of the probability of Weibull distribution.The average dynamic mechanical parameters used in the granite model has been listed in Table1 [31], and its density, bulk modulus, and dynamic shear modulus are 2.66 g/cm3, 25.59 GPa, and 21.83 GPa, respectively. Tensile strength and shear strength are set as 30 MPa and 75 MPa. In establishing models, it is considered that the distributions of the above parameters conform to Weibull distribution, and the vertical axis in Figure 4 is equally divided into 100 intervals. Then, the corresponding element numbers to one interval can be obtained by products of total elements number in the model times the probability value of the interval, and the average mechanical value at this interval was assigned to the corresponding elements. For example, in the n-th interval, the two boundary values of probability are Pn+1 and Pn, and the number of the corresponding element Nn can be calculated by(12)Nn=M⋅Pn+1−Pn,where M is the total elements number of the granite model. After determining the number of elements corresponding to every interval, they will be randomly assigned to the granite numerical model. As shown in Figure 3, elements in different colours belong to different intervals, with different dynamic mechanical parameters. In this way, the inhomogeneity of the granite model has been corrected by the Weibull distribution.Table 1
Average dynamic mechanical parameters of granite.
Densityρ (g/cm3)Dynamic shear modulusGd (GPa)Bulk modulusK (GPa)Tensile strengthσT (MPa)Shear strengthτs (MPa)2.6621.8325.5930.075.0As mentioned above, there is little plastic strain when fracture happens for granite due to brittleness characteristics. Therefore, a linear equation of state was used to express the relation of its deformation and pressure, which can be written as(13)P=k⋅ρρ0−1,where P was the hydrostatic pressure, k is the bulk modulus, and ρ/ρ0 is the ratio of the density of the deformed state and initial state. A modified principal stress failure criterion was used to evaluate the failure state of granite elements [29], which can be expressed as(14)σ1≥σTorτmax=σ1−σ32≥τS,where σT and τS are tensile strength and shear strength of granite, respectively, and σ1 and σ3 are the maximum and minimum principal stress. The modified principal stress failure criterion means that when maximum shear stress τmax exceeds the shear strength τS or maximum principal stress σ1 exceeds the tensile strength σT, the element will fail, and it will not stand any shear or tensile anymore.


### 2.2.2. JWL Equation for PETN

The military explosive, PETN, was used to offer blast loads in this study, which can set off about 10 GPa pressure on the borehole wall and cause rock breakage. The explosion is a very complex chemical reaction, with extremely high energy release during several microseconds, and the energy will make explosion products expand rapidly. In general, the JWL equation of state is used to model the relation of blast pressure and the volume of explosion products, which can be expressed by(15)P=A1−ωR1Ve−R1V+B1−ωR2Ve−R2V+ωEV,where P is the explosive pressure, V is the special volume of explosion products, E is the initial total energy, and A, B, R1, R2, and ω are the constants of explosions. For PETN[35], A = 348.62 GPa, B = 11.29 GPa, R1 = 7.0, R2 = 2.0, and ω = 0.247.


### 2.2.3. Shock EOS

As the copper and polyethene are directly subjected to strong blast shock waves, Mie-Gruneisenshock EOS is used to simulate their states in which the relationship of hydrostatic pressure and material density can be expressed by(16)p=ρ0c02μ1+μ1−s−1μ21−γ0μ2μ+1+γ0E,where c0 and s are the material constants, μ=ρ/ρ0−1, ρ is the current density, ρ0 is the initial density, γ0 is the Gruneisen constant γ0≈2s1−1, and E is the initial internal energy per unit volume. The material parameters of copper and polyethene are listed in Table 2.Table 2
Parameters of shock EOS of copper and polyethene.
Materialc0 (m/s)sCopper39401.489Polyethene29011.481


### 2.2.4. Verification of Numerical Model

All the models in the article were solved in AUTODTN code, which is a highly nonlinear explicit dynamic solver and does well in modelling rock fracture under blast loads [36–39]. To verify the accuracy of the model proposed in the study, a numerical result without in situ stresses was compared with the research results [19, 40], as shown in Figures 5(a) and 5(b). According to their results, there is a very serious damage zone (called the crushed zone) caused by blasts, and some radial cracks emanating from the crushed zone can be observed. The same phenomena were acquired in our model, as shown in Figure 5(c), so the granite numerical model proposed in the article can well simulate rock fracture characteristics under blast load.Figure 5
Comparison of blast-induced fractures. (a) Numerical result by Zhu et al. [40]. (b) Numerical result by Xie et al. [19]. (c) Numerical result in this study. (d) Damage details near the borehole.
(a)(b)(c)(d)As shown in Figure5(d), it is the details near the borehole. It was divided into three zones according to the damage degree and failure types of rocks. Zone I is the pure shear failure zone caused by the severe compressive function of blast shock waves. With blast shock waves declining into stress waves of high amplitude, it just leads to some element shear failure. Due to the hoop tensile stress of blast rarefaction waves, many radial cracks appear around the shear failure elements. Therefore, zone II called the severe damage zone, includes two failure types. In zone III, blast stress waves cannot cause elements shear failure due to lower compressive stresses. Some radial crack tips in zone II will keep propagating under hoop tensile stress, so zone III is also called radial cracks zone.


## 2.2.1. Granite Numerical Model

A suitable material model is critical for numerical simulations. Many scholars have applied the RHT model [32] and JHR model [33] in simulating rock dynamic mechanic behaviour, the excellent applicability of which has been verified by many cases. However, too many material parameters need to be considered, and their measurements are very complicated. A simple and practical way is to apply a linear elasticity model in rock material [34], which is based on the fact that rock is of high brittleness so that little plastic strain will happen when it breaks. Therefore, in this study, a linear elasticity model was used to model the relation of stress and strain, which can bring much convenience.In the granite interior, lots of microcracks are irregularly widely distributed, which will influence the homogeneity of rock mechanics properties. Because the microcracks conform toWeibull distribution in the rock mass [28], the rock mechanics properties are also considered to conform to this distribution, which can be expressed by(10)fx=mx0⋅xx0m−1⋅e−x/x0m,where fx is the probability density function of Weibull distribution, x0 is the statistical average of rock mechanical characteristics, such as the measured density and the measured elastic modulus, and m is the homogeneity coefficient, which represents inhomogeneity degree of microcracks distributions in the rock masses or dispersion degree of material properties. According to equation (10), the probability function of Weibull distribution can be obtained as follows:(11)Px=1−e−x/x0m,where Px is the probability value when the certain material property is less than the value of x. As shown in Figure 4, there are curves of the probability function with m of 20, 50, and 100.0. It can be found that the distribution of material mechanical property tends to close to its statistical average x0 with the increase of m.Figure 4
The curves of the probability of Weibull distribution.The average dynamic mechanical parameters used in the granite model has been listed in Table1 [31], and its density, bulk modulus, and dynamic shear modulus are 2.66 g/cm3, 25.59 GPa, and 21.83 GPa, respectively. Tensile strength and shear strength are set as 30 MPa and 75 MPa. In establishing models, it is considered that the distributions of the above parameters conform to Weibull distribution, and the vertical axis in Figure 4 is equally divided into 100 intervals. Then, the corresponding element numbers to one interval can be obtained by products of total elements number in the model times the probability value of the interval, and the average mechanical value at this interval was assigned to the corresponding elements. For example, in the n-th interval, the two boundary values of probability are Pn+1 and Pn, and the number of the corresponding element Nn can be calculated by(12)Nn=M⋅Pn+1−Pn,where M is the total elements number of the granite model. After determining the number of elements corresponding to every interval, they will be randomly assigned to the granite numerical model. As shown in Figure 3, elements in different colours belong to different intervals, with different dynamic mechanical parameters. In this way, the inhomogeneity of the granite model has been corrected by the Weibull distribution.Table 1
Average dynamic mechanical parameters of granite.
Densityρ (g/cm3)Dynamic shear modulusGd (GPa)Bulk modulusK (GPa)Tensile strengthσT (MPa)Shear strengthτs (MPa)2.6621.8325.5930.075.0As mentioned above, there is little plastic strain when fracture happens for granite due to brittleness characteristics. Therefore, a linear equation of state was used to express the relation of its deformation and pressure, which can be written as(13)P=k⋅ρρ0−1,where P was the hydrostatic pressure, k is the bulk modulus, and ρ/ρ0 is the ratio of the density of the deformed state and initial state. A modified principal stress failure criterion was used to evaluate the failure state of granite elements [29], which can be expressed as(14)σ1≥σTorτmax=σ1−σ32≥τS,where σT and τS are tensile strength and shear strength of granite, respectively, and σ1 and σ3 are the maximum and minimum principal stress. The modified principal stress failure criterion means that when maximum shear stress τmax exceeds the shear strength τS or maximum principal stress σ1 exceeds the tensile strength σT, the element will fail, and it will not stand any shear or tensile anymore.


## 2.2.2. JWL Equation for PETN

The military explosive, PETN, was used to offer blast loads in this study, which can set off about 10 GPa pressure on the borehole wall and cause rock breakage. The explosion is a very complex chemical reaction, with extremely high energy release during several microseconds, and the energy will make explosion products expand rapidly. In general, the JWL equation of state is used to model the relation of blast pressure and the volume of explosion products, which can be expressed by(15)P=A1−ωR1Ve−R1V+B1−ωR2Ve−R2V+ωEV,where P is the explosive pressure, V is the special volume of explosion products, E is the initial total energy, and A, B, R1, R2, and ω are the constants of explosions. For PETN[35], A = 348.62 GPa, B = 11.29 GPa, R1 = 7.0, R2 = 2.0, and ω = 0.247.


## 2.2.3. Shock EOS

As the copper and polyethene are directly subjected to strong blast shock waves, Mie-Gruneisenshock EOS is used to simulate their states in which the relationship of hydrostatic pressure and material density can be expressed by(16)p=ρ0c02μ1+μ1−s−1μ21−γ0μ2μ+1+γ0E,where c0 and s are the material constants, μ=ρ/ρ0−1, ρ is the current density, ρ0 is the initial density, γ0 is the Gruneisen constant γ0≈2s1−1, and E is the initial internal energy per unit volume. The material parameters of copper and polyethene are listed in Table 2.Table 2
Parameters of shock EOS of copper and polyethene.
Materialc0 (m/s)sCopper39401.489Polyethene29011.481


## 2.2.4. Verification of Numerical Model

All the models in the article were solved in AUTODTN code, which is a highly nonlinear explicit dynamic solver and does well in modelling rock fracture under blast loads [36–39]. To verify the accuracy of the model proposed in the study, a numerical result without in situ stresses was compared with the research results [19, 40], as shown in Figures 5(a) and 5(b). According to their results, there is a very serious damage zone (called the crushed zone) caused by blasts, and some radial cracks emanating from the crushed zone can be observed. The same phenomena were acquired in our model, as shown in Figure 5(c), so the granite numerical model proposed in the article can well simulate rock fracture characteristics under blast load.Figure 5
Comparison of blast-induced fractures. (a) Numerical result by Zhu et al. [40]. (b) Numerical result by Xie et al. [19]. (c) Numerical result in this study. (d) Damage details near the borehole.
(a)(b)(c)(d)As shown in Figure5(d), it is the details near the borehole. It was divided into three zones according to the damage degree and failure types of rocks. Zone I is the pure shear failure zone caused by the severe compressive function of blast shock waves. With blast shock waves declining into stress waves of high amplitude, it just leads to some element shear failure. Due to the hoop tensile stress of blast rarefaction waves, many radial cracks appear around the shear failure elements. Therefore, zone II called the severe damage zone, includes two failure types. In zone III, blast stress waves cannot cause elements shear failure due to lower compressive stresses. Some radial crack tips in zone II will keep propagating under hoop tensile stress, so zone III is also called radial cracks zone.


## 2.3. Application of In Situ Stresses

When considering in situ stress conditions, the practical engineering problem will be simplified into the model, as shown in Figure6. The lateral pressure Px and the vertical pressure Py just need to be considered. According to Saint Venant’s Principle, there is strong stress concentration around the borehole under in situ stresses, affecting blast results. The solutions of stress distributions near the borehole can be deduced based on elasticity theory, as expressed in equations (17)–(20):(17)σr=−Py+Px21−a2r2+Py−Px21+3a4r4−4a2r2cos2θ,(18)σθ=−Py+Px21+a2r2−Py−Px21+3a4r4cos2θ,(19)σz=vσr+σθ,where a is the radius of the borehole, r is the distance from the center of the borehole, σr is initial radial stress, σθ is initial hoop stress, and σz is principal stress along the z-direction.Figure 6
Sketch of applying in situ stresses.It is significant to realize the same initial stress distributions in the numerical model as the actual state. The traditional method is to conduct stresses initialization by importing the results of the implicit dynamic calculation or static calculation to an explicit dynamic model, which is a complex operation. In this study, a new method entirely based on explicit dynamic calculation was developed to acquire actual stress distributions near the borehole. First, the four edges of the model were set as the transmit boundary to simulate infinite rock media. Then, the corresponding constant stress boundaries, such asPx = 20 MPa and Py = 20 MPa, were applied on upper-lower boundaries and right-left boundaries. Finally, the explicit dynamic solving will be carried out for the model until it reaches an equilibrium state.Figure 7
σy distributions of the granite model as a function of time during stress initialization. (a) t = 8.27 μs. (b) t = 16.57 μs. (c) t = 24.87 μs. (d) t = 33.18 μs. (e) t = 41.48 μs. (f) t = 49.78 μs.In Figure7, there are the contour plots of σy at different times during stress initialization. Figures 7(a) and 7(b) show that stress waves set off by the constant stress boundaries are propagating in the model, and in Figures 7(c) and 7(d), it can be found that stress waves have passed the borehole and caused stress concentrations near it. With stress waves going forward, they will not be reflected due to the transmit boundaries’ existence when stress waves encounter the edges of the model, as shown in Figure 7(e). Therefore, the numerical model will easily reach the statics equilibrium state, as Figure 7(f) shows. To study the change of stresses during the stress initialization, four stress gauges were set in the locations of r=a,2a,3aand4a distance from the borehole center to record histories of hoop stress, as shown in Figure 8(a). It can be found that the stress value at r=a is much higher than others, which means that initial pressures caused strongly compressive stress concentration. With time increasing, all the curves of hoop stress gradually tend to be stable, and although there still exits some vibration, it is too small to need to be considered. According to the curves in Figure 8(a), it was thought that the numerical model had reached the equilibrium state at 90 μs, so it was set as the detonation time of the explosive.Figure 8
Stress details near the borehole during applying in situ stresses. (a) History ofσθ at gauge points. (b) Distributions of σθ and σr near the borehole.
(a)(b)To verify the accuracy of the method of stress initialization developed in the article, hoop stresses and radial stresses in the range of 25 mm distance from the borehole center were acquired from numerical results corresponding tot = 90 μs, which were compared with the theoretical solutions according to equations (18) and (19), as shown in Figure 8(b). It can be observed that there is a pretty good consistency between numerical solutions and theoretical solutions, which evidences that the method in this study can do well in simulating initial stress distributions near the borehole.


## 3. Effect of In Situ Stresses

In this study, rock fracture characteristics induced by blast loads were investigated under two types of initial in situ stress conditions. The first is to study the effect of the magnitude of initial pressure, in which initial lateral pressurePx is equal to initial vertical pressure Py and the magnitude of initial pressure are set as 10 MPa, 20 MPa, 30 MPa, and 40 MPa, respectively. The second is to explore the influence of pressure ratio (Px/Py) on rock breakage. Initial lateral pressure Px was kept as 5 MPa, and initial vertical pressure Py is 5 MPa, 15 MPa, 25 MPa, and 35 MPa, respectively. A relatively homogeneous granite was considered, so the homogeneity coefficient m was set as 100.0 in Section 3.1 and 3.2. The solving time for all the models is 200 μs to ensure the complete action of the explosion on rock media.


### 3.1. Numerical Results under First In Situ Stress Conditions

In this section, fracture characteristics of granite under the first type of in situ stress conditions were investigated. As shown in Figure9, there are pressure contour plots under different initial pressure magnitude versus time. As we can see, the four contour plots at 96.0 μs show that with initial pressure increasing, the pressure around the borehole at the early blast stage increase. With blast stress waves propagating, many initial radial cracks appeared around the borehole as shown in plots at 101.0 μs and 108.0 μs. With initial pressure increasing, fewer radial cracks can be observed around the borehole. This is because that higher initial pressure leads to higher inner pressure, which suppresses radial cracks growing and propagating. These phenomena can be seen clearly in the plots at 116.0 μs and 132.0 μs. There are more radial cracks that propagate faster in numerical models under lower initial pressure. High initial pressure almost arrests many radial cracks propagating owing to strong compression.Figure 9
Pressure contours versus time under different pressures.The final fracture pattern is shown in Figure10, in which blue zones are explosion products, yellow zones and purple zones mean shear failure and tensile failure, respectively, and green zones represent granite elements with no damage. It can be found that the size of cracks zones of granite models is decreasing with the magnitude of initial pressure according to the numerical results, and there are big differences among the pure shear failure zone (zone I), the severe damage zone (zone II), and the radial cracks zones (zone III) among for different model. To study the differences among different models in quantity, the numbers of failure elements in zone I, zone II, and zone III were counted, and the ratio of numbers of failure elements in a different zone and total elements number of the model under different initial pressure were calculated, as plotted in Figure 11.Figure 10
Comparison of blast-induced fractures under different pressures. (a)Px = Py = 10 MPa. (b) Px = Py = 20 MPa. (c) Px = Py = 30 MPa. (d) Px = Py = 40 MPa.
(a)(b)(c)(d)Figure 11
The ratio of failure elements as functions of in situ pressures.It can be found that the ratio of the number of failure elements in the radial crack zone is higher than that in other zones whenPx = Py = 10 MPa. The ratio of pure shear failure zone is increasing with initial pressure, and others gradually decrease. When initial pressure is 30 MPa or 40 MPa, the ratio of the number of shear failure elements in zone I have exceeded the ratio of other zones. Therefore, it can be concluded that the increase of initial pressure will promote shear failure near the borehole such as enlarging the size of the pure shear zone and suppress radial cracks propagating such as shorten the length of the radial cracks and reduced their number, which can be observed in Figure 10.As shown in Figure12, the magnitude of hoop compressive around the borehole increases with initial pressure, which will suppress radial crack propagating. Similarly, higher hoop compressive stress will not be conducive to in-plane shear failure causedbyτrθ; the shear stress τrθ under blast can be calculated by(20)τrθ=σre−σθe2+σrP−σθP2,where σθe and σre are the hoop and the radial compressive stresses induced by blast stress waves, and σrPandσθP are the initial hoop and radial stresses around the borehole induced by in situ stresses. According to dynamic elasticity theory, σre is much higher than σθe in magnitude, so σre−σθe is negative. σθP is higher than σrP as shown in Figure 8(b), so σrP−σθP is positive and increasing with initial pressure; consequently τrθ will gradually decrease in value, which makes it harder to cause shear failure in the in-plane. However, the numerical results show that a larger area of pure shear failure zone appears with initial pressure increasing, which may be resulted from the transformation of the mechanism of shear failure near borehole from in-plane shear failure mode to out-plane shear failure causedbyτrz.Figure 12
Hoop stress distributions along the edge of the borehole under different values ofPx, Py.To verify the inference above, two stress gauges were set in the same locations of the pure shear zones of models ofPx = Py = 10 MPa and Px = Py = 40 MPa to record the histories of hoop stress, radial stress, and stress in the z-direction, as shown in Figure 13. According to the curves, the increase of radial stresses at gauge points is much higher than hoop stresses when the blast stress wave passes. For the model of Px = Py = 10 MPa, it can be found that σr>σz>σθ at failure time; therefore, τrθ is the max shear stress and in-plane shear failure will happen. However, for the model of Px = Py = 40 MPa, σr>σθ>σz; therefore, the max shear stress is τrz so that out-plane shear failure will occur.Figure 13
History curves of three principal stressσr, σθ, and σz in shear failure locations. (a) Px = Py = 10 MPa. (b) Px = Py = 40 MPa.
(a)(b)To study the difference in failure mechanism among four models, the three principal stresses and two shear stresses,τrθandτrz at the same locations in four models corresponding to shear failure time were obtained as plotted in Figure 14. It can be observed that three principal stresses at failure time are increasing in value with the initial pressure, and hoop stresses have higher rising amplitude than the other two, which leads to in-plane shear stress τrθ gradually reduce, but out-plane shear stress τrz gradually increase in magnitude. When the initial pressure is bigger than 20 MPa, τrz has exceeded τrθ in magnitude, as shown in Figure 14, and consequently, the failure mode around the borehole has transformed from in-plane shear failure to out-plane shear failure.Figure 14
Three principal stresses and max shear stresses at the failure time as a function of in situ pressure.In this section, rock fracture characteristics under blast loads when initial lateral pressurePx is equal to initial vertical pressure Py were explored, and three typical damage zones are analyzed in the different magnitude of initial pressures. Based on the numerical results and analysis, it can be concluded that the initial hoop compressive stress around the borehole remarkably increase with the magnitude of initial pressure, which not only suppresses radial cracks propagation but also changes the mechanism of failure mode near the borehole from in-plane shear failure mode to out-plane shear failure mode and leads to larger pure shear failure zone.


### 3.2. Numerical Results under Second In Situ Stress Conditions

In this section, rock fracture characteristics subjected to blast load under the second type of in situ stresses were investigated. The initial lateral pressurePx is 5 MPa, and the initial vertical pressure is set as 5 MPa, 15 MPa, 25 MPa, and 35 MPa, respectively for four models. The pressure plots versus time are shown in Figure 15. It can be seen in the plots at 96.0 μs that there exists uneven pressure distribution around the borehole, and with Py increasing, higher pressure magnitude appears in the right part and left part of the borehole. In the pressure contour plots at 101.0 μs and 108.0 μs, many radial cracks grow and propagate under blast stress waves, and furthermore, more radial cracks can be seen in the upper and lower part of the borehole than its right and left part with Py increasing; this is because stronger compression suppresses radial cracks propagating in the right and left parts. Besides, there is almost the same radial cracks distribution in the upper and lower part of the borehole under different Py. The same phenomena can be observed in the contour plots at 116.0 μs and 132.0 μs.Figure 15
Pressure contours versus time under differentPy when Px = 5 MPa.The numerical fracture plots are shown in Figure16. There are some clear differences for the models under the different ratios of initial lateral pressure and initial vertical pressure. The shape of pure shear failure tends to change into an ellipse with initial vertical pressure increasing. The radial cracks are not uniformly distributed along the borehole anymore. In general, the distribution of radial cracks will have a great effect on blast engineering. Therefore, the main attention will focus on the influence of different pressure ratios on radial crack distributions.Figure 16
Comparison of blast-induced fractures under differentPy when Px = 5 MPa. (a) Px = 5 MPa, Py = 5 MPa. (b) Px = 5 MPa, Py = 15 MPa. (c) Px = 5 MPa, Py = 25 MPa. (d) Px = 5 MPa, Py = 35 MPa.
(a)(b)(c)(d)To analyze the effect of different pressure ratios on radial cracks distribution, every model in Figure16 is partitioned into six sector parts, and everyone accounts for 60 degrees. The ratio of the number of tensile failure elements in every sector and total element number of the model was acquired, which was used to characterize radial cracks density in different zones, as shown in Figure 17. When the initial vertical pressure Py is equal to 5 MPa, namely, the pressure ratio is 1.0; there is a relatively uniform radial crack distribution around the borehole. With the vertical pressure increasing, the density of radial cracks in zones of the sector from 60° to 120° and from 240° to 300° is larger than other zones, which means that it is easier for radial cracks to propagate along the acting direction of vertical pressure. In other directions, radial cracks density is gradually decreasing.Figure 17
Distributions of the ratio of tensile failure elements along the edge of the borehole under the different values ofPy when Px = 5 MPa. (a) Px = 5 MPa, Py = 5 MPa. (b) Px = 5 MPa, Py = 15 MPa. (c) Px = 5 MPa, Py = 25 MPa. (d) Px = 5 MPa, Py = 35 MPa.
(a)(b)(c)(d)It can be explained easily according to hoop stress distribution around the borehole, as shown in Figure18. It can be found that the hoop compressive stresses along the borehole in the sectors from 330° to 30° and from 150° to 210° have a remarkable increase with the initial vertical pressure, which will weaken the strength of hoop tensile stresses of blast stress waves. Therefore, radial crack initiation and propagation in these zones will be prevented. On the contrary, hoop compressive stresses along the borehole in the sectors from 60° to 120° and from 240° to 300° is slightly decreasing with the increase in initial vertical pressure, which leads to the correspondingly higher radial cracks density in these zones.Figure 18
Hoop stress distributions along the edge of the borehole under the different values ofPy when Px = 5 MPa.


## 3.1. Numerical Results under First In Situ Stress Conditions

In this section, fracture characteristics of granite under the first type of in situ stress conditions were investigated. As shown in Figure9, there are pressure contour plots under different initial pressure magnitude versus time. As we can see, the four contour plots at 96.0 μs show that with initial pressure increasing, the pressure around the borehole at the early blast stage increase. With blast stress waves propagating, many initial radial cracks appeared around the borehole as shown in plots at 101.0 μs and 108.0 μs. With initial pressure increasing, fewer radial cracks can be observed around the borehole. This is because that higher initial pressure leads to higher inner pressure, which suppresses radial cracks growing and propagating. These phenomena can be seen clearly in the plots at 116.0 μs and 132.0 μs. There are more radial cracks that propagate faster in numerical models under lower initial pressure. High initial pressure almost arrests many radial cracks propagating owing to strong compression.Figure 9
Pressure contours versus time under different pressures.The final fracture pattern is shown in Figure10, in which blue zones are explosion products, yellow zones and purple zones mean shear failure and tensile failure, respectively, and green zones represent granite elements with no damage. It can be found that the size of cracks zones of granite models is decreasing with the magnitude of initial pressure according to the numerical results, and there are big differences among the pure shear failure zone (zone I), the severe damage zone (zone II), and the radial cracks zones (zone III) among for different model. To study the differences among different models in quantity, the numbers of failure elements in zone I, zone II, and zone III were counted, and the ratio of numbers of failure elements in a different zone and total elements number of the model under different initial pressure were calculated, as plotted in Figure 11.Figure 10
Comparison of blast-induced fractures under different pressures. (a)Px = Py = 10 MPa. (b) Px = Py = 20 MPa. (c) Px = Py = 30 MPa. (d) Px = Py = 40 MPa.
(a)(b)(c)(d)Figure 11
The ratio of failure elements as functions of in situ pressures.It can be found that the ratio of the number of failure elements in the radial crack zone is higher than that in other zones whenPx = Py = 10 MPa. The ratio of pure shear failure zone is increasing with initial pressure, and others gradually decrease. When initial pressure is 30 MPa or 40 MPa, the ratio of the number of shear failure elements in zone I have exceeded the ratio of other zones. Therefore, it can be concluded that the increase of initial pressure will promote shear failure near the borehole such as enlarging the size of the pure shear zone and suppress radial cracks propagating such as shorten the length of the radial cracks and reduced their number, which can be observed in Figure 10.As shown in Figure12, the magnitude of hoop compressive around the borehole increases with initial pressure, which will suppress radial crack propagating. Similarly, higher hoop compressive stress will not be conducive to in-plane shear failure causedbyτrθ; the shear stress τrθ under blast can be calculated by(20)τrθ=σre−σθe2+σrP−σθP2,where σθe and σre are the hoop and the radial compressive stresses induced by blast stress waves, and σrPandσθP are the initial hoop and radial stresses around the borehole induced by in situ stresses. According to dynamic elasticity theory, σre is much higher than σθe in magnitude, so σre−σθe is negative. σθP is higher than σrP as shown in Figure 8(b), so σrP−σθP is positive and increasing with initial pressure; consequently τrθ will gradually decrease in value, which makes it harder to cause shear failure in the in-plane. However, the numerical results show that a larger area of pure shear failure zone appears with initial pressure increasing, which may be resulted from the transformation of the mechanism of shear failure near borehole from in-plane shear failure mode to out-plane shear failure causedbyτrz.Figure 12
Hoop stress distributions along the edge of the borehole under different values ofPx, Py.To verify the inference above, two stress gauges were set in the same locations of the pure shear zones of models ofPx = Py = 10 MPa and Px = Py = 40 MPa to record the histories of hoop stress, radial stress, and stress in the z-direction, as shown in Figure 13. According to the curves, the increase of radial stresses at gauge points is much higher than hoop stresses when the blast stress wave passes. For the model of Px = Py = 10 MPa, it can be found that σr>σz>σθ at failure time; therefore, τrθ is the max shear stress and in-plane shear failure will happen. However, for the model of Px = Py = 40 MPa, σr>σθ>σz; therefore, the max shear stress is τrz so that out-plane shear failure will occur.Figure 13
History curves of three principal stressσr, σθ, and σz in shear failure locations. (a) Px = Py = 10 MPa. (b) Px = Py = 40 MPa.
(a)(b)To study the difference in failure mechanism among four models, the three principal stresses and two shear stresses,τrθandτrz at the same locations in four models corresponding to shear failure time were obtained as plotted in Figure 14. It can be observed that three principal stresses at failure time are increasing in value with the initial pressure, and hoop stresses have higher rising amplitude than the other two, which leads to in-plane shear stress τrθ gradually reduce, but out-plane shear stress τrz gradually increase in magnitude. When the initial pressure is bigger than 20 MPa, τrz has exceeded τrθ in magnitude, as shown in Figure 14, and consequently, the failure mode around the borehole has transformed from in-plane shear failure to out-plane shear failure.Figure 14
Three principal stresses and max shear stresses at the failure time as a function of in situ pressure.In this section, rock fracture characteristics under blast loads when initial lateral pressurePx is equal to initial vertical pressure Py were explored, and three typical damage zones are analyzed in the different magnitude of initial pressures. Based on the numerical results and analysis, it can be concluded that the initial hoop compressive stress around the borehole remarkably increase with the magnitude of initial pressure, which not only suppresses radial cracks propagation but also changes the mechanism of failure mode near the borehole from in-plane shear failure mode to out-plane shear failure mode and leads to larger pure shear failure zone.


## 3.2. Numerical Results under Second In Situ Stress Conditions

In this section, rock fracture characteristics subjected to blast load under the second type of in situ stresses were investigated. The initial lateral pressurePx is 5 MPa, and the initial vertical pressure is set as 5 MPa, 15 MPa, 25 MPa, and 35 MPa, respectively for four models. The pressure plots versus time are shown in Figure 15. It can be seen in the plots at 96.0 μs that there exists uneven pressure distribution around the borehole, and with Py increasing, higher pressure magnitude appears in the right part and left part of the borehole. In the pressure contour plots at 101.0 μs and 108.0 μs, many radial cracks grow and propagate under blast stress waves, and furthermore, more radial cracks can be seen in the upper and lower part of the borehole than its right and left part with Py increasing; this is because stronger compression suppresses radial cracks propagating in the right and left parts. Besides, there is almost the same radial cracks distribution in the upper and lower part of the borehole under different Py. The same phenomena can be observed in the contour plots at 116.0 μs and 132.0 μs.Figure 15
Pressure contours versus time under differentPy when Px = 5 MPa.The numerical fracture plots are shown in Figure16. There are some clear differences for the models under the different ratios of initial lateral pressure and initial vertical pressure. The shape of pure shear failure tends to change into an ellipse with initial vertical pressure increasing. The radial cracks are not uniformly distributed along the borehole anymore. In general, the distribution of radial cracks will have a great effect on blast engineering. Therefore, the main attention will focus on the influence of different pressure ratios on radial crack distributions.Figure 16
Comparison of blast-induced fractures under differentPy when Px = 5 MPa. (a) Px = 5 MPa, Py = 5 MPa. (b) Px = 5 MPa, Py = 15 MPa. (c) Px = 5 MPa, Py = 25 MPa. (d) Px = 5 MPa, Py = 35 MPa.
(a)(b)(c)(d)To analyze the effect of different pressure ratios on radial cracks distribution, every model in Figure16 is partitioned into six sector parts, and everyone accounts for 60 degrees. The ratio of the number of tensile failure elements in every sector and total element number of the model was acquired, which was used to characterize radial cracks density in different zones, as shown in Figure 17. When the initial vertical pressure Py is equal to 5 MPa, namely, the pressure ratio is 1.0; there is a relatively uniform radial crack distribution around the borehole. With the vertical pressure increasing, the density of radial cracks in zones of the sector from 60° to 120° and from 240° to 300° is larger than other zones, which means that it is easier for radial cracks to propagate along the acting direction of vertical pressure. In other directions, radial cracks density is gradually decreasing.Figure 17
Distributions of the ratio of tensile failure elements along the edge of the borehole under the different values ofPy when Px = 5 MPa. (a) Px = 5 MPa, Py = 5 MPa. (b) Px = 5 MPa, Py = 15 MPa. (c) Px = 5 MPa, Py = 25 MPa. (d) Px = 5 MPa, Py = 35 MPa.
(a)(b)(c)(d)It can be explained easily according to hoop stress distribution around the borehole, as shown in Figure18. It can be found that the hoop compressive stresses along the borehole in the sectors from 330° to 30° and from 150° to 210° have a remarkable increase with the initial vertical pressure, which will weaken the strength of hoop tensile stresses of blast stress waves. Therefore, radial crack initiation and propagation in these zones will be prevented. On the contrary, hoop compressive stresses along the borehole in the sectors from 60° to 120° and from 240° to 300° is slightly decreasing with the increase in initial vertical pressure, which leads to the correspondingly higher radial cracks density in these zones.Figure 18
Hoop stress distributions along the edge of the borehole under the different values ofPy when Px = 5 MPa.


## 4. Discussion

In this section, we simply discuss the effect of the material heterogeneous characteristics on numerical results and differences between the three-dimensional case and the two-dimensional case. The heterogeneous characteristics of material usually have an influence on the results of numerical simulation. In this section, two kinds of homogeneity coefficientm = 20 and 100 were set to discuss how the heterogeneous characteristics of material affect granite fracture characteristics under blast stress wave. The numerical results are shown in Figures 19 and 20, respectively, for two kinds of in situ stress conditions. According to the numerical results, it can be found that main effect focuses on the size of shear failure zones, and there are lager shear zones in the mode with m in 20 than m in 100 for all the case with different in situ stress conditions. For influence on tensile failure zones, some obvious differences can be seen between Figures 19(c) and 19(f), and for other cases, there exist the same tensile failure distributions. So, it can be inferred that with the homogeneity coefficient m decreasing, a larger shear failure zone will appear around the borehole.Figure 19
Comparison of blast-induced fractures under differentPy when Px = 5 MPa in m = 20, 100. (a) Px = 10 MPa, Py = 10 MPa. (b) Px = 20 MPa, Py = 20 MPa. (c) Px = 30 MPa, Py = 30 MPa. (d) Px = 10 MPa, Py = 10 MPa. (e) Px = 20 MPa, Py = 20 MPa. (f) Px = 30 MPa, Py = 30 MPa.
(a)(b)(c)(d)(e)(f)Figure 20
Comparison of blast-induced fractures underPy = Px in m = 20, 100. (a) Px = 5 MPa, Py = 15 MPa. (b) Px = 5 MPa, Py = 25 MPa. (c) Px = 5 MPa, Py = 35 MPa. (d) Px = 5 MPa, Py = 15 MPa. (e) Px = 5 MPa, Py = 25 MPa. (f) Px = 5 MPa, Py = 35 MPa.
(a)(b)(c)(d)(e)(f)An actual blasting engineering is very complex compared with the plane models proposed in this article because it is a three-dimensional problem. As is well known that there is a big difference between wave propagating in two-dimensional space and three-dimensional space, so it is necessary to discuss the difference between the two-dimensional plane blast model and the three-dimensional blast model. A three-dimensional numerical model was established as shown in AUTODYN code as shown in Figure21, and its height is 100 mm, and its cross section is a square of 100 mm in length. The borehole depth is 50 mm, and the coupling method is similar to the plane model. For 3D models, in situ stress conditions of Px = Py were just considered. The explosive was detonated at its top.Figure 21
3D numerical model.The blast results are shown in Figure22. Clear shear failure zones can be observed around the borehole; it has continuous distributions along the borehole, and it is different for the tensile failure zones. To better understand the effect of in situ pressure on blast results, the ratio of shear failure zone and cross section area and the ratio of the number of tensile failure elements and total element number of models were calculated, as shown in Figures 23(a) and 23(b). It can be found that with the depth of cross section increasing, the ratio of the shear failure zone increases firstly and then decreases. During depth from 0 mm to 20 mm, initial pressure has little effect. During depth from 20 mm to 50 mm, the ratio of shear failure zone is increasing with the increment in pressure, which has the same conclusion as the plane strain models, so it can be inferred that with depth increasing, the stress conditions of the cross section of the borehole is tending to the plane strain conditions. According to Figure 23(b), it can be found that the number of tensile failure elements is negatively proportional to the pressure value, so the tensile failure in the models is suppressed by initial pressure, and it is similar to the result of plane strain models. In summary, there exists some differences between the 2D model and the 3D model, but some characteristics of the 3D blast results can be reflected by the 2D model.Figure 22
Blast-induced damage on the symmetrical plane under different pressure. (a)P = 10 MPa. (b) P = 20 MPa. (c) P = 30 MPa. (d) P = 40 MPa.
(a)(b)(c)(d)Figure 23
The ratio of shear failure zone at different cross sections and function of the ratio of tensile failure zones of models and in situ pressure.
(a)(b)


## 5. Conclusion

In this article, deep granite fracture characteristics under blast load were investigated by numerical simulations, and two types of in situ stress conditions were considered. A linear EOS was used to describe the relation of pressure and the volume of granite elements, andWeibull distributions were applied in numerical models to characterize the inhomogeneity of rocks. The modified max principal stress failure criterion was applied in evaluating the failure state of granite. A new stress initialization method entirely based on explicit dynamic calculation was developed, and its results agree well with the theoretical solution.Numerical simulations in the article focus on discussing three aspects. First, for a granite model with high homogeneity coefficient, two kinds of in situ stresses conditions are considered to investigate blast-induced fracture characteristics. It can be found that when initial lateral pressure is equal to initial vertical pressure, with the magnitude of initial pressure increasing, the size of pure shear failure gradually increases, and the length and number of radial cracks gradually decrease. Besides, a very interesting phenomenon was found that shear failure near the borehole will be transformed into the out-plane shear mode from the in-plane shear failure mode when initial pressure increases. When initial lateral pressure is invariable, with initial vertical pressure rising, radial cracks along the acting direction of vertical pressure will be promoted, and radial cracks in other directions will be prevented. Second, the effect of heterogeneous characteristics of material on granite fracture was discussed under different in situ stresses conditions, and numerical simulation results show that with the homogeneity coefficient decreasing, the shear failure zones around the borehole increase while the little effect can be observed about tensile failure characteristics. Thirdly, a simple three-dimensional numerical model was set to be compared with the results of two-dimensional cases. The same phenomena can be obtained that higher initial pressure leads to lower ratio of tensile failure elements, and in addition, it can be found that with the depth of cross section increasing, the ratio of the shear failure zone increases firstly and then decreases for three-dimensional cases.



---

*Source: 2579569-2021-09-15.xml*