Combined Finite-Discrete Element Modelling of Dynamic Rock Fracture and Fragmentation during Mining Production Process by Blast

A combined finite-discrete element method (FDEM) is proposed to model the dynamic fracture, fragmentation, and resultant muck-piling process during mining production by blast in underground mine. The key component of the proposed method, that is, transition from continuum to discontinuum through fracture and fragmentation, is introduced in detail, which makes the proposed method superior to the continuum-based finite element method and discontinuum-based discrete element method. The FDEM is calibrated by modelling the crater formation process by blast. The FDEM has well modelled the stress and fracture propagation and resultant fragmentation process. In addition, the proposed method has well captured the crushed zone, cracked zone, and the radial long crack zone. After that, the FDEM is employed to model the dynamic fracture and resultant fragmentation process by blast during sublevel caving process in an underground mine. Then the FDEM has well modelled the stress propagation process, as well as the fracture initiation and fragmenting process. Finally, the effects of borehole spacing and initial gas pressure are discussed. It is concluded that the FDEM is a value numerical approach to study the dynamic rock fracture process by blast.


Introduction
Blasting is widely used in geoengineering projects for rock fragmentation, hard rock tunnelling, and structure demolition [1]. To easily employ the blasting technology, empirical models or equations are proposed and implemented in the civil and mining industries [2][3][4]. Rosin and Rammler proposed the Rosin-Rammler equation to characterize the partial size distribution of material [5]. Kuznetsov developed a semiempirical equation to estimate size distribution of rock fragments [4]. Empirical models and equations are successfully implemented in the routine blasting designs due to the fact that only a few parameters are required. However, limited input parameters may result in an inaccurate prediction. In addition, it is time-and money-consuming to take several blasting tests during the blasting design using empirical models. As the blasting process is an extremely complex process, further studies are needed for benefiting the geoengineering in breaking geomaterials, for example, tunnelling, mining, and demolition.
For better understanding the rock fracture process, many experimental research studies have been done. Figure 1 illustrates the fracture patterns of fully contained explosive (Figure 1(a)) and the fracture pattern with one free surface (Figure 1(b)). As can be seen in Figure 1(a), the rock fracture pattern for a fully contained explosive with unlimited boundaries includes a crushed zone, cracked zone, and radial long crack zone. If there is a free surface (Figure 1(b)), the produced fragments will be pushed out by the high-pressure gas through the free surface. e crushed zone is mainly produced by the shock wave and compressive wave while the explosive is denoted. e cracked zone is formed due to the combination of the compressive wave and high-pressure gas. e long cracks are induced by the highpressure gas expansion. In addition, the compressive stress will be reflected at the free surface of the existing cracks or joints in the rock mass, and the tensile stresses are induced by the reflection. Due to the fact that the tensile strength is much smaller than the compressive strength for rock mass, the free surface will be easier to break and the blasting efficiency is improved.
Nowadays, instead of studying the rock blasting experimentally, various numerical studies also have implemented the rock fracture by blast due to the fact that the computing technologies are capable of carrying out largescale numerical calculation in a short time. In general, the numerical methods can be classified into three types, that is, continuum method, discontinuum method, and coupled or combined method. Many continuum methods have been employed in studying the rock blasting process, for example, ANSYS-LSDYNA [12], ABAQUS [13,14], LS-DYNA [15][16][17], and AUTOYN [18]. e continuum method can well model the stress wave propagation and stress reflection at the free surface. However, it faces difficulties in modelling the rock fragmentation and resultant muck-piling process. e discontinuum methods are also implemented in the rock blasting simulation, for example, UDEC [19], DDA [20], and DEM [21]. e discontinuum method can well model the rock fragmentation process, but it has a limitation in modelling the transition from continua to discontinua through fracture and fragmentation [22]. e emerging combined continuum-discontinuum method may be one of the best methods for modelling the entire rock blasting process since it combines the advantages of both continuum and discontinuum methods, and it overcomes the shortcomings of the both methods [23]. ere are several different types of combinations of continuum with discontinuum methods, such as the combined boundary element methodfinite element method (BEM-FEM), discrete element method-finite element method (DEM-FEM), discrete element method-boundary element method (DEM-BEM), and finite-discrete element method (FDEM) [24]. Among them, the FDEM might be the most widely used and further developed combined continuum-disctontinuum method. e FDEM method is initially developed by Munjiza for modelling the tensile fracture for concrete [25]. He then explained each component of the method and gave the open source, i.e., Y-code, for other researchers further developing the method conveniently [23,25,26]. After that, the FDEM method is further developed and implemented in modelling the fracture process of geomaterial under various loading and environmental conditions. Mohammadnejad et al. used a GPGPU-parallelized 3D FDEM method to model rock fracture process and analysed the effect of the meshes, loading rates, and specimen sizes [27]. Research at Alamos National Laboratory developed the Hybrid Optimization Software Suite (HOSS) as a hybrid multiphysics platform for the simulation of solid material behavior complemented with the latest technological enhancements for full fluidsolid interaction [28]. Yan et al. proposed a three-dimensional coupled thermomechanical model based on FDEM to model the thermal cracking [29]. Sun et al. proposed a novel thermal-mechanical coupling model which consists of a thermal part for the temperature field computation and the FDEM part for the crack evolution modelling [30]. e model is applied to investigate the thermal cracking process in functionally graded materials (FGMs) under different kinds of thermal loads. e modelled results indicates that the proposed method is useful for the fracture mechanics analysis and design of the FGMs structures. Farsi et al. used the FDEM method to simulate fracture propagation in fibrereinforced-concrete-(FRC-) lined tunnels [31]. e modelled results show the capabilities of FDEM method for better understanding of the fracture mechanics and crack propagation in FRC tunnels. Boyce et al. investigated the formation of mixed-mode fractures using the FDEM method, and the modelled results demonstrate that mixedmode fracture can be captured via numerical simulations [32].
In this research, a combined finite-discrete element method (FDEM) is proposed to model the rock fracture, fragment, and resultant muck-piling process in mining 2 Shock and Vibration production. e crater formation process by blast is modelled first to calibrate the proposed method. en the mining production process is modelled to demonstrate the potential application of the proposed method in the mining industry.

