# Rock Crushing Analysis of TBM Disc Cutter Assisted by Ultra-High-Frequency Loading

Authors: Dajun Zhao; Junpeng Han; Yu Zhou; Xiaoshu Lv; Shulei Zhang; Ke Gao; Yan Zhao
Journal: Shock and Vibration
 (2022)

---

## Abstract

In the construction of mountain tunnels, the increasing proportion of hard rock stratum enhances the difficulty of rock breaking by tunnel boring machine (TBM). As a type of cyclic load with ultrahigh frequency, ultrasonic vibration has the advantages of energy concentration and strong penetration. In the process of rolling and extrusion, the weakening of the rock by ultrahigh-frequency loading can improve the rock crushing efficiency of the TBM disc cutter. In this study, we established a physical model of rock cutting using a disc cutter assisted by ultrasonic vibration and obtained a motion equation. The discrete element software particle flow code (PFC) was selected to construct a heterogeneous granite model to demonstrate the mechanism and development of cracks inside the hard rock under a rolling disc cutter assisted by ultrasonic vibration. The results demonstrate that ultrasonic vibration helps the disc cutter construct a stronger stress field in the shallow layer, which promotes tensile damage of the surface rock. Vibration promotes the development of cracks inside the rock and accelerates the penetration of transverse cracks, which is conducive to the stripping of rock fragments. The addition of ultrahigh-frequency loading also reduces the fluctuation increase of cracks, which makes the operation state of the disc cutter more stable and avoids abnormal damage to tools.


---

## Body

## 1. Introduction

In recent years, the continuous in-depth development of underground spaces has gradually increased the amount of tunnel engineering, and the impact of various complex strata on construction has gradually become prominent. In the process of building tunnels in mountainous areas, the high-strength hard rock formation is a common and essential factor leading to slow construction progress, and efficient hard rock fragmentation has become a common concern [1, 2]. The tunnel boring machine (TBM) is commonly employed for rock excavation in tunneling construction. The working principle is that the disc cutter rolls and crushes the rock on the surface of the rock mass under the action of thrust and torque exerted on the cutter head. The main rock-breaking process is expressed as follows: after the tool is pressed into the rock, microcracks are produced inside the rock. The tool was gradually wedged into the rock with an increase in pressure, which caused cracks to scatter and expand. This approach eventually causes cracks between two adjacent tools to penetrate each other, and the resulting fragments and blocks of rock tend to fall off [3]. The pressure required to destroy the hard rock stratum is extremely high because these formations usually have high strength, and the rock can be very abrasive when the quartz content of the formation is high [4]. When the tool works under the conditions of strong extrusion and high grinding, its failure rate will be greatly accelerated. Therefore, the speed of tunneling under hard rock conditions is slower, and the tool wears faster, which greatly increases the construction cost [5].To improve the efficiency of TBM excavation and to reduce construction costs, many scholars have conducted studies on improving the efficiency and service life of disc cutters. Some scholars have investigated the tool wear mechanism and analyzed the factors affecting tool efficacy to extend the tool life per unit feed length [6]. Wang et al. [7–9] proposed the use of a high-pressure water jet to assist the rock breaking of the disc cutter, which reduces the stress of the disc cutter to reduce its abrasion. Liu and Wang [10] discovered that under the same conditions, two disc cutters sequentially crushing the rock form a larger breaking groove than double-disc cutters that are simultaneously working. Numerical analysis shows that the crack penetration effect is better under sequential breaking conditions. Tan et al. [11] compared rock-breaking using a disc cutter under impact dynamic load and static load conditions and discovered that cyclic impact loading can promote crack generation inside the rock and its extension in the horizontal direction. Zhao et al. [12] investigated the intrusion process of an indenter into a granite rock and concluded that coupled dynamic and static loading can significantly improve the rock-breaking effect.The vibration of the load may promote the growth of primary cracks and create new cracks, which is beneficial for the destruction of the rock structure [13]. Bagde and Petroš [14–16] tested the fatigue properties of rocks under dynamic loading and determined that the frequency and amplitude of the load affect the dynamic strength, Young's modulus, and dynamic axial stiffness of rocks. Ultrasonic vibration is a type of cyclic loading action at ultrahigh frequencies that is capable of generating the same or similar frequencies as the inherent frequency of hard rock [17]. Under the influence of certain high-frequency cyclic loads, the rock resonates, with fatigue damage occurring inside it [18]. Zhao et al. [19] used nondestructive thermal imaging technology to monitor the temperature climbing phenomenon of rocks before damage under ultrasonic vibration loads and concluded that the causal factors of rock damage are fatigue damage and thermal damage produced by alternating loads on rocks. The authors divided the damage process into three stages: elastic deformation, microcracking and yielding, and macroscopic cracking and failure [20]. Yin et al. [21] verified through uniaxial compression experiments that ultrasonic vibration effectively reduced the strength of granite. Fernando et al. [22] selected three types of rocks for rotary ultrasonic processing experiments and determined that ultrasonic vibration can reduce the rock cutting force and increase the cutting speed by three times.The detection of microscopic damage is essential in the study of rock damage processes; therefore, it is necessary to adopt a means to record the behavior of rock crack extension development. Numerical simulation methods are widely utilized in theoretical studies of rocks because they can demonstrate microscopic damage and macroscopic mechanical parameter changes during the operation of the rock model using computer software, which is more intuitive and easier to analyze than physical tests. The discrete element simulation method has attracted much attention because of its significant advantages in simulating crack evolution in rocks and the heterogeneous properties of rocks [23, 24]. Moon and Oh [25] analyzed the rock force state corresponding to the disc cutter spacing and penetration based on the discrete element numerical simulation method to obtain the optimal rock-breaking conditions. Jiang et al. [26] simulated the process of cutting rock using double-disc cutters and divided the rock fragmentation into three stages: loading, unloading, and residual leap fragmentation. Li et al. [27] established a wedge-shaped disc cutter crushing model with a new particle swarm component to depict the microcrack propagation and damage evolution and to analyze the crushing efficiency with different parameters.The existing impact dynamic load frequency applied to the disc cutter to improve the crushing efficiency is low, and a single impact is accompanied by a large load change, which is not conducive to the maintenance of the disc cutter life. There is still a lack of research on the application of ultrahigh-frequency loading to assist in disc cutter loading. Discrete element simulation technology can be employed to probe the internal damage evolution process of rocks from a microscopic perspective, which is important for revealing the mechanism of ultrasonic vibration-assisted rock breaking. In this study, a heterogeneous granite model was constructed using Particle Flow Code software. The development of rock crack extension under the normal rolling mode and ultrasonic vibration-assisted mode is analyzed by simulating the process of the disc cutter into granite rock. This research compares the rock-breaking effect of a disc cutter with the assistance of ultrasonic vibration, which provides a new idea for improving the efficiency of tunneling. Research on rock crushing of the TBM disc cutter assisted by ultrasonic vibration has also provided a theoretical basis.


## 2. Discrete Element Model and Parameter Calibration

The particle flow code is a discrete element model (DEM) based fine-scale analysis software that has a wide range of applications for simulating dynamic damage processes in heterogeneous rock materials.


### 2.1. Theory of the Discrete Element Method

#### 2.1.1. Fundamental Principle

Each particle in the model is a rigid body with a mass and surface that independently move and can generate displacement and rotation. The contact generated by the pairing of internal forces and torques defines the mode of interaction between two particles and is updated in real time with the identification and generation of new contacts during the calculation. The contact between two particles is achieved as soft contact, where rigid particles can overlap at the contact point, with a small overlap area, usually considered as a point, and the magnitude of the overlap is related to the contact force through the force-displacement law. The motion of the particles conforms to Newton's law of motion, and the motion of the wall is specified by the user and is independent of the contact force required to achieve the loading or restraint of the particles.The computation process takes the form of an iterative loop, in which the position and velocity are updated based on Newton's law of motion in a determined time step. Then, the model time is updated, the contact is created or deleted according to the particle position, and the forces and moments at the contact are updated according to the force-displacement law to start the next computation loop.


#### 2.1.2. Contact Model

The defined law of the interparticle interaction is a contact model. The mechanical behavior of the contact determines the accuracy of the discrete element calculation results. The micromechanical behavior between individual particles and walls is generated through contact and organized in series with each other, which in turn affects the whole by the local. For the different contact types, there were significant differences in the calculation of the contact forces. The research team selected the appropriate contact models to set the contact properties between particles and walls according to the characteristics of each model, and the types of contact models included the linear model, linear parallel bond model, and smooth-joint model. The operating rules for the three contact models in the 2D model are described.The linear model provides the force transfer behavior of an infinitesimal interface. The contact forceFcL includes linear component Fl and damping component Fd. Linear force Fl controls the linear elastic friction, whereas damping force Fd controls the bonding. The contact torque McL remains zero, allowing the particles to rotate relative to each other.The update of the forces and torque at the contact in the linear model satisfies:(1)FcL=Fl+Fd,McL=0.The forces decomposed in the tangential and normal directions:(2)Fl=−Fnln^c+Fsl,Fd=−Fndn^c+Fsd,where n^c is the contact unit-normal vector, which is directed from the center of piece 1 to the center of piece 2. For the 2D model, only when the surface gap gs<0 can the linear normal force be updated; otherwise, the force remains zero:(3)Fnl=kngs,where kn is the normal stiffness. To update the shear force, it is necessary to determine whether it is in a sliding state. The shear strength of the contact is calculated by the normal linear force:(4)Fsμ=−μFnl.The linear shear force value that should be obtained at the current moment is calculated through the initial shear force:(5)Fs∗=Fsl0−ksΔδs,where μ denotes the friction coefficient, ks is the shear stiffness, and ∆δs is the relative shear‒displacement increment.The final result of linear shear force:(6)Fsl=Fs∗,Fs∗≤Fsμ,FsμFs∗/Fs∗,Fs∗>Fsμ.In the full normal and full shear states, the update of the damping force is(7)Fnd=2βnmcknδ˙n,Fsd=2βsmcksδ˙s,where(8)mc=m1m2m1+m2,where βn is the normal critical damping, βs is the tangential critical damping, ks is the tangential stiffness, δ˙n is the normal relative displacement velocity, δ˙s is the tangential relative displacement velocity, and m1 and m2 are the masses of the particles at both ends of contact. In the slipping state, Fsd is zero.The linear parallel bond model can be used to describe the mechanical behavior of the bonded material. The bond is considered a set of springs with constant normal and shear stiffnesses, which are uniformly distributed in a rectangular plane within a certain range of contact between two particles. The model produces two bonding components: linear and parallel bonding, which is equivalent to the linear model. Parallel bonding works in parallel with the former. The linear components do not resist relative rotation, and parallel bonding generates forces (Fp) and torque (Mp). Thus, the component resists relative rotation. The contact force and torque between two particles using linear parallel bonding are(9)FcLPB=Fl+Fd+Fp,McLPB=Mp.In the bonded state or when the surface gap does not exceed zero, the model will update the contact through the force-displacement law.(10)Fp=−Fnpn^c+Fsp,Fnp=Fnp0+knSΔδn,Fsp=Fsp0−ksSΔδs,where ∆δn is the relative normal displacement increment, ∆δs is the relative shear displacement increment, and S is the cross-sectional area of contact. In the 2D model, the bond cross-section is rectangular.(11)Mp=Mp0−knIΔθb,where ∆θb is the relative bend-rotation increment, and I is the moment of inertia of the parallel bond cross-section.The smooth-joint model can simulate the shear-expansion mechanical behavior of a planar interface without considering the contact direction of the local particles distributed along the interface. Before the bond exceeds the strength limit or breaks, the bonded part has a linear elastic property, and the unbonded part is linearly elastic with expansive friction, which adapts to slip by imposing a Coulomb limit on the shear force. A frictional or bonded connection can be simulated by setting two sets of particles in contact with each other in this model. The contact force(12)FcSJ=Fs,McSJ=0,Fs=−Fnsn^j+Fss,where n^j is the unit-normal vector, which points into surface 2.In the bonded state,(13)Fns=Fns0+knSΔδne⌢,Fs∗=Fss0−ksSΔδse⌢,where Fns0 and Fss0 are the smooth-joint normal force and shear force at the beginning of the timestep, respectively, and ∆δne^ is the elastic portions of the normal displacement increments. ∆δse^ is the elastic portions of the shear displacement increments.The shear strength can be computed through the normal force:(14)Fsμ=−μFns.The shear force at this time is determined based on the shear strength:(15)Fss=Fs∗,Fs∗≤Fsμ,FsμFs∗/Fs∗,Fs∗>Fsμ.In the sliding state(Fs∗=Fnμ), the slip-change callback event will cause the normal force to increase with the expansion of the shear displacement:(16)Fns=Fns+Fs∗−Fsμkskntanψ,where ψ is the dilation angle.Figure1 shows the theoretical surface and particle states for the three contact types.Figure 1
Contact theory surface and particle state: (a) linear mode, (b) linear parallel bond model, and (c) smooth-joint model.
(a)
(b)
(c)


