# Efficient Investigation of Rock Crack Propagation and Fracture Behaviors during Impact Fragmentation in Rockfalls Using Parallel DDA

Authors: Lu Zheng; Yihan Wu; Wei Wu; Hong Zhang; Xinyan Peng; Xuelue Zhang; Xuezhen Wu
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

The study of the rock crack propagation and fracture behaviors during impact fragmentation is important and necessary for disaster evaluation of rockfalls. Discontinuous Deformation Analysis (DDA) incorporating virtual joints can offer a powerful tool to solve such a problem. In the analysis process, the computational efficiency is critical because the mesh must be very dense to make crack propagation more realistic. Thus, parallel DDA using OpenMP is applied. The flattened and precrack Brazilian disc tests are first reproduced, respectively, to verify the accuracy and efficiency of the parallel DDA with virtual joints. Then, the impact fragmentation process is simulated and validated with corresponding laboratory experiments in terms of crack propagation results. Furthermore, the effects of joint-slope angle, joint connectivity rate, and impact velocity on rock fracture behaviors are investigated. It is concluded that the peak number of cracks occurs when the joint-slope angle ranges between 30° and 45°; the higher impact velocity and joint connectivity rate tend to cause more cracks and larger damages to the specimen.


---

## Body

## 1. Introduction

Rockfall is one of the major geohazards in mountainous regions which is capable of threatening human life and properties [1, 2]. Since many infrastructures such as high-speed roads and railways are constructed inevitably through the area that is susceptible to rockfall, the investigation of development and mobility of rockfall is important [3–7].Rockfall has been studied for decades [8, 9]. The researches of rockfall mechanics can be classified into two categories [1]: cause and motion. Rockfall can be triggered by weathering or fracturing of its surrounding [10, 11] and other factors like earthquakes [12]. When triggered, the rock will move in several modes of motion [13–15]. In addition, researches about rockfall hazard assessment were also performed [16, 17]. Experiments are also served as an effective research approach about rockfall [18–20].Compared to those aspects mentioned above, less attention has been paid to the fragmentation during rockfall [4]. Rockfall fragmentation is commonly observed in both in situ investigations and experimental tests [4, 21]. Fragmentation shows a great influence on the mechanical behaviors of rockfalls, which may greatly alter the trajectory and increase the probability of impact [22–24]. Because the physical fragmentation progress in rockfalls is complex, many computational codes adopted empirical or semiempirical models for the simulation [25]. Various fragmentation consideration approaches have been applied in rockfall risk assessment, and most of them incorporate breakage models based on In Situ Block Size Distribution (IBSD) and Rockfall Block Size Distribution (RBSD) [26–29]. In recent years, some scholars tried to reveal the fragmentation mechanism in a physical way. The authors of [30] used the Discrete Element Method (DEM) to study the impact-induced fragmentation in rockfalls and claimed that large fragments are generated only when there are open preexisting fractures or when there are fully persistent closed fractures. De Blasio and Crosta [31] used DEM to study the fragmentation and boosting of rockfalls and rock avalanches. The authors of [32–34] used bond DEM particles to form an intact rock, and the simulation suggested that the collisions with the bottom floor produce fragmentations. Those studies yielded important and inspiring conclusions.Nevertheless, the physical simulation of fragmentation in rockfalls still needs to be investigated. The rocks used are usually bonded by sphere particles, which can be different from real situations. The preexisting defects of the rock may affect the crack propagation and dynamic behavior of the fragments. Moreover, many studies were focused on the effects of impact loads on the fragmentation, but the internal factors such as preexisting crack have not been thoroughly investigated in mechanic manners. It is important because the existing cracks in rocks may significantly weaken the internal structure and dominate the fragmentation behavior.To study the effect of preexisting cracks on rockfall fragmentation, two key issues need to be solved. First, a proper numerical tool needs to be selected, which can simulate the whole process of crack propagation, i.e., the continuity to discontinuity of medium, and the large deformation and displacement. Second, the computational efficiency should be sufficiently high because a careful investigation of crack propagation requires dense mesh. Based on those conditions, a parallel two-dimensional Discontinuous Deformation Analysis (2D DDA) on the OpenMP platform is a good choice. DDA was proposed by Shi [35]; it aims to analyze the evolution of blocky systems. Many studies involving DDA have been conducted since it has been proposed [36–43]. DDA has a complete theory about contacts and mechanics involving polygonal blocks, which is clearly suitable for crack propagation study. The virtual joint technique can be adopted to address the propagation of cracks. It refers to the joints that are not present in reality but do exist in models between blocks. The medium can be treated as continuous before cracking; when stress within the medium is great enough to create cracks, the blocks can detach each other through the virtual joints, and the virtual joints then become real joints. It has been introduced into DDA to simulate crack propagation [44–48]. When implementing virtual joints, the cracks extend through joints and the individual blocks keep intact. Apparently, to better describe the propagation of cracks in rockfalls, an intact rock should be meshed into blocks that are small enough. This brings the second issue, i.e., the computational efficiency. The authors have studied OpenMP-based parallel DDA, and the calculation can be significantly accelerated [49–51]. The parallel DDA is adopted in this study to facilitate the simulation of crack propagation in rockfalls.In this paper, the effect of the preexisting crack in rocks on rockfall impact fragmentation is studied by using parallel DDA based on OpenMP. First, the theory of DDA and the involved methods of virtual joints and parallel techniques are briefly introduced. The implementation of the methods into DDA is also presented. Then, validation examples about flattened and precrack Brazilian disc tests are provided, respectively. After that, the effects of joint-slope angle, joint connectivity rate, and impact velocity on rock crack propagation in rockfall impact fragmentation are discussed in detail. Finally, conclusions are drawn.


## 2. Theory of DDA

