Discrete element simulation for investigating fragmentation mechanism of hard rock under ultrasonic vibration loading

Assisted ultrasonic vibration technique can significantly improve the efficiency of hard rock drilling in petroleum and mineral engineering. In this study, to determine the fragmentation mechanism of rocks under ultrasonic vibration, numerical simulations using the discrete element method (DEM) were performed. A novel flat‐joint model (FJM), combined with an ultra‐high‐frequency loading boundary condition, was used to model the damage process of the hard rock under ultrasonic vibration loading. The numerical results demonstrated that the evolution of local strain and fragmentation were in good agreement with the experimental results. Based on the established model, the influence of loading parameters was investigated. Furthermore, by analyzing the development of the full strain field, crack orientations, and crack distribution, the fragmentation mechanism was revealed for the rock under ultrasonic vibration. Under ultra‐high‐frequency loading, the rock deformed in a heterogeneous manner and produced both compressive and tensile strain zones. The compressive zones were mainly distributed in the fringe and tensile zones in the top center. The generated tensile cracks caused by compression and tension in these two strain zones led to the rock failure.


| INTRODUCTION
With the gradual depletion of shallow resources, the energy exploration depth has been increasing and condition becomes complex. 1 Consequently, the explorations face increasingly harder rocks with higher strengths. Breaking hard rock rapidly in petroleum and mineral exploration is a popular and significant topic. [2][3][4][5][6] However, the current technologies cannot meet the demands of quick drilling. It is necessary to develop a more efficient technology for hard rock exploration. To achieve this goal, the assisted ultrasonic vibration technique (ultra-high-frequency cyclic loading), combined with the traditional drilling method, was first introduced by researchers at Aberdeen University. 7 To implement this technology more effectively, studies were conducted on rock-breaking characteristics under ultrasonic vibrations and optimization of the vibration parameters. The studies conducted by Wiercigroch et al 8 showed that the high-frequency axial vibration significantly enhanced the breaking rates compared to the traditional method. Subsequently, uniaxial compressive strength experiments were conducted on fine granite by Yin et al 9 to demonstrate that there exists a static force threshold for this technology. They found that the optimal value range of static loading was 200-300 N. An 3806 | TANG eT Al. ultrasonic vibration cutting process was investigated by Zhao et al 10 who found that the friction between the rock and the bit would be reduced under ultrasonic vibration loading. Based on the finite element method, Zhou et al 11 concluded that for the ultrasonic vibration technology, combining the ultrasonic dynamic loading and static loading, the static force applied on the top of the rock sample would increase its natural frequency. The above studies showed that compared with the traditional methods, assisted ultrasonic vibration technology appears to have unexpected effects on rock breaking. A thorough understanding of the fragmentation mechanism of rock subjected to ultrasonic vibration can have practical implications for its further engineering applications. However, there has been very little research in the investigation of fragmentation mechanism of rocks subjected to ultrasonic vibration.
Incidentally, in recent years, much work has been conducted on the damage process of rocks under cyclic loading. For example, Chen et al 12 investigated the characteristics of microcracks in granite during each stage of cyclic loading using a fluorescent method. Based on cyclic true triaxial loading experiments, Gao et al 13 found that the irreversible strains in σ1 and σ3 (principal stresses) directions changed linearly with cumulative damage, whereas the irreversible strains in σ 2 direction decreased nonlinearly. Yang et al 14 studied the influence of cracks on the mechanical properties of rocks under cyclic loading and found that the strength decreased with a decrease in the crack dip angle and an increase in the crack density. Liu et al 15 conducted constant-amplitude cyclic loading experiments and found that fractures preferentially propagated along the end or edge of the crack with strong stress sensitivity to form a shear failure plane at a high rate. Wang et al 16 found that the axial fatigue deformation could be divided into three stages, the initial stage, constant velocity stage, and accelerating stage. The irreversible deformation could be used to describe the damage accumulation process. Peng et al 17 obtained the energy evolution law during the cyclic loading process. They concluded that the larger the crack angle, the larger the energy storage coefficient of the rocks, and the smaller the energy dissipation coefficient. These works provide great help in exploring the fragmentation mechanism of rocks under ultrasonic vibration, although different loading frequencies may lead to different rock mechanical behaviors. 18,19 As the above studies were all related to low-frequency cyclic loading conditions, it is necessary to investigate the damage process and the development of cracks to analyze the failure mechanism of rocks under ultrasonic vibration loading. However, research on rock fragmentation mechanism has been limited because of the difficulty in measuring microcracks in real time during high-frequency vibration laboratory tests. 20 With the rapid development of computer technology and advanced numerical techniques, computational simulations help us overcome this limitation and determine the characterization of crack initiation and growth in more detail. [21][22][23][24][25][26][27] In recent years, the DEM has received considerable attention for its use in modeling the dynamic process of rock failure owing to its three advantages. [28][29][30][31] Firstly, considering its micromechanical foundation, the DEM is a useful simulation tool for understanding the dynamic mechanism of rock materials at a microscale directly. Secondly, real-time contact searching makes the DEM convenient for observing the dynamic process. Thirdly, in the DEM, the bond breakage between the grain elements could be useful in effectively simulating the dynamic crushing failure of the rock. By introducing a damage variable into the DEM model, Li et al 32 studied the failure mechanism of a rock containing two parallel flaws under uniaxial compression and found that the type of macrocracks did not depend on the dominant type of microcracks. Nhu et al 33 successfully reproduced the heterogeneous microstructure of cementitious materials and the actual crack development process by establishing a DEM model of the fatigue damage. Based on a DEM simulation, Ma et al 34 verified that shear failure was the primary failure pattern of the surrounding rock mass during tunneling. Thus, in our work, DEM simulations were used to investigate the fragmentation mechanism of rocks under ultrasonic vibration loading.
The following sections present the numerical simulations performed using particle flow code with two dimensions (PFC2D), which is a DEM-based simulation software, to investigate the fragmentation mechanism of hard rock under ultrasonic vibration and experimentally verify the applicability of the DEM simulation to high-frequency loading conditions. The remainder of this paper is organized as follows. The detailed features of the established rock model and methodology to simulate the ultrasonic vibration loading are described in Section 2. The applicability of the model is verified by comparing the numerical and experimental results. Section 3 describes the development of the strain field and crack evolution law and presents an analysis of the influence of the loading parameters. Section 4 presents a discussion of the fragmentation mechanism of rock under ultrasonic vibration according to the obtained results.