### 2.2. Microparameter Setting and Calibration

To make the physical and mechanical properties of the model more similar to the actual granite and to achieve a reliable simulation effect, the research team conducted indoor uniaxial compression tests on granite to obtain the macroscopic mechanical parameters. A discrete element model of granite was also constructed for the uniaxial compressive numerical simulation. Based on the indoor test and simulation results, the microparameters of the model were calibrated using the trial-and-error method.


#### 2.2.1. Rock Sample

The selected granite sample was medium-grained granite quarried from a mine in Jilin Province, China, and the main mineral composition of the sample included feldspar, quartz, and mica. The physical and mechanical parameters of granite were experimentally measured: the density is 2.79 g/cm3, the elastic modulus is 29.8 GPa, Poisson's ratio is 0.22, and the uniaxial compressive strength is 98.75 MPa. The mechanical parameters of the main rock-forming mineral components are listed in Table 1 [28].Table 1
Mechanical parameters of rock-forming minerals.
Rock-forming minerals
Quartz
Mica
Plagioclase
Alkali feldspar
Elastic modulusE (GPa)
90
40
81
67
Poisson’s ratioμ
0.3
0.28
0.28
0.27
Tensile strengthft (MPa)
100
60
90
90
Cohesionc (MPa)
55
30
45
45
Friction angleφ (°)
60
30
50
50
Crystal size (mm)
0.5–2.5
0.5–2.0
0.6–3.0
0.6–3.0
Proportion
30.50%
7.00%
45.30%
17.20%Figure2(a) shows the granite sample used for the indoor uniaxial compressive experiments; Figure 2(b) shows the granite particle model constructed during the numerical simulation, and Figure 2(c) shows the contact between two particles in the model. The constructed granite model contained 8818 particles, which were divided into four groups according to the mineral crystal type. The parameters for each group of particles were separately set according to the mineral mechanical parameters.Figure 2
Granite sample and numerical simulation model: (a) granite specimens, (b) granite heterogeneous model, and (c) granite contact force model.The model was built using the cluster command, which is a group of particles bonded with a certain strength to form a cluster. The bond between two particles within a cluster can be broken when the external force is sufficiently large to fail. Irregular clusters can closely reproduce the uncertain shapes of natural minerals. The difference in the parameter settings between clusters better describes the friction and interlocking behavior between grain boundaries, based on which the damage phenomena occurring inside the rock material is better demonstrated. Therefore, the mineral crystals within the rock are simulated by a cluster. The model generated 2203 particle clusters with a maximum of seven particles in each cluster. The clusters of particles in the model were randomly assigned to quartz, mica, plagioclase, and alkali feldspar, using Gaussian functions.A single-model setting approach produces higher compressive strength in the model. To more accurately restore the rock mechanical properties, the contact model was set in a hybrid manner. The contact between the wall and the particles was set to the linear model, the linear parallel bond model was used for the bond between two particles inside the particle cluster, and the smooth-joint model was chosen to construct the bonding relationship between two particle clusters. The bond strength between mineral grain clusters was set to 60% of the minimum bond strength of two adjacent minerals at that intersection [29]. The mechanical parameters applied in the model are presented in Table 2.Table 2
Model micromechanical parameters.
Contact model
Parameters
Quartz
Mica
Plagioclase
Alkali feldspar
Linear model
Effective modulusEL∗ (GPa)
10
10
10
10
Normal-to-shear stiffness ratiokL∗
0
0
0
0
Normal critical damping ratioβnL
0.5
0.5
0.5
0.5
Linear parallel bond model
Linear group
Effective modulusELP∗ (GPa)
60
20
40
35
Normal-to-shear stiffness ratiokLP∗
2.5
2.5
2.5
2.5
Friction coefficientμLP
0.6
0.6
0.6
0.6
Parallel-bond group
Bond effective modulusE¯LP∗ (GPa)
60
20
40
35
Bond normal-to-shear stiffness ratiok¯LP∗
2.5
2.5
2.5
2.5
Tensile strengthσ¯cLP (MPa)
500
300
400
400
Cohesionc¯LP (kN)
1000
600
800
800
Friction angleφ¯LP (°)
35
25
30
30
Smooth joint model
Normal stiffnessknSJ (GPa)
0.6 ×knLP
0.6 ×knLP
0.6 ×knLP
0.6 ×knLP
Shear stiffnessksSJ (GPa)
0.6 ×ksLP
0.6 ×ksLP
0.6 ×ksLP
0.6 ×ksLP
Tensile strengthσSJ (MPa)
12
12
12
12
CohesioncSJ (kN)
120
120
120
120
Joint friction angleφSJ (°)
30
30
30
30


#### 2.2.2. Uniaxial Compression Experiment and Simulation

According to the ISMR standard of the uniaxial compression test [30], the ratio of the test sample diameter to height was set to 2–2.5. We selected a Φ35 mm × 70 mm granite model as the object of study in both the uniaxial compression indoor test and discrete element simulation. In the simulation test, the top and bottom walls were set at a loading speed of 0.1 m/s. Figure 3 shows the uniaxial compression test and numerical simulation results for the rock sample.Figure 3
Uniaxial compressive test and simulation results: (a) uniaxial compressive test, (b) ball fragment in uniaxial compressive simulation, (c) DFN fractures in uniaxial compressive simulation, and (d) stress vs. strain curve.
(a)
(b)
(c)
(d)The variation in the particle structure damage and discrete fracture network (DFN) distribution during the numerical simulation coincided with the indoor test. The stress-strain curves obtained during the experiments and simulations are shown in Figure3, and the macroscopic mechanical data obtained from the simulations show agreement with the experimental data.


## 2.1. Theory of the Discrete Element Method

### 2.1.1. Fundamental Principle

Each particle in the model is a rigid body with a mass and surface that independently move and can generate displacement and rotation. The contact generated by the pairing of internal forces and torques defines the mode of interaction between two particles and is updated in real time with the identification and generation of new contacts during the calculation. The contact between two particles is achieved as soft contact, where rigid particles can overlap at the contact point, with a small overlap area, usually considered as a point, and the magnitude of the overlap is related to the contact force through the force-displacement law. The motion of the particles conforms to Newton's law of motion, and the motion of the wall is specified by the user and is independent of the contact force required to achieve the loading or restraint of the particles.The computation process takes the form of an iterative loop, in which the position and velocity are updated based on Newton's law of motion in a determined time step. Then, the model time is updated, the contact is created or deleted according to the particle position, and the forces and moments at the contact are updated according to the force-displacement law to start the next computation loop.


### 2.1.2. Contact Model

The defined law of the interparticle interaction is a contact model. The mechanical behavior of the contact determines the accuracy of the discrete element calculation results. The micromechanical behavior between individual particles and walls is generated through contact and organized in series with each other, which in turn affects the whole by the local. For the different contact types, there were significant differences in the calculation of the contact forces. The research team selected the appropriate contact models to set the contact properties between particles and walls according to the characteristics of each model, and the types of contact models included the linear model, linear parallel bond model, and smooth-joint model. The operating rules for the three contact models in the 2D model are described.The linear model provides the force transfer behavior of an infinitesimal interface. The contact forceFcL includes linear component Fl and damping component Fd. Linear force Fl controls the linear elastic friction, whereas damping force Fd controls the bonding. The contact torque McL remains zero, allowing the particles to rotate relative to each other.The update of the forces and torque at the contact in the linear model satisfies:(1)FcL=Fl+Fd,McL=0.The forces decomposed in the tangential and normal directions:(2)Fl=−Fnln^c+Fsl,Fd=−Fndn^c+Fsd,where n^c is the contact unit-normal vector, which is directed from the center of piece 1 to the center of piece 2. For the 2D model, only when the surface gap gs<0 can the linear normal force be updated; otherwise, the force remains zero:(3)Fnl=kngs,where kn is the normal stiffness. To update the shear force, it is necessary to determine whether it is in a sliding state. The shear strength of the contact is calculated by the normal linear force:(4)Fsμ=−μFnl.The linear shear force value that should be obtained at the current moment is calculated through the initial shear force:(5)Fs∗=Fsl0−ksΔδs,where μ denotes the friction coefficient, ks is the shear stiffness, and ∆δs is the relative shear‒displacement increment.The final result of linear shear force:(6)Fsl=Fs∗,Fs∗≤Fsμ,FsμFs∗/Fs∗,Fs∗>Fsμ.In the full normal and full shear states, the update of the damping force is(7)Fnd=2βnmcknδ˙n,Fsd=2βsmcksδ˙s,where(8)mc=m1m2m1+m2,where βn is the normal critical damping, βs is the tangential critical damping, ks is the tangential stiffness, δ˙n is the normal relative displacement velocity, δ˙s is the tangential relative displacement velocity, and m1 and m2 are the masses of the particles at both ends of contact. In the slipping state, Fsd is zero.The linear parallel bond model can be used to describe the mechanical behavior of the bonded material. The bond is considered a set of springs with constant normal and shear stiffnesses, which are uniformly distributed in a rectangular plane within a certain range of contact between two particles. The model produces two bonding components: linear and parallel bonding, which is equivalent to the linear model. Parallel bonding works in parallel with the former. The linear components do not resist relative rotation, and parallel bonding generates forces (Fp) and torque (Mp). Thus, the component resists relative rotation. The contact force and torque between two particles using linear parallel bonding are(9)FcLPB=Fl+Fd+Fp,McLPB=Mp.In the bonded state or when the surface gap does not exceed zero, the model will update the contact through the force-displacement law.(10)Fp=−Fnpn^c+Fsp,Fnp=Fnp0+knSΔδn,Fsp=Fsp0−ksSΔδs,where ∆δn is the relative normal displacement increment, ∆δs is the relative shear displacement increment, and S is the cross-sectional area of contact. In the 2D model, the bond cross-section is rectangular.(11)Mp=Mp0−knIΔθb,where ∆θb is the relative bend-rotation increment, and I is the moment of inertia of the parallel bond cross-section.The smooth-joint model can simulate the shear-expansion mechanical behavior of a planar interface without considering the contact direction of the local particles distributed along the interface. Before the bond exceeds the strength limit or breaks, the bonded part has a linear elastic property, and the unbonded part is linearly elastic with expansive friction, which adapts to slip by imposing a Coulomb limit on the shear force. A frictional or bonded connection can be simulated by setting two sets of particles in contact with each other in this model. The contact force(12)FcSJ=Fs,McSJ=0,Fs=−Fnsn^j+Fss,where n^j is the unit-normal vector, which points into surface 2.In the bonded state,(13)Fns=Fns0+knSΔδne⌢,Fs∗=Fss0−ksSΔδse⌢,where Fns0 and Fss0 are the smooth-joint normal force and shear force at the beginning of the timestep, respectively, and ∆δne^ is the elastic portions of the normal displacement increments. ∆δse^ is the elastic portions of the shear displacement increments.The shear strength can be computed through the normal force:(14)Fsμ=−μFns.The shear force at this time is determined based on the shear strength:(15)Fss=Fs∗,Fs∗≤Fsμ,FsμFs∗/Fs∗,Fs∗>Fsμ.In the sliding state(Fs∗=Fnμ), the slip-change callback event will cause the normal force to increase with the expansion of the shear displacement:(16)Fns=Fns+Fs∗−Fsμkskntanψ,where ψ is the dilation angle.Figure1 shows the theoretical surface and particle states for the three contact types.Figure 1
Contact theory surface and particle state: (a) linear mode, (b) linear parallel bond model, and (c) smooth-joint model.
(a)
(b)
(c)


