INFLUENCE OF BOUNDARY CONDITIONS IN DEM MODELS OF SUBLEVEL CAVING ON DILUTION AND RECOVERY

Discrete element modelling is common approach for analysis of dynamical processes that include granular materials. Method is commonly applied for the modeling of the gravity flow of the broken rock in caving operations. Herein, comparison of two different approaches in such modelling was presented. Two models, ring wide and symmetrically divided, were used for comparison. Ore dilution and recovery were monitored as main comparison parameters, but also the processing time of such models and spatial shape of extraction zones. Ring wide model show greater ore dilution and lower recovery, while half wide models show the opposite results. Cumulative dilution differs for 10%, while ore recoveries differ in approximately 15% for same extraction levels. Also, dilution entrance points are different, where ring wide model shows earlier dilution entrance at about 15% of extraction, while dilution enters at 25% in half wide models. Processing speed is increased in half wide model due to the reduction in number of elements, but overall processing time is nearly the same due to the frequent hang ups inside the half wide model. Shape and size of the extraction zones are very different. Half wide model shows much higher extraction zone in contrary with ring wide, while depth of extraction zone is greater in ring wide models which corresponds with higher and earlier entry of dilution.


INTRODUCTION
Sublevel caving (SLC) is underground mining method in which ore is blasted for the sublevel drifts using explosive charges in overhead ring (fan) disposition while hanging wall is subjected to spontaneous caving under the influence of gravity.Loading of ore causes that both ore and caved waste move towards the sublevel drift which is referred as gravity flow.At the beginning of the process pure ore is loaded and as extraction of one ring continues waste entry occurs.At certain point of time dilution of the ore/metal becomes too high that further extraction is not economically justified.High dilution is main disadvantage of the SLC mining method while ore recovery is directly related to the acceptable dilution level which usually varies between 20-40%.
There are many factors that influence the formation of the gravity flow in SLC and therefore dilution and recovery, of those, most dominant are blasting and hanging wall rock mass properties.
Fragmentation of blasted ore has high impact on SLC gravity flow and it is often identified as main influencing factor for ore recovery and dilution.Specific blasting pattern (fan) influences uneven spatial fragments size distribution and its well known that finer particles are moving more easily than coarser ones (Bridgwater, et al., 1978).In places where explosive charges are closer to each other, in plane of the boreholes, finer fragments are generated and vice versa.Another size differentiation occurs along the axis of sublevel drift, where finer fragments are generated closer to the boreholes and their size is increased with the distance from boreholes.Blasting and fragmentation impact has been extensively investigated (Kvapil, 1982;Stazhevskii, 1996;Selldén & Pierce, 2004;DeGagne & others, 2005;Rustan, 2000).Rustan (2000) has identified fragments size distribution, fragment shape and friction, density, cohesion and shear strength as main influencing parameters on gravity flow.
Beside ore properties, caved waste material has high impact on dilution of the ore.Influence of strength, friction and fragmentation change in caved waste material has been numerically investigated (DeGagne, et al., 2006;Lapčević & Torbica, 2017) while there is lack of in-situ reports on this subject.In general, caved material with lower friction angle will move faster, with lower frictional kinematic resistance, and will influence the greater dilution of ore.This is directly related to the roughness of waste rock blocks that depends on the actual geological setup and distance that caved material has traveled.During its gravity flow, waste material is being subjected to the change of its size and roughness of its fragments, while blasted ore has more or less constant properties everywhere in the mine.This change influences increasing of dilution as SLC production advances to the deeper parts of the mine.
Modelling of SLC gravity flow has been performed by many researchers over the years in order to investigate how different parameters influence its formation and behavior.Kvapil (1965) proposed the well-known ellipsoid theory that was mathematical formulation based on practical experiments of granular material flow.Kvapil conducted his research mainly on bunkers, but without the question this had major impact on SLC design and investigation of gravity flow.Chen (1997) and Gustafsson (1998) used stochastic method approach for gravity flow modelling.This method assumes stochastic nature of gravity flow where downward propagation of rock particles and upward propagation of void have certain probabilities.Plastic theory (Nedderman, 2005) has shown to be more applicable for analysis of material flow in hoppers.It is based on determination of particle velocities from stress distribution and provides best results in evenly grained materials.Discrete element method (DEM) has been successfully used for numerous simulations of granular material flow and as such for sublevel caving (Selldén & Pierce, 2004;DeGagne & others, 2005;DeGagne, et al., 2006;Lapčević & Torbica, 2017).DEM models provide very useful tool for gravity flow modelling since it is possible to investigate influence of different parameters that are normally hard to analyze in physical tests.Cellular automata (Castro, et al., 2009) is simplified approach similar to DEM and was created mainly to overcome high computational demands.Having all this considered, several modifications of SLC mining method were suggested (Milić, et al., 2014;Kacmar, 2008).