SIMULATION
In numerical simulations, compared to a 3D software, a 2D software has the advantage of providing insights into critical phenomena and mechanisms with less computational cost. Thus, 2D software programs have been widely used to investigate geomechanical problems. 35 In the present work, the PFC2D software was used to investigate the fragmentation mechanism of a hard rock under ultra-high-frequency loading.
PFC2D, which is based on the discrete element method, can be used to study the propagation of cracks in brittle rocks and rock-like materials (eg, granite and concrete). 33,36 In PFC2D, the material is modeled as a collection of particles (rigid circular disks). Each particle is in contact with the neighboring particles via contact elements. The parallel-bond model (PBM) has been widely used to simulate crack initiation and propagation of rocks or rock-like materials. 37,38 However, there are two main problems with regard to this model when it is used for simulating hard rock: (a) It yields an unrealistically low unconfined compressive strength-totensile strength ratio, and (b) it generates an excessively low internal friction angle. To solve these problems, Potyondy 39 proposed a new bond model called the FJM, which consists of common particle elements and flat-joint contacts between them. Wu et al 40

| Fundamentals of FJM
As shown in Figure 1, there are two notional surfaces, between which there is an interface divided into several parts in the FJM. The breakage of each of these parts contributes to the partial damage to the interface. The flat-joint contacts are rigidly connected to the grains and make the grain appear as a skirted particle. Thus, it could provide grain interlocking and rotational resistance even after the interface breaking (unbonded state), and the sample could be characterized by a realistic unconfined compressive strength-to-tensile strength ratio and internal friction angle.
After flat-jointed contacts are installed between the particles, forces are generated following the detection of a crack generation. The normal force is computed in a direct fashion and the shear force in an incremental fashion. 42 The normal stress (e) and shear stress (e) acting on an element are given by where F n e and F s e are the normal and shear forces acting on the element, respectively, and A (e) is the area of the element.
The normal interface stress is given by where k n is the normal stiffness, and g is the interface gap. The normal force-update law for a flat-joint contact is given by where is the normal interface stress in Equation (2), and the integration is performed over the element. An analytical expression for this integral is used to update F n e , which is then used to update (e) in Equation (1).
The tensile strength is specified as b . If the normal stress (e) > b , then the element breaks when a tensile crack is generated, and the element bond state is modified to an unbonded state.
The shear force-update law for a flat-joint contact is given by where k s is the tangential stiffness, and s is the tangential displacement. The shear force is then used to update (e) in Equation (1).
(2) (r) = 0, unbonded and g(r) ≥ 0 k n g(r), otherwise , The procedure considering the shear strength based on the bond state is as follows:

Unbonded state
The shear strength obeys the Coulomb sliding criterion, given by where r is the residual friction strength, and r is the residual friction angle. If (e) ≤ r , the shear stress remains as (e) ; otherwise, the shear-strength limit is exceeded, which implies that slip would occur.

Bonded state
The shear strength follows the Coulomb criterion with a tension cutoff, given by where c b is the bond cohesion, and b is the local friction angle. If (e) ≤ c , the shear stress remains as (e) ; otherwise, the shearstrength limit is exceeded. If the bond breaks in shear, and a shear crack is generated, the residual friction takes effect.

| Particle modeling and calibration
As is well known, there are many contact microparameters that need to be calibrated from the material's macroproperties before further simulations. The microparameters include the effective moduli, particle normal/shear stiffness ratio, flat-joint tensile strength, and cohesion. The corresponding macroproperties include Young's modulus, Poisson's ratio, UCS, and tensile strength. In addition to the microparameters, the microstructure (eg, slit contact, which is used to model the initial crack) also influences the macroproperties of the established model, such as the initial damage of the rock sample. 43 The fraction of slit contacts has been set to 10% in this study based on a previous study. 40 The properties of the flat-joint contact model and a full description of the microparameters can be found in Ref. 39,40 To comprehensively understand the influence of the microparameters on the numerical results and calibrate them effectively, numerous studies on the relationship between the microparameters and properties have been carried out, 40,41,44,45 and the following conclusions have been arrived at: The stiffness ratio is the most influential parameter. The macro-Poisson's ratio increases as the stiffness ratio increases. The other aforementioned macroproperties decrease as the stiffness ratio increases. The effective moduli appear to only influence the macro-Young's modulus and not any other macroproperties. The contact tensile strength and cohesion have a direct impact on the macro-tensile strength and UCS of the rock, respectively. The contact tensile strength also has a slight influence on the macro-Young's modulus, while the contact cohesion has no impact on the macro-tensile strength. These aforementioned macroproperties could be obtained from the UCS and Brazilian tests. In this study, the microparameters were calibrated using the procedure proposed by Li et al 46 based on these two test simulations with a sample with a height and diameter of 100 and 50 mm, respectively. In the UCS test simulations, the loading velocity at the upper and bottom walls was 0.1 m/s. On the other hand, for the Brazilian test simulations, it was 0.0075 m/s, under which the specimen could be regarded to be in a quasistatic state. 40 According to Li et al, 46 the ratio between the sample diameter and maximum particle diameter should be greater than or equal to 20, and the larger the ratio, the more stable the result. Therefore, the number of particles was determined based on this standard and model run time. Finally, the particle sizes were controlled by following a uniform distribution, with maximum and minimum radii of 0.24 and 0.32 mm, respectively. The detailed calibration procedure is as follows: The microparameters used in this study are listed in Table 1. Comparison of the macroproperties and failure patterns between the experimental and numerical results is provided in Table 2 and Figure 2. It can be found that the model based on the above microparameters could reflect the macromechanical behavior of the rock. Hence, this model could be used in further simulations.

