The Use of Simulation Software using the Discrete Element Method (DEM) for the Process of Materials Comminution

. Comminution processes are one of the most common processes for processing energy materials, e.g. coal, biomass, and post-recycling elements. The hitherto unsolved problem is the high energy consumption of machines and the lack of precise descriptions of the phenomenon of comminution in terms of the relationship between the design features of mills and the properties of comminuted materials. The dynamic development of simulation techniques based on advanced models and the method of discrete elements allows for a certain mapping of occurring phenomena. The purpose of the work is to illustrate the possibility of using simulation software based on the discrete element method to model the grinding processes in the shredders grinding assemblies. The paper presents aspects of modeling the shape and size of particles, their interactions and contacts with mills structural elements, as well as aspects of crushing modeling in RockyDem software.


Introduction
The need for modeling of shredding and flow of bulk materials is mainly due to machine problems occurring during their processing: loss of ability to move (transport) the bulk material by blockage of dosing channels, transporting in feed supply systems, processing processes, too high intermolecular friction, which causes the effect of material clogging, which in consequence leads to material losses or deterioration of the quality of processed, transported material, energy start, decrease in efficiency and effectiveness of processing lines, and thus increase costs [1][2][3][4][5]. Therefore, it seems rational to model and search for structural solutions of machines working with loose materials, taking into account the modeling of the behaviour of these materials in contact with the working surface of these devices. The issues of modelling in mechanical engineering have been reflected in many works, e.g. in the field of surface condition modelling, thermal stability and surface properties and fatigue damage analysis [6][7][8][9].
The main goal of optimization of machinery and equipment design for grinding of bulk materials is to minimize costs, reduce energy consumption, while maintaining and even improving the performance and functionality of these units [10][11][12][13][14][15][16][17][18][19]. Improvement of these parameters is connected with appropriate selection of constructional features. For this purpose computer support methods are used, e.g. [20][21][22][23]. The commonly used analytical methods make it possible to determine the geometric and strength characteristics of machines, however, they do not take into account the behaviour of loose materials during shredding [24][25][26][27].The appropriate characteristics of the materials, especially in terms of grain size, have a great influence on the quality of the final products in many branches, e.g. in the medical industry, in the manufacture of artificial organs or devices supporting basic activities [28][29][30][31][32]. Computer applications based on solving problems with the use of numerical methods become helpful in this respect. One of the methods used to solve calculation problems of grinding and flow of bulk materials is the Discrete Elements Method (DEM), which allows to simulate the behavior of materials under machine conditions [1,2,33]. The Discrete Elements Method has recently been one of the most frequently used methods in the analysis and simulation of the flow of bulk, granular materials, supporting the optimization of machine and equipment design. Much attention is paid to issues related to modeling of mixing processes [34,35], storage [1] and transport [36] of granular materials and grinding, mainly of brittle materials, e.g. rocks and minerals [37][38][39]. The topic of modelling granular materials, e.g. cereal grains, is much less frequently addressed [37]. Zeng at al. [37] performed a simulation of rice grain cracking using the method of discrete elements, examining the influence of impact velocity and humidity on the crack course. Their considerations concerned single grains, which in industrial applications does not take into account the interactions between particles during crushing. Karwat et al. [2] indicated that in modeling the flow of bulk materials it is important to correctly calibrate input parameters on a macro scale, i.e. it is more important to map the behaviour of the entire deposit than the dependence of individual particles. Boikov et al. [39] also pointed out the importance of identifying the material behavior in the macro scale. They noted that in the case of spherical particles it is necessary to determine 7 input parameters of the simulation, which, if they change in a linear way, can be described by means of an approximation function containing these parameters. The development of a reliable simulation model, as close as possible to the actual behavior of bulk materials, requires correct calibration of the model input parameters, e.g. particle shape and size, coefficients of friction and restitution or bulk density, which are determined experimentally in laboratory tests [2,35,38,39]. The purpose of the work is to illustrate the possibility of using simulation soft-ware based on the discrete element method to model the grinding processes in the shredders grinding assemblies. The paper presents aspects of modeling the shape and size of particles, their interactions and contacts with mills structural elements, as well as aspects of crushing modeling in RockyDem software.