## 2.1.1. Fundamental Principle

Each particle in the model is a rigid body with a mass and surface that independently move and can generate displacement and rotation. The contact generated by the pairing of internal forces and torques defines the mode of interaction between two particles and is updated in real time with the identification and generation of new contacts during the calculation. The contact between two particles is achieved as soft contact, where rigid particles can overlap at the contact point, with a small overlap area, usually considered as a point, and the magnitude of the overlap is related to the contact force through the force-displacement law. The motion of the particles conforms to Newton's law of motion, and the motion of the wall is specified by the user and is independent of the contact force required to achieve the loading or restraint of the particles.The computation process takes the form of an iterative loop, in which the position and velocity are updated based on Newton's law of motion in a determined time step. Then, the model time is updated, the contact is created or deleted according to the particle position, and the forces and moments at the contact are updated according to the force-displacement law to start the next computation loop.


## 2.1.2. Contact Model

The defined law of the interparticle interaction is a contact model. The mechanical behavior of the contact determines the accuracy of the discrete element calculation results. The micromechanical behavior between individual particles and walls is generated through contact and organized in series with each other, which in turn affects the whole by the local. For the different contact types, there were significant differences in the calculation of the contact forces. The research team selected the appropriate contact models to set the contact properties between particles and walls according to the characteristics of each model, and the types of contact models included the linear model, linear parallel bond model, and smooth-joint model. The operating rules for the three contact models in the 2D model are described.The linear model provides the force transfer behavior of an infinitesimal interface. The contact forceFcL includes linear component Fl and damping component Fd. Linear force Fl controls the linear elastic friction, whereas damping force Fd controls the bonding. The contact torque McL remains zero, allowing the particles to rotate relative to each other.The update of the forces and torque at the contact in the linear model satisfies:(1)FcL=Fl+Fd,McL=0.The forces decomposed in the tangential and normal directions:(2)Fl=−Fnln^c+Fsl,Fd=−Fndn^c+Fsd,where n^c is the contact unit-normal vector, which is directed from the center of piece 1 to the center of piece 2. For the 2D model, only when the surface gap gs<0 can the linear normal force be updated; otherwise, the force remains zero:(3)Fnl=kngs,where kn is the normal stiffness. To update the shear force, it is necessary to determine whether it is in a sliding state. The shear strength of the contact is calculated by the normal linear force:(4)Fsμ=−μFnl.The linear shear force value that should be obtained at the current moment is calculated through the initial shear force:(5)Fs∗=Fsl0−ksΔδs,where μ denotes the friction coefficient, ks is the shear stiffness, and ∆δs is the relative shear‒displacement increment.The final result of linear shear force:(6)Fsl=Fs∗,Fs∗≤Fsμ,FsμFs∗/Fs∗,Fs∗>Fsμ.In the full normal and full shear states, the update of the damping force is(7)Fnd=2βnmcknδ˙n,Fsd=2βsmcksδ˙s,where(8)mc=m1m2m1+m2,where βn is the normal critical damping, βs is the tangential critical damping, ks is the tangential stiffness, δ˙n is the normal relative displacement velocity, δ˙s is the tangential relative displacement velocity, and m1 and m2 are the masses of the particles at both ends of contact. In the slipping state, Fsd is zero.The linear parallel bond model can be used to describe the mechanical behavior of the bonded material. The bond is considered a set of springs with constant normal and shear stiffnesses, which are uniformly distributed in a rectangular plane within a certain range of contact between two particles. The model produces two bonding components: linear and parallel bonding, which is equivalent to the linear model. Parallel bonding works in parallel with the former. The linear components do not resist relative rotation, and parallel bonding generates forces (Fp) and torque (Mp). Thus, the component resists relative rotation. The contact force and torque between two particles using linear parallel bonding are(9)FcLPB=Fl+Fd+Fp,McLPB=Mp.In the bonded state or when the surface gap does not exceed zero, the model will update the contact through the force-displacement law.(10)Fp=−Fnpn^c+Fsp,Fnp=Fnp0+knSΔδn,Fsp=Fsp0−ksSΔδs,where ∆δn is the relative normal displacement increment, ∆δs is the relative shear displacement increment, and S is the cross-sectional area of contact. In the 2D model, the bond cross-section is rectangular.(11)Mp=Mp0−knIΔθb,where ∆θb is the relative bend-rotation increment, and I is the moment of inertia of the parallel bond cross-section.The smooth-joint model can simulate the shear-expansion mechanical behavior of a planar interface without considering the contact direction of the local particles distributed along the interface. Before the bond exceeds the strength limit or breaks, the bonded part has a linear elastic property, and the unbonded part is linearly elastic with expansive friction, which adapts to slip by imposing a Coulomb limit on the shear force. A frictional or bonded connection can be simulated by setting two sets of particles in contact with each other in this model. The contact force(12)FcSJ=Fs,McSJ=0,Fs=−Fnsn^j+Fss,where n^j is the unit-normal vector, which points into surface 2.In the bonded state,(13)Fns=Fns0+knSΔδne⌢,Fs∗=Fss0−ksSΔδse⌢,where Fns0 and Fss0 are the smooth-joint normal force and shear force at the beginning of the timestep, respectively, and ∆δne^ is the elastic portions of the normal displacement increments. ∆δse^ is the elastic portions of the shear displacement increments.The shear strength can be computed through the normal force:(14)Fsμ=−μFns.The shear force at this time is determined based on the shear strength:(15)Fss=Fs∗,Fs∗≤Fsμ,FsμFs∗/Fs∗,Fs∗>Fsμ.In the sliding state(Fs∗=Fnμ), the slip-change callback event will cause the normal force to increase with the expansion of the shear displacement:(16)Fns=Fns+Fs∗−Fsμkskntanψ,where ψ is the dilation angle.Figure1 shows the theoretical surface and particle states for the three contact types.Figure 1
Contact theory surface and particle state: (a) linear mode, (b) linear parallel bond model, and (c) smooth-joint model.
(a)
(b)
(c)


## 2.2. Microparameter Setting and Calibration

To make the physical and mechanical properties of the model more similar to the actual granite and to achieve a reliable simulation effect, the research team conducted indoor uniaxial compression tests on granite to obtain the macroscopic mechanical parameters. A discrete element model of granite was also constructed for the uniaxial compressive numerical simulation. Based on the indoor test and simulation results, the microparameters of the model were calibrated using the trial-and-error method.


### 2.2.1. Rock Sample

The selected granite sample was medium-grained granite quarried from a mine in Jilin Province, China, and the main mineral composition of the sample included feldspar, quartz, and mica. The physical and mechanical parameters of granite were experimentally measured: the density is 2.79 g/cm3, the elastic modulus is 29.8 GPa, Poisson's ratio is 0.22, and the uniaxial compressive strength is 98.75 MPa. The mechanical parameters of the main rock-forming mineral components are listed in Table 1 [28].Table 1
Mechanical parameters of rock-forming minerals.
Rock-forming minerals
Quartz
Mica
Plagioclase
Alkali feldspar
Elastic modulusE (GPa)
90
40
81
67
Poisson’s ratioμ
0.3
0.28
0.28
0.27
Tensile strengthft (MPa)
100
60
90
90
Cohesionc (MPa)
55
30
45
45
Friction angleφ (°)
60
30
50
50
Crystal size (mm)
0.5–2.5
0.5–2.0
0.6–3.0
0.6–3.0
Proportion
30.50%
7.00%
45.30%
17.20%Figure2(a) shows the granite sample used for the indoor uniaxial compressive experiments; Figure 2(b) shows the granite particle model constructed during the numerical simulation, and Figure 2(c) shows the contact between two particles in the model. The constructed granite model contained 8818 particles, which were divided into four groups according to the mineral crystal type. The parameters for each group of particles were separately set according to the mineral mechanical parameters.Figure 2
Granite sample and numerical simulation model: (a) granite specimens, (b) granite heterogeneous model, and (c) granite contact force model.The model was built using the cluster command, which is a group of particles bonded with a certain strength to form a cluster. The bond between two particles within a cluster can be broken when the external force is sufficiently large to fail. Irregular clusters can closely reproduce the uncertain shapes of natural minerals. The difference in the parameter settings between clusters better describes the friction and interlocking behavior between grain boundaries, based on which the damage phenomena occurring inside the rock material is better demonstrated. Therefore, the mineral crystals within the rock are simulated by a cluster. The model generated 2203 particle clusters with a maximum of seven particles in each cluster. The clusters of particles in the model were randomly assigned to quartz, mica, plagioclase, and alkali feldspar, using Gaussian functions.A single-model setting approach produces higher compressive strength in the model. To more accurately restore the rock mechanical properties, the contact model was set in a hybrid manner. The contact between the wall and the particles was set to the linear model, the linear parallel bond model was used for the bond between two particles inside the particle cluster, and the smooth-joint model was chosen to construct the bonding relationship between two particle clusters. The bond strength between mineral grain clusters was set to 60% of the minimum bond strength of two adjacent minerals at that intersection [29]. The mechanical parameters applied in the model are presented in Table 2.Table 2
Model micromechanical parameters.
Contact model
Parameters
Quartz
Mica
Plagioclase
Alkali feldspar
Linear model
Effective modulusEL∗ (GPa)
10
10
10
10
Normal-to-shear stiffness ratiokL∗
0
0
0
0
Normal critical damping ratioβnL
0.5
0.5
0.5
0.5
Linear parallel bond model
Linear group
Effective modulusELP∗ (GPa)
60
20
40
35
Normal-to-shear stiffness ratiokLP∗
2.5
2.5
2.5
2.5
Friction coefficientμLP
0.6
0.6
0.6
0.6
Parallel-bond group
Bond effective modulusE¯LP∗ (GPa)
60
20
40
35
Bond normal-to-shear stiffness ratiok¯LP∗
2.5
2.5
2.5
2.5
Tensile strengthσ¯cLP (MPa)
500
300
400
400
Cohesionc¯LP (kN)
1000
600
800
800
Friction angleφ¯LP (°)
35
25
30
30
Smooth joint model
Normal stiffnessknSJ (GPa)
0.6 ×knLP
0.6 ×knLP
0.6 ×knLP
0.6 ×knLP
Shear stiffnessksSJ (GPa)
0.6 ×ksLP
0.6 ×ksLP
0.6 ×ksLP
0.6 ×ksLP
Tensile strengthσSJ (MPa)
12
12
12
12
CohesioncSJ (kN)
120
120
120
120
Joint friction angleφSJ (°)
30
30
30
30


