Simulation of FLY-ROCK DISTANCe as a Function of Blast Conditions: a Case Study in Vietnam

The prediction of fly-rock distances is a big problem in the blasting areas of open-pit mines. The fly-rock distance plays a crucial role in the reduction and control of blasting accidents in quarries. This paper proposes the modelling of the contact dynamics as a non-smooth discrete element method (NSCD). Then, the fly-rock phenomenon is simulated using this NSCD method. This approach was to develop a model that correlates to blast conditions, initial fly-rock/rock-fall velocities and permits the computation of fly-rock range from randomization orbits. The results showed that the NSCD method is a good means for the simulation with the variability of blasting parameters. This method is to relate the initial fly-rock velocity to blast conditions and then uses ballistic trajectories to compute the maximum fly-rock distance. Finally, it should be noted that the proposed simulation of fly-rock trajectories which shows the distribution of fly-rock bounce heights above the ground level with the different coefficient of restitution range as a function of blast conditions. It should be used in the Ninh Dan limestone quarry belonging to the Song Thao Cement, Phu Tho province (Vietnam), and it should be directly used in the same other quarries.


Introduction
Blasting is a common rock fragmentation technique utilized in mining operations such as quarrying. When a blasting agent is initiated, the chemical energy of is converted into mechanical energy, it is used in the actual breakage and displacement of the rock mass, and the rest of energy is spent in the spreading like ground vibration, the seismic effect of the blasting, fly-rocks, back breaks, etc. (Kutuzov, 1992). The blasting process, some stone pieces travels to a certain distance away from the blast area. Massive mine explosion sends rocks flying everywhere. This undesirable projection of stones is termed as fly-rock. Fly-rock is defined as the rock propelled beyond the blast area by the force of an explosion. According to previous researchers, controllable factors (burden, stemming, drilling parameters, and powder factor) and uncontrollable factors (geotechnical and geological conditions) should be used in predicting fly-rock since their effects on the occurrence of fly-rock (Monjezi, 2011). Fly-rock is the source of most of the injuries and property damage in the majority of blasting accidents in quarries (Bajpayee, 2004).
Rock-fall is a fragment of rock detached by sliding, topping, or talling, that fall along a vertical or sub-vertical cliff. Rock-fall is the natural downward motion, the movement of rock from the slope that is steep the rock continues to move down the slope which may be by free falling, bouncing, rolling or sliding depending on the slope angle. Ritchie (1963) identified the characteristics of rock-fall motion relative to a slope's configuration and height, to determine the expected impact distance of a rock-fall from the base of open-pit benches. Irrespective of a rock's shape, the rock's mode of travel down the slope is a function of the slope angle. Fig. 1 shows an empirical design table of the recommended minimum rock catch-ment area with and depth, based on the slope height and slope angle. Fig. 1, fly-rocks or rock-fall in trajectory seldom give a high bounce after impact. The degree of rock-fall promotion depends on the environmental factors causing physical and chemical weathering, and on the bedrock type (Dorren 2003).
There are many different models for calculating runout zones of rock-fall events such as empirical models, process-based models and GIS-based models (Geophysical Information System) (Dorren 2003;Pierson et al. 2001). Empirical rock-fall models are generally based on relationships between topographical factors and the length of the runout zone of one or more rock-fall events (Tianchi, 1983;Keylock and Dommas, 1999). Rock-falls can be a hazard for many facilities in mountainous terrain, particularly in areas with high rainfall and freezing temperatures, and seismic events occur (Wyllie 2014). The randomness of the rock-fall process has many sources. The stochastic model is capable of the simulation both the average behaviour and the variability of the process. (Preh et al. 2015).
Fly-rocks or rock-fall is one of the most hazardous events in blasting operations of quarries. The fly-rock range is as a function of blast conditions. After the rock has been detached and starts to move, it descends the slope in different modes of motion such as figure 1b, which depends on the mean slope gradient. In this paper, we have theoretical research to determine the parameters and influencing factors to simulate the orbit of fly-rock and falling rocks. By simulation model of fly-rock and rock-fall are used as a basis for of forecasting and proposed solutions to protected areas surrounding quarries. Then, we will discuss the models of fly-rock simulation. These models fly-fall simulations can be observed the geometric modelling of the boulder. A boulder is considered to be a Submission date: 06-03-2020 | Review date: 22-09-2020 point-like particle or a rigid homogenous sphere with a certain radius. A boulder may be a rectangular block with three parameters (height, width, depth) or be an arbitrary polyhedron. Finally, fly-rock simulation the models for rigid bodies with the contact interaction (rock and terrain with Signorini's law, Coulomb's law of dry friction). It is possible to classify according to the modelling of the contact interaction such as rebound impact models; non-smooth contact dynamics method (NSCD). The results of the application to investigate the verification examples present the simulation of fly-rock trajectories with the different coefficient of restitution in some quarries Vietnam.