| Crack mechanisms in the three mechanical simulation tests
It is widely accepted that tensile cracking is the primary damage mechanism of brittle rock in compression tests, such as the Brazilian test and the UCS test. Most of the tensile cracks formed during these two tests are caused by compression. [47][48][49][50][51][52][53] On the other hand, in the direct tension test, most of the tensile cracks are induced by tension. 54 These are the two main failure modes for rock materials. Exploring one method which could reflect the cracking mechanism based on PFC2D would be useful to expand the usage of this kind of numerical method.
This section presents the investigations on the crack orientation of rock samples from the UCS test, Brazilian test, and direct tension test, and the PFC2D simulations. The UCS and Brazilian test simulation are the same as presented in Section 2.2. In the direct tension test simulations, the tension force was applied to the two ends of the rock sample and the loading rate was 0.0075 m/s. 40 Figure 3 shows the tensile crack orientations from these three tests. In the UCS and Brazilian tests, although the number of cracks varied greatly, the orientations of most tensile cracks were between 60 and 120°; these were parallel or subparallel to the loading direction. In the direct tension test, the crack orientation was quite different from that of the other two tests; here, the mean orientations of the tensile cracks were in the ranges of 0-30° and 150-180°, which were perpendicular or subperpendicular to the loading direction. Based on these results, it can be concluded that the crack orientation obtained by PFC2D adequately reflected the mechanism of crack evolution.

| Modeling the rock breaking by ultrasonic vibration
The full-scale ultrasonic vibration test system used in this study is shown in Figure 4 (1-ultrasonic power, 2-ultrasonic transducer, 3-static loading device, 4-dynamic strain data collection device, 5-computer, 6-strain gauge, and 7-rock sample with a diameter and height of 36 and 72 mm, respectively). An ultrasonic dynamic load combined with a vertical static load was applied to the top surface of the samples. To monitor the strain in real time, a dynamic strain data acquisition device was used, and the strain gauge was attached to the top surface of the sample. During the tests, the static loading F staic and loading frequency f were chosen as 200 N and 30 000 Hz, respectively. The experiments were conducted at ambient temperature.
To simulate the experimental conditions, that is, to model the breaking of the rock by ultrasonic vibration, an intact numerical specimen with a width and height of 36 and 72 mm, respectively (ie, with same size as that of the rock used in the test). was constructed first ( Figure 5A). The microparameters were ensured to be consistent with those shown in Table 1. Each numerically intact specimen was discretized into 10 118 particles with 27 022 contacts. The bonded and slit contact distributions are shown in Figure 5A, and the slit contacts were used to simulate the initial cracks.
The role of static loading in the ultrasonic vibration system is to enlarge the loading amplitude, 55 which means that only the dynamic loading should be simulated. It can be verified from related studies that the velocity boundary condition could closely replicate the dynamic behavior of ultrasonic vibration. [56][57][58] Thus, in the present work, the wall element was used, which could be applied with velocity boundary conditions to simulate the ultrasonic vibration loading. The velocity boundary condition was set according to the procedure described below.
As can be found from Zhou, 59 the displacement amplitude of the test part obtained from the strain gauge of the rock was approximately 10 μm under ultrasonic vibration with the selected loading parameters (ie, f = 30 000 Hz and F static = 200 N). Thus, the motion equation of the wall element is as follows.
where y is the displacement in the y-direction (m); f is the loading frequency (Hz); and t is the simulation time (s).
The velocity boundary condition applied to the wall element was where v y is the velocity of the wall in the y-direction (m∕s).
If the simulation boundary was set the same as that in the experiments shown above, the computational time was approximately 1 month per sample with an Intel ® Xeon ® CPU E3-1231 v3 3.40 GHz, 16 GB RAM processor. To save the calculation time, an "amplitude enlargement" assumption was employed to reduce the computational cost. For this assumption, it was based on the relationship between the ratio of the maximum loading stress and UCS of rock material obtained from Cerfontaine's research. 60 A linear regression is adopted to describe S-N curve in log 10 basis. The S-N rule is shown in Equation (10), with the ratio increases, the failure cycle number will decrease. Thus, in the present work, the loading amplitude of the wall element was enlarged with an upscaling factor f E ; in this case, we set it as 5.3. Although the loading amplitude is still relatively low, it would greatly decrease the simulation time. After determining the enlargement factor of the loading amplitude, the factor between the simulation time and real time should be obtained. The stress curve before and after the increase in loading amplitude is shown in Figure 6. It can be found from the Figure 6, before increasing the loading amplitude, the ratio between the maximum loading stress and UCS is 0.0434. The ratio increases to 0.1458 after enlargement. Then, these two ratios are substituted into Equation (10). The decreasing factor for the loading cycle number is 10 3.089 .
The enlarged displacement and velocity boundary are shown in Equation (11).
where y′ and v ′ y are the enlarged displacement and velocity of the wall in the y-direction (m∕s).
In addition to the crack evolution process, 13 × 25 measurement circles were installed in this model to record the strain development, as shown in Figure 5B (only some of the measurement circles are labeled in this figure).