### 2.2.2. Uniaxial Compression Experiment and Simulation

According to the ISMR standard of the uniaxial compression test [30], the ratio of the test sample diameter to height was set to 2–2.5. We selected a Φ35 mm × 70 mm granite model as the object of study in both the uniaxial compression indoor test and discrete element simulation. In the simulation test, the top and bottom walls were set at a loading speed of 0.1 m/s. Figure 3 shows the uniaxial compression test and numerical simulation results for the rock sample.Figure 3
Uniaxial compressive test and simulation results: (a) uniaxial compressive test, (b) ball fragment in uniaxial compressive simulation, (c) DFN fractures in uniaxial compressive simulation, and (d) stress vs. strain curve.
(a)
(b)
(c)
(d)The variation in the particle structure damage and discrete fracture network (DFN) distribution during the numerical simulation coincided with the indoor test. The stress-strain curves obtained during the experiments and simulations are shown in Figure3, and the macroscopic mechanical data obtained from the simulations show agreement with the experimental data.


## 2.2.1. Rock Sample

The selected granite sample was medium-grained granite quarried from a mine in Jilin Province, China, and the main mineral composition of the sample included feldspar, quartz, and mica. The physical and mechanical parameters of granite were experimentally measured: the density is 2.79 g/cm3, the elastic modulus is 29.8 GPa, Poisson's ratio is 0.22, and the uniaxial compressive strength is 98.75 MPa. The mechanical parameters of the main rock-forming mineral components are listed in Table 1 [28].Table 1
Mechanical parameters of rock-forming minerals.
Rock-forming minerals
Quartz
Mica
Plagioclase
Alkali feldspar
Elastic modulusE (GPa)
90
40
81
67
Poisson’s ratioμ
0.3
0.28
0.28
0.27
Tensile strengthft (MPa)
100
60
90
90
Cohesionc (MPa)
55
30
45
45
Friction angleφ (°)
60
30
50
50
Crystal size (mm)
0.5–2.5
0.5–2.0
0.6–3.0
0.6–3.0
Proportion
30.50%
7.00%
45.30%
17.20%Figure2(a) shows the granite sample used for the indoor uniaxial compressive experiments; Figure 2(b) shows the granite particle model constructed during the numerical simulation, and Figure 2(c) shows the contact between two particles in the model. The constructed granite model contained 8818 particles, which were divided into four groups according to the mineral crystal type. The parameters for each group of particles were separately set according to the mineral mechanical parameters.Figure 2
Granite sample and numerical simulation model: (a) granite specimens, (b) granite heterogeneous model, and (c) granite contact force model.The model was built using the cluster command, which is a group of particles bonded with a certain strength to form a cluster. The bond between two particles within a cluster can be broken when the external force is sufficiently large to fail. Irregular clusters can closely reproduce the uncertain shapes of natural minerals. The difference in the parameter settings between clusters better describes the friction and interlocking behavior between grain boundaries, based on which the damage phenomena occurring inside the rock material is better demonstrated. Therefore, the mineral crystals within the rock are simulated by a cluster. The model generated 2203 particle clusters with a maximum of seven particles in each cluster. The clusters of particles in the model were randomly assigned to quartz, mica, plagioclase, and alkali feldspar, using Gaussian functions.A single-model setting approach produces higher compressive strength in the model. To more accurately restore the rock mechanical properties, the contact model was set in a hybrid manner. The contact between the wall and the particles was set to the linear model, the linear parallel bond model was used for the bond between two particles inside the particle cluster, and the smooth-joint model was chosen to construct the bonding relationship between two particle clusters. The bond strength between mineral grain clusters was set to 60% of the minimum bond strength of two adjacent minerals at that intersection [29]. The mechanical parameters applied in the model are presented in Table 2.Table 2
Model micromechanical parameters.
Contact model
Parameters
Quartz
Mica
Plagioclase
Alkali feldspar
Linear model
Effective modulusEL∗ (GPa)
10
10
10
10
Normal-to-shear stiffness ratiokL∗
0
0
0
0
Normal critical damping ratioβnL
0.5
0.5
0.5
0.5
Linear parallel bond model
Linear group
Effective modulusELP∗ (GPa)
60
20
40
35
Normal-to-shear stiffness ratiokLP∗
2.5
2.5
2.5
2.5
Friction coefficientμLP
0.6
0.6
0.6
0.6
Parallel-bond group
Bond effective modulusE¯LP∗ (GPa)
60
20
40
35
Bond normal-to-shear stiffness ratiok¯LP∗
2.5
2.5
2.5
2.5
Tensile strengthσ¯cLP (MPa)
500
300
400
400
Cohesionc¯LP (kN)
1000
600
800
800
Friction angleφ¯LP (°)
35
25
30
30
Smooth joint model
Normal stiffnessknSJ (GPa)
0.6 ×knLP
0.6 ×knLP
0.6 ×knLP
0.6 ×knLP
Shear stiffnessksSJ (GPa)
0.6 ×ksLP
0.6 ×ksLP
0.6 ×ksLP
0.6 ×ksLP
Tensile strengthσSJ (MPa)
12
12
12
12
CohesioncSJ (kN)
120
120
120
120
Joint friction angleφSJ (°)
30
30
30
30


## 2.2.2. Uniaxial Compression Experiment and Simulation

According to the ISMR standard of the uniaxial compression test [30], the ratio of the test sample diameter to height was set to 2–2.5. We selected a Φ35 mm × 70 mm granite model as the object of study in both the uniaxial compression indoor test and discrete element simulation. In the simulation test, the top and bottom walls were set at a loading speed of 0.1 m/s. Figure 3 shows the uniaxial compression test and numerical simulation results for the rock sample.Figure 3
Uniaxial compressive test and simulation results: (a) uniaxial compressive test, (b) ball fragment in uniaxial compressive simulation, (c) DFN fractures in uniaxial compressive simulation, and (d) stress vs. strain curve.
(a)
(b)
(c)
(d)The variation in the particle structure damage and discrete fracture network (DFN) distribution during the numerical simulation coincided with the indoor test. The stress-strain curves obtained during the experiments and simulations are shown in Figure3, and the macroscopic mechanical data obtained from the simulations show agreement with the experimental data.


## 3. Disc Cutter Model

Scholars have shown that crack propagation in rocks during disc cutter rock breaking is mainly caused by normal thrust, and the tangential force mainly squeezes the crushed rock with a relatively small value [31, 32]. Normal thrust has a major role in the rock-breaking process of the disc cutter, and the entire process can be simplified to the process of disc cutter intrusion into the rock mass [33, 34]. Therefore, the two-dimensional model can reflect the process of crack initiation and propagation. Considering the computational efficiency limitations, a two-dimensional rock crushing model is selected to simulate the working process.


### 3.1. Tool and Rock Model

In actual working conditions, the rolling of the disc cutter along the rock surface is provided by forceFr, and the normal force Fv exerts a squeezing effect on the rock. A groove with depth h is created in the rock during rolling. The reaction force exerted by the rock on the disc cutter is parabolically distributed along the interface between them. By virtue of the feature that ultrasonic vibration can weaken the rock strength, adding ultrasonic vibration force Fu to the cutter can reduce force Fv that needs to be exerted on the rock. The service life of the disc cutter may be extended. The ultrasonic vibration force Fu was applied in the normal direction, exerting an auxiliary ultrahigh-frequency vibration. The stress distribution is shown in Figure 4(a). By building a physical model, the interaction between the cutter shaft and the rock surface components is analyzed. Ultrasonic vibration is added to the hob crushing process and finally applied to the rock surface. The process of ultrasonic waves being added to the rock crushing process and applied to the rock surface can be determined. The physical model constructed based on this structure is shown in Figure 4(b). The motion equation for this system can be established as follows:(17)m1Z¨1=−k1Z1−c1Z˙1+k12Z2−Z1+c12Z˙2−Z˙1+k13Z3−Z1+c13Z˙3−Z˙1,m2Z¨2=k12Z1−Z2+c12Z˙1−Z˙2+k3Z3−Z2+c3Z˙3−Z˙2+k2qt−Z2+c2q˙t−Z˙2+m2g,m3Z¨3=k13Z1−Z3+c13Z˙1−Z˙3+k3Z2−Z3+c3Z˙2−Z˙3+m3g,where m1, m2, and m3 are the equivalent mass of the cutter shaft, disc cutter, and ultrasonic vibrator; k1, k2, and k3, are the equivalent stiffness of the cutter shaft, disc cutter, and ultrasonic vibrator; k12 is the equivalent stiffness of oil film on the clearance fit surface between the disc cutter and the cutter shaft; k13 is the equivalent stiffness of oil film on the clearance fit surface between the disc cutter and the shaft ultrasonic vibrator and cutter shaft; Z1, Z2, and Z3, are the displacement of the cutter shaft, disc cutter, and ultrasonic vibrator; qt is the displacement of the effective rock crushing surface; c1, c2, and c3, are the damping coefficients of the cutter shaft, disc cutter, and ultrasonic vibrator, respectively; c12 is the damping coefficient of the mating surface between the disc cutter and the shaft; c13 is the damping coefficient of the mating surface between the ultrasonic vibrator and the cutter shaft.Figure 4
Rock breaking model of the disc cutter: (a) stress distribution between the rock and disc cutter, (b) physical model, (c) simplified physical model, and (d) discrete element model.A complex model has a significant impact on calculation difficulty; therefore, reasonable simplification is necessary. Assuming that the influence of the fitting surface is disregarded, the model is further simplified (Figure4(c)). The motion equation of the simplified model is(18)m1+m3Z¨3=−k1Z3−c1Z˙3+k3Z2−Z3+c3Z˙2−Z˙3+m1+m3g,m2Z¨2=k3Z3−Z2+c3Z˙3−Z˙2+k2qt−Z2+c2q˙t−Z˙2+m2g.At the initial moment, the values of both displacementZ0 and velocity Z˙0 is 0. Z3 is the displacement of the ultrasonic vibration, and its value is(19)Z3=Asin2πft,where A is the amplitude and f is the frequency. Then Z˙3 and Z¨3 are obtained with Z3.Ideally, the crushing pit of the rock is changed to a conical curve. The coordinate system is established based on the starting point of the crushing pit, and its displacement is expressed as follows:(20)qx=ax2+bx.For the crushing pit with lengthw and depth h:(21)qx=−4hw2x2+4hwx.Substitutingx=vt into the equation, we obtain:(22)qt=−4hv2w2t2+4hvwt,where v is the linear velocity of the rolling disc cutter.Substituting equations (19) and (22) into (18), the relationship between the disc cutter and the rock crushing pit is obtained. In the solution of the equation, q is positively related to A, that is, the addition of ultrasonic vibration can increase depth h of the groove.Based on the analysis of the rock-breaking process of the disc cutter in Figures4(a)–4(c), we discovered the basic prototype of the numerical simulation model. A heterogeneous granite particle flow model with a size of 300 mm × 150 mm was established. The model contains a total of 13768 particles and 3584 particle clusters, in which the minimum particle radius is 0.9 mm and the ratio of the maximum particle size to the minimum particle size is 1.22. The particle size values were randomly specified using Gaussian functions. Considering that the strength and stiffness of the disc cutter rim are high relative to the rock material in practical engineering applications, the disc cutter is characterized by a rigid circular wall with a diameter of 432 mm. The disc cutter model is shown in Figure 4(d).