Proposed by Dr. Shi [35], DDA has become a widely used tool in simulations involving discrete geological models such as landslides and rockfalls. The basic element in a DDA model is called block, which is an arbitrarily shaped polygon with constant stress and strain. The displacement variables are written into a vector form:(1)Di=u0,v0,r0,εx,εy,γxy,where Di is the deformation matrix of block i in a blocky system; u0 and v0 are the x- and y-wise displacements at the centroid x0,y0, respectively; r0 is the rotation of the block i; εxand εy represent the x- and y-wise strains of the block, respectively; and γxy is the shear strain of the block i. The displacement of a point x,y within block i can be calculated by Di using the following equation:(2)uv=Tix,y⋅Di,where Tix,y denotes displacement transformation matrix at Px,y, and it is calculated by(3)Tix,y=10−y−y001x−x0x−x00y−y0/20y−y0x−x0/2.The dynamic behavior of a system can be described by(4)MD¨+CD˙+KD=F,where D, D˙, and D¨ denote the matrices of displacement, velocity, and acceleration, respectively, and M, C, and K represent the matrices of mass, damping, and stiffness, respectively, of a system subject to the forcing matrix F. The damping matrix C takes the form of(5)C=ηM,where η is the viscosity. In this study, the viscosity η equals zero and thus no viscous damping is introduced. The energy loss is caused by the friction between blocks because the Mohr–Coulomb yield criterion is adopted to control the block sliding.Equation (4) is solved by Newmark’s β and γ method, with parameters γ=1.0 for velocity weighting and β=0.5 for acceleration weighting:(6)Dn+1=Dn+ΔtD˙n+Δt21−2βD¨n+2βD¨n+12,(7)D˙n+1=D˙n+Δt1−γD¨n+γD¨n+1,where superscript n denotes the calculation step. Taking Dn=0 as the displacement at the start of the current calculation step and solving for the acceleration D¨n+1 from equation (6) result in(8)D¨n+1=2Δt2Dn+1−2ΔtD˙n.Substituting equations (8) into (7), the velocities at the end of the current calculation step are obtained as(9)D˙n+1=2ΔtDn+1−D˙n.Substituting equations (8) and (9) into (4) forms the global form(10)K˜Dn+1=F˜,where K˜ and F˜ are generalized stiffness matrix and force matrix, respectively. Assuming that a blocky system contains n blocks, equation (10) can be rewritten as(11)K11K12K21K22⋯K1n⋯K2n⋮⋮Kn1Kn2⋮⋮⋯KnnD1D2⋮Dn=F1F2⋮Fn,where Kij is the stiffness submatrices and Fi is the loading submatrices of block n.


## 3. Parallel DDA with Virtual Joints and Its Verification

To model rock crack propagation during impact fragmentation in rockfalls, the domain of interest is generally discretized into a large number of triangular blocks using a triangular mesh generation method. Hence, a large-scale blocky system must be involved. To ensure efficiency, the parallel DDA with virtual joints is adopted to solve such problems.


### 3.1. Parallel DDA with Virtual Joints

The virtual joint in DDA is a general technique to study rock crack propagation, as shown in Figure1. The essential of modeling rock crack propagation lies in the representation of an intact rock by gluing adjacent blocks through virtual joints and in the representation of the cracking by removing the linkage between the glued blocks. More specifically, the gluing of blocks can be realized by specifying strong strength parameters (friction, cohesion, and tensile strength) for virtual joints; meanwhile, the rock cracking is implemented by converting the virtual joints into real joints (with weak strength parameters). There are thus two kinds of interfaces between blocks: virtual and real interfaces. The real interfaces represent the real discontinuities (including the existing joints or fractures, the primary cracks, and the newly propagated cracks) whereas the virtual interfaces indicate the block boundaries that are artificially cut in the continuous domain [45]. The crack propagating paths in an intact rock are predetermined by the virtual joint system, and any crack propagates along the virtual interfaces.Figure 1
Schematic diagram of virtual joints.With an aim to simulate continuous rock, a cohesive algorithm is applied to glue the adjacent blocks together through virtual interfaces for preventing their detachment. In DDA, when the contact between two blocks is in a locked contact state, both normal and shear contact springs are applied to prevent normal penetration and relative shear displacement. Once the contact force exceeds the joint strength, the corresponding contact spring is removed, and the contact state is changed. To ensure the bonding between adjacent blocks, the strength parameters of the virtual joints take the same values as those of the rock material. When the failure criterion is satisfied, the bonding between the adjacent blocks fails and the virtual interface becomes a real interface, which results in the crack initiation or propagation occurred at this interface. There are two types of failure: the tensile failure and the shear failure that are, respectively, along the normal and tangential directions. The tensile failure adopts the maximum tensile stress criterion while the shear failure employs the Mohr–Coulomb criterion. The maximum tensile and shear forces are, respectively, computed by(12)ft,max=σtl,(13)fs,max=cl+fntanφ,where φ, c, and σt are, respectively, the internal friction angle, cohesion, and tensile strength of the rock material; l is the contact length of the virtual joint; and fn is the normal contact force. To determine the failure mode, two ratios (r1 and r2) are used and calculated as(14)r1=fnft,max,(15)r2=fsfs,max,where fs is the shear contact force. If r1>1 or r2>1, the tensile failure or the shear failure occurs, leading to crack initiation or propagation along the virtual joint plane. When both r1 and r2 exceed 1, their values should be compared further: if r1>r2, the tensile failure occurs; otherwise, the shear failure occurs. The flowchart of the cracking judgment process for virtual joints is presented in Figure 2.Figure 2
Flowchart of the cracking judgment process for virtual joints.Nevertheless, with the incorporation of virtual joints into DDA, the problem of computational ability becomes prominent because the mesh must be very dense to make crack propagation more realistic. In this case, a parallel technique based on OpenMP is adopted to accelerate the simulation. The previous work of [51] is followed in which a parallel framework for 2D DDA was developed and the parallel 2D DDA was applied to study a large-scale rockslide (containing numerous blocks) induced by earthquake. The parallelization for 2D DDA is briefly introduced here, and the interested readers are suggested to refer to [51].A 2D DDA program can be regarded as an assembly of subroutines, and each of them has certain purposes and tasks. The main subroutines are contact detection, matrix assembly, equation solver, and contact post judgment. Those subroutines can consume as much as 96% of the total calculation time according to our previous research [51]. Therefore, by parallelizing the major subroutines, computational efficiency can be significantly improved, and the modeling of crack propagation will be easier. In this parallel strategy, independent tasks, for example, contact detection for every pair of blocks, are distributed to several threads that are forked by a master thread. These slave threads are then executed separately. However, for the tasks with data dependency or data race, some additional measures have to be implemented in advance, such as modifications of code structure or using atomic operation. Based on the previous work, the subroutine to complete the cracking judgment process for virtual joints is also executed in parallel. This is implemented by placing the compiler direction “#pragma omp parallel for” before the for loop that iterates every virtual joint, as shown in the upper shaded box in Figure 2.


### 3.2. Verification for Accuracy and Efficiency

