A Hybrid FEM/DEM Approach to Estimating Rock Block Breakage due to Gravity Free-Fall

An understanding of the fragmentation of rock blocks from gravity free-fall has applications to rock fall analysis as well as block cave mining. In the case of a rock fall, blocks splitting mechanics can reduce the overall block energy and therefore reduce the travel distance of the blocks. Conversely, splitting may result in a larger number of smaller fragments that may still pose a risk despite their reduced dimensions. In the case of block caving, splitting of blocks it is believed to aid the flow of the ore column, to reduce the potential of hang-ups, and to allow for easier extraction and processing of the ore. Traditional methods for block fragmentation analysis use empirical relationships and rule-based approaches that heavily rely on the block geometry, rather than block-to-block interactions. In this paper we present a study of fragmentation processes using a hybrid finite-element/discrete-element method (FEM/DEM). The approach is capable to account for the numerical instability generally associated with the simulation of high velocity surface interactions and subsequent fracturing. To date the analysis has focused on the simulation of free-falling blocks onto a fixed surface. The initial and final block breakage has been compared against parameters including the roughness and curvature of the impacted surface, and the rock block orientation in space during free fall. We believe the results could provide useful insight into designing catchments for rock fall, as well as an increased understanding of the complexities of secondary fragmentation estimates and the range of fragmentation that could be expected in block caving.


Introduction
The breakage of rock blocks, also called fragmentation, is an important concept in many aspects of rock engineering including rock fall analysis and block cave mining. The fragmentation of a rock block during rock fall can reduce the overall block energy and reduce the travel distance of resultant subblocks; however, splitting blocks results in more fragments with trajectories that are difficult to predict [1]. The expected fragmentation of a block cave mine is a key decision parameter when determining the feasibility of such projects [2]. The reduced block size due to block splitting may aid in the flow of the ore column and reduce the potential for rock blocks to obstruct drawpoints (hang-ups) and smaller rock blocks are also easier and more cost-effective to transport and process. [3].
Typically, the estimation of rock block fragmentation often heavily relies on empirical relationships and rule-based approaches that depend on initial block geometry rather than block-to-block interactions. For example, in the mining industry, Block Cave Fragmentation (BCF) software [4] is often used to determine the fragmentation of a rock block in the ore column. This software is based on the principle that blocks with high aspect ratios are more prone to splitting. The software works in cycles to repeatedly In rock fall projects, fragmentation of rock blocks has only recently started to become a topic of interest. A fractal model of block size resulting from rock fall has been suggested [5] and discrete element modelling of gravity free fall of a spherical rock has been explored by others [6].
In order to understand the fragmentation of a rock block due to gravity, a hybrid finite element/discrete element model (FEM/DEM) has been developed. Hybrid FEM/DEM analysis has been used to simulate a variety of rock engineering problems [7], including laboratory strength tests [8][9][10][11], failure of hard rock pillars [12], brittle failure of rock slopes [13], and block cave development and subsidence analysis [14]. FEM/DEM models use a combination continuum and discontinuum; deformation of discrete elements is determined by continuum mechanics while the interaction and motion of the discrete elements is determined by discontinuum mechanics [15]. In other words, each body is represented by a single discrete element that can interact with other discrete elements, and within each discrete element is a finite element mesh that is used to model the deformation of the body. The transition from a continuum to discontinuum is determined through fracture mechanics constitutive models which split discrete elements into multiple discrete elements each with a respective finite element mesh. There are many hybrid FEM/DEM codes available, for this study ELFEN [16] was utilized. It includes a coupled, elastoplastic fracture-mechanics criterion that models the continuum/discontinuum transition realistically.
This model can be used to estimate the breakage of a rock block due to gravity free-fall in a variety of different applications. The remainder of the paper outlines the creation of this model and the steps taken to manage model instability, a sensitivity analysis of some key modelling parameters, and preliminary results of a geometric study.   3 towards the minor principal stress and propagating towards the major principal stress when differential stress exceeds the tensile strength.

