Numerical Study on Crack Propagation by Using Softening Model under Blasting

Amixed failure criterion, which combined themodifiedmaximumprincipal stress criterion with the damagemodel of tensile crack softening, was developed to simulate crack propagation of rock under blasting loads. In order to validate the proposed model, a set of blasting models with a crack and a borehole with different incident angles with the crack were established. By using this model, the property of crack propagation was investigated. The linear equation of state (EOS) was used for rock, and the JWL EOS was applied to the explosive. In order to validate the numerical simulation results, experiments by using PMMA (polymethyl methacrylate) with a crack and a borehole were carried out. The charge structure and incident angle of the blasting experimental model were the same as those in the numerical models. The experiment results agree with the numerical simulation results.


Introduction
Blasting is an economic and efficient excavation method and is widely used in engineering of quarrying, mining, and tunnel excavation.The propagation of blasting-induced stress wave in cracked rock mass would introduce crack expansion and failure, which would possibly induce destructively geologic hazards, such as rock burst and slope slump.The propagation property of stress wave in cracked rock mass and its dynamic mechanical response have drawn massive attention since plentiful joints and cracks exist within rock mass.It is an important topic in rock dynamics and relevant disciplines to study the way cracks propagate under dynamic loading so as to know how the wing cracks develop at crack tips [1] as the case in static compressive loading, which has been well investigated by using complex function method [2][3][4][5].This research helps with predicting the dynamic strength and structural stability of cracked rock mass and provides theoretical basis for blasting design, improving blasting efficiency.
The time taken by the explosive loading is extremely short, ranging, usually, from a couple of microseconds to dozens of microseconds.That renders the recording of the whole physical field impossible experimentally, whereas the state of the field can be described by means of numerical simulation, which is less expensive, and can be easily realized.However, the numerical simulation accuracy depends largely on the quality of the numerical model employed.Hence, a combined use of numerical simulation and experiment was carried out in this paper.In the numerical study, Preece and Thorne [6] used 3D finite element techniques and a damage constitutive model to study the detonation timing and fragmentation.Donzé et al. [7] applied a model based on the discrete element to investigate the importance of stress waves on the initiation and propagation of radial fractures during the dynamic loading phase of explosion.By using UEDC code, Chen and Zhao [8] have simulated blasting wave propagation in joint rock mass.Using AUTODYN code and the modified maximum principal stress, Zhu et al. [9,10] established dynamic computation model for rock under dynamic loading, and by using this model, the crater blasting [11] and the zonal disintegration phenomenon [12] have been investigated.However, brittle materials, such as rock, act in inelastic way, which is caused by the nucleation, propagation, and connection of microcracks.Softening models are often used to describe the inelasticity of brittle material [13,14].Feenstra and De Borst [15,16] used the softening model based on 2 Advances in Materials Science and Engineering Rankine plasticity criterion to simulate the dynamic fracture process of concrete.In this paper, the softening model applied is a function of fracture energy, plastic strain, and the grid size, and it can limit the results of the calculation of grid sensitivity to a certain extent.In experimental study, Rossmanith et al. [17] proved that the failure mechanism of PMMA is similar to that of rock under dynamic loading.The PMMA possesses high strength and transmittance.Its fracture modes can be easily observed.The dynamic fracture mechanism of it has been widely studied by researchers [18][19][20].
A numerical code, AUTODYN [21,22], is applied in this study.It is an explicit finite difference code for solving a variety of nonlinear problems in solid, fluid, and gas dynamics.AUTODYN code has been successfully applied in the study of rock fracturing by Zhu et al. [9][10][11][12].Based on the computational model proposed by Zhu et al. [9,10] and the softening damage model for tensile cracking based on Rankine plasticity criterion, a crack propagation model of crack propagation under blasting loads was established to simulate the crack initiation and propagation and the impact of the orientation angle between the borehole and the preexisting crack.PMMA plates with an existing preset crack were used in blasting experiment, where the charge structures and the distance between the borehole and the crack tip were kept as constants, while the orientation angle between the borehole and the crack varied.The results from blasting experiment and those from the simulation were compared.

