Numerical and experimental simulation of damaged rock with randomly oriented cracks by shock disturbance

The aim of this study is to investigate the effect of shock-disturbed cracks on the dynamic fragmentation of granite. Considering the complex behavior of rock materials, the Walsh’s model was revisited and extended by including the stress effect required to close an initially open crack and examining the unloading process in detail. This analysis leads to closed-form expressions for loading and unloading portions of the effective Young’s modulus, as functions of the crack density, characteristic aspect ratio, and crack friction coefficient. Subsequently, the effective Young’s modulus and cutting force are simulated and the influence of cracks is studied. The analysis results with different crack density and disturbed frequency are compared in terms of effective Young’s modulus and cutting force. Finally, the tool and damaged rock model with randomly oriented cracks by shock disturbed at a different frequency was demonstrated by the test. The good agreement between the simulation results and experimental data demonstrates the validity of the simulation method.


Introduction
A large diameter rotary drilling rig is often used in mining of large-scale surface mines.However, this rig is inefficient, and its cone bits experience serious wear when used in hard rock metal mines.Therefore, the new technology of shock-induced drilling is presented, and the principle is shown in Fig. 1.The strength of the rock mass is weakened by the redistribution of stress around the induced hole, so cutting a rock mass in the damaged zone becomes easy [1][2][3].
To study the cutting features of the damaged rock, constructing the tool-damaged rock model and the damaged rock-crack model is necessary to predict the cutting force on the damaged rock by shock disturbance.
A considerable effort has been made recently by researchers to numerically simulate the tool-rock interaction and rock fracture response during percussive drilling.Forquin and Hild [4] used a continuum approach to study the dynamic fragmentation in a wide range of brittle materials because of impact loadings.A plasticity model introduced by Swenson [5] is employed together with an anisotropic damage model developed by Denoual [6], Forquin and Hild [4].This model was used to describe the fragmentation of two grades of limestone subjected to edge-on impact tests [7] and the strength improvement and failure pattern of microconcrete [8].Saadati et al. [9] applied the model on granite and investigated the rock fragmentation process and the force-penetration response at percussive drilling.It was shown that the plasticity model is an appropriate tool for this purpose as it deals with both the dynamic fragmentation as a result of tensile stress and the plasticity-like deformation in compression.However, the finite element method based on the continuum mechanics theory of material simulation has serious problems in representing properly discontinuities of the material occurring during rock cutting.
The discrete element method can take into account most of kinds of discontinuities and material failure characterized with multiple fractures which make it a suitable tool to study rock cutting.Among the numerical studies that are performed in 2D, a work of Liu [10] and Wang et al. [11] can be mentioned.Both of these studies are performed using in-house numerical codes and do not account for inelastic strains in the rock.Another example of 2D numerical simulation of percussive drilling is the work by Saksala [12] where a viscoplasticity consistency model and an isotropic damage description are coupled to take both the compressive and tensile responses of a rock into account.The model has been improved recently to deal with 3D simulations [13].Furthermore, Chiang et al. [14] also used a 3D finite element (FE) approach to simulate the rock fragmentation during the drilling process by using a linear Mohr envelope together with a tension cutoff plane.Among the studies using discontinuous mechanics approach, Thuro et al. [15] used a Particle Flow Code (PFC2D) to investigate the crack pattern in drilling and its correlation with existing foliations.However, this research did not consider the effect of cracks on the cutting force.
Various research studies of a cracks analysis for the damage rock have been reported.In 1965, Walsh [16,17] published a set of three papers that provide the conceptual basis of the influence of cracks on the elastic rock deformation.Under hydrostatic loading [16], open cracks (i.e., at low stresses) initially contribute an excess compliance to the rock.As the stress increased, each crack closes up at a stress that is roughly equal to , where is the Young's modulus of the uncracked rock, and is the initial aspect ratio of the crack.After a crack is closed, it is assumed not to contribute to rock behavior under hydrostatic loading, because such loading is assumed to cause no resolved shear stress on a closed crack face.This model has been extended [18] by adding a cohesion term to the frictional constitutive law for the cracks, thereby allowing the model to predict yield-type behavior, but again retained the simplification of assuming that the compressive stress required to initially close a crack is negligible.However, this research did not consider the effect of disturbance frequency on the cracks.
In the present investigation, the effect of the presence of cracks, which are produced by the former impact of the drill bit in the damage rock, is studied in detail and the simulations are compared with experimental observations.The different densities and orientations of the cracks are introduced in a numerical model, and their effect on the Young's modulus is studied.Furthermore, the numerical tool is used to explore the influence of cracks on cutting force at different disturbance frequency.