RESULTS
The results from the numerical simulations and the physical experiments for the hard rock subjected to ultrasonic vibration are presented in the following subsections.  with that obtained by Wang et al 61 The numerically simulated strain curve agreed well with the experimental curve, except for the curve smoothness. The experimental strain curve rose and fell nonlinearly. In addition to the smoothness, compared with the experimental data, the strain curve in the simulation decreased more slowly in the initial stage. These differences were possibly owing to the difference between the FJM and the actual rock. 62 The initial cracks in the actual rock were open, which could influence the system response and deformation behavior of the rock samples under ultrasonic vibration. 23,63 It was difficult to faithfully replicate the features of the cracks in the FJM. It would be interesting to simulate a more realistic behavior of cracks, and this would constitute our future work. The fragmenting process in the numerical modeling is shown in Figure 8. The irregular blocks with different colors, except those in azure, represent fragments, and the azure ones are the matrix of the sample. The fragmentation process predicted by the modeling was similar to that observed in the experiments (Figure 9), with fragments first occurring at the corner and then propagating to the center. In this process, with an increase in the number of cycles, the rate of fragment generation increased. This result was consistent with that of the strength tests by Yin et al. 9 The size distributions of the fragments in the simulations and laboratory tests are shown in Figures 10 and 11, respectively. The green line in Figure 10B represents a fractured contact, and the black line represents a bonded contact. The fragments predicted by the modeling agree well with those generated in the test, with several large fragments produced along with smaller fragments and some fines. Figure 10C illustrates the development of the fragment size distributions F I G U R E 6 Stress curves before and after the increase in the loading amplitude with the factor of 5.3

F I G U R E 7 Comparison between experimental and numerical strain curve
with different vibration durations in this simulation. It shows that in the initial stage (cycle num = 5 000 000), most fragments had a size of 0.04 m, and more than 50% of the total fragments had sizes between 1.33 and 1.51 mm. After 10 000 000 cycles, the fragments became relatively uniform; in addition, some large fragments were generated. After this point, the proportion of small fragments with a size of approximately 0.86 mm increased dramatically. This was because the vibration would cause the recrushing of large fragments. With further vibration, this phenomenon became more apparent.
In general, according to the above numerical and experimental results, it can be concluded that PFC2D could faithfully replicate the dynamic process of the hard rock specimens under ultrasonic vibration loading condition.