### 3.2. Loading Conditions

Because the disc cutter does not have the power to actively generate rolling in actual operation, it is subjected to the combined effect of thex-direction force parallel to the free face of rocks and the y-direction force perpendicular to the rock face. The rolling motion state of the disc cutter is produced by the contact friction between the rock and the disc cutter. Based on the above-given considerations, the velocity boundary conditions in the x- and y-directions were separately set for the disc cutter in this simulation. The rolling center of the disc cutter is determined by real-time detection of the position coordinates of the updated disc cutter wall, and the linear and angular velocities of the rolling are determined by the calculation of the displacement change amount, by which the definition of the rolling state is completed. Based on the above-given analysis, the equations of motion of the wall unit were set as follows:(23)x=vxt,y=vy0t+A⋅sin2πft,where x is the displacement of the wall along the horizontal direction (m), y is the displacement of the wall along the vertical direction (m), t is the physical time of the model (s), A is the amplitude of the ultrasonic vibration load (m), and f is the frequency of the ultrasonic vibration load (Hz).According to the actual working conditions, the horizontal rolling velocity (vx) of the disc cutter in the model was set to 5 mm/s, the disc cutter penetration velocity into the rock vy0 was 3 mm/s, and the ultrasonic vibration load applied in the vertical direction had a frequency of 20 kHz and an amplitude of 20 μm. The velocity conditions applied to the wall unit are presented as follows:(24)vy=vy0+A⋅2πfcos2πft.Figure5 shows the velocity curve of the motion of the disc cutter along the vertical direction and the displacement variation curve with time under the action of ultrasonic vibration assistance. Points a-e are the five points of one complete cycle of ultrasonic vibration after the stable operation of the model, corresponding to 0, T/4, T/2, 3T/4, and T moments in the ultrasonic vibration cycle, respectively, where T indicates the ultrasonic vibration period.Figure 5
Disc cutter velocity and displacement curve along the vertical direction: (a) velocity; (b) displacement.
(a)
(b)


## 3.1. Tool and Rock Model

In actual working conditions, the rolling of the disc cutter along the rock surface is provided by forceFr, and the normal force Fv exerts a squeezing effect on the rock. A groove with depth h is created in the rock during rolling. The reaction force exerted by the rock on the disc cutter is parabolically distributed along the interface between them. By virtue of the feature that ultrasonic vibration can weaken the rock strength, adding ultrasonic vibration force Fu to the cutter can reduce force Fv that needs to be exerted on the rock. The service life of the disc cutter may be extended. The ultrasonic vibration force Fu was applied in the normal direction, exerting an auxiliary ultrahigh-frequency vibration. The stress distribution is shown in Figure 4(a). By building a physical model, the interaction between the cutter shaft and the rock surface components is analyzed. Ultrasonic vibration is added to the hob crushing process and finally applied to the rock surface. The process of ultrasonic waves being added to the rock crushing process and applied to the rock surface can be determined. The physical model constructed based on this structure is shown in Figure 4(b). The motion equation for this system can be established as follows:(17)m1Z¨1=−k1Z1−c1Z˙1+k12Z2−Z1+c12Z˙2−Z˙1+k13Z3−Z1+c13Z˙3−Z˙1,m2Z¨2=k12Z1−Z2+c12Z˙1−Z˙2+k3Z3−Z2+c3Z˙3−Z˙2+k2qt−Z2+c2q˙t−Z˙2+m2g,m3Z¨3=k13Z1−Z3+c13Z˙1−Z˙3+k3Z2−Z3+c3Z˙2−Z˙3+m3g,where m1, m2, and m3 are the equivalent mass of the cutter shaft, disc cutter, and ultrasonic vibrator; k1, k2, and k3, are the equivalent stiffness of the cutter shaft, disc cutter, and ultrasonic vibrator; k12 is the equivalent stiffness of oil film on the clearance fit surface between the disc cutter and the cutter shaft; k13 is the equivalent stiffness of oil film on the clearance fit surface between the disc cutter and the shaft ultrasonic vibrator and cutter shaft; Z1, Z2, and Z3, are the displacement of the cutter shaft, disc cutter, and ultrasonic vibrator; qt is the displacement of the effective rock crushing surface; c1, c2, and c3, are the damping coefficients of the cutter shaft, disc cutter, and ultrasonic vibrator, respectively; c12 is the damping coefficient of the mating surface between the disc cutter and the shaft; c13 is the damping coefficient of the mating surface between the ultrasonic vibrator and the cutter shaft.Figure 4
Rock breaking model of the disc cutter: (a) stress distribution between the rock and disc cutter, (b) physical model, (c) simplified physical model, and (d) discrete element model.A complex model has a significant impact on calculation difficulty; therefore, reasonable simplification is necessary. Assuming that the influence of the fitting surface is disregarded, the model is further simplified (Figure4(c)). The motion equation of the simplified model is(18)m1+m3Z¨3=−k1Z3−c1Z˙3+k3Z2−Z3+c3Z˙2−Z˙3+m1+m3g,m2Z¨2=k3Z3−Z2+c3Z˙3−Z˙2+k2qt−Z2+c2q˙t−Z˙2+m2g.At the initial moment, the values of both displacementZ0 and velocity Z˙0 is 0. Z3 is the displacement of the ultrasonic vibration, and its value is(19)Z3=Asin2πft,where A is the amplitude and f is the frequency. Then Z˙3 and Z¨3 are obtained with Z3.Ideally, the crushing pit of the rock is changed to a conical curve. The coordinate system is established based on the starting point of the crushing pit, and its displacement is expressed as follows:(20)qx=ax2+bx.For the crushing pit with lengthw and depth h:(21)qx=−4hw2x2+4hwx.Substitutingx=vt into the equation, we obtain:(22)qt=−4hv2w2t2+4hvwt,where v is the linear velocity of the rolling disc cutter.Substituting equations (19) and (22) into (18), the relationship between the disc cutter and the rock crushing pit is obtained. In the solution of the equation, q is positively related to A, that is, the addition of ultrasonic vibration can increase depth h of the groove.Based on the analysis of the rock-breaking process of the disc cutter in Figures4(a)–4(c), we discovered the basic prototype of the numerical simulation model. A heterogeneous granite particle flow model with a size of 300 mm × 150 mm was established. The model contains a total of 13768 particles and 3584 particle clusters, in which the minimum particle radius is 0.9 mm and the ratio of the maximum particle size to the minimum particle size is 1.22. The particle size values were randomly specified using Gaussian functions. Considering that the strength and stiffness of the disc cutter rim are high relative to the rock material in practical engineering applications, the disc cutter is characterized by a rigid circular wall with a diameter of 432 mm. The disc cutter model is shown in Figure 4(d).


## 3.2. Loading Conditions

Because the disc cutter does not have the power to actively generate rolling in actual operation, it is subjected to the combined effect of thex-direction force parallel to the free face of rocks and the y-direction force perpendicular to the rock face. The rolling motion state of the disc cutter is produced by the contact friction between the rock and the disc cutter. Based on the above-given considerations, the velocity boundary conditions in the x- and y-directions were separately set for the disc cutter in this simulation. The rolling center of the disc cutter is determined by real-time detection of the position coordinates of the updated disc cutter wall, and the linear and angular velocities of the rolling are determined by the calculation of the displacement change amount, by which the definition of the rolling state is completed. Based on the above-given analysis, the equations of motion of the wall unit were set as follows:(23)x=vxt,y=vy0t+A⋅sin2πft,where x is the displacement of the wall along the horizontal direction (m), y is the displacement of the wall along the vertical direction (m), t is the physical time of the model (s), A is the amplitude of the ultrasonic vibration load (m), and f is the frequency of the ultrasonic vibration load (Hz).According to the actual working conditions, the horizontal rolling velocity (vx) of the disc cutter in the model was set to 5 mm/s, the disc cutter penetration velocity into the rock vy0 was 3 mm/s, and the ultrasonic vibration load applied in the vertical direction had a frequency of 20 kHz and an amplitude of 20 μm. The velocity conditions applied to the wall unit are presented as follows:(24)vy=vy0+A⋅2πfcos2πft.Figure5 shows the velocity curve of the motion of the disc cutter along the vertical direction and the displacement variation curve with time under the action of ultrasonic vibration assistance. Points a-e are the five points of one complete cycle of ultrasonic vibration after the stable operation of the model, corresponding to 0, T/4, T/2, 3T/4, and T moments in the ultrasonic vibration cycle, respectively, where T indicates the ultrasonic vibration period.Figure 5
Disc cutter velocity and displacement curve along the vertical direction: (a) velocity; (b) displacement.
(a)
(b)


## 4. Mechanical Behavior and Fracture Evolution

This subsection presents the numerical simulation results of the damage to hard rock caused by disc cutter rolling when ultrasonic vibration is applied to assist the action. Discrete element simulation results characterize the micromechanical properties and damage development during the simulation based on the distribution characteristics of the contact force chain, discrete fracture network, etc.


### 4.1. Force Distribution Characteristics

The stresses within the medium in the particle flow model were specifically expressed. In the iterative calculation, based on the particle displacement obtained by Newton’s law of motion, the force and moment of contact between two particles or between the particle and the wall are calculated using the force-displacement law. At this point, the properties and magnitude of the force at each contact in the model are available. The force of each contact was combined to form the overall distribution of the model force chain, and the macroscopic stresses were characterized by the microscopic contact force.


#### 4.1.1. Contact Force Chain

Figure6 shows the variation in the contact force chain inside the granite model during the rolling loading of the disc cutter under a period of ultrasonic vibration after the model was stable. The green line in the figure indicates the tensile stress chain, and the blue line indicates the compressive stress chain. Figures 6(a)–6(e) shows the characteristics of the contact force chain distribution for each of the five points marked in Figure 5 for one ultrasonic loading cycle. Points a and c are the highest point and lowest point of the disc cutter position, respectively, in one ultrasonic vibration cycle.Figure 6
Distribution of the contact force chain: (a)—(e) change in the force chain in one period.In the center of the upper surface of the rock model, which is the application point of force by the disc cutter, the rock is always under compressive stress. Inside the rock, the contact force between two particles shows alternating transformations of tensile and compressive stresses. When the disc cutter started to move down from the highest point of the periodic motion (pointa), a fan-shaped compressive stress zone of approximately 120° was generated below the application point of the internal force of the rock. The area of the stress zone expanded with the downward movement of the disc cutter (Figures 6(b) and 6(c)). The tensile stress area symmetrically appears on both sides of the compressive stress area, which is particularly evident in the shallow region of the rock. The corresponding position of the disc cutter at point d is the median point in the cyclic motion, after which the disc cutter moves to the highest point position. At this time, the initial fan-shaped compressive stress zone transformed into a semicircular annular wave and propagated downward. When a new compressive stress zone was formed (Figure 6(e)), a ring-shaped tensile stress zone was interspersed between it and the previous compressive stress wave. The compressive stress distribution was stronger in the vertical direction, and the tensile stress was dominant on the lateral and shallow surfaces.