Combined Finite-Discrete Element Method
A combined finite-discrete element method (FDEM) used in the following sections has been developed by Liu et al. on the basis of their previous enriched finite element codes RT2D [33] and TunGeo3D [34] and the open-source combined finite-discrete element libraries Y2D and Y3D originally developed by Munjiza and Xiang et al., respectively. e combined method consists of the following components: contact detection, contact interaction between individual bodies, deformability and transition from continuum to discontinuum, temporal integration scheme, and computational fluid dynamics [23,[35][36][37][38]. Among them, the transition from continuum to discontinuum is the unique feature that makes the FDEM distinguish from continuum methods and discontinuum methods. us, the transition from continuum to discontinuum is introduced in detail. In addition, rock blasting process involves interaction between solid discrete elements and gas [23], that is, interaction between high-pressure gas and the blasting production. As it plays a significant role in modelling rock blasting, the gas pressure as surface load on the rock around the borehole and its flow into subsequently generated fractures are also introduced.

Transition from Continuum to Discontinuum.
e FDEM for modelling the sublevel caving induced by the blast can have one discrete element or a number of discrete elements with general shape and size. e interactions of the discrete elements are governed by contact law: a stiffness in the normal direction and a stiffness and friction angle in the tangential direction. For a single discrete element, it can be discretised into finite elements and simulated through the continuum law. e transition from continuum to discontinuum occurs through the separation of the finite elements according to the stress strengths of the material. For the mining production by blast, the entire rock mass is considered as a discrete element, which is then discretised into finite elements. e transitions from continuum to discontinuum occur through the fracture and fragmentation by blast, that is, separation of the finite elements.
In the FDEM, the finite elements in a single discrete element are bonded together through a bonding stress, and the fractures are assumed to coincide with element surface during fracture. e bonding stress is regarded as the function of the size of separation δ, and it can be divided into two components, that is, in the normal and tangential directions of the finite element surface, as follows: (1) where n and t are the unit vectors in the normal and tangential directions, respectively, of the surface at such a point and δ n and δ s are the magnitudes of the components of δ in the normal and tangential directions, respectively. Figure 2 demonstrates the relationship of the bonding stress and the separation of the finite element in pure normal direction for mode I fracture, that is, tensile fracture, and the relationship of the bonding stress and the slinging of the finite element along the element surface in pure tangential direction for mode II fracture, that is, shear fracture.
As illustrated in Figure 2(a), the bonding stress in the normal direction, that is, tensile stress, increases with the increase of the opening displacement of the bonded finite elements before it reaches the critical displacement, that is, δ np , prescribed by the tensile strength of the element, σ t . During this period, that is, δ n < δ np , no fracture occurs and the opening displacement of the bonded finite element will recover if the bonding stress disappears. en, as the opening displacement of finite element increases, that is, δ n ≥ δ np , the fracture is assumed to occur, and the bonding stress gradually decreases. When the opening displacement of the finite elements is beyond the ultimate opening displacement δ nu , tensile failure has completed. For the bonding stress in the pure mode fracture, it can be expressed by the opening displacements as follows: where D is a damage variable between 0 and 1, f(D) is damage function described in the mechanical damage model [39], and all other parameters have the same meanings as those introduced before. Figure 2(b) illustrates the relationship of the sliding displacement of the adjacent finite elements along the element surfaces with the bonding stress in the tangential direction, that is, shear stress. Before the sliding displacement of the adjacent element reaches a critical value, that is, δ sp , described by the shear strength of the element, σ s , no crack occurs, and the sliding of the adjacent element can recover if the tensile stress is removed. As the sliding displacement increases, that is, δ s ≥ δ sp , a shear crack occurs, and the tensile stress decreases gradually with the increase of the sliding displacement. As the sliding displacement of element continues to increase and eventually exceeds a residual opening displacement δ sr , the shear crack completes and bonding stress becomes a purely frictional resistance. e relationship of the sliding displacement and the boding stress in the tangential direction can be expressed as follows:

Shock and Vibration
where D is a damage variable between 0 and 1, g(D) is the damage function described in the mechanical damage model [39], ∅ f is the joint residual friction angle, and all other parameters have the same meanings as those introduced before.
In most of the cases, the rock fracture is produced by the combination of the shear stress and tensile stress. In that case, if the opening and sliding of the displacement of the adjacent finite element satisfy the following equation, the mixed-modes I-II fracture occurs.
In the FDEM, the rock fractures are governed by the energy release rate, instead of the rock strengths. For the tensile failure, that is, mode-I fracture, the fracture release rate G f I is equated to area under the curve of bonding stress and the opening displacement depicted by equation (5) and as illustrated in Figure 2(a). For the shear failure, that is, mode-II fracture, the fracture release rate G f II can be calculated using equation (6).
After the fracture and fragmentation, the integration of the motion of either the initial discrete elements or the newly formed discrete elements is conducted by central difference explicit time integration scheme. e scheme can be formulated as follows: where U, _ U, and € U are the displacement, velocity, and acceleration, respectively; F is the sum of the body forces, contact forces, and external loads together with any damping forces; M is the mass associated; and t is the current time, while Δt is the time increment.

Detonation-Induced Gas Expansion and Flow through the Fracturing Rock Mass.
In the sublevel caving, the detonation-induced gas plays a significant role for the rock fracture and fragmentation process, especially expanding the fractures and pushing the fragments around the boreholes. us, the detonation-induced gas should be taken into account while modelling the mining production process by blast. As the explosive is detonated, a shock wave coupled with highpressure gas is generated. en the crushed zone is produced, which lies in the vicinity around the borehole. e shock wave transmits into the rock and continues to break the rock mass. e detonation-induced gas involves interaction of the produced fractures and fragments and penetrates into the fracture voids to accelerate the fracturing process. In general, the detonation gas exerts pressure on the rock, causing the rock to damage, fracture, and fragment.
As most of the finite element methods, the Jones-Wilkins-Lee (JWL) equation of state (EOS) in equation (8) is used to model the interaction between detonation production and surrounding rocks:  Shock and Vibration where P is the instantaneous pressure at any time, A, B, R 1 , R 2 , and ω are the material constants, ρ 0 and ρ are the densities of the explosive and the detonation products, respectively, and E m0 is the specific energy.
After that, the gas flow through the fracturing model [23] is implemented in the FDEM to model the gas flow through the fractures in the rock by blast.
us, the spatial and temporal distribution of the gas pressure can be obtained according to the distance from the detonation point.

Calibration of the Combined Finite-Discrete
Element Method e FDEM has been calibrated by the uniaxial compressive strength and Brazilian tensile strength test under quasistatic and dynamic loading conditions in our previous research [40][41][42][43]. us, in this section, instead of modelling of the conventional compressive and tensile strength test used in laboratory, the FDEM is calibrated through modelling the crater formation process by blast. Test of the crater formation is widely used to study the blast phenomena, the behavior and destructive power of different explosives, and the response of the soils and rocks under this type of load [44]. us, the crater formation process is modelled using the FDEM and the modelled results are compared with those well documented in the literature to calibrate the proposed method. Figure 3 shows the geometrical model and the numerical model for the crater blasting. As illustrated in Figure 3(a), the crater blasting is simplified as a two-dimensional plane strain problem, and only the vertical direction is considered. A single borehole with a radius of 50 mm is prefabricated in the model and placed at 2 m far from the top surface. e boundaries of the other three surfaces are put far away from the borehole to eliminate the influence due to the reflected waves from the boundaries. As shown in Figure 3(b), the numerical model is discretised using triangular elements. e dense elements are employed in the interested area, that is, the vicinity of the borehole, to analyse the stress prorogation and crack initiation and propagation. e material properties of the rock specimen are Young's modulus E � 60 GPa, Poisson's ratio v � 0.26, density ρ � 2600 kgm − 3 , tensile strength t � 20 MPa, compressive strength c � 200 MPa, internal friction angle V � 30 o , surface friction coefficient u � 0.1, mode-I fracture energy release rate G fI � 50 Nm − 1 , and mode-II fracture energy relate rate G fII � 250 Nm − 1 . Figure 4 illustrates the temporal distribution of the minor principal stress wave propagation and the fracture initiation and propagation of the fracture during the crater blasting process. Following the solid mechanics regulations, the compressive stresses are taken as negative, while the tensile stresses are regarded as negative. In additional, in Figure 4(b), the tensile failure is marked using red color, while the compressive failure is marked in blue color. While the explosive in the borehole is detonated, a very rapid chemical reaction occurs and results in great heat and high-density gas [45]. en, the gas rapidly expands the surrounding rock mass and emits an intense pressure wave, that is, compressive stress, to the rock. As the pressure wave propagates radially out of the borehole (Figure 4 en the compressive stress continues to propagate and reaches the top free surface (Figure 4(b), at 0.55 ms). e compressive stress can be reflected and turn to be tensile stress. As the tensile strength of the rock is much smaller than the shear strength, the rock is easy to break by tensile stress. As can be seen, in the free surface, tensile failure is produced (Figure 4(b), at 0.55 mm). Meanwhile, the highpressure gas in the contained crushed zone penetrates the cracks and accelerates the crack propagation. us, long cracks are produced around the borehole (Figure 4(b), at 0.75 ms), which eventually reaches top surface (Figure 4(b), at 0.75 and 1 ms). e coalescence of the long radial cracks and the tensile cracks at free surface causes more fragments produced at the free surface (Figure 4(b), at 1.5 and 3 ms). Finally, the fragments in the crushed zone and the free surface are thrown out due to the high-pressure gas expansion (Figures 4(a) and 4(b), at 4.5 ms), and the blasting crater forms. e FDEM reproduces the rock fracture and fragmentation during the crater formation process by rock blast. e FDEM can well model the three zones, that is, crushed zone, cracked zone, and long radial cracks zone, as illustrated in Figure 5. Comparing modelled results with Figure 1, the modelled results agree well with those well documented in the literature [45].

Stress and Fracture Propagation and the Resultant Fragment Casting.
is research is trying to model a mining production process by blast to show the capabilities and potential application for the FDEM in mining production. e example is taken for a metal mine in Yunnan Province of China. e in situ leaching or solution mining technology is employed in the mine production. e mining technology can help to recover minerals, for example, copper and uranium from the fractured ore.
us, the ore does not need to be moved outside of the underground mine, which makes the mining production process more efficient and more economical. e initial process for the in situ leaching involves the drilling of holes into the ore deposit. Explosive or hydraulic fracturing Shock and Vibration can be used to create open pathway for the penetration for solution. In this mine, the blasting is employed to break the ore and finally result in fragments pilling. us, the solution efficiency can be improved for metal recovery. Figure 6 depicts the layout of a stope in underground mining in Yunnan Province of China. As illustrated in Figure 6, the dashed line indicates the current excavation area by blast. e drilling tunnel has been excavated in advance for providing working space for works. In addition, the open-off cuts have been excavated for every 9 m in length for providing swelling space for the blasted rock mass. As mentioned before, the in situ leaching or solution mining technology is employed in this mine. In addition, the rock fragments will be collected first, and then the solution technology for mine piles is taken for obtaining the minerals. us, the mining technology, that is, sublevel caving, is implemented to collect the mining fragments. Figure 7 gives a 3D model for better underlining the mining method of sublevel caving and provides the layout of the boreholes. e dimensional parameters for the boreholes correspond to those listed in Table 1. e combined finite element modelling of the dynamic rock fracture and fragmentation processes is considered as a plane strain problem, and B-B and C-C profiles are taken into account. Figure 8(a) depicts the sublevel caving used by underground mines for the C-C profile. e boreholes are placed in the C-C profile.
ere are a few C-C profiles depending on borehole spacing. e explosives in the borehole of the C-C profiles are detonated simultaneously, while they are detonated in sequence to recover the ore for boreholes in different profiles. e rock parameters used in this section are the same as those in Section 3. However, the acceleration in the vertical direction is set as 10 times gravity to save the computational time. Figure 8(b) shows the corresponding numerical model. It can be seen that the model is discretised using the triangle elements, while the interested area employs the dense meshes. Figure 9 shows the stress propagation, while Figure 10 illustrates the rock fracture and fragmentation processes. e explosives in the borehole in C-C profile are detonated from the bottom and prorogate upwards. As illustrated in Figure 9, the explosives are detonated at 0 ms; then the stress waves are produced and propagate from the bottom to the top of the boreholes (Figure 9, from 0 m to 0.5 m). At 1 ms, the stress waves reach the top of the boreholes. e stress waves continue to propagate outwards and interact with the rock mass ( Figure 9, from 1.5 ms to 15 ms). After 18 ms, the stress waves in the rock have attenuated and are eventually not able to break the rock anymore. However, the area initially with boreholes has been broken into fragments (Figure 9, at 18 ms). As time goes by, the fragments finally fall down. As illustrated in Figure 10, the fractures are produced immediately after the explosives are detonated, and the rock mass turns into fragments almost at the same time ( Figure 10, at 0.5 ms). With the explosions propagating upwards, the rock is broken into fragments, while the explosives pass way (Figure 10, from 0.5 to 15 ms). After 15 ms, the excavation area is full of fragments, which then fall down to the drilling tunnel ( Figure 10, from 18 to 100 ms). Figure 11 depicts the geometrical and numerical models for the sublevel caving used by underground mines for the B-B profile. e models only selected the section between two open-off cuts, which is used for providing swelling space for the excavated rock mass. e spaces between boreholes are 1 m. e explosives are detonated from the left borehole to the right borehole. e first four rings are detonated simultaneously, while the last four rings are detonated with 0.25 ms delay. Figure 12 illustrates the stress propagation, while Figure 13 shows the fragmentation process for the sublevel caving by blast in the B-B profile. e first four rings are detonated simultaneously, and the stress waves propagate upwards, as shown in Figure 12, at 0.2 ms. en the stress waves reach the top of the rings and begin to attenuate (Figure 12, from 0.4 to 1.2 ms). e rock mass with first four rings turned into fragments by blasts. en stress waves initiation of the last four rings is not very obvious, but the stress waves can be seen at 24 ms and 26 ms. Finally, the excavation area is completely broken into fragments. Figure 13 shows the same phenomenon as that in Figure 12 in terms of rock fracture process. e fractures first initiate from the bottom of the first four rings as they are detonated first. en the fractures propagate upwards and finally reach the top of the rings. us, the rock area with the first four rings turns into fragments ( Figure 13, at 1.2 ms). e last four rings experience the same process, and the rock mass is excavated finally (Figure 13, at 50 and 400 ms).

e Effect of the Row Spacing and the Initial Gas Pressure.
It is essential to choose the appropriate row spacing and amount of the explosives for the excavation by blast. As this is concerned to the economical and efficient blasting design in the mining industry, in this section, the modelled results  Shock and Vibration with different row spaces and initial gas pressures are compared to illustrate their relationships with blasting results. Figure 14(a) shows the modelled result with row spacing of 1 m and initial gas pressure of 0.5 GPa, while Figure 14(b) shows the modelled result with row spacing of 0.75 m and initial gas pressure of 0.5 GPa. It can be seen that, with the lower spacing, the rock is broken into smaller fragments. Figure 14(c) shows the modelled results with row spacing of 1 m and initial gas pressure of    0.5 GPa. Compared with Figure 14(a), it can be seen that, with the same row spacing, the lower initial gas pressure leads to much bigger fragments. us, the FDEM can well indicate the effect of the initial gas pressure for gas blasting results.  modelling the dynamic rock fracture and fragmentation process, the crater blasting process is modelled here using the finite element method. e finite element software ABAQUS is employed to model the process, while the same sizes of the model and the meshes for the numerical model are employed, as shown in Figure 3. All the rock parameters are the same as those mentioned in Section 3.

Comparing the
As illustrated in Figure 15, at 0.15 ms, after the explosive is detonated, the crushed zone is produced immediately, while the compressive stress propagates outward the borehole.
en, as the compressive stress continues to propagate, no more cracks can be observed obviously ( Figure 15, at 0.25 ms). As the compressive reaches the top free surface (Figure 15, at 0.35 ms), the stress is reflected to tensile strength. However, only cracks at the very top of the model are produced. Although the compressive stress and the reflected tensile strength continue to propagate, no more fractures are produced and the rock mass at the free surface area is not broken into fragments.
us, the throwing process from the free surface and the final muck-piling process are not observed.
Compared with the combined finite-discrete element modelled result (Figure 4), the finite element method can only well model the stress propagation process. Although it can model the rock fracture to some extent by removing the finite elements when stress meets the rock strengths, it is hard to model the rock fragmentation process and the resultant fragment casting. us, as the FDEM can model the entire blasting process, the potential application in mining industries is enormous.

Conclusion
In this study, the FDEM is implemented to model the rock fracture process and the resultant fragment muck-piling process by blast in an underground mine. To model entire rock fracture and fragmentation process by blast, the FDEM develops the transition from continuum to discontinuum thorough fracture and fragmentation, the detonation-induced gas expansion, and flow through the fracturing rock. en, the FDEM is implemented to model the crater formation process by blast to calibrate the proposed method. e proposed method not only can model the rock fracture process and resultant fragmentation but also can well capture the main characteristics of rock blasting in the single borehole, that is, formation of the crushed zone, cracked zone, and the radial long crack zone. e three modelled zones are compared with those well documented in literature and show good agreements. After that, the FDEM is employed, modelling the rock fracture process and the resultant fragment muck-piling process by blast during the mining production.
e proposed method has well modelled the stress and fracture propagation processes and the resultant fragmentation muck-piling process. Finally, the effects of the borehole row spacing and the initial gas pressure are discussed. e advantages of the FDEM in modelling the dynamic rock fracture and fragmentation by blast are discussed by comparison with the finite element modelled result. e results of this study conclude that e FDEM has well modelled the fracture and fragmentation process by blast in mining production process due to the implementation of the key component, that is, transition from continuum to discontinuum, in the proposed method.
e key component makes the FDEM take the advantage of both the continuum-based finite element method and discontinuum-based discrete element method and can free transition from the continuum modelling to the discontinuum modelling. e combined method is a valuable numerical approach to further study the dynamic behavior of rock by blast as it has well demonstrated that the proposed method not only can well model the entire rock blasting process, that is, fracture initiation, propagation, and resultant fragment casting process, but also can capture main blast-induced zones, that is, crushed zone, cracked zone, and the radial long crack zone.

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

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.