| Crack evolution under ultrasonic vibration
In the simulation process, most of the generated cracks were tensile cracks. Figure 12 shows the results for acoustic emission (AE) count and damage variable from the ultrasonic vibration simulations. This figure clearly shows that the ringing response changed from sparse to dense with an increase in the number of cycles and that AE events occurred throughout the stage. In other words, as the number of cycles increased, the rate of microcrack generation gradually accelerated, until failure occurred eventually. At the beginning of the simulation, the AE signal was weak, and only a few cracks were generated. As the loading continued, the signal became distinct, and some exploding points occurred, which represented the generation of macroscopic cracks and fragments. During the last stage, the number of AE events further increased when the released energy was significant, indicating that the cracks in the sample experienced an unsteady propagation stage.
Next, the distribution of crack orientation was analyzed to understand the fragmentation mechanism, and the results are shown in Figure 13. The angle shown in the rose diagram is the angle between the crack plane and the horizontal axis, and the value is the number of cracks. Based on the conclusion given in Section 2.3, the tensile cracks, whose orientations were in the range of 0-30° or 150-180°, were caused by tension, while those in the range of 60-120° were caused by compression. Figure 13 shows that in the initial stage, microcracks were mainly compression-induced tensile cracks with a mean orientation of 90°. As the number of cycles increased, tension-induced tensile cracks with orientations perpendicular to the loading direction were gradually generated. Both types of cracks continued to grow with the calculating cycles; however, the latter grew faster, and eventually, they were close in number.
A quantitative analysis of these two types of cracks is shown in Figure 14. It can be seen that the total number of these two types of cracks accounted for more than 93% of all the cracks. During the whole simulation, the growth rate of compression-induced tensile cracks decreased proportionately from 75% to 50%. The tension-induced tensile cracks grew at a faster rate, resulting in a 50% increase from their original proportion of 25%. Finally, the number of these two types of cracks tended to be the same. Figure 15 intuitively shows the spatial distribution of these two types of cracks over time. The red points represent the tension-induced tensile cracks, while the green points represent compression-induced tensile cracks. In the early F I G U R E 8 Fragment process of rock in simulation stages, more compression-induced cracks occurred at the edge of the sample. Subsequently, a large number of tension-induced cracks were generated, and these extended toward the center. The compression-induced cracks were more likely to be distributed around the periphery of the sample, while tension-induced cracks were distributed in the center.

| Distribution of the time-dependent strain field
The Equation (12) is the interaction law in the particle system which is one of the basic theories of PFC2D. 64 (12)  where p =u∕L p and L = 2R, in which R is an intrinsic particle length scale; u represents the overlap/separation/sliding distance of the two contacting/neighboring particles; L p is the representative length of the particle; p represents the particle strain; and k is the stiffness. Thus, not only did a stress limitation exist between the two bonded particles, but also a strain limitation. On the other hand, Wang et al 61 concluded that the axial strain could be a better option to describe the mechanical behavior of rock under cyclic loading conditions. Thus, the axial strain was analyzed in this section.
The strain field in the y-direction (y-strain field), based on the PFC2D measurement cycles, is shown in Figures 16 and  17. It can be seen that the strain field was not uniform. There were two different kinds of strain zones in the upper part of the sample: One was the compressive strain zone appearing in the upper left corner, and the other was the tensile strain zone distributed in the center. During cycles 5 000 000-20 000 000, the compressive strain first increased and then decreased, which was caused by the damage of the rock at the upper left corner of the sample. The tensile strain in the original tensile zone gradually increased, and a new and distinct tensile strain zone appeared at the lower right corner of the sample. Combined with the crack distribution, it can be observed that the development of the strain field corresponded to the spatial distribution of the cracks. Tension-induced tensile cracks mainly occurred in the tensile strain region and were mostly distributed in the center of the sample. In contrast, the compression-induced tensile cracks were mainly generated in the compressive strain region and were mostly distributed in the periphery.

| Influence of loading parameters on fragmentation efficiency
This section presents an analysis of three loading parameter combinations. The chosen parameters were determined by the design principles of the ultrasonic transducer (ie, under the same power, the higher the vibration frequency, the lower the amplitude 65 ). The loading parameter combinations are listed in Table 3. The vibration duration was the same as that described in Section 2.4.  Figure 18A illustrates the crack evolution of the rock samples as a function of the three vibration loading parameters. Specifically, as the vibration frequency increased from 20 to 40 Hz, the total crack number first increased to 3430 and then sharply decreased to 780. Such a phenomenon could be explained by considering the stress level of the rock sample under vibration and the duration of the vibration loadings. Compared to the loading frequency of 20 kHz, the rock sample was closer to the resonance condition with a frequency of 30 kHz. Thus, the stress level was higher for the latter case even though its amplitude was lower than that of the former ( Figure 18B,C). To minimize the influence of cracks on the system response, these contact force data points were obtained in the same crack condition as the red frame shown in Figure 17A. For Case 3, on the one hand, as the frequency was higher than 30 kHz, the rock sample was also far from the resonance condition. On the other hand, the amplitude in Case 3 was the lowest. Thus, its stress level was the lowest ( Figure 18D). Furthermore, in Case 3, the load varied more rapidly, and thus, the duration of the load acting on the tested specimens was shorter. The microcracks had insufficient time to develop and coalesce, and more cycles were required to damage the specimen. 66 Thus, the key to use this technology for a high efficiency was selecting a vibration frequency that is close to the resonance frequency of the rock sample and improving the vibration amplitude. Figure 18E,F depicts the failure scenarios of the tested specimens under different vibration loadings. It is evident that all the failed specimens exhibited similar failure modes, namely the splitting failure mode. It can thus be concluded that the failure mode of the hard rock under ultrasonic vibration was independent of the loading parameters (ie, loading frequency and amplitude).