Methodology
The algorithm of the general resolution used in the Non-Smooth Contact Dynamics method was developed in the laboratory LMGC with the works of Jean 1999 (Dubois & Jean 2006). Non-Smooth dynamics in the contact frame is based on a node by node resolution of the contact. The laws of Signorini for unilateral contact and of Coulomb for friction are not regularized (Radjai & Richefeu 2009). The NSCD method is distinguishable from the original smooth Discrete Elements method (DEM) (Cundall, 1980) due to the following features: an implicit scheme for integrating the time discretized dynamical equation; a non-regularized interaction law (Signorini unilateral contact and Coulomb dry friction); the possibility for finite element discretization in order to take into account the mechanical behaviour of rock blocks. This approach used was to develop a model that correlates blast conditions, initial fly-rock/rock-fall velocities, and permits the computation of fly-rock rang from randomization orbit. This method used in the present study is to relate initial fly-rock velocity to blast conditions and then use ballistic trajectories to compute the maximum fly-rock range.

The NSCD method
The NSCD method relies on a unique formulation of the mathematical and mechanical background, allowing us to deal with some extended kinds of laws. For the non-smooth in time, the occurrence of velocity jumps is a well know feature of the second-order dynamics with unilateral constraints on the position even with continuous media.
The contact forces between two bodies are bound by the principle of mutual actions. At any time of the evaluation of the system, one needs to define the interaction locus and an associated local frame in order to describe the interaction behaviour. It is assumed that one is able to define for each point (C) of the candidate boundary (unique) nearest point (A) on the antagonist boundary. It allows defining for each couple of points a local frame (t, n, s): with n the normal vector of the antagonist boundary and (s, t) two vectors of its tangential plane. The calculation of contact forces in the NSCD method is performed in two steps. The result of the interaction of the antagonist body (A) on the candidate body (C) can be considered equal to the force rα acting at the contacts point between these two bodies. The normal vector nα pointing from antagonist to candidate body two tangential vectors and, which define the tangential space by respecting this convention sα x tα = nα. On the other hand, we denote the gap distance between bodies in the normal direction. This value will be negative if there is interpenetration between the bodies.
By the definition of a linear mapping Hα that relates the local forces to the global forces, verifying the lolling equation: where Hα(q) is a mapping which contains the local information about contactors.
Finally, the global contact forces can be obtained by the relation: The same procedure is employed for the velocity calculation, and the velocity of the bodies can be expressed in the local frame.
The contact conditions are solved at the local level, and the impulse force is R and the relative velocity U. The dynamic equations are solved at the global level. The global impulse force is noted r, and the global velocity vector is noted. H and H* are the mapping operators to pass from local to global unknowns (Jean, 1999; Mereau, 1998, Dubois, 2011), Fig. 2.
In these many choices: punctual contact with extended law (transmission of torque); multi-punctual contacts with classical interaction laws; continuous surfaced description; etc. Less trivial in usual cases: not strictly convex (cubes, bricks, etc.); only locally convex (general polyhedron, triangulated surface); not convex at all, it may be decomposed in not strictly convex shapes. Only in simplest cases (rigid body with strictly convex boundary), the interaction locus may be considered as punctual.

Contact laws
Collision detection is a vast problem. The most popular approach is defining a contact force R = k*g at each contact point where k is stiffness and g is a measure of the interpenetration between a pair of objects. The value k depends on the nature of the objects, the type of interaction, and other elements of the simulation. No matter how k is chosen. Hence, when g is negative, which generate force. g is positive when the objects are no longer in contact.