Fundamentals for DEM simulation of comminution
Over the years, the method of discrete elements originally developed by Cundall and Strack [40,41] has developed dynamically and is used not only for modeling bulk material flows, but also for shredding modeling [35,37,42]. The first DEM models for rock crushing were developed in the early 21st century. Potyondy and Cundall [43] developed and introduced a particle-bonded model for crushing processes, which assumed that single particles consist of spheres glued together. This model was developed and modified by Cho et al. [44], Tan et al. [45], Ivars et al. [46]. A different strategy was introduced by Cleary et al. [47], according to which spherical particles are replaced by densely packed smaller spheres (constituting offspring) after exceeding the value of energy corresponding to the fracture, which made it possible to map crushing processes [47,48]. This method has been extended so that it is possible to model non-spherical particles [49]. DEM is a numerical method that treats particles individually and motion equations are solved at the level of single particles, which has made DEM applicable to modeling discontinuous centers [48]. In order to predict the behavior of a particle, Newtonian equations of motion and inter-particle contact models are solved [48]. Some of them will be presented in the next section.

Possibilities of simulation software
Many commercial applications using DEM are currently available, e.g. RockyDEM, EDEM, Newton, ThreeParticle/CAE, PASSAGE®/DEM Software, Bulk Flow Analyst or open-open source codes such as LIGGGHTS (LAMMPS improved for general granular and granular heat transfer simulations). The possibilities of computer modelling of granular materials using the Discrete Elements Method (DEM) will be presented on the example of the commercial RockyDEM software, which enables to model the movement of bulk materials in machines and devices such as: belt conveyors, shredders, silos. The software also allows to model particles of any shape and model their motion on the basis of determining the type of contacts, coefficients of friction and restitution and to shred non-spherical particles. In RockyDEM, the calculations are mainly based on Newton's laws and gravity ( Fig. 1) including intermolecular contacts (normal forces, tangential forces and adhesive forces), with information at the particle scale ( Fig. 2) [35,50]. Motion equations are solved using the finite difference method.  To perform the simulation it is possible to import your own geometry or to select conveyor and feeder models from the application library. For moving parts, a suitable type of motion is assigned: translational, rotational, vibrating, force-fed, periodic. It is then necessary to prepare your own particle shape or to choose the shape and size of the particle from the available shapes ( fig. 3 Particle shapes: a) sphere, b) straight fiber, c) polyhedron, d) sphero-cylinder, e) spheropolygon, f) sphero-polyhedron, g) briquette, h) faceted cylinder.
Simultaneous simulation of particles of different shapes and dimensions is possible. The dimensions can be mapped on the scale of the actual grain dimensions, based on screen dimensions or calculated equivalent sphere diameter. The next step in preparing the simulation is to determine material properties such as density or bulk density, Yung's modulus or Poisson's coefficient, and then contact parameters (e.g. static and dynamic friction coefficient, tangential stiffness ratio, restitution coefficient) resulting from selected contact models (Fig. 2). For non-spherical particles two shredding models are available: Ab-T10 model and Tavares breakage model. Particle decay is performed using the Voronoi fracture algorithm [51]. Ab-T10 is a probability model of particle decay calculated from the value of accumulated contact energy [52]: where:e cumcumulative specific energy, S -selection function coefficient, L -particle size,L refreference size.
Tavares breakage model is based on the probability of upper-truncated log-normal distribution of unit crack energy [53,54]: where:e* -relative particle specific fracture energy,e 50median particle specific fracture energy, 2variance of the log-normal distribution of fracture energy. For mapping the creation of particles of a given size after grinding, an appropriate particle size distribution based on Gaudin-Schumann or Incomplete Beta can be adjusted. In order to determine the grinding parameters it is necessary to perform one of the following tests: drop weight testing, impact load target testing or repeated impact testing to obtain the minimum particle size distribution energy for a specific reference dimension. The analysis module allows to determine particle velocity, collision energy, particle size and mass after grinding as well as to determine the loads exerted by the material on machine components, including wear modelling.