SDMTC (Softening Damage Model for Tensile Cracking)
Before the description of SDMTC (softening damage model for tensile cracking), we need to introduce JMFDFV (the joint method of finite difference and finite volume) for the SDMTC employed in this paper.

Brief Introduction of JMFDFV.
The JMFDFV introduced in this paper is a numerical method based on the Lagrangian description and can be used to solve solid kinetic problems.Figure 1 is a set of meshes for rectangular element cell under dynamic loading.1∼4 denote cell number;  ∼  denote node number.The quantities, pressure (), density (), stress (), strain (), strain rate ( ε ), temperature (), and mass (), are defined within the element cell and remain unchanged in a time cycle.The quantities, coordinate (, ), displacement (, ), speed ( u , V ), acceleration ( ü , V ), and node force (  ,   ), are defined on nodes.According to Newton's second law of motion, the acceleration components of node , ü and V , can be written as where  denotes the total mass of domain .The node speed at time  + 1/2 can be written as where  denotes the number of time step and Δ denotes time step.The node displacement at time  + 1 can be written as By making use of Green's integral formula and the four node speeds of element , the strain rates of element , ε  , ε  , and γ  , can be written as where   denotes the area of element .Using the strain rates in (4), the stress deviators,   ,   , and   , in step  + 1 can be obtained as where  denotes shear modulus and ė volume strain rate.In 2D cases, ė = ε  + ε  .
Advances in Materials Science and Engineering 3 The stresses of element  can be written as where  denotes mean pressure and is determined by the material EOS in the model.The node force is determined by integral on both ends of the equations of wave propagation: According to Green's integral formula, ( 7) can be rewritten as where  denotes the area of domain  and  is the closed curve .The following equations can be obtained using ( 7) and ( 8): Taking into consideration (9) and the mesh in Figure 1, the following finite difference equations for 2D can be obtained as where superscript 1, 2, 3, and 4 denote elements 1, 2, 3, and 4 in Figure 1, respectively.In practice, besides the type of forces in (10), there are other nodal forces that should be considered.They are pseudoviscous forces, hourglass damping, and antitangle forces, and if the node is located on the boundaries of the load, external forces should be incorporated.The computational method adopted in this paper used the numerical method for pseudoviscous force proposed by Wilkins [23].Adopt the method proposed by Century Dynamics Inc. [22] and the hourglass damping is determined, namely, adding a set of corrective forces in the four nodes of the cell.The external forces are divided equally between the two nodes on the boundary, and the boundary forces due to external forces are therefore determined.Above all, the nodal force is a resultant force.The process starting from (1) to ( 10) is a cycle of the computation.The nodal acceleration is calculated after calculating the nodal force.By iteration of this calculation process, the quantities, node speed, displacement, acceleration, element stress, and strain rate of each node and element in different times can be obtained.