Signorini condition
The NSCD of objects in contact has been extensively studied in the mechanical simulation of a jointed rock mass in contact. A contradiction with Signorini's condition of contact (Fig. 3a) which states that there exists a complementarity relation between the interpenetration distance g and the normal contact force RN at the point of contact.
Where g is the algebraic distance between two bodies at the point of contact and R N is the amplitude of normal force needed to solve the contacts. In the case of frictional sliding, a tangential component R T is introduced, leading to a contact force R = (R N , R T ). For the dynamical problem, it is more natural to formulate the unilateral contact in term of velocities with assuming g(t o ) ≥ 0 then t > t o if g(t) ≤ 0 then U N ≥0, RN ≥ 0, U N *R N =0 else R N =0.

Contact law of Coulomb
Following the contact law of Coulomb (Fig. 3b), the sliding at the contact point only appears when the tangential component of the contact forces between two objects is larger than the sliding threshold. This condition is given as: (4) With U = (U N , U T ) is the sliding velocity vector between one object to another. μ = tan (φ) is friction coefficient and φ is the internal friction angle. It is existed a parameter ρ > 0 that U T = -ρ.R T with frictional contact (Signorini-Coulomb).

Impact laws (Newton restitution)
Signorini unilateral and Coulomb dry friction are described and proposed the force laws for the non-impulsive contact force. The laws mean a non-regularized interaction law.
When shocks occur in a rigid body collection, the equation of motion and the interaction law are not sufficient to describe all the physics of the problem properly. Therefore, the local phenomena as the inelastic behaviour of materials at the interface and the global phenomena as the wave propagation in the body bulk should be taken into account long-distance effects due to simultaneous impact.
In the case of binary shocks, three kinds of restitution can be used as law Newton restitution, Poisson restitution or Energy or Strong restitution. Newton restitution relates the velocity after (U+) to the velocity before (U-) impact. Poisson restitution, which relates restitution impulsion (Rr) to the compression impulsion (Rc) according to the decomposition of the shocks in a compression and restitution phase. Energy or Strong restitution, which relates restitution energy to compression energy. With hard unilateral constraints, impulsive contact forces arise whenever contact closes with a negative relative velocity U-< 0. There is a jump in the generalized velocities U from the impulsive contact forces with that ve-locity in the normal direction is nonnegative U N ≥0. The velocity jump is called an impact law. We use the kinematic of coefficients of restitution (en) such as Newton restitution, is defined as: The kinematic of coefficients of restitution was analyzed to the normal and tangential components with respect to the plan surface, defining the normal (e n ) and the tangential (e t ) coefficients of restitution. Impact law for the normal component is well understood; it is not the case of the tangential component (for the rock terrain interaction as such elastic effects are absent). e n [0, 1] is the coefficient of restitution for the normal direction. The case with e n = 1 corresponds to a reflection of the normal velocity, whereas the smaller than one additionally dissipates energy. The choice of the restitution coefficient is a difficult task for complex structures.

The rock geometry
The rock geometry is considered to be 2D/3D convex polyhedron and can, therefore, be defined as the convex hull of a finite point sets in the space. Having many choices such as punctual contact with extended law (transmission of torque); multi-punctual contacts with classical interaction laws; continuous surfaced description; etc. Less trivial in usual cases: not strictly convex (cubes, bricks, etc.); only locally convex (general polyhedron, triangulated surface); not convex at all, it may be decomposed in not strictly convex shapes (Dubois 2012). Only in simplest cases (rigid body with strictly convex boundary), the interaction locus may be considered as punctual, table 1.