FRAGMENTATION MECHANISM
Based on the numerical results, there were mainly two types of cracks generated under ultrasonic vibration: compression-induced F I G U R E 1 5 Time and spatial distribution of the two types of tensile cracks F I G U R E 1 6 Strain field in the ydirection: A, Cycle num = 5 000 000; B, Cycle num = 100 000 (the positive value represents the tensile strain, and the negative value represents the compressive strain) cracks and tension-induced cracks. These two types of cracks were caused by the inhomogeneous deformation of the rock samples. The mechanism for the formation of compression-induced cracks was determined, as shown in Figure 19(A). The translucent circles with dotted lines represent the original location of the particles, and the opaque circles represent the locations of the particles after the deformation. Particles 2 and 3 were forced apart by the radial load owing to the displacement difference between particles 1 and 4, causing the restraining bond between particles 2 and 3 to produce a tensile deformation.
Notably, the tensile strength of the rock was considerably less than its compressive strength and the tensile limitation strain of the rock was accordingly less than its compressive limitation strain. Thus, the subaxially aligned tensile microcrack was formed when the tensile strain was higher than the tensile limitation strain of the bonded assemblies of the circular or spherical particles, and it was formed from the compression-induced tensile cracks.
The mechanism for the formation of tension-induced tensile cracks is depicted in Figure 19(B). In the tensile strain zone, the particles were stretched in the y-direction, which was parallel to the loading direction. Thus, particle 1 was apart from particles 2 and 3. Accordingly, a relative displacement occurred between particle 1 and particles 2 and 3, which induced a tensile strain inside the contact between these particles. With a gradual increase in the relative displacement and tensile strain, subradial cracks were eventually generated, when the tensile strain was higher than the tensile limitation strain of the bonded contact.

| CONCLUSIONS
To determine the fragmentation mechanism of hard rock under ultrasonic vibration, numerical simulations based on PFC2D were performed with a flat-joint contact model. In order to accelerate the simulation, the "amplitude enlargement" assumption has been employed. The numerical simulation results for the rock specimens subjected to ultrasonic vibration showed good agreement with the experimental results in terms of the damage velocity, fragment size distribution, and local strain development. Upon comparing three loading parameter combinations, it was observed that the most effective was the case with a frequency of 30 kHz and an amplitude of 10 μm. The factors responsible for this result included a high-stress level owing to the resonance condition, a high amplitude, and a relatively long loading duration. Thus, for engineering practice, the key for using this technology for a high efficiency is selecting a vibration frequency that is close to the resonance frequency of the rock sample and improving the vibration amplitude.
By analyzing the development of the full-field strain, crack orientation, and crack distribution, the fragmentation mechanism in the hard rock under ultrasonic vibration was revealed. Cracking is mainly caused by uneven deformation of the sample resulting from the ultra-high-frequency loading. There are two different strain zones inside the sample: the compressive strain zone and the tensile strain zone. In the compressive strain zone, most cracks are in the orientations of 0-30° and 150-180°, which is caused by compressive deformation. On the other hand, in the tensile strain zone, the cracks were primarily orientated in a range of 60-120° and are mainly caused by tensile deformation. Eventually, the proportions of these two cracks in the specimen become almost equal. The damage gradually accumulates and causes the final failure of the sample. In summary, these numerical results are expected to improve the understanding of the fragmentation mechanism of hard rock under ultrasonic vibration and broaden the application scope of PFC2D and FJM to ultra-high-frequency loading conditions.