SDMTC.
Basic physical quantities can be obtained by the computational method introduced above.To describe the failure process in reality, however, an effective failure criterion is needed.Thus MMTS (the modified maximum tensile stress) is adopted in this paper.The MMTS is an instantaneous failure law.The meaning can be expressed as follows: element loses its ability to bear tensile and shear stresses, when the principal stress of the calculating element reaches the tensile strength or when the maximum shear stress reaches the shear strength of the material, but its ability to bear compression is maintained.The law can be written as where  1 ,  3 denote two principal stresses,  max denotes the maximum shear stress, and   ,   denote the tensile strength and shear strength of material, respectively.In this paper, stress, strain, and strain rate are defined and kept unchanged within element in a time cycle.The numerical method using MMTS can easily lead to the dependence of calculation results on grid.In fact, element loses its bearing capacity gradually as crack propagates through it.Figure 2 sheds light on the gradual failure process of an element under static loading; that is,  max increases to a certain value that is larger than or equal to the tensile strength   before it decreases to a value smaller than   and reaches 0. According to the SIF (stress intensity factor) of edge crack and the judgment formula for critical stress in LEFM (linear elastic fracture mechanics), the relationship between critical stress of element and crack length  can be listed as where  IC denotes fracture toughness,   fracture energy,  elastic modulus, and  crack length.In order to illustrate the process of gradual failure for an element, the author adopted the linear SDMTC algorithm and Figure 3 shows the stress-strain curve of this method.
In this model, the MMTS determines the initiation stage of element failure, after which the maximum tensile principal   stress  max , crack strain  cr , and the crack strain   for the totally fractured element are recorded, respectively.  is defined as a function of   , dynamic tensile strength   , and the critical length  in the direction of maximum tensile principal stress: The damage variable Dam is used to describe the damage level, defined as a function of  cr and   , where  cr is obtained by the Backward-Euler method: In damage process, the maximum tensile principal stress is defined as a function of dynamic tensile strength   and damage variable Dam:

The Calculation of Crack Strain.
The calculation of crack strain is the key to the realization of SDMTC.The crack strain calculated by Backward-Euler method must reflect and incorporate the essence of material failure.Brittle damage is tensile failure developed along the direction of principal stress.Hence, Rankine's plastic yield criterion is chosen: where  1 ,  3 denote the principal stresses and σ (  ) denotes the equivalent stress.This paper uses Backward-Euler method to obtain the crack strain.When the trial stress surpasses the yield surface, associated Backward-Euler method is employed to back-flow the stress to the yield surface.This method is similar to that described by Crisfield [24] showed in Figure 4. Associated flow is used in the meridian plane and -plane.

Numerical Blasting Model
Using the numerical method introduced above and PMMA parameters obtained by experiment, the 2D numerical blasting model shown in Figure 5 is established.

Blasting Model.
In the numerical blasting model shown in Figure 5, the decoupling charge is adopted, decoupling coefficient being set as 2. Dimension of the model is 400 mm × 400 mm, with one borehole and one thorough centric crack.To delve into the impact of the orientation angle between crack and borehole on crack propagation, a set of models, with orientation angle being 0 ∘ , 15 ∘ , 30 ∘ , 45 ∘ , 60 ∘ , or 75 ∘ , are established.Free boundary is chosen as the boundary condition.The charge and PMMA are set with different subgrids, and the mesh near the borehole is shown in Figure 6.The mesh size determines the computational efficiency and the accuracy of result.However, it is difficult to estimate the sensitivity of our approach.Usually, we determine the best mesh size by trying different mesh sizes, like the method used by Y. Hao and H. Hao [25].The total number of elements is more than 160000 and the maximum mesh size is less  than 2 mm.These meshes are seen to be fine enough for the required precision of calculation.With additional details on the interactions between subgrids, such as gap, one type of Lagrange-Lagrange coupling can be found in [21,22].No mesh is set for air for the inaccessibility of computation caused by the large blasting-induced deformation.

Material Model.
PMMA falls into the catalog of isotropic material, has been widely studied, shows similar characteristics to rock under dynamic loading, machines easily, and has excellent transparency for crack observation.In accordance with the charge amount, the linear EOS is adopted for PMMA [22]: where  denotes pressure,  bulk modulus, and / 0 the ratio of present density and initial density.The JMFDFV and MMTS are used as the failure models for PMMA with the parameters listed in Table 1.
The dynamic strength of PMMA is difficult to determine, because it changes largely at different environment temperature and strain rate [27].The shear strength is determined from the charge amount and experiment results.The static values [28] of tensile strength and GC were adopted and the tensile strength was modified experientially because later experiment was taken at temperature 313 K.The dynamic strength parameters of PMMA are listed in Table 2.
The JWL (Johns-Wilkins-Lee) EOS (equation of state) has been adopted for the ANFO charge: where  denotes pressure,  denotes the initial internal energy of denotation product,  denotes relative volume, and , ,  1 ,  2 , and  denote material constants.With strength model and failure criterion ignored, the charge will turn into ideal gas after blasting.The material parameters for ANFO (ammonium nitrate/fuel oil) are listed in Table 3.