FEM/DEM Model
The hopper was fixed in the x and y directions and a global gravity value of -9.81 m/s 2 was applied to the system.
A uniform triangular finite element mesh was used with a target side length of 0.1 m. Meaning the smallest equilateral triangle formed has an area of approximately 4.3e-3 m 2 . The models were not permitted to fracture below the minimum mesh size; further discussion on the impacts of this are included below. Only the rock block was permitted to fracture, and the hopper was modelled as a finite element continuum.
The material parameters used in the models are shown in table 1 and the modelling parameters are shown in table 2.  Determination of these parameters stemmed from similar FEM/DEM models as well as guidance on typical values provided from the literature [15,16].
Of particular note, the contact field parameter is typically chosen at 10-20% of the side length of the mesh and represents the thickness of the contact layer corresponding to the maximum permissible penetration. In this case, a larger contact field was chosen due to the high velocities seen in the model. Small contact fields allowed for extremely high stresses to develop in the model leading to model instability.
Similarly, the search zone parameter is typically chosen to be equal to the side length of the mesh elements and represents the size of the buffer zone that will be used to search for contacts. In this case, a slightly larger search zone was used at the expense of computational efficiency to ensure reasonable element penetrations. Since elements cannot fracture smaller than the mesh size, the energy that would be absorbed through the fracturing process is either damped, converted to kinetic energy, or converted to elastic energy. When the would-be fracturing energy is converted to kinetic energy, elements exhibit considerably large velocities. Similarly, when converted to elastic energy, the elements exhibit high degrees of deformation. These high degrees of deformation often result in model instability. The deformations and velocities are considered modelling artifacts and not representative of reality, therefore, a deletion sequence was added to the model that removes elements with velocities greater than 20 m/s and elements with timesteps less than 1e-10 seconds. The elements that are removed are assumed to fragment less than the model mesh size and are accounted for in the results processing. Figure 2 presents the results of the fragmentation analysis on a percent passing chart for the first output time step for which the block has contacted the hopper (initial) and for the output time step for which motion has stopped (final). It can be observed that the largest fragment formed during the initial contact does not reduce in size considerable as the model progresses towards equilibrium resulting in a rotation of the curve around the largest fragment from initial to final rather than a translation. Similar fragmentation changes are noted in literature during secondary fragmentation assessments of block cave mines [19].

Model Sensitivity
Damping is a key factor in fracturing models and is difficult to approximate as there is no real-world equivalent. Therefore, a sensitivity analysis was conducted on the base case model to see the effects of different values of contact damping, top frequency damping, and rigid body motion damping. The contact damping parameter, which reduces oscillations of edges between contact and non-contact states, was tested using additional values of 0.1, 0.5, and 0.7. The top frequency damping parameter, which represents the proportion of high frequencies that will be critically damped, was tested using additional values of 1, 5, and 7. The rigid body motion damping, which is momentum dependent damping, was also tested using values of 1 and 3 in addition to the base value of 2.
Higher values of contact damping appear to reduce the number of deleted particles in the model due to the kinetic and elastic energy artifacts. It is unclear how many elements will break below the mesh size, in reality, so this would be a key parameter to adjust during model calibration. Broadly, the top frequency damping parameter had minimal effect on the resulting block size distributions of the models and therefore the highest damping frequency was used for the remainder of the analysis to improve model stability.
The rigid body motion damping has little effect on the models when reduced to 1, however when increased to 3, the models slow in velocity considerably and the fragmentation decreases significantly. This was considered not representative of reality even though it would increase model stability and the base case of 2 was used for rigid body damping for the remainder of the assessment.

