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

,e 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. ,us, parallel DDA using OpenMP is applied. ,e flattened and precrack Brazilian disc tests are first reproduced, respectively, to verify the accuracy and efficiency of the parallel DDA with virtual joints. ,en, 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.


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 highspeed 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][4][5][6][7].
Rockfall has been studied for decades [8,9]. e 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][14][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][19][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][23][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][27][28][29]. In recent years, some scholars tried to reveal the fragmentation mechanism in a physical way. e 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. e authors of [32][33][34] used bond DEM particles to form an intact rock, and the simulation suggested that the collisions with the bottom floor produce fragmentations.
ose studies yielded important and inspiring conclusions. Nevertheless, the physical simulation of fragmentation in rockfalls still needs to be investigated. e rocks used are usually bonded by sphere particles, which can be different from real situations. e 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][37][38][39][40][41][42][43]. DDA has a complete theory about contacts and mechanics involving polygonal blocks, which is clearly suitable for crack propagation study. e 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. e 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][45][46][47][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. is brings the second issue, i.e., the computational efficiency. e authors have studied OpenMP-based parallel DDA, and the calculation can be significantly accelerated [49][50][51]. e 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. e implementation of the methods into DDA is also presented. en, 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.

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. e basic element in a DDA model is called block, which is an arbitrarily shaped polygon with constant stress and strain. e displacement variables are written into a vector form: where D i is the deformation matrix of block i in a blocky system; u 0 and v 0 are the xand y-wise displacements at the centroid (x 0 , y 0 ), respectively; r 0 is the rotation of the block i; ε x and ε y represent the xand y-wise strains of the block, respectively; and c xy is the shear strain of the block i. e displacement of a point (x, y) within block i can be calculated by D i using the following equation: where T i (x, y) denotes displacement transformation matrix at P(x, y), and it is calculated by e dynamic behavior of a system can be described by 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. e damping matrix C takes the form of where η is the viscosity. In this study, the viscosity η equals zero and thus no viscous damping is introduced. e 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 c method, with parameters c � 1.0 for velocity weighting and β � 0.5 for acceleration weighting: 2 Advances in Civil Engineering where superscript n denotes the calculation step. Taking D n � 0 as the displacement at the start of the current calculation step and solving for the acceleration € D n+1 from Substituting equations (8) into (7), the velocities at the end of the current calculation step are obtained as Substituting equations (8) and (9) into (4) forms the global form 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 where K ij is the stiffness submatrices and F i is the loading submatrices of block n.

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.

Parallel DDA with Virtual
Joints. e virtual joint in DDA is a general technique to study rock crack propagation, as shown in Figure 1. e 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). ere are thus two kinds of interfaces between blocks: virtual and real interfaces. e 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]. e crack propagating paths in an intact rock are predetermined by the virtual joint system, and any crack propagates along the virtual interfaces.
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. ere are two types of failure: the tensile failure and the shear failure that are, respectively, along the normal and tangential directions. e tensile failure adopts the maximum tensile stress criterion while the shear failure employs the Mohr-Coulomb criterion. e maximum tensile and shear forces are, respectively, computed by 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 f n is the normal contact force. To determine the failure mode, two ratios (r 1 and r 2 ) are used and calculated as Advances in Civil Engineering where f s is the shear contact force. If r 1 > 1 or r 2 > 1, the tensile failure or the shear failure occurs, leading to crack initiation or propagation along the virtual joint plane. When both r 1 and r 2 exceed 1, their values should be compared further: if r 1 > r 2 , the tensile failure occurs; otherwise, the shear failure occurs. e flowchart of the cracking judgment process for virtual joints is presented in Figure 2. 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. e 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. e 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. e main subroutines are contact detection, matrix assembly, equation solver, and contact post judgment. ose subroutines can consume as much as 96% of the total calculation time according to our previous research [51]. erefore, 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. ese 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. is 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.