#### 4.1.2. Contact Force

The observation of the contact forces inside the model (Figure7) indicates that the stress values inside the rock exhibit cyclic variation, with the contact point always being subjected to larger stress values. The disc cutter is in a state of downward penetration into the rock at all times during reciprocating motion under the action of ultrasonic vibration force. The stress distribution inside the rock at point c (Figure 7(c)) shows a significant difference. At this moment, the disc cutter is at the lowest point of one periodic motion and accelerates upward. The force exerted by the disc cutter on the rock surface decreased with upward movement, whereas the stress concentration generated during the previous downward movement propagated downward in the form of waves.Figure 7
Distribution of contact force (a)–(e) change of contact force in one period.


### 4.2. Fracture Evolution Characteristics

#### 4.2.1. Distribution of the Discrete Fracture Network

The discrete fracture network is a commonly used method of fracture construction and measurement in the numerical simulation of rock masses [35–37]. Figure 8 shows the type and evolution distribution with the loading time of the discrete fracture network of the rock in the pure rolling condition compared with the ultrasonic vibration-assisted loading condition. The blue lines in the two figures indicate tensile cracks caused by compression, and the red lines indicate tensile cracks caused by stretching.Figure 8
Distribution of the discrete fracture network with time: (a) DFN-type distribution under a single disc cutter, (b) DFN-type distribution assisted by ultrasonic vibration, blue lines represent tensile cracks caused by compression and red lines by tension, (c) DFN changes with time under a single disc cutter, and (d) DFN changes with time assisted by ultrasonic vibration.The results show that the trend of fracture development does not significantly differ between the two loading states, and both show the generation trend of spreading from shallow to deep and from the middle to all around with an increase in loading time. The ultrasonic vibration-assisted loading mode tends to generate more fractures in the surface layer of the rock than the conventional loading mode, in which rock cracks tend to initiate at a certain depth below the rock surface. At the early stage of loading, the ultrasonic vibration generated more cracks caused by stretching in the shallow layer of the rock. The fracture zone was generated at the same defect location in both modes, the width of the fracture zone generated by ultrasonic vibration-assisted loading was larger, and the generation time was relatively earlier.


#### 4.2.2. Microcrack

The number of cracks was counted, and the curve of the cracks inside the rock with loading time was obtained (Figure9). It is widely recognized that tensile failure characterizes the brittle failure of materials. The number of tensile cracks accounted for more than 60% in both loading modes, confirming that tensile damage is the main form of cracking within the rock during rock crushing by the disc cutter [38]. In the conventional loading mode, a crack surge occurs inside the rock, which is also referred to as the leaping rock-breaking phenomenon and verified by the crack number statistics. In comparison, the auxiliary effect of the ultrahigh circumferential alternating cyclic loading caused smoother crack expansion, with only one crack surge. Moreover, this surge phenomenon has characteristics similar to those of a surge in the traditional loading mode. Presumably, it is damaged at the same defect inside the rock. A comparison of the initial lifting time of the curves revealed that the application of an ultrahigh circumferential alternating load caused a significant increase in the onset time and number scale of the internal cracking expansion of the rock, with a 56.3% advance in time and a 30.3% increase in the total number of cracks.Figure 9
Crack number variation curve with time.


#### 4.2.3. Crack Azimuth Angle

During the whole simulation process, the iterative computation lasted 1.11 × 107 times. In the entire iterative operation process, we sampled the whole model 10 times in total, and adjacent sampling nodes were separated by 1.11 × 106 iterations. At each sampling node, we scanned the entire granite model and recorded the azimuth information of all cracks inside the rock at the current timestep. By comparing the propagation characteristics of cracks inside rock under two loading modes, the difference in action mechanisms can be discovered.The azimuthal evolution characteristics of the cracks inside the rock at different sampling moments in the conventional single rolling loading mode are shown in Figure10, and those of the ultrahigh circumferential cyclic load-assisted loading mode are shown in Figure 11. The crack azimuth information was recorded in the form of a rose diagram, in which the peripheral circumferential coordinate values 0–180 represent the azimuth of the crack, the inner circular contours represent the intensity of the crack, and the petal length of the rose diagram characterizes the frequency of the crack at a specific orientation.Figure 10
Evolution characteristics of the rock crack azimuth angle under pure rolling of the disc cutter: (a)–(j) crack azimuth distribution at 10 sampling points.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Figure 11
Evolution characteristics of the rock crack azimuth of the ultrasonic vibration assisting disc cutter: (a)–(j) crack azimuth distribution at 10 sampling points.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)As shown in Figure10, in the single-rolling mode, the number of cracks generated at the beginning of loading is small, the azimuth angle of the cracks is also maintained in a small range, and the angle of the crack in the second half of the loading is stable to maintain the vertical bias toward the forward direction. After adding the auxiliary effect of ultrahigh circumferential cyclic loading (Figure 11), more cracks were generated in the direction of the auxiliary force (90°) at the beginning of loading. The azimuth of the cracks had a larger range, although the main azimuth angle development trend was similar to that shown in Figure 10 in the later period. However, it is evident that the addition of the auxiliary load causes more cracks in the lateral direction, and there is a significant increase in the frequency of cracks in the direction of the main crack generation.


## 4.1. Force Distribution Characteristics

The stresses within the medium in the particle flow model were specifically expressed. In the iterative calculation, based on the particle displacement obtained by Newton’s law of motion, the force and moment of contact between two particles or between the particle and the wall are calculated using the force-displacement law. At this point, the properties and magnitude of the force at each contact in the model are available. The force of each contact was combined to form the overall distribution of the model force chain, and the macroscopic stresses were characterized by the microscopic contact force.


### 4.1.1. Contact Force Chain

Figure6 shows the variation in the contact force chain inside the granite model during the rolling loading of the disc cutter under a period of ultrasonic vibration after the model was stable. The green line in the figure indicates the tensile stress chain, and the blue line indicates the compressive stress chain. Figures 6(a)–6(e) shows the characteristics of the contact force chain distribution for each of the five points marked in Figure 5 for one ultrasonic loading cycle. Points a and c are the highest point and lowest point of the disc cutter position, respectively, in one ultrasonic vibration cycle.Figure 6
Distribution of the contact force chain: (a)—(e) change in the force chain in one period.In the center of the upper surface of the rock model, which is the application point of force by the disc cutter, the rock is always under compressive stress. Inside the rock, the contact force between two particles shows alternating transformations of tensile and compressive stresses. When the disc cutter started to move down from the highest point of the periodic motion (pointa), a fan-shaped compressive stress zone of approximately 120° was generated below the application point of the internal force of the rock. The area of the stress zone expanded with the downward movement of the disc cutter (Figures 6(b) and 6(c)). The tensile stress area symmetrically appears on both sides of the compressive stress area, which is particularly evident in the shallow region of the rock. The corresponding position of the disc cutter at point d is the median point in the cyclic motion, after which the disc cutter moves to the highest point position. At this time, the initial fan-shaped compressive stress zone transformed into a semicircular annular wave and propagated downward. When a new compressive stress zone was formed (Figure 6(e)), a ring-shaped tensile stress zone was interspersed between it and the previous compressive stress wave. The compressive stress distribution was stronger in the vertical direction, and the tensile stress was dominant on the lateral and shallow surfaces.


### 4.1.2. Contact Force

The observation of the contact forces inside the model (Figure7) indicates that the stress values inside the rock exhibit cyclic variation, with the contact point always being subjected to larger stress values. The disc cutter is in a state of downward penetration into the rock at all times during reciprocating motion under the action of ultrasonic vibration force. The stress distribution inside the rock at point c (Figure 7(c)) shows a significant difference. At this moment, the disc cutter is at the lowest point of one periodic motion and accelerates upward. The force exerted by the disc cutter on the rock surface decreased with upward movement, whereas the stress concentration generated during the previous downward movement propagated downward in the form of waves.Figure 7
Distribution of contact force (a)–(e) change of contact force in one period.


## 4.1.1. Contact Force Chain

Figure6 shows the variation in the contact force chain inside the granite model during the rolling loading of the disc cutter under a period of ultrasonic vibration after the model was stable. The green line in the figure indicates the tensile stress chain, and the blue line indicates the compressive stress chain. Figures 6(a)–6(e) shows the characteristics of the contact force chain distribution for each of the five points marked in Figure 5 for one ultrasonic loading cycle. Points a and c are the highest point and lowest point of the disc cutter position, respectively, in one ultrasonic vibration cycle.Figure 6
Distribution of the contact force chain: (a)—(e) change in the force chain in one period.In the center of the upper surface of the rock model, which is the application point of force by the disc cutter, the rock is always under compressive stress. Inside the rock, the contact force between two particles shows alternating transformations of tensile and compressive stresses. When the disc cutter started to move down from the highest point of the periodic motion (pointa), a fan-shaped compressive stress zone of approximately 120° was generated below the application point of the internal force of the rock. The area of the stress zone expanded with the downward movement of the disc cutter (Figures 6(b) and 6(c)). The tensile stress area symmetrically appears on both sides of the compressive stress area, which is particularly evident in the shallow region of the rock. The corresponding position of the disc cutter at point d is the median point in the cyclic motion, after which the disc cutter moves to the highest point position. At this time, the initial fan-shaped compressive stress zone transformed into a semicircular annular wave and propagated downward. When a new compressive stress zone was formed (Figure 6(e)), a ring-shaped tensile stress zone was interspersed between it and the previous compressive stress wave. The compressive stress distribution was stronger in the vertical direction, and the tensile stress was dominant on the lateral and shallow surfaces.


## 4.1.2. Contact Force

The observation of the contact forces inside the model (Figure7) indicates that the stress values inside the rock exhibit cyclic variation, with the contact point always being subjected to larger stress values. The disc cutter is in a state of downward penetration into the rock at all times during reciprocating motion under the action of ultrasonic vibration force. The stress distribution inside the rock at point c (Figure 7(c)) shows a significant difference. At this moment, the disc cutter is at the lowest point of one periodic motion and accelerates upward. The force exerted by the disc cutter on the rock surface decreased with upward movement, whereas the stress concentration generated during the previous downward movement propagated downward in the form of waves.Figure 7
Distribution of contact force (a)–(e) change of contact force in one period.


## 4.2. Fracture Evolution Characteristics

### 4.2.1. Distribution of the Discrete Fracture Network