Discrete element method
Discrete element method is numerical method commonly used for modeling of discrete dynamic processes.It was originally developed for problems in rock mechanics by Cundall (1971).This method assumes discretization of space domain into spherical, polyhedral or arbitrary shaped discrete elements and their movement by the Newton's laws of motion.Discrete elements can be bonded together to represent complex shapes and their behavior in dynamic processes.Cundall proposed the concept of the method which is also referred as Distinct Element Method, while Williams proposed general formulation for 3D DEM (Williams, 1985) and Shi (1992) came up with Discontinuous Deformation Analysis.Munjiza (2004) proposed the combined Finite-Discrete element method that couples significant advantages of both methods.
Application of the DEM is very wide and beside the rock engineering it is used for modelling of problems that include granular materials and fluids, powders and even food, chemical reactions and astrophysical processes.Many commercial and open sources DEM codes are available like PFC (Itasca, 2008), YADE (Šmilauer, 2015), LIGGGHTS (Kloss, et al., 2012) and many others.Each of the codes has its own specifics and decision on which one should be used depends on the problem that is being modelled.

DEM model setup
Herein, YADE code has been used for the creation of the DEM models.This code is open source and has large community of developers and users, while its application has been shown successfully in numerous examples.Some of applications include slope stability (Bonilla-Sierra, et al., 2015;Scholtès & Donzé, 2012), soil mechanics (Tran, et al., 2014), rockfall protection systems (Thoeni, et al., 2014), stress distribution around underground openings (Tran, et al., 2012) and many other.Results that are obtained from DEM simulations can be much more reliable than those obtained from other methods.Beside its numerous advantages, there are also disadvantages of the method, and one of main disadvantages is long computational time.Due to this problem, it is common practice that models are reduced so they consist of smaller number of elements.Improper reductions may have impact on the outcome, and this is going to be explained in this paper for the case of modeling of sublevel caving.
Modeling procedure consists of creation of external boundaries of model inside of which discrete elements will be distributed.External boundary can be created in CAD software and exported in STL or GTS format (other formats are also supported) and then imported into the model.External boundaries are represented as rigid surfaces with mechanical properties that define their interaction with the rest of the model segments.In this case, external boundary represents the end on one sublevel drift where blasted ore is loaded, and it is large enough to encompass the close perimeter of the rock mass (caved rock) that has the influence on the mining process.Two different models are created for the purpose of this research, one that is wide as one ring and the second that is symmetrically divided along the axis of the sublevel drift using the surface elements.Void space inside of the external boundary is filled with two types of spherical discrete elements, one that represent the ore, and other that represent the waste material.Half-wide model contains only the half of the number of elements present in ring wide model.Such reduction of the half-wide model should reduce the time of the model convergence.Layout of models can be seen in Figure 1.Ore and waste elements have different mechanical properties (Table 1), as well as different spatial size distributions.Sublevel caving mining method is well known by application of fan blasting.In this case blasted ore has uneven fragment size distribution in space.At places where explosive charges are closer to each other fragments are finer, while at places where explosive charges are more distant fragments are coarser.This is true for the plane in which boreholes are, but also there is fragments size differentiation along the axis of the sublevel drift.Finer fragments are closer to the explosive charge and coarser fragments are formed away from the boreholes.This is illustrated in Figure 2 for case of models used in this research.Draw of two subsequent fans is simulated.
Waste size distribution depends on the size of primary blocks in the rock mass and size distribution used in this case is shown in Figure 3. where D and R are ore dilution and recovery respectively, mw and mo are extracted masses of waste and ore respectively and Mo is the total mass of the blasted ore.
In order to obtain comparable results from different models it is necessary to have same convergence criteria for each of them.In this case, models are processed until 120% of extraction is reached.This means that total amount of extracted material is equal to 120% of total blasted ore mass.
During the execution of the model hangup or jamming of the draw point may occur as it is also seen in practice.In order to avoid jamming of the model additional routine is called at each 10000 iterations of the model.This routine applies zero friction to the spherical elements that are found in close perimeter of the draw point.After 100 iterations frictional properties are restored to original values.In this manner, only small number of elements are being frictionless for very small time period so consistency of the models is preserved.Jamming of draw points in simulation influences large time wasting that may occur between two calls of routine.More frequent calling of such routine would have possibly negative impact on simulation outcome since frictionless elements in gravity flow would move faster than those with friction.This would inevitably disturb the gravity flow and the results.
All DEM simulations require definition of time step which represents time increment used in each iteration.Time step needs to be sufficiently small in order to keep simulation stable and its value depends on number and radii of particles.In some cases, time step has constant value since nature of simulation is that no significant change to the stable time step occur.However, simulations where number of particles is changed during the simulation may suffer from such constant value.Decreasing number of particles may allow time step increasing and vice versa.YADE code has built-in routine that can be called after certain number of iterations to adjust the time step to the current state of the simulation based on its global stiffness.This routine is explained YADE documentation (Šmilauer, 2015).Since sublevel caving simulation assumes constant change of number of particles this method for time step calculation is adopted.This provides same conditions for each simulation and makes their execution optimized.