Shock-induced drilling principle
The principle of shock-induced drilling is shown in Fig. 1.First, the impactor is used for drilling the center hole; the strength of rock mass is weakened by the redistribution of stress around the induced center hole.Second, the little cutting force on the damage zone can crush the rock. is the axial pressure, is the impact load, and is the rotary torque.Fig. 2 is the force analysis diagram of the damaged rock around the induced hole during tooth intrusion.is the radius of the induced hole, is the radius of the tooth intrusion, − is the radius of the weakened zone of rock around the induced hole.Fig. 3 is the enlarged image of region C in Fig. 2, Fig. 3(a) is the cracks' propagation conditions after the drill tooth penetration, Fig. 3(b) is the model of the cracks' propagation during the tooth intrusion.Fig. 4(a) is the structure of the crack 1 in Fig. 2(b), assuming that it is constituted by three different shapes of elliptical cracks (because of the elliptical crack has two surfaces and the crack tip, so it's similar to the actual crack), is the half-length of elliptical crack 1, is the half-length of the elliptical crack 1, is the angle of the elliptical crack 1, is the angle between elliptical crack 1 and elliptical crack 2, is the center distance between the elliptical crack 1 and elliptical crack 2, changing the parameters and numbers of the elliptical crack, different shapes' cracks can be obtained.Fig. 4 (b) is the force analysis of single elliptical crack around the induced hole during tooth intrusion.
The dimensionless functional relationships that link macroscopic and microscopic parameters can be postulated for a discrete element model [13]: = Φ , .
The effective Young's modulus of the specimen can be defined by: where is Young's modulus of the uncracked rock, is a certain length parameter, is a characteristic area related to the discrete element model, , are the contact stiffness in the normal and tangential direction respectively, , are the interface strength in the normal and tangential direction respectively, Δ is the excess energy stored in the rock due to the presence of a crack and ℎ is the volume of the plate.
If the crack is open and the plane, in which it lies (not its normal vector), is oriented at an angle to the direction of the applied stress, under the assumption of the plane stress, the excess energy term is equal to [17]: This result can be transformed into one appropriate for plane strain by replacing by /(1 − ), where is the Poisson's ratio of the uncracked rock.To simplify the notation, the case of plane stress is considered, where does not appear.However, the term (1 − ) is usually very close to the unity, in any event.Hence, the effective modulus of a rock containing this single crack will be: where = / ℎ is the crack density parameter.The case is considered when the crack is closed; the two opposing faces are thus in contact and sliding past each other.The energy, supplied by the externally applied stress that is required to cause the crack faces to slide, is given by [17]: where sin cos = is the resolved shear stress along the crack plane and is the frictional stress that resists the sliding.Adopting the simplest law of sliding friction, = is assumed, where is the friction coefficient, and is the resolved normal stress acting on the crack surface.Following the Walsh's argument, if is the normal stress necessary to cause the crack faces to close up and come into contact, then the actual effective value of the resolved normal stress acting on the crack will be = − .However, = sin ; thus, in combining all of these relations, the difference between the resolved shear stress and the frictional resistive shear stress is: In conjunction with the general Eq. ( 4), the effective elastic modulus for the rock containing a single sliding crack will be: = 1 + 2 sin cos − sin + sin cos . (10)

Open or closed crack criterion
Although the relations of the effects on the elastic modulus of an open crack or a closed sliding crack are obtained, a criterion is derived to decide whether a given crack is open or closed, or to choose if it is sliding when closed.The normal stress is required to close an elliptical crack whose initial aspect ratio is /2, where is defined to be the ratio of the minor axis to major axis, so that < 1.Furthermore, is defined as the initial aspect ratio of the crack at zero stress.The original Walsh's model and other aforementioned models correspond to the idealized case of = 0.
The resolved normal stress on the face of the crack is = sin ; thus, the criterion for a crack to be open at an applied compressive stress is < = /2.Hence, at a given stress , a crack will be open if < , and closed if > , where the critical angle is given by: Cracks at all orientations will be open at stresses below , which is at normalized stresses below 1.When = 1, ( ) = /2, only the cracks that are precisely normal to the loading direction will close.As increases, cracks that are oriented at smaller angles to the loading axis will close, because Eq. (11) shows that the critical angle ( ) decreases as increases.

Sliding crack criterion
For a crack to be sliding, it must already be closed.Hence, a necessary condition for sliding is > .However, the sliding can occur only if the resolved shear stress exceeds the frictional resistive stress, which is from Eq. ( 8): Using standard trigonometric identities, Eq. ( 12) can be rearranged as follows: Eq. ( 14) can be rearranged as follows: A manipulation of Eq. ( 13) shows that, when = 1 , = /2 , and when → ∞ , → arctan(1/ ).Moreover, for all normalized stresses > 1, the critical angle for sliding, , as given by Eq. ( 13), exceeds the critical angle for closure, , as given by Eq. (11).
To summarize, if < 1, a crack at any orientation to the direction of externally applied stress will be open.If > 1 and if < , cracks will be open, and closed if > , where the critical angle for the closure is given by Eq. (11).A closed crack will be sliding if < < , and will not be sliding if < < /2, where the critical angle for sliding is given by Eq. ( 13).The situation is illustrated in Fig. 5.

Loading and unloading of rock containing randomly oriented cracks
The rock that contains a distribution of cracks is supposed, each with the same size and with their orientation angles uniformly distributed.By symmetry, only the range of values 0 ≤ ≤ /2 is considered.To simplify the calculations, the no-interaction approximation is invoked, in which the excess energy because of each crack is computed as if it were an isolated crack in an infinite body [17].This approximation is very accurate at low crack densities, but becomes gradually inaccurate as crack density increases.
In the regime < 1, all cracks remain open.If the interactions between cracks is continually ignored, the contributions of each crack to the compliance will be additive, and thus the right-hand side of Eq. ( 6) over the range 0 ≤ ≤ /2, is averaged to find: where = / is the standard two-dimensional crack density parameter.If a normalized elastic modulus is defined as = / , then the normalized modulus for a body containing randomly distributed open cracks will be = 1/(1 + ).
For > 1, the right-hand term in Eq. ( 6) must be integrated for those cracks that are open, < , and the right-hand side of Eq. ( 10) for those closed cracks that are sliding, < < .As mentioned above, the cracks, for which > will be closed and not sliding do not contribute to the energy budget.The temporary notation 1/ = 1 + + is introduced, where the terms are compliances due to the open and sliding cracks.
From Eq. ( 6), the compliance because of open cracks when > 1 is: represents the total crack density; the proportions of cracks that are open or are sliding are accounted for by the stress-dependent limits of integration.
From Eq. ( 10), the compliance because of the sliding cracks is given by: Thus, the full expression for the elastic modulus is given by: The applied compressive axial stress is assumed to be increased to a certain value , after which it begins to decrease.A decrease in the compressive stress can be considered as the super-position of a tensile stress or a stress of the opposite sign.If the applied stress is reduced from to a certain value ; where = − Δ , an associated decrease Δ will occur in the shear stress acting on the crack as well as a decrease Δ will occur in the frictional resistive stress.Reverse sliding will occur if [17]: The energy of the reverse sliding is given by: where: So Eq. ( 18) can be rearranged in terms of the normalized stresses as follows: Using standard trigonometric identities, Eq. ( 18) can be rearranged as follows: Therefore, the condition for reverse sliding can be summarized as follows: where the critical (maximum) angle for reverse sliding is given by: The excess compliance contributed by the cracks that are undergoing reverse sliding will be given by: The full expression for the elastic modulus is now given by:

Simulation and test results
The bonded particle model of the rock material was built in the particle flow code software (PFC2D).Because the particle flow samples could reflect the macroscopic mechanical properties of real rock materials, a calibration of the microscopic parameters of particle model was necessary.The simulated system is shown in Fig. 6.Fig. 6(a) is the cutting simulation of the whole cone bit; Fig. 6(b) is the cutting simulation of a single drill tooth.In Fig. 6(a), the right and bottom boundary is fixed in the calculation process, while the left boundary is servo-controlled to provide periodic disturbance force in simulations.The rock material in this study was granite purchased from the Yuelu Mountain in the Hunan province.Considering great differences in stiffness between particles and the cutter, the cutter with the same size as that in experimental tests is simplified into rigid walls.The set of microscopic parameters of the simulation is: = 150 N, = 121 N, = = 3.6×10 6 N/m, = 0.5.Table 1 shows the comparison of the simulated values and the macro parameters of real rock material.The error between the simulated values and the real values is very negligible.The experimental system is shown in Fig. 7. Fig. 7(a) and 7(b) shows the experimental principle and the experimental system, respectively.The specimen size is 400 mm×400 mm×500 mm (length×width×height).Before the test, the specimen surface is polished by a grinder.The induced hole is then drilled with an impactor at different frequencies of 0, 10, and 15 Hz.The diameter of the center hole is 100 mm, and the depth of the hole is 400 mm.Then the specimen is pushed to the right place by an oil cylinder.Finally, the cutting force of a single tooth around the center hole was tested under the same axial loads.The simulated and experimental results are shown in Figs.8- (1) Fig. 8 represents, as a function of the normalized stress, the evolution of the critical angles for closure, sliding, and reverse sliding, as given by Eqs. ( 11), (13), and (19), respectively.The values of these critical angles are not given for < 1, because = 1 is the stress required to close the first crack.The friction coefficient is considered as = 0.6.The cracks, whose orientation is below the closure criteria, can be seen as open cracks.Closure proceeds at a rate such that, for example, 70 % of the cracks are already closed when = 4. Nevertheless, a significant portion of the cracks remain open even during high stresses.The cracks mobilized by frictional sliding during loading are those cracks lying between the critical angles for closure and sliding (or reverse sliding during unloading).Reverse sliding does not occur at the beginning of unloading.Moreover, Eq. ( 19) can easily show that as → , → ( )and it does not depend on the friction coefficient.Fewer cracks are mobilized more often during unloading than during loading.
(2) When the disturbance force is constant, the cracks of the region A in Fig. 6 under the condition of different disturbance frequency, shown in Fig. 9. Fig. 9 (3) Fig. 10 shows the stress dependence of the effective modulus and the resulting stress-strain curves for the three values of the crack density, = 0.2, 0.5, and 0.8, for a fixed value of the friction coefficient, = 0.6.Fig. 10 shows the equivalent plots for three values of the friction coefficient = 0.3, 0.5, and 0.8, for a fixed value of the crack density = 0.5.
Considering that no-interaction approximation is used for calculating the effective modulus, the compliance is a linear function of crack density.Hence, if all other parameters are held constant, the compliance (1/ ) will be a linear function of .The resulting effect of increasing the crack density is more pronounced compliance, more non-linearity, and a larger hysteresis loop, as illustrated in Fig. 10.The only parameter whose influence cannot be seen clearly from the equations is the friction coefficient .Physically, the effect of sliding is expected to decrease with increasing values of , thus resulting in increased stiffness and less pronounced non-linearity and hysteresis, as illustrated in Fig. 11.The extreme case of = ∞, although not physically realistic, is not shown in Fig. 11, but would result in a perfectly reversible stress-strain curve, without any hysteresis.
Research studies commonly assert that after a phase of non-linearity caused by the closure of the cracks, the subsequent nearly linear regime is representative of the compression of an uncracked elastic body.This result is a good approximation under a hydrostatic stress, but not under impact loading.For realistic values of the friction coefficient, as shown in Figs. 10 and 11, the frictional sliding of the cracks provides the excess compliance and causes a non-linear behavior, although the nonlinearity is less pronounced due to the closure of the cracks.For high values of stress, even when the cracks are mainly closed, the effective modulus during loading is never equal to the uncracked modulus, and has a value of approximately 70 %-90 % of .The uncracked rock modulus can only be inferred from the slope of the stress-strain curve at the beginning of unloading, because this slope does not depend on any of the microstructural parameters with the cracks being closed and stuck.
When the invasion load is constant, the drill tooth cutting phenomenon under different disturbance frequencies is shown in Fig. 12. Fig. 12( In Fig. 13, Fracture condition of the granite under different disturbance frequencies of the impactor are shown.It can be obtained that the volume and the scope of the fractured rock are increasing with the increase of the disturbance frequency.This shows that the crack density is increasing with the increase of the disturbance frequency, and the elastic modulus decreases with the increase of the crack density, so the drill teeth intrusion more easily and the volume of the fractured rock under the same intrusive force is bigger.
The cutting forces between the simulation and experiment under different disturbance frequencies are shown in Fig. 14.In the experiment, the cutting force decreased when the frequency increase was greater than the values in the simulation results.The difference occurs due to the nonlinearity of rock materials.The peak values of the different frequencies in the simulation were 60, 30, and 25 KN, respectively, whereas in the test, these were 80, 40, and 20 KN, respectively.The fluctuation of the cutting force is caused by shock disturbance.There are nearly 25 % errors between the simulation and experimental data.These errors are caused by the differences between the calculated values and the actual values.
From the aforementioned analysis, the error of the simulation is acceptable for the application in this present paper.As such, the model can be used to analyze and predict the cutting force of the damaged hard rock.

Conclusions
A new technology of shock-induced drilling was presented, and a crack model of the damaged rock around the induced hole was built.Subsequently, the effective elastic modulus of a damaged rock under the condition of the impact loads was analyzed, and the results showed that the effective elastic modulus significantly decreased with the increase of the crack density.
A numerical simulation and laboratory results, disturbance frequency of the cutting force of the damaged rock showed that the cutting force decreased with the increase of the disturbance frequency.When the disturbance frequency was 15 Hz, the cutting force was about one third of the force without disturbance.
The numerical simulation could reveal the cutting phenomenon.However, an error between numerical results and the experimental measurements still existed.Thus, the developed model is required to be further verified.More comprehensive experimental testing procedure is needed to calibrate the model, including a triaxial compression test, a direct shear test and a fracture mechanics test.The use of these test results for a model calibration could allow better model tuning and the subsequent improvement of the qualitative results during rock cutting.

Fig. 1 . 2 .
Fig. 1.Shock induced drilling principle of Fig. 2. Force analysis of rock around induced hole during tooth intrusion

Fig. 5 .
Fig. 5. Cracks at different orientations, under given compressive stress and impact loads

Fig. 6 .
Numerical simulation of rock cutting test a) numerical simulation principle b) cutting phenomenon of single drill tooth