The mechanical parameters such as tensile strength and fracture toughness of rock material can be determined by using the flattened Brazilian disc test. To verify the accuracy of the parallel DDA, the splitting process of the flattened Brazilian disc test is numerically studied. The Brazilian disc has a diameter of100mm and flats at each end with a central angle of 10° according to [52]. Three numerical models with different block quantities (2045, 4107, 6081) are constructed, as shown in Figure 3. The disc is placed between two platens each with a thickness of 20mm. The top platen serves as the loading end and the bottom platen is fixed. The mechanical properties for blocks to form the flattened Brazilian disc are the density of 2050kg/m3, Young’s modulus of 10.53GPa, and Poisson’s ratio of 0.12.Figure 3
Three numerical models with different block quantities for the flattened Brazilian disc test. (a) 2045 blocks. (b) 4107 blocks. (c) 6081 blocks.
(a)(b)(c)The strength parameters for virtual joints and real joints in the calculation are listed in Table1. In the loading process, a time-dependent displacement constraint is applied at the loading platen at a loading rate of 2mm/s. The calculation parameters in DDA are normal spring stiffness kn of 10GPa, maximum time step size Δt of 1×10−7s, and maximum displacement ratio of 1×10−4. In the simulation, gravity is not considered.Table 1
Strength parameters for virtual joints and real joints when DDA simulating the flattened Brazilian disc test.
Virtual jointsReal jointsFriction angleφ (°)38.520Cohesionc (MPa)8.291Tensile strengthσt (MPa)7.861×10−5Under the diametrical compression, the crack initiation and propagation are reproduced. The cracks first initiate from the loading regions under the top platen and above the bottom platen. With the increasing loading, the cracks propagate along the vertical middle line, and finally, they coalesce to form a penetrative crack, dividing the specimen into two parts. The splitting process is similar to that in the experiment. The final splitting failure for the three numerical models is shown in Figure4, compared with the experimental result [52]. The simulated results are in good agreement with the experiment one, suggesting the accuracy of the parallel DDA with virtual joints. In addition, a large number of microcracks appear in the top and bottom ends of the Brazilian disc, which may be attributed to the excessive stiffness of the platens.Figure 4
The final splitting failure for the three numerical models with different block quantities of blocks, compared with the experiment result [52]. (a) 2045 blocks. (b) 4107 blocks. (c) 6081 blocks. (d) experiment result.
(a)(b)(c)(d)On the other hand, the efficiency of the parallel DDA to analyze the rock crack propagation is investigated. The speedup ratio, defined as the ratio between serial and parallel CPU times, is calculated. The CPU times and the calculated speedup ratios for the three numerical models are compared in Figure5. One can observe that the parallel DDA has approximately 5 times enhanced efficiency with 6 OpenMP threads used, compared with the serial computing.Figure 5
The CPU times and the calculated speedup ratios for the three numerical models of the flattened Brazilian disc test. (a) CPU times. (b) Calculated speedup ratios.
(a)(b)In addition, the crack propagation of the precracked Brazilian disc (with a diameter of100mm) is also numerically studied. The single crack in the middle part of the disc has a length of 25mm, and different values of the crack inclination (30°, 45°, 60°, 75°) are considered. By making a compromise between efficiency and accuracy, the domain of interest is discretized into around 4000 smaller blocks for each case [53]. The numerical model with a crack inclination of 45° is depicted in Figure 6.Figure 6
Numerical model of the precracked Brazilian disc with45° crack.The same mechanical properties and calculation parameters in DDA are adopted. The simulated final failure for the precracked Brazilian disc with different crack inclination angles is presented in Figure7, compared with the experimental results [54]. It can be seen that both the crack propagation path and the fracture feature obtained by the parallel DDA accord well with the experiment results.Figure 7
Comparison of the crack propagation paths from (a) the simulated and (b) the experimental results [54] for the precracked Brazilian disc with different crack inclination angles. a1 and b1 are 30 degrees, a2 and b2 are 45 degrees, a3 and b3 are 60 degrees, and a4 and b4 are 75 degrees.
(a)(b)


## 3.1. Parallel DDA with Virtual Joints

The virtual joint in DDA is a general technique to study rock crack propagation, as shown in Figure1. The essential of modeling rock crack propagation lies in the representation of an intact rock by gluing adjacent blocks through virtual joints and in the representation of the cracking by removing the linkage between the glued blocks. More specifically, the gluing of blocks can be realized by specifying strong strength parameters (friction, cohesion, and tensile strength) for virtual joints; meanwhile, the rock cracking is implemented by converting the virtual joints into real joints (with weak strength parameters). There are thus two kinds of interfaces between blocks: virtual and real interfaces. The real interfaces represent the real discontinuities (including the existing joints or fractures, the primary cracks, and the newly propagated cracks) whereas the virtual interfaces indicate the block boundaries that are artificially cut in the continuous domain [45]. The crack propagating paths in an intact rock are predetermined by the virtual joint system, and any crack propagates along the virtual interfaces.Figure 1
Schematic diagram of virtual joints.With an aim to simulate continuous rock, a cohesive algorithm is applied to glue the adjacent blocks together through virtual interfaces for preventing their detachment. In DDA, when the contact between two blocks is in a locked contact state, both normal and shear contact springs are applied to prevent normal penetration and relative shear displacement. Once the contact force exceeds the joint strength, the corresponding contact spring is removed, and the contact state is changed. To ensure the bonding between adjacent blocks, the strength parameters of the virtual joints take the same values as those of the rock material. When the failure criterion is satisfied, the bonding between the adjacent blocks fails and the virtual interface becomes a real interface, which results in the crack initiation or propagation occurred at this interface. There are two types of failure: the tensile failure and the shear failure that are, respectively, along the normal and tangential directions. The tensile failure adopts the maximum tensile stress criterion while the shear failure employs the Mohr–Coulomb criterion. The maximum tensile and shear forces are, respectively, computed by(12)ft,max=σtl,(13)fs,max=cl+fntanφ,where φ, c, and σt are, respectively, the internal friction angle, cohesion, and tensile strength of the rock material; l is the contact length of the virtual joint; and fn is the normal contact force. To determine the failure mode, two ratios (r1 and r2) are used and calculated as(14)r1=fnft,max,(15)r2=fsfs,max,where fs is the shear contact force. If r1>1 or r2>1, the tensile failure or the shear failure occurs, leading to crack initiation or propagation along the virtual joint plane. When both r1 and r2 exceed 1, their values should be compared further: if r1>r2, the tensile failure occurs; otherwise, the shear failure occurs. The flowchart of the cracking judgment process for virtual joints is presented in Figure 2.Figure 2
Flowchart of the cracking judgment process for virtual joints.Nevertheless, with the incorporation of virtual joints into DDA, the problem of computational ability becomes prominent because the mesh must be very dense to make crack propagation more realistic. In this case, a parallel technique based on OpenMP is adopted to accelerate the simulation. The previous work of [51] is followed in which a parallel framework for 2D DDA was developed and the parallel 2D DDA was applied to study a large-scale rockslide (containing numerous blocks) induced by earthquake. The parallelization for 2D DDA is briefly introduced here, and the interested readers are suggested to refer to [51].A 2D DDA program can be regarded as an assembly of subroutines, and each of them has certain purposes and tasks. The main subroutines are contact detection, matrix assembly, equation solver, and contact post judgment. Those subroutines can consume as much as 96% of the total calculation time according to our previous research [51]. Therefore, by parallelizing the major subroutines, computational efficiency can be significantly improved, and the modeling of crack propagation will be easier. In this parallel strategy, independent tasks, for example, contact detection for every pair of blocks, are distributed to several threads that are forked by a master thread. These slave threads are then executed separately. However, for the tasks with data dependency or data race, some additional measures have to be implemented in advance, such as modifications of code structure or using atomic operation. Based on the previous work, the subroutine to complete the cracking judgment process for virtual joints is also executed in parallel. This is implemented by placing the compiler direction “#pragma omp parallel for” before the for loop that iterates every virtual joint, as shown in the upper shaded box in Figure 2.