Calibration and simulation
The calibration and validation of the model is a key element in creating a simulation using the discrete elements method, which determines the reliability and credibility of the results obtained [55]. When preparing a simulation and deciding on the choice of the calibration method, it is important to be aware of the difference between the concepts of material properties and model parameters. Material properties are essentially physical characteristics of the material, which can be determined on a micro scale (relations between single particles) or on a macro scale (properties determined for a set of particles) [55]. On the other hand, the model parameters refer to selected physical properties of the material defined in the micro scale and may vary depending on the selected contact model or used in DEM code modeling [55,56]. Depending on the research possibilities and approaches to determine material properties (at the particle or particle group level) two approaches to the calibration process are distinguished: the Bulk Calibration Approach (BCA) and the Direct Measuring Approach.(DMA) [55,57] (Fig. 4). The Bulk Calibration Approach. BCA is based on the calibration and mapping of the bulk material behavior in the DEM model based on results of bulk material properties tests. In this case, the bulk density, bulk stiffness, angle of repose and others are determined in the laboratory tests, after which the input parameters are changed interactively in the DEM model until the material response is consistent with the experimental tests [55,57]. Unfortunately, the disadvantage of this approach is that once determined parameters for one sy-mulation do not have to be suitable for the same material in another type of simulation [57]. Particle size and shape determination. In preparing the DEM model it is important to determine the parameters of particle shape and size. The shape and size of particles for different moisture levels of biomass grains is determined by means of granulometric analysis, thus it is possible to obtain the distribution of particle size distribution, characteristic dimensions and grain shape parameters [58][59][60].
Determination of bulk density. The bulk density of a material is a parameter that determines the volume of the material in the loose state, taking into account the pores and space between the grains. The value of the bulk density is influenced by, among others, humidity, specific density as well as the degree of filling of the preset volume. The bulk density is determined from the relation [2,36,58]: where: ρ B -bulk density in loose state, g•cm -3 , m 2weight of measuring cylinder and sample in loose state, g, m 1weight of measuring cylinder, g, Vcapacity of measuring cylinder, cm 3 .
Determination of the angle of repose. The angle of natural embankment is called the angle that is formed between the plane and the cone formed during free pouring of material on the plane. This angle may be measured directly or by reference to a relationship [2,35].