Simulation Results and Analysis
It is well known that the blasting load in the borehole is p wave in 2D case and p wave scattering near the crack is shown in Figure 7. Crack blocks the spread of stress wave, diffracted wave generates at the crack tip, and reflected wave generates on the crack surface.As you can image, when the incident wave impinges the crack with a certain angle, the crack tip nearer to the borehole cannot completely block the incident wave.At the crack tip on the far end, there exists only the diffraction wave on the side of the crack surface due to the crack surface blocking.

Crack Initiation at the Far End for Incident Angular
Degree 0 ∘ .Only a single gauge point is set in the vicinity of crack tip, gauge #3 shown in Figure 5, to analyze the cause for the initiation at far end, because the loading is symmetrical for the case when the incident angle is 0 ∘ .Figure 8(a) shows the time history of the major principal stress.It can be seen that the major principal stress turns from compression to tensile stress, and as it reaches the tensile strength, material is damaged and no tensile stress is sustained.This is because when the incident angle is 0 ∘ , neither reflection nor diffraction is aroused by stress wave.The driving force for crack initiation and propagation is the same as the case without a crack, that is, the circular tensile stress, which is generated by the radial movement of materials, result of compression pressure on the borehole.The circular tensile stress gives rise to crack initiation and propagation.compression stress encounters the crack tip, diffraction wave generates.The general wave field consists of both incident wave and diffraction wave.The domain in the crack back is in a stress-free state; the only acting wave is diffraction wave composed of tensile wave and shear wave.Hence, the crack initiation at the far end begins in the direction away from the borehole.

The Impact of Crack Orientation on Crack Initiation at Far
End.Based on the initial analysis of the crack tip at far end for incident angular degrees 0 ∘ and 45 ∘ , two conclusions are reached.First, when the incident angular degree is relatively small, the cause for the crack initiation at far end is mainly because of the circular tensile stress.While the incident angular degree grows, the diffraction of stress wave plays the major role.Figure 9 compares the damage status of the models, in angular degree range 0-75 ∘ , 80 s after detonation.As the incident angle increases, the initiation angle at far end increases.When the angular degree reaches 75 ∘ , crack initiation takes form at a distance from the preset crack tip.

Comparison of Experiment and Simulation Results.
In order to verify the validation of the results from numerical simulation, corresponding blasting experiment was carried out in the blasting laboratory, Southwest University of Science and Technology, using PMMA.The environment temperature is 313 K. Specimen size is 400 * 400 * 8 mm and ANFO is employed as the explosive.The experiment results were listed in Figure 10.This figure is part of the specimen and it was taken after explosion.Dark paper was used as the background paper for photographing (we put the specimen on a dark paper, and then we take it using camera, Canon power shot A3300IS, vertical to specimen.).It can be seen that the angular degree for crack initiation at far end is less than 90 ∘ when crack angle falls in 0 ∘ -30 ∘ range.When crack angle stays in range 45 ∘ -75 ∘ , the angular degree for crack initiation is a bit larger than 90 ∘ .The experiment results illustrated in Figure 9 are analogous to the simulation results in Figure 10, meaning the simulation is reliable.

Figure 2 :
Figure 2: The maximum tensile stress varies with crack propagation.

Figure 6 :
Figure 6: Mesh of the explosive and the PMMA.

Figure 7 :
Figure 7: p wave scattering near the crack.

Figure 8 :
Figure 8: Time history curves of major principal stress at gauges for incident angle of (a) 0 ∘ and (b) 45 ∘ .

Figure 9 :
Figure 9: Material status of rock samples as a function of degree  at t = 80 us after detonation in the borehole.

Table 2 :
Dynamic strength parameters of PMMA.