## 3.2. Verification for Accuracy and Efficiency

The mechanical parameters such as tensile strength and fracture toughness of rock material can be determined by using the flattened Brazilian disc test. To verify the accuracy of the parallel DDA, the splitting process of the flattened Brazilian disc test is numerically studied. The Brazilian disc has a diameter of100mm and flats at each end with a central angle of 10° according to [52]. Three numerical models with different block quantities (2045, 4107, 6081) are constructed, as shown in Figure 3. The disc is placed between two platens each with a thickness of 20mm. The top platen serves as the loading end and the bottom platen is fixed. The mechanical properties for blocks to form the flattened Brazilian disc are the density of 2050kg/m3, Young’s modulus of 10.53GPa, and Poisson’s ratio of 0.12.Figure 3
Three numerical models with different block quantities for the flattened Brazilian disc test. (a) 2045 blocks. (b) 4107 blocks. (c) 6081 blocks.
(a)(b)(c)The strength parameters for virtual joints and real joints in the calculation are listed in Table1. In the loading process, a time-dependent displacement constraint is applied at the loading platen at a loading rate of 2mm/s. The calculation parameters in DDA are normal spring stiffness kn of 10GPa, maximum time step size Δt of 1×10−7s, and maximum displacement ratio of 1×10−4. In the simulation, gravity is not considered.Table 1
Strength parameters for virtual joints and real joints when DDA simulating the flattened Brazilian disc test.
Virtual jointsReal jointsFriction angleφ (°)38.520Cohesionc (MPa)8.291Tensile strengthσt (MPa)7.861×10−5Under the diametrical compression, the crack initiation and propagation are reproduced. The cracks first initiate from the loading regions under the top platen and above the bottom platen. With the increasing loading, the cracks propagate along the vertical middle line, and finally, they coalesce to form a penetrative crack, dividing the specimen into two parts. The splitting process is similar to that in the experiment. The final splitting failure for the three numerical models is shown in Figure4, compared with the experimental result [52]. The simulated results are in good agreement with the experiment one, suggesting the accuracy of the parallel DDA with virtual joints. In addition, a large number of microcracks appear in the top and bottom ends of the Brazilian disc, which may be attributed to the excessive stiffness of the platens.Figure 4
The final splitting failure for the three numerical models with different block quantities of blocks, compared with the experiment result [52]. (a) 2045 blocks. (b) 4107 blocks. (c) 6081 blocks. (d) experiment result.
(a)(b)(c)(d)On the other hand, the efficiency of the parallel DDA to analyze the rock crack propagation is investigated. The speedup ratio, defined as the ratio between serial and parallel CPU times, is calculated. The CPU times and the calculated speedup ratios for the three numerical models are compared in Figure5. One can observe that the parallel DDA has approximately 5 times enhanced efficiency with 6 OpenMP threads used, compared with the serial computing.Figure 5
The CPU times and the calculated speedup ratios for the three numerical models of the flattened Brazilian disc test. (a) CPU times. (b) Calculated speedup ratios.
(a)(b)In addition, the crack propagation of the precracked Brazilian disc (with a diameter of100mm) is also numerically studied. The single crack in the middle part of the disc has a length of 25mm, and different values of the crack inclination (30°, 45°, 60°, 75°) are considered. By making a compromise between efficiency and accuracy, the domain of interest is discretized into around 4000 smaller blocks for each case [53]. The numerical model with a crack inclination of 45° is depicted in Figure 6.Figure 6
Numerical model of the precracked Brazilian disc with45° crack.The same mechanical properties and calculation parameters in DDA are adopted. The simulated final failure for the precracked Brazilian disc with different crack inclination angles is presented in Figure7, compared with the experimental results [54]. It can be seen that both the crack propagation path and the fracture feature obtained by the parallel DDA accord well with the experiment results.Figure 7
Comparison of the crack propagation paths from (a) the simulated and (b) the experimental results [54] for the precracked Brazilian disc with different crack inclination angles. a1 and b1 are 30 degrees, a2 and b2 are 45 degrees, a3 and b3 are 60 degrees, and a4 and b4 are 75 degrees.
(a)(b)


## 4. Study on Crack Propagation and Fracture Behaviors in Impact Fragmentation

In this section, the crack propagation of jointed rock mass under impact loads will be explored. The disc models of single jointed rock mass under different conditions are numerically analyzed. The fracturing behaviors during impact fragmentation as well as the effects of different factors on the impact fragmentation are studied.The model for rockfall impact fragmentation comprises a rectangular collision plate (100mm×20mm) and a disc (with a diameter of 100mm) formed by numerous smaller blocks (around 5500∼5600). The collision plate is placed at the bottom and keeps fixed in the whole simulation. The disc vertically impacts the collision plate. The disc contains a single crack with a width of 0.4mm in the middle part. Different joint-slope angles θ (0°, 15°, 30°, 45°, 60°, 75°, and 90°) as well as different joint connectivity rates γ (0.25, 0.50, and 0.75) are considered. A representative numerical model is shown in Figure 8. To save time, the free fall process is not calculated, and instead, the disc with an initial velocity of v0 impacts the collision plate as the simulation begins. The same mechanical properties listed in Section 3 are employed for both blocks and joints. The calculation parameters in DDA for dynamic simulation are set as normal spring stiffness kn of 0.35GPa, maximum time step size Δt of 1×10−7s, and maximum displacement ratio of 1×10−4. In the simulation, gravity is considered.Figure 8
A representative numerical model (θ=45° and γ=0.50) of rockfall impact fragmentation simulation.Considering different values of the joint-slope angleθ, the joint connectivity rate γ, and the initial impact velocity v0, various cases (as listed in Table 2) are studied.Table 2
Different values of the joint-slope angleθ, the joint connectivity rate γ, and the initial impact velocity v0 for different cases.
Case no.θ (°)γv0 (m/s)Case no.θ (°)γv0 (m/s)100.253.58450.253.52150.253.59450.503.53300.253.510450.753.54450.253.511450.253.05600.253.512450.253.56750.253.513450.254.07900.253.514450.255.015450.256.0