arctan
h AoR D  (4) where: AoRnatural angle of the repose, hheight of the cone, mm, Dlength of the cone base, mm. In the secondary case, image analysis and a quick-frame camera are used to record the arrangement of material. The rotating /drum method is used to determine the dynamic AoR, while others in Fig. 5 allows to determine the static AoR. The Direct Measuring Approach.(DMA). In DMA, the material parameters for individual particles are determined, which then serve as the corresponding input data in the DEM model. The problem is the determination of material properties of very small sizes [55,57]. There is also no certainty that after entering the obtained parameters into the DEM model the material behavior will correspond to the actual one [55,61,62]. The appropriate precision of the model in this case depends on the accuracy of the particle shape mapping and selection of the appropriate contact model [55,63]. Determination of particle-surface coefficient of friction. Determination of the coefficient of friction can be carried out using the test of the slope equilibrium with a variable slope angle and interchangeable substrates. The material is placed on a movable plate and starts to lift the plate until the material starts to slip evenly. In such an arrangement the tilt angle of the equatorial or the length of the equatorial throw together with its height can be directly measured (Fig. 2a) [35]. Based on the results obtained, the coefficient of friction can be calculated from the formula:  6. Example of determining the coefficient of friction by means of the inclined test a) for the particle-plate pair, b) for the particle-pair -particle.
When considering the movement of bulk material similar to a spherical or cylindrical shape, enabling it to roll over an inclined surface, the rolling friction moments and the values of the rolling friction coefficient (RFC) should be determined. For spheres for viscoelastic materials, the RFC without slip is [64] where: ƞ 1 , ƞ 2dissipative constants. In work [64] it was shown that knowing the RC for a given material, the value of RFC can be calculated. Determination of the internal friction coefficient (particle-particle). The coefficient of internal friction between particles can also be carried out using the experiment on a sloping equilibrium. In this case, particles of the material are additionally permanently placed on the equilibrium, after which the loose particles will slide (Fig. 2b) [12]. Determination of the restitution coefficient. The restitution coefficient is a measure of particle inflexibility during a collision. It is extremely important for determining and predicting the trajectory of materials after the impact. In simple terms, it describes the inflexible collision of a particle with a plane, and can be defined as [65,66]: where: v -particle velocity before collision, m•s -2 , v'velocity of the particle after collision, m•s -2 . For an ideal fully elastic collision, the coefficient of restitution is 1. In practice, this coefficient is within the range from 0-1. Simplifying with the omission of friction forces, this coefficient can be determined on the basis of the knowledge of changes of potential energy of the body E p1 falling from the height h 1 and potential energy after the reflection of E p2 from the surface to the height h 2 according to formula [35]: The RC is a measure of the inelasticity of particles during a collision. In simplified terms, it describes the collision of an inflexible particle with a plane, and can be defined as [65,66] 1 1 1 ' where: v -particle velocity before collision, m•s -2 , v'particle velocity after collision in time t i , m•s -2 . It can be measured as the time delay between successive impacts of particles on a hard surface by piezoelectric force sensors or accelerometers mounted in the plate. Delay, and thus RC, can also be measured by recording the sound of the particles hitting the disc [67]. RC is a function of impact velocity, particles masses, radii and materials parameters [67]. Depending on interaction force law, the RC can be determined by solving the equations of motion. For linear-dashpot force, the RC is constant and independent of collision speed [68] In fact, RC is dependent on speed and its fluctuations are observed, which result from the microscopic roughness of the surface of the spheres hitting the plate, causing the transfer of energy between the translational and rotational degree of freedom. RC increases as the speed of impact decreases. At constant speed, the fluctuation probability distribution of RC is non-Gaussian [58]. Calibration of shredding parameters. Shredding modelling requires, in addition to calibrating the contact parameters, the determination of the input parameters for particle division. The literature describes a number of tests that allow the determination of the grinding energy and particle size after grinding, including Bond Ball Mill Grindability test, Bond Rod Mill Grindability test, Bond Low-Energy Impact test, JKTech drop-weight test, SAG mill comminution test [69]. The brief characteristics of these tests are presented in Table 1. The Ab-T10 is one of the shredding models in the RockyDEM software and its parameters can be determined with the JKTech drop-weight test, so it will be discussed in more detail. The JKTech drop-weight test is divided into three parts: the first determination of resistance to impact cracking (material size 63-13.2 mm), the second determination of resistance to abrasion, breaking (particles size 53-37.5 mm), the third determination of mean density and dispersion for 30 particles [69].
In the Impact test, five dimensional fractions of the material are tested in three series differing in the shredding energy. In each series, 10 to 30 particles are shredded. In the tests, the level of shredding energy is known, which is determined by the weight and height from which it is lowered onto the material. The samples are then subjected to a sieve analysis after shredding. As a result of the test, a crack curve t 10 characteristic of the material is obtained, as determined by equation [69]: where: t 10 -the percent weight of fragments thatpasses 0,1 of its original size A, bmodel parameters, E csunit shredding energy(kWh/t).