The discrete fracture network is a commonly used method of fracture construction and measurement in the numerical simulation of rock masses [35–37]. Figure 8 shows the type and evolution distribution with the loading time of the discrete fracture network of the rock in the pure rolling condition compared with the ultrasonic vibration-assisted loading condition. The blue lines in the two figures indicate tensile cracks caused by compression, and the red lines indicate tensile cracks caused by stretching.Figure 8
Distribution of the discrete fracture network with time: (a) DFN-type distribution under a single disc cutter, (b) DFN-type distribution assisted by ultrasonic vibration, blue lines represent tensile cracks caused by compression and red lines by tension, (c) DFN changes with time under a single disc cutter, and (d) DFN changes with time assisted by ultrasonic vibration.The results show that the trend of fracture development does not significantly differ between the two loading states, and both show the generation trend of spreading from shallow to deep and from the middle to all around with an increase in loading time. The ultrasonic vibration-assisted loading mode tends to generate more fractures in the surface layer of the rock than the conventional loading mode, in which rock cracks tend to initiate at a certain depth below the rock surface. At the early stage of loading, the ultrasonic vibration generated more cracks caused by stretching in the shallow layer of the rock. The fracture zone was generated at the same defect location in both modes, the width of the fracture zone generated by ultrasonic vibration-assisted loading was larger, and the generation time was relatively earlier.


### 4.2.2. Microcrack

The number of cracks was counted, and the curve of the cracks inside the rock with loading time was obtained (Figure9). It is widely recognized that tensile failure characterizes the brittle failure of materials. The number of tensile cracks accounted for more than 60% in both loading modes, confirming that tensile damage is the main form of cracking within the rock during rock crushing by the disc cutter [38]. In the conventional loading mode, a crack surge occurs inside the rock, which is also referred to as the leaping rock-breaking phenomenon and verified by the crack number statistics. In comparison, the auxiliary effect of the ultrahigh circumferential alternating cyclic loading caused smoother crack expansion, with only one crack surge. Moreover, this surge phenomenon has characteristics similar to those of a surge in the traditional loading mode. Presumably, it is damaged at the same defect inside the rock. A comparison of the initial lifting time of the curves revealed that the application of an ultrahigh circumferential alternating load caused a significant increase in the onset time and number scale of the internal cracking expansion of the rock, with a 56.3% advance in time and a 30.3% increase in the total number of cracks.Figure 9
Crack number variation curve with time.


### 4.2.3. Crack Azimuth Angle

During the whole simulation process, the iterative computation lasted 1.11 × 107 times. In the entire iterative operation process, we sampled the whole model 10 times in total, and adjacent sampling nodes were separated by 1.11 × 106 iterations. At each sampling node, we scanned the entire granite model and recorded the azimuth information of all cracks inside the rock at the current timestep. By comparing the propagation characteristics of cracks inside rock under two loading modes, the difference in action mechanisms can be discovered.The azimuthal evolution characteristics of the cracks inside the rock at different sampling moments in the conventional single rolling loading mode are shown in Figure10, and those of the ultrahigh circumferential cyclic load-assisted loading mode are shown in Figure 11. The crack azimuth information was recorded in the form of a rose diagram, in which the peripheral circumferential coordinate values 0–180 represent the azimuth of the crack, the inner circular contours represent the intensity of the crack, and the petal length of the rose diagram characterizes the frequency of the crack at a specific orientation.Figure 10
Evolution characteristics of the rock crack azimuth angle under pure rolling of the disc cutter: (a)–(j) crack azimuth distribution at 10 sampling points.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Figure 11
Evolution characteristics of the rock crack azimuth of the ultrasonic vibration assisting disc cutter: (a)–(j) crack azimuth distribution at 10 sampling points.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)As shown in Figure10, in the single-rolling mode, the number of cracks generated at the beginning of loading is small, the azimuth angle of the cracks is also maintained in a small range, and the angle of the crack in the second half of the loading is stable to maintain the vertical bias toward the forward direction. After adding the auxiliary effect of ultrahigh circumferential cyclic loading (Figure 11), more cracks were generated in the direction of the auxiliary force (90°) at the beginning of loading. The azimuth of the cracks had a larger range, although the main azimuth angle development trend was similar to that shown in Figure 10 in the later period. However, it is evident that the addition of the auxiliary load causes more cracks in the lateral direction, and there is a significant increase in the frequency of cracks in the direction of the main crack generation.


## 4.2.1. Distribution of the Discrete Fracture Network

The discrete fracture network is a commonly used method of fracture construction and measurement in the numerical simulation of rock masses [35–37]. Figure 8 shows the type and evolution distribution with the loading time of the discrete fracture network of the rock in the pure rolling condition compared with the ultrasonic vibration-assisted loading condition. The blue lines in the two figures indicate tensile cracks caused by compression, and the red lines indicate tensile cracks caused by stretching.Figure 8
Distribution of the discrete fracture network with time: (a) DFN-type distribution under a single disc cutter, (b) DFN-type distribution assisted by ultrasonic vibration, blue lines represent tensile cracks caused by compression and red lines by tension, (c) DFN changes with time under a single disc cutter, and (d) DFN changes with time assisted by ultrasonic vibration.The results show that the trend of fracture development does not significantly differ between the two loading states, and both show the generation trend of spreading from shallow to deep and from the middle to all around with an increase in loading time. The ultrasonic vibration-assisted loading mode tends to generate more fractures in the surface layer of the rock than the conventional loading mode, in which rock cracks tend to initiate at a certain depth below the rock surface. At the early stage of loading, the ultrasonic vibration generated more cracks caused by stretching in the shallow layer of the rock. The fracture zone was generated at the same defect location in both modes, the width of the fracture zone generated by ultrasonic vibration-assisted loading was larger, and the generation time was relatively earlier.


## 4.2.2. Microcrack

The number of cracks was counted, and the curve of the cracks inside the rock with loading time was obtained (Figure9). It is widely recognized that tensile failure characterizes the brittle failure of materials. The number of tensile cracks accounted for more than 60% in both loading modes, confirming that tensile damage is the main form of cracking within the rock during rock crushing by the disc cutter [38]. In the conventional loading mode, a crack surge occurs inside the rock, which is also referred to as the leaping rock-breaking phenomenon and verified by the crack number statistics. In comparison, the auxiliary effect of the ultrahigh circumferential alternating cyclic loading caused smoother crack expansion, with only one crack surge. Moreover, this surge phenomenon has characteristics similar to those of a surge in the traditional loading mode. Presumably, it is damaged at the same defect inside the rock. A comparison of the initial lifting time of the curves revealed that the application of an ultrahigh circumferential alternating load caused a significant increase in the onset time and number scale of the internal cracking expansion of the rock, with a 56.3% advance in time and a 30.3% increase in the total number of cracks.Figure 9
Crack number variation curve with time.


## 4.2.3. Crack Azimuth Angle

During the whole simulation process, the iterative computation lasted 1.11 × 107 times. In the entire iterative operation process, we sampled the whole model 10 times in total, and adjacent sampling nodes were separated by 1.11 × 106 iterations. At each sampling node, we scanned the entire granite model and recorded the azimuth information of all cracks inside the rock at the current timestep. By comparing the propagation characteristics of cracks inside rock under two loading modes, the difference in action mechanisms can be discovered.The azimuthal evolution characteristics of the cracks inside the rock at different sampling moments in the conventional single rolling loading mode are shown in Figure10, and those of the ultrahigh circumferential cyclic load-assisted loading mode are shown in Figure 11. The crack azimuth information was recorded in the form of a rose diagram, in which the peripheral circumferential coordinate values 0–180 represent the azimuth of the crack, the inner circular contours represent the intensity of the crack, and the petal length of the rose diagram characterizes the frequency of the crack at a specific orientation.Figure 10
Evolution characteristics of the rock crack azimuth angle under pure rolling of the disc cutter: (a)–(j) crack azimuth distribution at 10 sampling points.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Figure 11
Evolution characteristics of the rock crack azimuth of the ultrasonic vibration assisting disc cutter: (a)–(j) crack azimuth distribution at 10 sampling points.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)As shown in Figure10, in the single-rolling mode, the number of cracks generated at the beginning of loading is small, the azimuth angle of the cracks is also maintained in a small range, and the angle of the crack in the second half of the loading is stable to maintain the vertical bias toward the forward direction. After adding the auxiliary effect of ultrahigh circumferential cyclic loading (Figure 11), more cracks were generated in the direction of the auxiliary force (90°) at the beginning of loading. The azimuth of the cracks had a larger range, although the main azimuth angle development trend was similar to that shown in Figure 10 in the later period. However, it is evident that the addition of the auxiliary load causes more cracks in the lateral direction, and there is a significant increase in the frequency of cracks in the direction of the main crack generation.


## 5. Discussion

### 5.1. Fracture Morphology of Load Invading Rock

According to previous studies, the mechanism of rock fragmentation by a disc cutter is that a hemispherical dense nucleus is generated after the disc cutter squeezes the rock [39]. The stress passes through the dense area to the deep area of the rock and cracks gradually form. This finding was verified by the simulation process. In Figure 8, we observe that during rock fragmentation, the starting position of crack initiation is located at a small distance from the free surface of the rock. The semicircular region within this distance is expected to be a dense nucleus.Based on the analysis of the cavity expansion theory [40], the interior of the rock during cutter indentation can be divided into three regions: dense core, plastic zone, and elastic zone (Figure 12(a)). The dense nucleus generated in the contact area between the disc cutter and the rock transferred the pressure load outward and generated a fan-shaped plastic zone. Radial compressive stress Fr, and circumferential tensile stress Ft are applied in the junction region between the elastic zone and the plastic zone. Shear deformation at the bottom of the plastic zone is the cause of the lateral crack. Combining the DFN-type distribution shown in Figures 8(a) and 8(b) with the particle displacement, the approximate relationship of the partition in the rock can be obtained. Figure 12(b) shows the particle displacement and DFN distribution under rolling load, and Figure 12(c) shows the distribution under ultrasonic assisted load. From the distribution of tension cracks caused by tension in the red-boxed region in Figure 12(b), combined with the distribution of force at the junction of the elastic zone and the plastic zone in Figure 12(a), we know that Ft is the causative factor for the tensile crack (red crack) in this region. The distribution of tension cracks also represents the boundary between the elastic zone and the plastic zone. Figure 12(c) shows the distribution of cracks and particle displacements after the addition of ultrasonic vibration assistance. The particle displacement near the cutter is larger than that in Figure 12(b), and the number and distribution of cracks caused by tension are significantly larger than the rolling load. This finding illustrates that the addition of ultrasonic vibration widens the area of the plastic zone, and strengthens area Ft between the plastic zone and the elastic zone. Based on the previous conclusions of a research group on ultrahigh cyclic loading, the generation of cracks inside the rock under this loading condition is caused by inhomogeneous tensile and compressive deformation inside the rock [17]. Therefore, the addition of ultrahigh circumferential cyclic loading produces more tensile and compressive deformation in the shallow layers of the rock, which also reduces the area of the dense core and enhances tensile stresses at the junction of the elastic and plastic zones. Finally, we observed more cracks caused by tension.Figure 12
Internal stress and partition of rock: (a) internal partition of rock; (b) fractures and partition under single disc cutter; (c) fractures and partition assisted by ultrasonic vibration. Blue lines represent tensile cracks caused by compression and red lines represent tensile cracks by tension.


### 5.2. Crack Propagation under Ultrasonic Vibration

