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

State Key Laboratory of Hydraulics and Mountain River Engineering, School of Architecture and Environment, Sichuan University, Chengdu 610065, China Failure Mechanics & Engineering Disaster Prevention and Mitigation, Key Laboratory of Sichuan Province, Sichuan University, Chengdu 610065, China MOE Key Laboratory of Deep Earth Science and Engineering, College of Architecture and Environment, Sichuan University, Chengdu 610065, China


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. erefore, it is very significant to investigate rock fracture characteristics under in situ stresses [4], which can optimize blasting design. 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][9][10][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 to Weibull distribution in rock masses [28]. erefore, it will be used to correct the numerical models in this study.
In another aspect, the Damage index D, 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. erefore, 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.
e 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. e 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.

Dynamic Finite Volume Method.
e 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 Figure 2.
For node A in Figure 2, its accelerations € u in x-direction and € v in y-direction can be calculated as follows: where F x and F y 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 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: Based on the velocities of the nodes A, B, C, and D, the strain rate of the element ABCD can be obtained by According to the linear elastic constitutive, the components of the deviatoric stress tensor for the two-dimensional case can be calculated by

Shock and Vibration
S n+1 where G is the shear modulus, _ ε x , _ ε y , and _ ε xy are the components of Cauchy strain tensor, and _ e is the volume strain rate. erefore, the components of Cauchy stress tensor can be obtained as follows: 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: Applying Green's formula in the above integrating equations, the following equation can be obtained: 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 F x and F y can be expressed as follows: e 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.

Numerical Model.
In deep underground blasting engineering, the displacement of rocks along the direction of a borehole is usually strongly constrained. erefore, a numerical model based on the plane strain condition was established in this study, as shown in Figure 3, which consists of four parts: granite, copper, polyethene, and PETN. e granite model is 200 mm in width and height, and the radius of the borehole is 3.5 mm. e charge method is similar to research of Banadaki and Mohanty [31]. e 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. e 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.

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. erefore, 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 to Weibull distribution in the rock mass [28], the rock mechanics properties are also considered to conform to this distribution, which can be expressed by where f(x) is the probability density function of Weibull distribution, x 0 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: where P(x) 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 x 0 with the increase of m. e average dynamic mechanical parameters used in the granite model has been listed in Table 1 [31], and its density, bulk modulus, and dynamic shear modulus are 2.66 g/cm 3 , 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. en, 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 P n+1 and P n , and the number of the corresponding element N n can be calculated by 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.

Shock and Vibration
As mentioned above, there is little plastic strain when fracture happens for granite due to brittleness characteristics. erefore, a linear equation of state was used to express the relation of its deformation and pressure, which can be written as 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 where σ T and τ S are tensile strength and shear strength of granite, respectively, and σ 1 and σ 3 are the maximum and minimum principal stress. e 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.

JWL Equation for PETN.
e 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. e 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 where P is the explosive pressure, V is the special volume of explosion products, E is the initial total energy, and A, B, R 1 , R 2 , and ω are the constants of explosions. For PETN [35],

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 where c 0 and s are the material constants, μ � (ρ/ρ 0 ) − 1, ρ is the current density, ρ 0 is the initial density, c 0 is the Gruneisen constant (c 0 ≈ 2s 1 − 1), and E is the initial internal energy per unit volume. e material parameters of copper and polyethene are listed in Table 2.

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][37][38][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

Shock and Vibration
As shown in Figure 5(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. erefore, 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.

Application of In Situ Stresses.
When considering in situ stress conditions, the practical engineering problem will be simplified into the model, as shown in Figure 6. e lateral pressure P x and the vertical pressure P y 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. e solutions of stress distributions near the borehole can be deduced based on elasticity theory, as expressed in equations (17)- (20): 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.
It is significant to realize the same initial stress distributions in the numerical model as the actual state. e 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. en, the corresponding constant stress boundaries, such as P x � 20 MPa and P y � 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.
In Figure 7, 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). erefore, 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, 3a and 4a distance from the borehole center to record histories of hoop stress, as shown in Shock and Vibration 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.
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 to t � 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.

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. e first is to study the effect of the magnitude of initial pressure, in which initial lateral pressure P x is equal to initial vertical pressure P y and the magnitude of initial pressure are set as 10 MPa, 20 MPa, 30 MPa, and 40 MPa, respectively. e second is to explore the influence of pressure ratio (P x /P y ) on rock breakage. Initial lateral pressure P x was kept as 5 MPa, and initial vertical pressure P y is 5 MPa, 15 MPa, 25 MPa, and 35 MPa, respectively. A relatively homogeneous granite was considered, so the

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 Figure 9, 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. is is because that higher initial pressure leads to higher inner pressure, which suppresses radial cracks growing and propagating. ese phenomena can be seen clearly in the plots at 116.0 μs and 132.0 μs. ere 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. e final fracture pattern is shown in Figure 10, 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.
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 when P x � P y � 10 MPa. e 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. erefore, 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 Figure 12, 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 inplane shear failure (caused by τ rθ ); the shear stress τ rθ under blast can be calculated by   where σ e θ and σ e r are the hoop and the radial compressive stresses induced by blast stress waves, and σ P r and σ P θ are the initial hoop and radial stresses around the borehole induced by in situ stresses. According to dynamic elasticity theory, σ e r is much higher than σ e θ in magnitude, so (σ e r − σ e θ ) is negative. σ P θ is higher than σ P r as shown in Figure 8(b), so (σ P r − σ 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 (caused by τ rz ).
To verify the inference above, two stress gauges were set in the same locations of the pure shear zones of models of P x � P y � 10 MPa and P x � P y � 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 P x � P y � 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 P x � P y � 40 MPa, σ r > σ θ > σ z ; therefore, the max shear stress is τ rz so that out-plane shear failure will occur. 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.
In this section, rock fracture characteristics under blast loads when initial lateral pressure P x is equal to initial vertical pressure P y 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.

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. e initial lateral pressure P x is 5 MPa, and the initial vertical pressure is set as 5 MPa, 15 MPa, 25 MPa, and 35 MPa, respectively for four models. e 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 P y 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 P y 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 P y . e same phenomena can be observed in the contour plots at 116.0 μs and 132.0 μs. e numerical fracture plots are shown in Figure 16. ere are some clear differences for the models under the different ratios of initial lateral pressure and initial vertical pressure. e shape of pure shear failure tends to change into an ellipse with initial vertical pressure increasing. e 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. erefore, the main attention will focus on the influence of different pressure ratios on radial crack distributions.
To analyze the effect of different pressure ratios on radial cracks distribution, every model in Figure 16 is partitioned into six sector parts, and everyone accounts for 60 degrees. e 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 P y 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. It can be explained easily according to hoop stress distribution around the borehole, as shown in Figure 18. It can be found that the hoop compressive stresses along the borehole in the sectors from 330°to 30°and from 150°to 210°h ave a remarkable increase with the initial vertical pressure, which will weaken the strength of hoop tensile stresses of blast stress waves. erefore, 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.

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 twodimensional case. e heterogeneous characteristics of material usually have an influence on the results of numerical simulation. In this section, two kinds of homogeneity coefficient m � 20 and 100 were set to discuss how the heterogeneous characteristics of material affect granite fracture characteristics under blast stress wave. e 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.
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 twodimensional 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 Figure 21, and its height is 100 mm, and its cross section is a square of 100 mm in length. e borehole depth is 50 mm, and the coupling method is similar to the plane model. For 3D models, in situ stress conditions of P x � P y were just considered. e explosive was detonated at its top.
e blast results are shown in Figure 22. 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.

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, and Weibull distributions were applied in numerical models to characterize the inhomogeneity of rocks. e 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 inplane 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. irdly, a simple threedimensional numerical model was set to be compared with the results of two-dimensional cases. e 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.

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

Conflicts of Interest
e authors declare that they have no conflicts of interest.

Authors' Contributions
D. W., Z. Z., and W. G. conceptualized the study; D. W. and M. W. performed data curation and formal analysis; Z. Z. and M. W. were responsible for funding acquisition; D. W. and W. G. investigated the study and were responsible for methodology and software; D. W. prepared the original draft; Z. Z. and W. G. reviewed and edited the manuscript; D. W., Y. S., and M. W. visualized the study; M. W. performed study supervision; Y. S. and W. G. validated the data.