Example of simulating particle crushing using DEM
The shredding simulation model was developed in the Rocky DEM environment. The crushing of grain was modelled, based on the results of grain size, force, stress and compression work of rice in a static compression test presented in [70]. First, the spatial form of rice grain was modeled and the working arrangement of the piston and plate for particle compression was mapped on the Instron 5966 machine (Fig. 7a). The simulation assumed a mass of the piston equal to 100 kg, and took into account the interaction of gravity force and earth acceleration (g=9, 81 m/s 2 ). The piston descending speed was set at 0.005 m/s. Table 2 shows the dimensional parameters of the modelled rice grain and Figure  7b shows the modelled grain shape.  . 7. a) Modelled piston-plate compression system, b) Modelled shape of rice grain.
Then contacts and friction coefficients were determined and verified. Mindlin-Deresiewicz contact model was used, adhesion forces were omitted. Other parameters used in the contact model are shown in Table 3. The Ab-T10 and Gaudin-Shumann models were chosen as the grinding model for both grains for the reproduction of particle size distribution after grinding. The grinding modeling parameters for rice are shown in Table 4.  Figure 8 shows the squeezed grain of rice and corn, generated. It should be noted that particles with different geometries were obtained. A dependency was used to assess the compatibility of the simulation model with the real one: where: ∆prelative error of the designated simulation parameter of the batch material motion, x 1sthe value of the designated parameter based on the simulation, x 1ethe value of the parameter designated based on the experiment.
The particle size distributions obtained as a result of simulation in a computer program and as a result of an experiment on a testing machine are similar to each other, although in the case of simulation finer particles were obtained ( Fig. 9 and Fig. 10). It should be pointed out that the simulation concerned one grain, while the results from the strength tests concern 100 grains, so discrepancies may result. The largest relative error of the share of particles after grinding for the simulation results was recorded for particles in size classes below 0.68 mm, for these particles the error was over 50%. As the particle size increased, the error value decreased logarithmically. The constructed simulation model best represented the percentage share of particles with dimensions above 1.46 mm, where the relative error was below 11%. The average relative error of the cumulative percentage of particles after grinding was 28.43%. The developed simulation model will be well suited to predict the size distribution of rice grinding particles to larger fractions, as indicated by the obtained errors of the cumulative percentage share of particles depending on the fraction size. It should be noted that in the real experiment the biological material -rice grains, which are characterized by great dimensional and structural diversity, was crushed, which causes a significant scattering of the results of physical experiments and determination of their properties (see: [70]). The method of discrete elements in case of prediction of granulometric composition takes the best results of compliance for the grinding of brittle materials, rocks for which structural diversity is smaller within one type of material, and calibration methods give more reliable results. Given the large variety of rice grains and the fact that the models available in the simulation software are rather adapted to the shredding of brittle materials, it is necessary to conduct research and develop calibration methods that can be successfully used in modelling the shredding of biological materials.

Summary
The method of discrete elements has recently been one of the most frequently used methods in the analysis and simulation of the flow of bulk and granular materials, supporting the optimization of machine and equipment design. Much attention is paid to issues related to the modelling of mixing processes, storage and transport of granular materials and crushing, mainly of brittle materials such as rocks and minerals. Calibrating and developing an adequate model using the discrete elements method is a time-consuming process that requires a number of experiments to determine the input parameters that determine the behavior of the bulk material. Unfortunately, sometimes Reliative error, %

Particle size, mm
Reliative error Average relative error correct determination of input materials will not guarantee that the simulation model will reflect reality. Simulation programs for flow modeling and grinding of materials by means of discrete elements facilitate the process of modeling, as they do not require entering own code sequences and the user can use ready-made solutions which shorten the time of simulation and calculations. Sometimes this can be a limitation, although in most programs it is also possible to enter your own macros and numerical codes. The process of shredding biological materials takes place in an uncontrolled and difficult to predict manner. In the case of these materials, the development of simulation models reflecting the actual behavior is problematic, especially as the properties of the materials change within one species, variety and also with the change in humidity. This paper presents the possibilities of simulation software and draws attention to the most important elements when determining the input parameters of the simulation model and the calibration process.