From the perspective of microscopic mineral crystals, there are three types of cracks in the rock [41]: intragranular cracks, grain-boundary cracks, and transgranular cracks, as shown in Figure 13(c). Figures 13(a) and 13(b) show the distribution of cracks inside granite at the same position after the action of two loading modes. Intergranular cracks and grain-boundary cracks are observed in both loading modes. However, ultrasonic vibration promotes the generation of transgranular cracks and induces more intragranular cracks.Figure 13
Crack types in granite after loading: (a) cracks in granite under ultrasonic vibration-assisted loading and single rolling loading (b) local magnification of mineral crystal cracks, and (c) types of mineral crystal cracks.Ultrasonic vibration is a process of cyclic loading with an extremely small period. In the first half of one cycle, the impact body strikes the rock with kinetic energy, causing the displacement of shallow mineral crystals. The impact between two crystals generates compression and tensile force at crystal junctions, and the original cracks in the rock will be closed. During the second half cycle, the return stroke of the impactor returns the crystal to its original position. The pressure exerted on the crystal is released, and the change in displacement causes tension to appear between the two crystals. Tensile stress is also generated inside the crystal under the tension of the surrounding crystals. The ultrahigh frequency loading causes mineral crystals to be stretched and compressed at high speed. The stress on the crystal alternates between tension and compression, and fatigue failure is more likely to occur in this case, so more intragranular cracks are generated. The collision between two crystals also intensifies. The stress concentration created with the collision from the crystal junction will cause failure at weak points and extend to the two crystals with further collision. Transgranular cracks are thus generated.


### 5.3. Influence Depth of Ultrasonic Vibration

Due to the crushing energy generated by the stress wave, as well as the transformation into friction heat of mineral crystals, the stress wave cannot indefinitely propagate inside the rock. The load has a certain influence depth on the action of the rock. Hakailehto [42] proposed an attenuation model of stress waves in rock. He considered that the propagation of stress waves in rock satisfies:(25)σ=σe+σee−al1−e−alσ0−σe/σe,where σe is the rock fracture strength, and a is the brittleness coefficient of rock. l is the propagation distance, and σ0 is the initial stress of the stress wave on the rock surface.This result shows that when the stress wave is incident to the rock, the portion exceeding the rock's fracture strength exponentially decays with the propagation length.Figure8 shows that the DFN distributions obtained under the two loading modes are different to some extent. Comparing the crack distribution generated in the rock under rolling load and ultrasonic assisted load, many crack distributions under the two loading modes are coincident, with more cracks in the ultrasonic loading mode. For the cracks in the ultrasonic assisted loading that are the same as the rolling load mode, we judge them as the cracks generated by the rolling load. The remaining cracks are regarded as ultrasonic-induced. Figure 14(a) overlays the DFN distribution diagram under two loading modes into one diagram for display. The yellow lines represent the crack caused by the rolling load, and the green lines represent the crack caused by the ultrasonic vibration. With increasing depth, the granite is divided into a highly affected area with numerous cracks, a slightly affected area, and an unaffected area without obvious cracks. The division of regions only considers the promoting effect of ultrasonic vibration on cracks without considering the effect of the rolling disc cutter on rock cracks. In the high-affected area near the load surface, the stress wave contains higher energy and produces more cracks. The loss of crushing energy and crystal friction heat attenuates the stress wave, and only a small amount of fatigue damage occurs in the direction of principal stress in the slightly affected zone. The stress wave energy of ultrasonic vibration is lower than the rock strength, which is not enough to produce fracture failure.Figure 14
Influence depth of ultrasonic vibration: (a) cracks generated by two loading modes and (b) position change in rock particles under ultrasonic vibration. Yellow lines represent cracks caused by rolling load, and green lines represent ultrasonic cracks.


## 5.1. Fracture Morphology of Load Invading Rock

According to previous studies, the mechanism of rock fragmentation by a disc cutter is that a hemispherical dense nucleus is generated after the disc cutter squeezes the rock [39]. The stress passes through the dense area to the deep area of the rock and cracks gradually form. This finding was verified by the simulation process. In Figure 8, we observe that during rock fragmentation, the starting position of crack initiation is located at a small distance from the free surface of the rock. The semicircular region within this distance is expected to be a dense nucleus.Based on the analysis of the cavity expansion theory [40], the interior of the rock during cutter indentation can be divided into three regions: dense core, plastic zone, and elastic zone (Figure 12(a)). The dense nucleus generated in the contact area between the disc cutter and the rock transferred the pressure load outward and generated a fan-shaped plastic zone. Radial compressive stress Fr, and circumferential tensile stress Ft are applied in the junction region between the elastic zone and the plastic zone. Shear deformation at the bottom of the plastic zone is the cause of the lateral crack. Combining the DFN-type distribution shown in Figures 8(a) and 8(b) with the particle displacement, the approximate relationship of the partition in the rock can be obtained. Figure 12(b) shows the particle displacement and DFN distribution under rolling load, and Figure 12(c) shows the distribution under ultrasonic assisted load. From the distribution of tension cracks caused by tension in the red-boxed region in Figure 12(b), combined with the distribution of force at the junction of the elastic zone and the plastic zone in Figure 12(a), we know that Ft is the causative factor for the tensile crack (red crack) in this region. The distribution of tension cracks also represents the boundary between the elastic zone and the plastic zone. Figure 12(c) shows the distribution of cracks and particle displacements after the addition of ultrasonic vibration assistance. The particle displacement near the cutter is larger than that in Figure 12(b), and the number and distribution of cracks caused by tension are significantly larger than the rolling load. This finding illustrates that the addition of ultrasonic vibration widens the area of the plastic zone, and strengthens area Ft between the plastic zone and the elastic zone. Based on the previous conclusions of a research group on ultrahigh cyclic loading, the generation of cracks inside the rock under this loading condition is caused by inhomogeneous tensile and compressive deformation inside the rock [17]. Therefore, the addition of ultrahigh circumferential cyclic loading produces more tensile and compressive deformation in the shallow layers of the rock, which also reduces the area of the dense core and enhances tensile stresses at the junction of the elastic and plastic zones. Finally, we observed more cracks caused by tension.Figure 12
Internal stress and partition of rock: (a) internal partition of rock; (b) fractures and partition under single disc cutter; (c) fractures and partition assisted by ultrasonic vibration. Blue lines represent tensile cracks caused by compression and red lines represent tensile cracks by tension.


## 5.2. Crack Propagation under Ultrasonic Vibration

From the perspective of microscopic mineral crystals, there are three types of cracks in the rock [41]: intragranular cracks, grain-boundary cracks, and transgranular cracks, as shown in Figure 13(c). Figures 13(a) and 13(b) show the distribution of cracks inside granite at the same position after the action of two loading modes. Intergranular cracks and grain-boundary cracks are observed in both loading modes. However, ultrasonic vibration promotes the generation of transgranular cracks and induces more intragranular cracks.Figure 13
Crack types in granite after loading: (a) cracks in granite under ultrasonic vibration-assisted loading and single rolling loading (b) local magnification of mineral crystal cracks, and (c) types of mineral crystal cracks.Ultrasonic vibration is a process of cyclic loading with an extremely small period. In the first half of one cycle, the impact body strikes the rock with kinetic energy, causing the displacement of shallow mineral crystals. The impact between two crystals generates compression and tensile force at crystal junctions, and the original cracks in the rock will be closed. During the second half cycle, the return stroke of the impactor returns the crystal to its original position. The pressure exerted on the crystal is released, and the change in displacement causes tension to appear between the two crystals. Tensile stress is also generated inside the crystal under the tension of the surrounding crystals. The ultrahigh frequency loading causes mineral crystals to be stretched and compressed at high speed. The stress on the crystal alternates between tension and compression, and fatigue failure is more likely to occur in this case, so more intragranular cracks are generated. The collision between two crystals also intensifies. The stress concentration created with the collision from the crystal junction will cause failure at weak points and extend to the two crystals with further collision. Transgranular cracks are thus generated.


## 5.3. Influence Depth of Ultrasonic Vibration

Due to the crushing energy generated by the stress wave, as well as the transformation into friction heat of mineral crystals, the stress wave cannot indefinitely propagate inside the rock. The load has a certain influence depth on the action of the rock. Hakailehto [42] proposed an attenuation model of stress waves in rock. He considered that the propagation of stress waves in rock satisfies:(25)σ=σe+σee−al1−e−alσ0−σe/σe,where σe is the rock fracture strength, and a is the brittleness coefficient of rock. l is the propagation distance, and σ0 is the initial stress of the stress wave on the rock surface.This result shows that when the stress wave is incident to the rock, the portion exceeding the rock's fracture strength exponentially decays with the propagation length.Figure8 shows that the DFN distributions obtained under the two loading modes are different to some extent. Comparing the crack distribution generated in the rock under rolling load and ultrasonic assisted load, many crack distributions under the two loading modes are coincident, with more cracks in the ultrasonic loading mode. For the cracks in the ultrasonic assisted loading that are the same as the rolling load mode, we judge them as the cracks generated by the rolling load. The remaining cracks are regarded as ultrasonic-induced. Figure 14(a) overlays the DFN distribution diagram under two loading modes into one diagram for display. The yellow lines represent the crack caused by the rolling load, and the green lines represent the crack caused by the ultrasonic vibration. With increasing depth, the granite is divided into a highly affected area with numerous cracks, a slightly affected area, and an unaffected area without obvious cracks. The division of regions only considers the promoting effect of ultrasonic vibration on cracks without considering the effect of the rolling disc cutter on rock cracks. In the high-affected area near the load surface, the stress wave contains higher energy and produces more cracks. The loss of crushing energy and crystal friction heat attenuates the stress wave, and only a small amount of fatigue damage occurs in the direction of principal stress in the slightly affected zone. The stress wave energy of ultrasonic vibration is lower than the rock strength, which is not enough to produce fracture failure.Figure 14
Influence depth of ultrasonic vibration: (a) cracks generated by two loading modes and (b) position change in rock particles under ultrasonic vibration. Yellow lines represent cracks caused by rolling load, and green lines represent ultrasonic cracks.


## 6. Conclusion

In this study, we propose a new type of rock-crushing method using disc cutters. With the aid of ultrahigh-frequency loading, alternating stress was generated inside the rock, inducing fatigue failure. The strength of the rock was weakened under the influence of resonance, which improved the efficiency of rock crushing by the disc cutter. Based on the calibration of rock model parameters through laboratory experiments, a heterogeneous granite model was constructed using the discrete element method. Numerical simulations of the rock fragmentation process under the two loading modes were performed to demonstrate the mechanism and development of cracks under the rolling disc cutter assisted by ultrasonic vibration. The following conclusions were drawn:(1)
With the aid of ultrahigh-frequency loading, the disc cutter produced periodic stress waves in the rock. The alternating transformation of tensile and compressive stresses significantly contributes to the generation of stress concentrations at the original defects inside the rock, which leads to the eventual fatigue failure of the rock material. More intense tensile stress appears in the shallow layer of the rock, and more microcracks are observed here, which is conducive to the tensile failure of shallow rock.(2)
Ultrahigh-frequency loading promotes the impact between two mineral crystals, which enhances the development of internal cracks, advances the generation time of fractures, and improves the development of fractures in the transverse direction. Under these conditions, the transverse penetration of fractures inside the rock and the stripping of rock fragments become more favorable.(3)
The addition of a superhigh cycle load also reduces the fluctuation increase in cracks, which alleviates the jump crushing phenomenon. Thus, the running state of the disc cutter is more stable, which is conducive to reducing uneven and intense impact loads applied to the tool. Thus, unconventional wear and damage to tools are avoided.(4)
The stress wave attenuates due to the transformation of rock crushing energy and the consumption of crystal friction heat. The promotion of ultrasonic vibration is not infinite. With increasing depth, the influence of ultrasonic vibration gradually weakened. [41].



---

*Source: 1177745-2022-11-03.xml*