### 4.1. Analysis of Fracture Behaviors in Impact Fragmentation

Figure9(a) shows the failure states for the impact fragmentation modeling when different joint-slope angles θ (0°, 30°, 45°, 60°, 90°) are considered. In the different cases, the joint connectivity rate γ is 0.25, and the initial impact velocity v0 is 3.5 m/s. When θ=0°, the fragmentation first occurs at the region near the impacting point, and individual blocks are cut out by the newly appeared cracks. The cracks then initiate from the middle of the joint and propagate radially toward the impact point and the top end of the specimen, and finally, they coalesce to form a penetrative crack. When θ=45°, two wing cracks perpendicular to the joint are first formed at the two ends of the prefabricated joint and gradually propagate toward the impact point and the top end of the specimen, and finally, two macrocracks are formed. When θ=30° or 60°, the crack initiation and propagation are similar to those in the case of θ=45°. When θ=90°, the stress wave is affected when it propagates from the impact point to the lower end of the joint, and it repeatedly oscillates between the impact point and the lower end of the joint, resulting in a large amount of debris near the impact point. At the same time, the cracks propagate from the upper and lower ends of the joint toward the ends of the specimen in the vertical direction, finally leading to the penetration failure of the specimen. By making a comparison between the simulated failure patterns presented in Figure 9(a) and the experiment results [55] shown in Figure 9(b), the fracturing behaviors are quite consistent.Figure 9
Comparison of failure states for the impact fragmentation modeling obtained from (a) the parallel DDA and (b) the experiment when the different joint-slope angles are considered (γ=0.25 and v0=3.5m/s). a1 and b1 are 0 degrees, a2 and b2 are 30 degrees, a3 and b3 are 45 degrees, a4 and b4 are 60 degrees, and a5 and b5 are 90 degrees.
(a)(b)In addition, the crack mechanism in impact fragmentation is clarified based on the calculation results in Case 4 (θ=45°, γ=0.25, and v0=3.5m/s). The penetrative crack formed in the specimen is the resultant failure mode of the stress distribution. In the simulation, the maximum principal stresses at different measurement points are recorded, as presented in Figure 10. The tensile stress is positive while the compressive stress is negative. It can be observed that after the specimen impacts the collision board, only the part near the impact point presents relatively large compressive stress whereas the tensile stress dominates the other parts. Measurement point A near the impact point shows gradually decreasing compressive stress after reaching the maximum compressive stress peak, suggesting the occurrence of the plastic deformation inside the specimen. The measurement points (B, C, and D) placed along the vertical direction of the specimen successively attain the peak values from the bottom to the top. At about 265μs, the measurement point E near the lower end of the joint first gets the peak stress and then rapidly decreases, which enables the stress redistribution. At about 285μs, the measurement point F near the upper end of the joint also reaches the peak stress. All the peak stresses of these two measurement points near the joint exceed the tensile strength of the specimen, and thus the tensile failure occurs. The peak values of stresses at the measurement points G and H placed near the middle of the joint are small so that the failure will not occur in the middle of the joint.Figure 10
Time histories of the maximum principal stresses at different measurement points in Case 4 (θ=45°, γ=0.25, and v0=3.5m/s).Since the tensile strength of the rock material is much smaller than its compressive strength, the tensile failure typically occurs during the impact failure process. In order to reflect the failure of the specimen more clearly, the simulated maximum principal stress contours are plotted in Figure11. One can also observe that relatively large compressive stresses appear in the area near the impact point of the specimen, but the compressive stress is only distributed in a small range near the impact point.Figure 11
The simulated maximum principal stress contours in Case 4 (θ=45°, γ=0.25, and v0=3.5m/s). (a) 100 μs, (b) 300 μs, (c) 500 μs, and (d) 800 μs.
(a)(b)(c)(d)At about300μs, the maximum tensile stress exceeds the tensile strength of the specimen, causing the virtual joints of these blocks near the two ends of the preset joint to fail, and the cracks are thus initiated. At about 500μs, the lower crack propagates toward the impact point, and the upper crack propagates toward the top end of the specimen. At about 800μs, the cracks coalesce to form a penetrative crack, and the specimen is completely penetrated and destroyed. The failure of the specimen can be divided into three stages, namely, crack initiation, crack propagation, and crack penetration.In addition, the simulated maximum principal stress contours in Case 1 (θ=0°, γ=0.25, and v0=3.5m/s) and Case 7 (θ=90°, γ=0.25, and v0=3.5m/s) are shown in Figure 12. In Case 1, many microcracks appear in the periphery of the preset joint at about 400μs, and macrocracks sprout upward in the middle of the joint. The crack propagation gradually stagnates for a while. At about 900μs, the radial cracks continue to propagate toward the impact point or the top end of the specimen. At 1400μs, the upper cracks expand to the top end and the lower cracks expand to the impact point, causing the specimen to be completely penetrated and destroyed. In Case 7, at about 300μs, the largest maximum principal stresses at the lower and upper ends of the joint, where the cracks initiate and propagate. At about 400μs, the radial cracks continue to expand; due to the collision between the specimen and the collision board, stress superimposition is generated at the impact point and the lower end of the joint, resulting in fragmentation near the impact point of the specimen. At about 1200μs, the cracks coalesce to form a penetrative crack, and the specimen is completely penetrated and destroyed.Figure 12
The simulated maximum principal tensile stress contours in (a) Case 1 and (b) Case 7 (γ=0.25 and v0=3.5m/s). a1 is 400 μs, a2 is 900 μs, a3 is 1400 μs, b1 is 300 μs, b2 is 400 μs, and b3 is 1200 μs.
(a)(b)


### 4.2. Effect of Joint-Slope Angles