Fly-rock simulations with blast conditions
The approach used in the present study is to relate initial fly-rock velocity to blast conditions and then use ballistic tra-  Tab. 2. Odpowiednie parametry strzałowe dla kamieniołomu wapienia Ninhdan jectories to compute the maximum fly-rock range. This approach his entirely justified because the effects of air friction are quite small for typical fly-rock sizes and velocities. The maximum fly-rock range can define a "blasting area" in which no personnel or equipment should be present during a blast. Thus, determination of initial fly-rock velocity completely determines the maximum fly-rock range.
We need to know the information about the true burden and knowledge of the joint sets before planning our blast. There are some parameters important for estimating the fly-rock such as geological conditions (rock types, joints sets, faults, crushed zones, points, filled joint) will with a great deal of confidence, give uncontrolled results like fly-rock; design and performance (burden, especially in the first row, collaring, orientation and drilling, borehole deviation measurements, total explosive load per borehole, borehole diameter, powder factors, timing-delays between rows, stemming, etc.). Breaking burden and breakout occurs only at the "vertical" free face in the region of length are caused the fly-rock due to blasting, Fig. 3.
The Gurney equation (1943) successfully predicts initial velocities (v 0 ) of metal plates and metal fragments propelled by explosives. Consequently, it is logical to attempt to adapt the Gurney approach to the determination of initial velocities of rocks propelled by explosives to fly-rock velocities (v 0 ) obtained in bench blasting. (6) where √2E is Gurney constant with the characteristic of the explosive used; a form of the function (f) depends on the geometry of the system; c and m are the masses of explosive and material, respectively; (c/m) such as powder factors.
Where: W/l is explosive weight per unit length of borehole; l is explosive column length; ρ m is the density of the rock; b is a burden to the free face; b 2 btan(α⁄2) is break- Fig. 6. Results of the simulation of rock-falls trajectories from the top crest source area in Mathematica which shows the distributions of rock-fall bounce heights above the ground level with the coefficient of restitution 0.6 (x-axis is measured fly-rock distance, m; y-axis is fly-rock heights above the ground level with shallow or deep benches) Rys. 6. Wyniki symulacji trajektorii obrywów skalnych z obszaru źródłowego szczytu wierzchołkowego w Mathematica, na której przedstawiono rozkłady wysokości odrzutów skał nad poziomem gruntu ze współczynnikiem restytucji 0,6 (oś x jest mierzona odległością latających skał, m; oś y to wysokości latających skał nad poziomem ziemi z płytkimi lub głębokimi wyrobiskami) out width at the free face; α is breakout angles in bench blasting.
For most explosives √(2E' )=D/3 where D is the detonation velocity of explosive. For ANFO, which is the explosive used in most surface mine blasts: √(2E' )= 0.44D.
We can predict initial velocities (v 0 ) of rocks propelled by explosives to fly-rock velocities (v 0 ) obtained in bench blasting. We can determiner or measured fly-rock velocities such as: Grannite Bench blast: Dolomite and Limestone Bench blast: The powder factor is estimated with: Powder Factor ≈ 0.5ρ m (c/m) is for shallow benches, Powder Factor ≈ 0.8ρ m (c/m) is for deep benches.

Study areas
The limestone quarry has been studied through this research. The Ninh Dan limestone quarry belongs to the Song Thao Cement. It locates in the Phu Tho province. The quarries are one of the large production cement companies of Vietnam with an output of over 2 million tons per year. Discontinuity surveys were performed on 5 benches of a limestone quarry situated in the Nord of Vietnam, 70 km North of Hanoi. 55 planes were collected, which is a small sample size, but sufficient to identify sets of discontinuities. The studied slope has benches from 10 to 15 meters height with a dip of 70° and an azimuth varying from N100 to N120 (Fig. 5). Mechanical and physical properties of the rock mass in the quarry Ninhdan: for the rock mass characteristics, the friction coefficient equals 0.5, and the rock mass density equals 2.7 g/cm3. The in-situ friction coefficient for this type of rock is approximately 0.65, but in the simulation presented in this study, the value 0.5 was considered for safety reasons. Rigid type elements were used for the boundaries of the models and the dynamic simulation.
Based on the geologic conditions, properties of rock, the current state of blasting operation, used equipment of the Ninh Dan limestone quarry: Suitable explosives consist of ANFO (water-resistant and no water-resistant types). They can be loaded under individual and combined forms; Appropriate blasting parameters are calculated for all kinds of rocks, different drill hole diameters, and bench heights, various types of explosives, different bucket capacities of excavators. The appropriate blasting parameters calculated for specific conditions at Ninh Dan limestone quarry as follows in Tab. 2.
Comparison of observed and computed fly-rock velocities in limestone bench blast with explosive consist ANFO, spall velocities in limestone, initial fly-rock velocities calculated are minimum 2 m/sec and maximum 3 m/sec. The model is given by the Mathematica code.