Block Rotation
The role of the rotation of the square rock block on initial and final fragmentation was explored. The zero-rotation condition was set to be when the flat edge of the block is parallel to the lower hopper surface. Models were created with 15-degree incremental rotations both clockwise (positive) and counterclockwise (negative). Therefore, a total of six models were used to determine the effect of rock block rotation at angles of -30-degrees, -15-degrees, 0-degrees, 15-degrees, 30-degrees, and 45-degrees (the base case) Broadly, the initial fragmentation curves of the rock blocks for all models except the 0-degree model are similar. The 0-degree model shows high levels of initial fragmentation due to line-contact of the rock block to the hopper. The final fragmentation curves also show the 0-degree model has the highest fragmentation, with both the positive and negative 15-degree models showing slightly lower fragmentation, similar to the base case, and the positive and negative 30-degree models showing the least fragmentation. The slight eccentricity of these models breaks the block unevenly such that the largest fragment is more than half the size of the original block, the other models tend to break the block such that the largest fragment is half the size of the original block or less.
In order to better understand the differences between these models, the block size percentiles of D80, D60, and D30 were plotted against the rotation angle ( figure 3). These plots show the same trends noted from the fragmentation curves, but also show the similarities and differences in the model fragmentation with block rotation angle.
Intuitively, it would be expected that the models would be perfectly symmetrical (i.e., the counterclockwise models would be the same as the clockwise models). Although they show similarities, the models are not identical indicating that the fracturing in these models is a random process and is slightly mesh dependent. The randomness and mesh dependency of these types of models is discussed in Munjiza [15], however; further study should be conducted to characterize these issues.

Surface Curvature
The curvature of the contact surface has also been explored in this study. Two models were made with the lower surface of the hopper being either concave or convex. The height of the apex of the concavity/convexity was modelled to be one metre.
The results of these models show that both the convex and concave models predict less fragmentation than the flat model ( figure 4). However, the mechanisms behind these reduced fragmentations appear to be different. In the concave model, the concavity reduces the ability of the rock blocks to travel laterally which results in less high velocity block-to-block interactions. In the convex model, after initial breakage, the rock blocks are immediately carried away from each other to the edges of the model limiting the block-to-block interaction. Overall, the concavity/convexity of the models reduces blockto-block interactions and therefore reduces the overall fragmentation.

Surface Roughness
An analysis of the roughness of the contact surface was conducted. To quantify the roughness of the surface, the JRC chart from Barton and Choubey [20] was digitized. The JRC profiles represent surfaces that are 10 cm long, therefore, the profile needs to be scaled in order to fit the 10 m wide base of the hopper. The JCR value is not scale independent. From the chart in Barton [21] the profile needed to be scale by a factor of 80 in the vertical direction and a factor of 100 in the horizontal direction in order to retain the original JRC value when scaled from 10 cm to 10 m. These profiles were used as the hopper base geometry for 10 models ranging from JRC 0 through 18 at steps of 2.
Block rotation angles of 0-degrees and 45-degrees were both modelled to characterize the difference in line-contact versus point-contact of the blocks. The initial fragmentation of the 0-degree block rotation models is much more variable compared to the 45-degree block rotation models, with a minimum block size observed at a JRC value of 6 and a maximum block size observed at a JRC value of 16.
Overall, it appears the JRC number does not have any particular control on the resulting fragmentation. It should be noted that a large change in the surface characteristics of the JRC profiles can be seen starting around JRC of 10. For lower JRC profiles, the surfaces show small scale roughness, for JRC 10 and higher, that roughness becomes more wavey in nature. This could result in the change of fragmentation character observed in these models. Another method of characterising the roughness of the hopper surface may show stronger fragmentation trends.

Discussion
The models presented above allow for simple variation of geometric and strength parameters to observe the change in block fragmentation due to gravity free-fall. The observations taken from these models can be used to inform calculations for block fragmentation in block caving and rock fall settings. As an example, BCF software [4] assumes that blocks will break into two equally sized portions and as shown in these models, that is not the case. Further study is ongoing to characterise more realistic block breakage trends that can help predict expected fragmentation using the models presented here as a starting point.