To study the effect of the joint-slope angle on the impact fragmentation, seven cases (Case 1–Case 7) with different joint-slope angles (θ=0°, 15°, 30°, 45°, 60°, 75°, and 90°) are considered. Some discussions have been conducted in Section 4.1. The cumulative number of cracks generated inside the specimen can reflect the internal damage: the more cracks in the specimen, the more serious the internal damage it suffers. Figure 13 shows the comparison of the cumulative crack number at about 1800 μs. It can be seen that the number of accumulated cracks in the specimen increases continuously in the range of 0°∼30° and reaches the maximum when the joint-slope angle θ is set to 30°∼45°. After that, the number of accumulated cracks decreases with an increase of the joint-slope angle θ. This suggests that the internal damage of the specimen after impacting the collision plate is relatively large in the range of 30°∼45°.Figure 13
Comparison of the cumulative crack number at about1800μs when different joint-slope angles are considered.Figure14 shows the comparison of the incremental crack number over time when different joint-slope angles are considered. The largest incremental number of cracks occurs at 150μs∼300μs for the specimens with the joint-slope angles of 0° and 15°, it occurs at 300μs∼450μs for the specimens with the joint-slope angles of 30°, 45°, and 60°, and it occurs at 450μs∼600μs for the specimens with the joint-slope angles of 70° and 90°. One can conclude that, with the decrease of the joint-slope angle, the largest incremental number of cracks occurs earlier. Possibly because when the joint-slope angle is smaller, after being reflected by the joint surface, the stress wave is superimposed at the middle or the end of the joint to exceed the tensile strength, thereby initiating cracks and making the incremental number of cracks earlier reach the peak. When the joint-slope angle equals 90°, the joint will just slightly hinder the stress wave from propagating upwards and most of the stress wave will reflect downwards after reaching the top end of the specimen. Therefore, the incremental number of cracks will reach the peak later.Figure 14
Comparison of the incremental crack number over time when different joint-slope angles are considered.


### 4.3. Effect of the Joint Connectivity Rate

To study the effect of the joint connectivity rate on the impact fragmentation, three cases (Case 8–Case 10) with different joint connectivity rates (γ=0.25,0.50,0.75) are considered. Figure 15 shows the simulated maximum principal stress contours in Case 8 (γ=0.25) and Case 10 (γ=0.75).Figure 15
The simulated maximum principal tensile stress contours in (a) Case 8 and (b) Case 10 (θ=45° and v0=3.5m/s). a1 and b1 are 200 μs, a2 and b2 are 400 μs, and a3 is 800 μs and b3 is 1000 μs.
(a)(b)When the joint connectivity rate is smaller, the smaller damage of the specimen appears. This is because when its length is shorter, the joint shows less influence on the upward propagation of the stress wave. When the joint length is infinitely small, the specimen can be regarded as an intact rock mass without joints. The specimen damage increases with an increase in the joint connectivity rate. It can be seen from Figure15(b) that the initial cracks appear at the upper position of the lower end of the preset joint, and the stress concentrates at the lower edge of the joint surface.Figure16 shows that, with the increasing joint connectivity rate, the cumulative crack number grows, leading to greater internal damage of the specimen. This result further confirms the analysis above. For the specimens with the joint connectivity rates of 0.25 and 0.50, the growth of the cumulative crack number is slower, and the incremental crack number (as shown in Figure 17) reaches the peak at about 450μs; after 1200μs, the crack propagation basically terminates. When the joint connectivity rate is 0.75, the joint has a greater influence on the stress wave. The cumulative number of cracks in the specimen increases sharply, reaching the peak value at 600μs. After 1500μs, the crack propagation basically terminates. At 1800μs, the cumulative crack number of the specimen with the joint connectivity rate of 0.75 is 3.05 and 5.06 times that with the joint connectivity rate of 0.50 and 0.25, respectively. Therefore, the greater the joint connectivity rate, the greater the damage to the specimen.Figure 16
Comparison of the cumulative crack number over time when different joint connectivity rates are considered.Figure 17
Comparison of the incremental crack number over time when different joint connectivity rates are considered.


### 4.4. Effect of Impact Velocity

To study the effect of the joint-slope angle on the impact fragmentation, five cases (Case 11–Case 15) with different impact velocities (v0=3.0m/s,3.5m/s,4.0m/s,5.0m/s,6.0m/s) are considered.Figure18 shows the simulated maximum principal stress contours in Case 11 (v0=3.0m/s), while Figure 19 shows those results in Case 14 (v0=5.0m/s). When the impact velocity is smaller, the damage to the specimen is lower. As shown in Figure 18, the initial crack stops after extending for a certain distance, and the energy is completely consumed in the process of crack propagation, making the cracks unable to penetrate the specimen. When the specimen rebounds upwards, the cracks shrink and become microcracks. As the impact velocity increases, the peak value of the stress also increases. As shown in Figure 19, affected by the impact stress wave, the specimen forms a high-stress zone at the joint end. The wing cracks are first generated at the joint ends, and then, secondary cracks also initiate at the joint ends. Specimen debris appears near the impact point after the collision. As the impact velocity increases, the specimen is more smashed.Figure 18
The simulated maximum principal stress contours in Case 11 (θ=45°, γ=0.25, and v0=3.0m/s). (a) 200 μs, (b) 700 μs, and (c) 900 μs.
(a)(b)(c)Figure 19
The simulated maximum principal stress contours in Case 14 (θ=45°, γ=0.25, and v0=5.0m/s). (a) 200 μs, (b) 1200 μs, (c) 1500 μs, and (d) 3000 μs
(a)(b)(c)(d)Figure20 compares the cumulative number of cracks in the specimen under different impact velocities. It can be seen that when the impact velocity is 3.0m/s or 3.5m/s, the cumulative number of cracks hardly increases before 300μs. After 300μs, it begins to grow slowly, and the growth stops after 900μs, which indicates that the stress wave in the specimen is completely absorbed after 900μs and the damage process ends. Meanwhile, when the impact velocity is greater than 4 m/s, the cumulative number of cracks grows slowly before 300μs but grows rapidly between 300μs and 900μs. The growth of cracks continues even after 900μs, suggesting that the energy generated by the impact is not completely consumed when the specimen is penetrated, the stress wave is still propagating in the specimen, and the damage process is still continuing.Figure 20
Comparison of the cumulative crack number over time when different impact velocities are considered.Figure21 shows the comparison of the incremental number of cracks in the specimen under different impact velocities. It can be seen that, within the range of 3.0∼5.0m/s, the number of cracks reaches the peak value at 300∼450μs, and the crack grows the fastest in the whole process; after that, the number of cracks in the 3.0m/s specimen drops rapidly, which indicates that the cracks no longer propagate. The greater the impact velocity is, the slower the decrease in the number of cracks is and the slower the energy dissipates. The comparison results suggest that the higher impact velocity leads to greater damage to the specimen.Figure 21
Comparison of the incremental crack number over time when different impact velocities are considered.


## 4.1. Analysis of Fracture Behaviors in Impact Fragmentation