Modelling of fly-rock
An assessment of rock-fall was undertaken to identify the potential fate of blocks that may be detached from unstable pit lope.
faces and therefore, to improve the safety on the openpit mines. The simulation method proposed in this paper with 2D and 3D simulation model for rock-fall based on the point-mass or the rigid homogenous sphere or the rectangular block or be an arbitrary polyhedron on the non-smooth contact dynamics method. The analysis was completed using Mathematica (Fig. 6). and NSCD method in LMGC90, program that simulates the trajectories of rocks falling from the slope (Nguyen et al. 2015). We have the fly-rock simulations with Mathematica code and the LMGC90 (general-purpose and open-source software).

Results and discussion
The first in the environment Mathematica (Verdel 1997), a trajectory is modelled as a two-dimensional rock-fall simulation largely based on the slope geometry (point-mass 200kg, velocity along with X 2 m/s, and velocity along with Y -1 m/s). By using statistical analyses, the method includes the computation of the probable trajectories, energy, velocity and bounce height envelopes for individual rock blocks. The sloping cut can be modelled by the program, that the ultimate resting locations of rock-falls can be determined and the results graphed with comprehensive statistics calculated. The program requires input data such as the slope roughness, restitution of normal and tangential rock energy, coefficient of rolling friction and rock mechanical parameters.
Besides the uncertainty in the definition of the input parameters, stochastic variability is used Figs. 7 and 8. The model is given by the mathematical code in Tab. 3.
Secondly, some models simulate the movement at the slope surface during rock-fall with detailed characterizations for bouncing, sliding and rolling. Code LMGC90, which is a general-purpose and open-source software, is capable of modelling a large collection of deformable or undeformable particles of various shapes with various interaction laws. Models have used Coulomb's law of friction for calculating rolling and sliding velocities. We developed the models using both a tangential and normal coefficient for the efficiency of collision.
Boulder rock-falls usually generate faster movement as longer runout distances. A boulder was recorded to have fallen from the potentially unstable slope blocks and travelled away from the toe of benches. It must be identified for the prevention of impacts on people and assets from rock-falls in open-pit mines. These include flexible rock-fall barriers, identifying effective berm widths, bunds constructed on production berms and using draped mesh. We can be observed in horizontal locations for maximum trajectories of falling blocks are shown in Fig. 8.
Fly-rocks are considered as significant risk security in quarries. Blocks falling from high up on benches can travel into the pit floor may destroy mining infrastructure as present a serious safety hazard for mine personnel. In the present paper, we have been used the simulation techniques of nonsmooth contact dynamics with stochastic approaches and adapted for the efficient simulation of rock-fall trajectories. The terrain geometry is constructed by section 2D and 3D.

Conclusion
Fly-rock prediction is a complex issue in the mining industry because at first, many parameters influence fly-rock phenomenon that can be divided into two categories: control-lable and uncontrollable. Most of these parameters accompany uncertainty due to variability in blasting parameters. In this study, the presented method allows for the 2D, 3D simulation of rock-fall events on arbitrary slopes and with arbitrary rock geometries. The influence of shape on the rolling behaviour of bodies can be studied with the presented simulation method in the environmental Mathematica and the code LMGC90. The models show that a barrier located on the pit floor and would be required to prevent falling rocks reaching the mining operation zone safely. The method with this model is constituted of major controllable blasting parameters such as burden, stemming, powder factor, delay interval between adjacent rows, and time delay respectively at the bottom and total explosive.
In this paper, the obtained results from these methods are the powder factor, mean charge per blast hole and blast hole length, whereas based correlation sensitivity these parameters are powder factor, stemming, and burden. The approach used was to develop a model that correlates blast conditions and initial fly-rock velocities, and permits the computation of flyrock rang from randomization orbit and radius as explosively-propelled fly-rock from vertical rock faces or benchtops; measured fly-rock velocities from mining and explosives; the model indicates that for fly-rock. Results of the simulation of fly-rock trajectories from the top crest source area, which shows the distributions of rock-fall bounce heights above the ground level with the coefficient of restitution 0.6. The simulation of fly-rock trajectories from the top crest source area, which shows the distributions of fly-rock bounce heights above the ground level with shallow or deep benches.
It means that in order to develop the fly-rock trajectories by the randomization orbit and radius based on the collected data in the Ninh Dan limestone quarry belongs to the Song Thao Cement, Phu Tho province. Consequently, this indirect relationship is observed in the results of sensitivity analysis.