Verification for Accuracy and Efficiency.
e 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. e Brazilian disc has a diameter of 100 mm and flats at each end with a central angle of 10°a ccording to [52]. ree numerical models with different block quantities (2045, 4107, 6081) are constructed, as shown in Figure 3. e disc is placed between two platens each with a thickness of 20 mm. e top platen serves as the loading end and the bottom platen is fixed. e mechanical properties for blocks to form the flattened Brazilian disc are the density of 2050 kg/m 3 , Young's modulus of 10.53 GPa, and Poisson's ratio of 0.12. e strength parameters for virtual joints and real joints in the calculation are listed in Table 1. In the loading process, a time-dependent displacement constraint is applied at the loading platen at a loading rate of 2 mm/s. e calculation parameters in DDA are normal spring stiffness k n of 10 GPa, maximum time step size Δt of 1 × 10 − 7 s, and maximum displacement ratio of 1 × 10 − 4 . In the simulation, gravity is not considered.
Under the diametrical compression, the crack initiation and propagation are reproduced. e 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. e splitting process is similar to that in the experiment. e final splitting failure for the three numerical models is shown in Figure 4, compared with the experimental result [52]. e 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.
On the other hand, the efficiency of the parallel DDA to analyze the rock crack propagation is investigated. e speedup ratio, defined as the ratio between serial and parallel CPU times, is calculated. e CPU times and the calculated speedup ratios for the three numerical models are compared in Figure 5. One can observe that the parallel DDA has approximately 5 times enhanced efficiency with 6 OpenMP threads used, compared with the serial computing.
In addition, the crack propagation of the precracked Brazilian disc (with a diameter of 100 mm ) is also numerically studied. e single crack in the middle part of the disc has a length of 25 mm, 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]. e numerical model with a crack inclination of 45°is depicted in Figure 6. e same mechanical properties and calculation parameters in DDA are adopted. e simulated final failure for the precracked Brazilian disc with different crack inclination angles is presented in Figure 7, 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.

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. e disc models of single jointed rock mass under different conditions are numerically analyzed. e fracturing behaviors during impact fragmentation as well as the effects of different factors on the impact fragmentation are studied. 4 Advances in Civil Engineering e model for rockfall impact fragmentation comprises a rectangular collision plate (100 mm × 20 mm) and a disc (with a diameter of 100 mm ) formed by numerous smaller blocks (around 5500 ∼ 5600). e collision plate is placed at the bottom and keeps fixed in the whole simulation. e disc vertically impacts the collision plate. e disc contains a single crack with a width of 0.4 mm in the middle part. Different joint-slope angles θ (0°, 15°, 30°, 45°, 60°, 75°, and 90°) as well as different joint connectivity rates c (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 v 0 impacts the collision plate as the simulation begins. e same mechanical properties listed in Section 3 are employed for both blocks and joints. e calculation parameters in DDA for dynamic simulation are set as normal spring stiffness k n of 0.35 GPa, maximum time step size Δt of 1 × 10 − 7 s, and maximum displacement ratio of 1 × 10 − 4 . In the simulation, gravity is considered.
Considering different values of the joint-slope angle θ, the joint connectivity rate c, and the initial impact velocity v 0 , various cases (as listed in Table 2) are studied. Figure 9(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    Figure 9(a) and the experiment results [55] shown in Figure 9(b), the fracturing behaviors are quite consistent. In addition, the crack mechanism in impact fragmentation is clarified based on the calculation results in Case 4 (θ � 45°, c � 0.25, and v 0 � 3.5m/s). e 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. e 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   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. e 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.  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 Figure 11. 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.

Analysis of Fracture Behaviors in Impact Fragmentation.
At about 300μ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. e 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°, c � 0.25, and v 0 � 3.5 m/s) and Case 7 (θ � 90°, c � 0.25, and v 0 � 3.5 m/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. e 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.

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. e 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 θ.
is suggests that the internal damage of the specimen after impacting the collision plate is relatively large in the range of 30°∼ 45°. Figure 14 shows the comparison of the incremental crack number over time when different joint-slope angles are considered. e 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. erefore, the incremental number of cracks will reach the peak later.

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 (c � 0.25, 0.50, 0.75) are considered. Figure 15 shows the simulated maximum principal stress contours in Case 8 (c � 0.25) and Case 10 (c � 0.75).
When the joint connectivity rate is smaller, the smaller damage of the specimen appears. is 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. e specimen damage increases with an increase in the joint connectivity rate. It can be seen from Figure 15(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. Figure 16 shows that, with the increasing joint connectivity rate, the cumulative crack number grows, leading to greater internal damage of the specimen.
is 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. e 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. erefore, the greater the joint connectivity rate, the greater the damage to the specimen.

Effect of Impact Velocity.
To study the effect of the jointslope angle on the impact fragmentation, five cases (Case 11-Case 15) with different impact velocities (v 0 � 3.0m/s, 3.5m/s, 4.0m/s, 5.0m/s, 6.0m/s) are considered. Figure 18 shows the simulated maximum principal stress contours in Case 11 (v 0 � 3.0 m/s), while Figure 19 shows those results in Case 14 (v 0 � 5.0 m/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. e 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 20 compares the cumulative number of cracks in the specimen under different impact velocities. It can be seen that when the impact velocity is 3.0 m/s or 3.5 m/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   Figure 12: e simulated maximum principal tensile stress contours in (a) Case 1 and (b) Case 7 (c � 0.25 and v 0 � 3.5 m/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.
slowly before 300 μs but grows rapidly between 300 μs and 900 μs. e 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 21 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.0 m/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.0 m/s specimen drops rapidly, which indicates that the cracks no longer propagate. e greater the impact velocity is, the slower the decrease in the number of cracks is and the slower the energy dissipates. e comparison results suggest that the higher impact velocity leads to greater damage to the specimen.

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. en, validation examples about both flattened and precrack Brazilian disc tests are provided. en, 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°. e study in this paper still has space to improve. e 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.

Data Availability
e data used to support the findings of this study are included within the article and available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.