Figure9(a) shows the failure states for the impact fragmentation modeling when different joint-slope angles θ (0°, 30°, 45°, 60°, 90°) are considered. In the different cases, the joint connectivity rate γ is 0.25, and the initial impact velocity v0 is 3.5 m/s. When θ=0°, the fragmentation first occurs at the region near the impacting point, and individual blocks are cut out by the newly appeared cracks. The cracks then initiate from the middle of the joint and propagate radially toward the impact point and the top end of the specimen, and finally, they coalesce to form a penetrative crack. When θ=45°, two wing cracks perpendicular to the joint are first formed at the two ends of the prefabricated joint and gradually propagate toward the impact point and the top end of the specimen, and finally, two macrocracks are formed. When θ=30° or 60°, the crack initiation and propagation are similar to those in the case of θ=45°. When θ=90°, the stress wave is affected when it propagates from the impact point to the lower end of the joint, and it repeatedly oscillates between the impact point and the lower end of the joint, resulting in a large amount of debris near the impact point. At the same time, the cracks propagate from the upper and lower ends of the joint toward the ends of the specimen in the vertical direction, finally leading to the penetration failure of the specimen. By making a comparison between the simulated failure patterns presented in Figure 9(a) and the experiment results [55] shown in Figure 9(b), the fracturing behaviors are quite consistent.Figure 9
Comparison of failure states for the impact fragmentation modeling obtained from (a) the parallel DDA and (b) the experiment when the different joint-slope angles are considered (γ=0.25 and v0=3.5m/s). a1 and b1 are 0 degrees, a2 and b2 are 30 degrees, a3 and b3 are 45 degrees, a4 and b4 are 60 degrees, and a5 and b5 are 90 degrees.
(a)(b)In addition, the crack mechanism in impact fragmentation is clarified based on the calculation results in Case 4 (θ=45°, γ=0.25, and v0=3.5m/s). The penetrative crack formed in the specimen is the resultant failure mode of the stress distribution. In the simulation, the maximum principal stresses at different measurement points are recorded, as presented in Figure 10. The tensile stress is positive while the compressive stress is negative. It can be observed that after the specimen impacts the collision board, only the part near the impact point presents relatively large compressive stress whereas the tensile stress dominates the other parts. Measurement point A near the impact point shows gradually decreasing compressive stress after reaching the maximum compressive stress peak, suggesting the occurrence of the plastic deformation inside the specimen. The measurement points (B, C, and D) placed along the vertical direction of the specimen successively attain the peak values from the bottom to the top. At about 265μs, the measurement point E near the lower end of the joint first gets the peak stress and then rapidly decreases, which enables the stress redistribution. At about 285μs, the measurement point F near the upper end of the joint also reaches the peak stress. All the peak stresses of these two measurement points near the joint exceed the tensile strength of the specimen, and thus the tensile failure occurs. The peak values of stresses at the measurement points G and H placed near the middle of the joint are small so that the failure will not occur in the middle of the joint.Figure 10
Time histories of the maximum principal stresses at different measurement points in Case 4 (θ=45°, γ=0.25, and v0=3.5m/s).Since the tensile strength of the rock material is much smaller than its compressive strength, the tensile failure typically occurs during the impact failure process. In order to reflect the failure of the specimen more clearly, the simulated maximum principal stress contours are plotted in Figure11. One can also observe that relatively large compressive stresses appear in the area near the impact point of the specimen, but the compressive stress is only distributed in a small range near the impact point.Figure 11
The simulated maximum principal stress contours in Case 4 (θ=45°, γ=0.25, and v0=3.5m/s). (a) 100 μs, (b) 300 μs, (c) 500 μs, and (d) 800 μs.
(a)(b)(c)(d)At about300μs, the maximum tensile stress exceeds the tensile strength of the specimen, causing the virtual joints of these blocks near the two ends of the preset joint to fail, and the cracks are thus initiated. At about 500μs, the lower crack propagates toward the impact point, and the upper crack propagates toward the top end of the specimen. At about 800μs, the cracks coalesce to form a penetrative crack, and the specimen is completely penetrated and destroyed. The failure of the specimen can be divided into three stages, namely, crack initiation, crack propagation, and crack penetration.In addition, the simulated maximum principal stress contours in Case 1 (θ=0°, γ=0.25, and v0=3.5m/s) and Case 7 (θ=90°, γ=0.25, and v0=3.5m/s) are shown in Figure 12. In Case 1, many microcracks appear in the periphery of the preset joint at about 400μs, and macrocracks sprout upward in the middle of the joint. The crack propagation gradually stagnates for a while. At about 900μs, the radial cracks continue to propagate toward the impact point or the top end of the specimen. At 1400μs, the upper cracks expand to the top end and the lower cracks expand to the impact point, causing the specimen to be completely penetrated and destroyed. In Case 7, at about 300μs, the largest maximum principal stresses at the lower and upper ends of the joint, where the cracks initiate and propagate. At about 400μs, the radial cracks continue to expand; due to the collision between the specimen and the collision board, stress superimposition is generated at the impact point and the lower end of the joint, resulting in fragmentation near the impact point of the specimen. At about 1200μs, the cracks coalesce to form a penetrative crack, and the specimen is completely penetrated and destroyed.Figure 12
The simulated maximum principal tensile stress contours in (a) Case 1 and (b) Case 7 (γ=0.25 and v0=3.5m/s). a1 is 400 μs, a2 is 900 μs, a3 is 1400 μs, b1 is 300 μs, b2 is 400 μs, and b3 is 1200 μs.
(a)(b)


## 4.2. Effect of Joint-Slope Angles

To study the effect of the joint-slope angle on the impact fragmentation, seven cases (Case 1–Case 7) with different joint-slope angles (θ=0°, 15°, 30°, 45°, 60°, 75°, and 90°) are considered. Some discussions have been conducted in Section 4.1. The cumulative number of cracks generated inside the specimen can reflect the internal damage: the more cracks in the specimen, the more serious the internal damage it suffers. Figure 13 shows the comparison of the cumulative crack number at about 1800 μs. It can be seen that the number of accumulated cracks in the specimen increases continuously in the range of 0°∼30° and reaches the maximum when the joint-slope angle θ is set to 30°∼45°. After that, the number of accumulated cracks decreases with an increase of the joint-slope angle θ. This suggests that the internal damage of the specimen after impacting the collision plate is relatively large in the range of 30°∼45°.Figure 13
Comparison of the cumulative crack number at about1800μs when different joint-slope angles are considered.Figure14 shows the comparison of the incremental crack number over time when different joint-slope angles are considered. The largest incremental number of cracks occurs at 150μs∼300μs for the specimens with the joint-slope angles of 0° and 15°, it occurs at 300μs∼450μs for the specimens with the joint-slope angles of 30°, 45°, and 60°, and it occurs at 450μs∼600μs for the specimens with the joint-slope angles of 70° and 90°. One can conclude that, with the decrease of the joint-slope angle, the largest incremental number of cracks occurs earlier. Possibly because when the joint-slope angle is smaller, after being reflected by the joint surface, the stress wave is superimposed at the middle or the end of the joint to exceed the tensile strength, thereby initiating cracks and making the incremental number of cracks earlier reach the peak. When the joint-slope angle equals 90°, the joint will just slightly hinder the stress wave from propagating upwards and most of the stress wave will reflect downwards after reaching the top end of the specimen. Therefore, the incremental number of cracks will reach the peak later.Figure 14
Comparison of the incremental crack number over time when different joint-slope angles are considered.