RESULTS
YADE DEM code provides options to save results in various formats so different post processing techniques can be applied.One way is to process results using Paraview software (Henderson, et al., 2004) that has been developed specially for post processing of different data sets.YADE has built-in routine that saves geometrical and physical data, so it can be post processed.Data is saved after defined number of iterations and if this number is larger than one, only one set of discrete stages of simulation are saved.However, having in mind that simulations may have more than several millions of iterations, saving data on each hundred or few hundred iteration have no significant impact on later post processing of results.Another way of representation is using plots of calculation parameters, which is done using GnuPlot (Williams, et al., 2017).
Ore recovery and dilution are main parameters of interest in this simulation and herein focus is given of difference between cumulative recovery and dilution obtained from two simulations.Ore dilution from both models is presented in Figure 4.As it can be seen there is obvious difference between cumulative dilution curves, where dilution rate in ring wide model reaches above 35% of cumulative dilution while half wide model reaches 25%.Another thing that may be noted is that shape of the dilution curves is different.In ring wide model dilution curve start at about 15% of extraction and has steep inclination up to the 60% of extraction, and from this point curve tends to converge at the 40% of dilution.On the other side, dilution curve from the half wide model start at about 25% of extraction and is almost linear until its end.However, difference between dilution is from 10-15% depending on the actual extraction rate that is observed, which is, in terms of ore and waste quantities, great difference.In terms of the achieved dilution rates in production both scenarios are possible depending on the actual mining conditions.rough and their gravitational movement would show somehow different behavior that in those models.It could be concluded that half-wide model may overcome or suspend the reduction of fragment shape, but this is not investigated in this research.Limitations that are related to the half-wide model used here is that symmetry plane reduces the draw point size and coarser waste particle could be impossible to pass.This arises additional question, should be particle size of waste modified and adjusted to the size of draw point size?
Figure 6 Comparison of extraction volumes in two models at 90% and 120% of extraction Recovery of ore is important factor in sublevel caving and its estimation is not easy task.
Figure 7 illustrates the recovery of ore for DEM models used in this research.From here it can be seen that curves are same until the extraction of 20% and from this point on difference between two curves is rising.At the final stage, when extraction rate reaches 120%, difference between two curves is around 15%. High recovery of ore in half-wide model corresponds to the height of the extraction zone.In practice, both scenarios are possible and can be acceptable if amount of dilution is in its acceptable limits.Main reason to use the half wide models is reduction of the processing time of such model.However, this expectation was not met in this case.In both cases models converged in about 80 hours, even if half wide model had higher average number of iterations per second.Reason for such timings is that in half wide model very frequent hang ups of particles occur after the entrance of the coarser waste material (Figure 8).Symmetrical division surface in half wide model fundamentally changes the boundary conditions and reduces the cross section of extraction point in that model.This means that such model could not simulated the caving of the coarser waste since it would be physically impossible for such particles to pass to the sublevel drift.

CONCLUSION
Numerical modelling is common approach in solving of engineering tasks.Depending on the scale of the phenomena that is modeled and the actual numerical method, depends the computational cost of such analysis.Using DEM for the modeling of sublevel caving it is possible to estimate the outcome of the production for the known mining conditions.Herein, comparison of two different modelling approaches is presented.
Comparing the two modeling approaches in DEM modelling of the sublevel caving process it was shown that fundamental differences exists.Different boundary conditions in half wide model are the main reasons for the existence of such differences.In this manner cross section of the drawing point is reduced and coarser waste material than material that is used in this example would physical be unable to pass to the sublevel drift.Also, this impacts the geometrical shape of the extraction zone making it possible to achieve the ore recovery rate of almost 90% in the model which is uncommon for the sublevel applications.On the other side, cumulative dilution curves are different in their shape and final resulting values show great difference (10-15%).
It was experienced that reduction in number of discrete elements increased the processing speed, but cumulative processing time remained the same due to the frequent hang ups in the half wide model.
As a conclusion it should be mentioned that all reductions of the real system will impact the discrepancies of the real world and model results.Modeling of the sublevel caving using DEM should be based on the reconstruction of the rock mass state prior to the caving.This means that initial model should be created and processed.State of the rock mass after the processing of such model should be used for the constitution of next (adjacent) model in order to recreate the representative initial conditions inside the model.Scale of the model should be as large as necessary to eliminate the boundary effect.Therefore, using the appropriate modelling approach one should tend to achieve more realistic estimate of the mining outcome.

Figure 1
Figure 1 Geometry used for two DEM models a) whole ring wide model b) symmetrically divided model containing only a half of elements

Figure 2 Figure 3
Figure 2 Particle size distribution for the ore packing

Figure 4 Figure 5
Figure 4 Cumulative dilution comparison

Figure 7
Figure 7 Ore recovery comparison for two models

Figure 8
Figure 8 Hang up of material in half wide model

Table 1
Material properties used in DEM models