Fig. 7 .Fig. 8 .
Fig. 8. Relationship between normalized stress and critical anglesThe comparisons of the experimental and simulation results are shown in Figs.8-14.The following conclusions can be made:(1) Fig.8represents, as a function of the normalized stress, the evolution of the critical angles for closure, sliding, and reverse sliding, as given by Eqs.(11),(13), and (19), respectively.The values of these critical angles are not given for < 1, because = 1 is the stress required to close the first crack.The friction coefficient is considered as = 0.6.The cracks, whose orientation is below the closure criteria, can be seen as open cracks.Closure proceeds at a rate such that, for example, 70 % of the cracks are already closed when = 4. Nevertheless, a significant portion of the cracks remain open even during high stresses.The cracks mobilized by frictional sliding during loading are those cracks lying between the critical angles for closure and sliding (or reverse sliding during unloading).Reverse sliding does not occur at the beginning of unloading.Moreover, Eq. (19) can easily show that as → , → ( )and it does not depend on the friction coefficient.Fewer cracks are mobilized more often during unloading than during loading.(2)When the disturbance force is constant, the cracks of the region A in Fig.6under the condition of different disturbance frequency, shown in Fig.9.Fig.9(a), Fig. 9(b) and Fig. 9(c), are the cracks' propagation conditions of the region A under the condition of 0, 10 and 15 Hz.In contrast to Fig. 9(a), Fig. 9(b), and Fig. 9(c), the cracks of the region A increased as the disturbance frequency increased.This shows that the elastic modulus of granite decreased with an increase of the disturbance frequency.

Fig. 9 .
(a), Fig. 9(b) and Fig. 9(c), are the cracks' propagation conditions of the region A under the condition of 0, 10 and 15 Hz.In contrast to Fig. 9(a), Fig. 9(b), and Fig. 9(c), the cracks of the region A increased as the disturbance frequency increased.This shows that the elastic modulus of granite decreased with an increase of the disturbance frequency.Cracks' propagation conditions of region A in Fig. 6 under different disturbance frequencies

Fig. 12 .Fig. 13 .
a), Fig. 12(b) and Fig. 12(c) are the drill tooth cutting phenomenon under the condition of 0, 10 and 15 Hz.In contrast to the Fig. 12(a), in the Fig. 12(b) and Fig. 12(c), it can be obtained that the cracks' density increases as the disturbance frequency increases, the volume and the scope of the fractured rock are also increasing with the increase of the disturbance frequency.It is shown that the efficiency of the shock-induced drilling could be improved by increasing the disturbance frequency.Cutting phenomenon of drill tooth under different disturbance frequencies aFracture condition of the granite under different disturbance frequencies a) Simulated results b) Experimental results Fig. 14.Comparison between simulation and experimental results of cutting force under different disturbance frequency

Table 1 .
Comparison of simulation values and macro parameters of real rock material