## 4.3. Effect of the Joint Connectivity Rate

To study the effect of the joint connectivity rate on the impact fragmentation, three cases (Case 8–Case 10) with different joint connectivity rates (γ=0.25,0.50,0.75) are considered. Figure 15 shows the simulated maximum principal stress contours in Case 8 (γ=0.25) and Case 10 (γ=0.75).Figure 15
The simulated maximum principal tensile stress contours in (a) Case 8 and (b) Case 10 (θ=45° and v0=3.5m/s). a1 and b1 are 200 μs, a2 and b2 are 400 μs, and a3 is 800 μs and b3 is 1000 μs.
(a)(b)When the joint connectivity rate is smaller, the smaller damage of the specimen appears. This is because when its length is shorter, the joint shows less influence on the upward propagation of the stress wave. When the joint length is infinitely small, the specimen can be regarded as an intact rock mass without joints. The specimen damage increases with an increase in the joint connectivity rate. It can be seen from Figure15(b) that the initial cracks appear at the upper position of the lower end of the preset joint, and the stress concentrates at the lower edge of the joint surface.Figure16 shows that, with the increasing joint connectivity rate, the cumulative crack number grows, leading to greater internal damage of the specimen. This result further confirms the analysis above. For the specimens with the joint connectivity rates of 0.25 and 0.50, the growth of the cumulative crack number is slower, and the incremental crack number (as shown in Figure 17) reaches the peak at about 450μs; after 1200μs, the crack propagation basically terminates. When the joint connectivity rate is 0.75, the joint has a greater influence on the stress wave. The cumulative number of cracks in the specimen increases sharply, reaching the peak value at 600μs. After 1500μs, the crack propagation basically terminates. At 1800μs, the cumulative crack number of the specimen with the joint connectivity rate of 0.75 is 3.05 and 5.06 times that with the joint connectivity rate of 0.50 and 0.25, respectively. Therefore, the greater the joint connectivity rate, the greater the damage to the specimen.Figure 16
Comparison of the cumulative crack number over time when different joint connectivity rates are considered.Figure 17
Comparison of the incremental crack number over time when different joint connectivity rates are considered.


## 4.4. Effect of Impact Velocity

To study the effect of the joint-slope angle on the impact fragmentation, five cases (Case 11–Case 15) with different impact velocities (v0=3.0m/s,3.5m/s,4.0m/s,5.0m/s,6.0m/s) are considered.Figure18 shows the simulated maximum principal stress contours in Case 11 (v0=3.0m/s), while Figure 19 shows those results in Case 14 (v0=5.0m/s). When the impact velocity is smaller, the damage to the specimen is lower. As shown in Figure 18, the initial crack stops after extending for a certain distance, and the energy is completely consumed in the process of crack propagation, making the cracks unable to penetrate the specimen. When the specimen rebounds upwards, the cracks shrink and become microcracks. As the impact velocity increases, the peak value of the stress also increases. As shown in Figure 19, affected by the impact stress wave, the specimen forms a high-stress zone at the joint end. The wing cracks are first generated at the joint ends, and then, secondary cracks also initiate at the joint ends. Specimen debris appears near the impact point after the collision. As the impact velocity increases, the specimen is more smashed.Figure 18
The simulated maximum principal stress contours in Case 11 (θ=45°, γ=0.25, and v0=3.0m/s). (a) 200 μs, (b) 700 μs, and (c) 900 μs.
(a)(b)(c)Figure 19
The simulated maximum principal stress contours in Case 14 (θ=45°, γ=0.25, and v0=5.0m/s). (a) 200 μs, (b) 1200 μs, (c) 1500 μs, and (d) 3000 μs
(a)(b)(c)(d)Figure20 compares the cumulative number of cracks in the specimen under different impact velocities. It can be seen that when the impact velocity is 3.0m/s or 3.5m/s, the cumulative number of cracks hardly increases before 300μs. After 300μs, it begins to grow slowly, and the growth stops after 900μs, which indicates that the stress wave in the specimen is completely absorbed after 900μs and the damage process ends. Meanwhile, when the impact velocity is greater than 4 m/s, the cumulative number of cracks grows slowly before 300μs but grows rapidly between 300μs and 900μs. The growth of cracks continues even after 900μs, suggesting that the energy generated by the impact is not completely consumed when the specimen is penetrated, the stress wave is still propagating in the specimen, and the damage process is still continuing.Figure 20
Comparison of the cumulative crack number over time when different impact velocities are considered.Figure21 shows the comparison of the incremental number of cracks in the specimen under different impact velocities. It can be seen that, within the range of 3.0∼5.0m/s, the number of cracks reaches the peak value at 300∼450μs, and the crack grows the fastest in the whole process; after that, the number of cracks in the 3.0m/s specimen drops rapidly, which indicates that the cracks no longer propagate. The greater the impact velocity is, the slower the decrease in the number of cracks is and the slower the energy dissipates. The comparison results suggest that the higher impact velocity leads to greater damage to the specimen.Figure 21
Comparison of the incremental crack number over time when different impact velocities are considered.


## 5. Conclusions

In this paper, the effect of the preexisting crack in rocks on rockfall fragmentation is studied by using parallel DDA with virtual joints based on OpenMP. Although fruitful research regarding rockfalls has been conducted, the influence of preexisting cracks in rocks on the fragmentation of rockfalls is still worth studying. First, the parallel DDA with virtual joints is presented. Then, validation examples about both flattened and precrack Brazilian disc tests are provided. Then, DDA simulation cases of the rockfall impact fragmentation process are performed and validated with corresponding laboratory experiments in terms of crack propagation results. Further, the effects of joint-slope angle, joint connectivity rate, and impact velocity on fracture behaviors during rockfall impact fragmentation are discussed in detail. It is concluded that the higher impact velocity and joint connectivity rate tend to cause larger damages to the specimen. When the joint-slope angle increases from 0° to 90°, the number of cracks tops in the range of 30° to 45°.The study in this paper still has space to improve. The cases used in this study are in a 2D situation. As a future study, the 3D analysis is desired since the 3D motion and rotation obviously affect the dynamic behavior of the rockfalls. In addition, the damping effect is another factor that should be addressed, especially in rockfalls. Further studies will focus on those mentioned aspects.



---

*Source: 5901561-2021-08-31.xml*