Semi-Autogenous Wet Grinding Modeling with CFD-DEM

: The paper highlights the features of constructing a model of a wet semi-autogenous grinding mill based on the discrete element method and computational ﬂuid dynamics. The model was built using Rocky DEM (v. 4.4.2, ESSS, Brazil) and Ansys Fluent (v. 2020 R2, Ansys, Inc., United States) software. A list of assumptions and boundary conditions necessary for modeling the process of wet semi-autogenous grinding by the ﬁnite element method is presented. The created model makes it possible to determine the energy-coarseness ratios of the semi-autogenous grinding (SAG) process under given conditions. To create the model in Rocky DEM the following models were used: The Linear Spring Rolling Limit rolling model, the Hysteretic Linear Spring model of the normal interaction forces and the Linear Spring Coulomb Limit for tangential forces. When constructing multiphase in Ansys Fluent, the Euler model was used with the primary phase in the form of a pulp with a given viscosity and density, and secondary phases in the form of air, crushing bodies and ore particles. The resistance of the solid phase to air and water was described by the Schiller–Naumann model, and viscosity by the realizable k-epsilon model with a dispersed multiphase turbulence model. The results of the work methods for material interaction coefﬁcients determination were developed. A method for calculating the efﬁciency of the semi-autogenous grinding process based on the results of numerical simulation by the discrete element method is proposed. The dependence of the probability of breakage of particles on speciﬁc impact energy.


Introduction
The Finite Element Method (FEM) and computational fluid dynamics (CFD) are branches of fundamental numerical modeling associated with the calculation of all contact interactions that occur between conditionally indivisible objects. Applied to the analysis of comminution processes, FEM began to be used with the development of computer technology at the end of the 20th century [1]. This method is especially widely used to optimize the design of comminution equipment [2,3].
The Discrete Element Method (DEM) is a numerical technique for predicting the behavior of bulk solids. DEM is a mesh-free method and does not solve the continuum equations of motion; this is how DEM differs from FEM. DEM modeling of tumbling mills reveals the characteristics of the mill media motion [4], predicting the wear of lifters, lining [5,6], grinding media [7] and energy consumption under various conditions [8,9]. However, the number of particles involved even in grinding processes largely limits the application of this method today. With the development of computer technology and software, the FEM finds more and more application in the field of comminution processes modeling and begins to compete with and complement software for steady-state modeling of comminution processes (JKSimMet, USIM PAC, etc.) [10][11][12].
The main advantage of fundamental modeling over empirical modeling is the ability to effectively analyze any factors affecting the comminution process. The combined use of these two directions of modeling allows engineers, designers, and scientists to Minerals 2021, 11,485 2 of 17 reduce the financial and labor resources required to justify the adopted technical and technological solutions [13][14][15].
In the coming decades, subject to the preservation of the dynamics of the development of computer technology, FEM modeling [16,17] will increasingly displace empirical modeling due to its universality.
The paper presents the concept of an approach to the design of a combined DEM-CFD model of semi-autogenous grinding, and the processing of the results obtained. The main motivation for this work is the need to study various technical and technological parameters, the influence of which on the SAG process cannot be investigated in static modeling.

Materials and Methods
The object of the study is a semi-autogenous grinding mill SAG 34 × 17 (the diameter of the mill is 34 feet, the length is 17 feet), which processes charges of copper-nickel ores. The ores that make up the charge are divided into 3 conventional types. For each type of ore, the strength characteristics were determined using the drop weight test methods (JK DWT [11,17,18], SMC Testing [19], Australia) and F. Bond test [20] for ball grinding. The results are presented in Table 1. The test results show that the studied ores M1 and M2 based on their strength properties belong to very soft types of ores, and ore T1 to medium ones. That is, the ores composing the charge are contrasting in their strength properties.
The mass fractions of each type of ore in the charge processed by the plant in 2019 are shown in Figure 1. financial and labor resources required to justify the adopted technical and technological solutions [13][14][15].
In the coming decades, subject to the preservation of the dynamics of the development of computer technology, FEM modeling [16,17] will increasingly displace empirical modeling due to its universality.
The paper presents the concept of an approach to the design of a combined DEM-CFD model of semi-autogenous grinding, and the processing of the results obtained. The main motivation for this work is the need to study various technical and technological parameters, the influence of which on the SAG process cannot be investigated in static modeling.

Materials and Methods
The object of the study is a semi-autogenous grinding mill SAG 34x17 (the diameter of the mill is 34 feet, the length is 17 feet), which processes charges of copper-nickel ores. The ores that make up the charge are divided into 3 conventional types. For each type of ore, the strength characteristics were determined using the drop weight test methods (JK DWT [11,17,18], SMC Testing [19], Australia) and F. Bond test [20] for ball grinding. The results are presented in Table 1.  18.60 The test results show that the studied ores M1 and M2 based on their strength properties belong to very soft types of ores, and ore T1 to medium ones. That is, the ores composing the charge are contrasting in their strength properties.
The mass fractions of each type of ore in the charge processed by the plant in 2019 are shown in Figure 1. The number of the SAG unit samplings for one year showed that the feed size of the F80 mill varied in the range from 47,603 to 103,993 μm, and for the P80 from 920 to 5267 μm.
The contrast in the strength properties of the processed ores and the instability of the mass fraction of each type of ore in the charge have a negative effect on the process of semi-autogenous grinding:

•
In the mill, overgrinding of soft types of ores occurs, leading to sliming and a decrease in the efficiency of subsequent beneficiation processes [22][23][24]; The number of the SAG unit samplings for one year showed that the feed size of the F80 mill varied in the range from 47,603 to 103,993 µm, and for the P80 from 920 to 5267 µm.
The contrast in the strength properties of the processed ores and the instability of the mass fraction of each type of ore in the charge have a negative effect on the process of semi-autogenous grinding:

•
In the mill, overgrinding of soft types of ores occurs, leading to sliming and a decrease in the efficiency of subsequent beneficiation processes [22][23][24]; • Accumulation of critical fraction ore in the mill is observed; • The instability of the composition ( Figure 1) and the ore size lead to the instability of the size of the mill product. The mass fractions of various types of ores are selected based on the useful component content in them.
The aim of the work was to create a stable FEM model of the wet semi-autogenous grinding process and to quantify the efficiency of processing copper-nickel ores under given conditions. The SAG 34 × 17 model was created by the discrete element method and computational fluid dynamics in Rocky DEM and Ansys Fluent software.

Setting the Modeling Boundary Conditions
DEM modeling of bulk materials is a special case of the FEM method. Ore particles represented in such a model are indivisible objects, between the surfaces of which interactions are calculated by various types of contact models. With an increase in the number of particles participating in the simulation, and with a decrease in their size, the total number of fixed contacts grows exponentially. In this regard, DEM modeling is inevitably associated with the assumptions inherent in the FEM method [25,26]. In this study, the following assumptions were made:

•
The shape of the ore particles is represented by polyhedrons with 10 vertices; • The minimum particle size shown in the model is 5 mm, particle size distribution is shown in Figure 2 and the grinding media is represented by 100% balls Ø125 mm; • Part of the solid fraction is modeled as a suspension in the pulp of a given density and viscosity.

•
The design area of the mill is represented by a linear section of the cylindrical part of the drum 500 mm long, the end walls of which are solid surfaces with zero coefficients of friction in relation to particles and pulp. • Accumulation of critical fraction ore in the mill is observed; • The instability of the composition ( Figure 1) and the ore size lead to the instability of the size of the mill product. The mass fractions of various types of ores are selected based on the useful component content in them.
The aim of the work was to create a stable FEM model of the wet semi-autogenous grinding process and to quantify the efficiency of processing copper-nickel ores under given conditions. The SAG 34x17 model was created by the discrete element method and computational fluid dynamics in Rocky DEM and Ansys Fluent software.

Setting the Modeling Boundary Conditions
DEM modeling of bulk materials is a special case of the FEM method. Ore particles represented in such a model are indivisible objects, between the surfaces of which interactions are calculated by various types of contact models. With an increase in the number of particles participating in the simulation, and with a decrease in their size, the total number of fixed contacts grows exponentially. In this regard, DEM modeling is inevitably associated with the assumptions inherent in the FEM method [25,26]. In this study, the following assumptions were made:

•
The shape of the ore particles is represented by polyhedrons with 10 vertices; • The minimum particle size shown in the model is 5 mm, particle size distribution is shown in Figure 2 and the grinding media is represented by 100% balls Ø125 mm; • Part of the solid fraction is modeled as a suspension in the pulp of a given density and viscosity.
• The design area of the mill is represented by a linear section of the cylindrical part of the drum 500 mm long, the end walls of which are solid surfaces with zero coefficients of friction in relation to particles and pulp.

Contact Models
The choice of contact interaction models was dictated by a decrease in the computational power consumed during modeling. The following models were used in the DEM simulation module: particle rolling resistance model-not chosen (due to the use of polyhedrons instead of spheres); normal contacts model-Linear Spring Dashpot [27,28]; adhesion model-disabled (due to the use of three-phase modeling, adhesion of particles inside the pulp should not be observed); and tangential contact model-Linear Spring Coulomb Dashpot.
The following models were used in the CFD module: Euler's multiphase model [29] with three phases (the primary phase is pulp, the secondary phases are ore particles, balls

Contact Models
The choice of contact interaction models was dictated by a decrease in the computational power consumed during modeling. The following models were used in the DEM simulation module: particle rolling resistance model-not chosen (due to the use of polyhedrons instead of spheres); normal contacts model-Linear Spring Dashpot [27,28]; adhesion model-disabled (due to the use of three-phase modeling, adhesion of particles inside the pulp should not be observed); and tangential contact model-Linear Spring Coulomb Dashpot.
The following models were used in the CFD module: Euler's multiphase model [29] with three phases (the primary phase is pulp, the secondary phases are ore particles, balls and air). In a pair of pulp-air phases, a symmetric interaction was used, in a pair of particles-pulp and particle-air phases-the Schiller-Naumann interaction model [30]. The orthogonal cell size of the computational grid was 60 mm, with a single cell for the entire width of the computational zone. The laminar viscosity model was used. The air had a Minerals 2021, 11, 485 4 of 17 density of 1.225 kg/m 3 and a viscosity of 1.789 × 10 −5 kg/m·s (reference data) [31]; the pulp had a density of 1992 kg/m 3 and a viscosity of 1.3 kg/m·s (according to the results of laboratory determination of viscosity with a disk viscometer and density by measuring with the Marcy system). The method for calculating the combined pressure and velocity fields was Phase Coupled SIMPLE [32]. Spatial discreteness of pressures was calculated by the PRESTO! [33].

Calibration of Material Interaction Coefficients
To achieve the best convergence of DEM simulation results with experimental data, it is necessary to carry out calibration tests [34,35]. Calibration tests are laboratory experiments designed to find coefficients used in models of contact interaction of particles with each other and with geometry surfaces in DEM simulations. These include the coefficients of friction and restitution between material pairs "ore-ore", "ore-steel" and "steel-steel".
The essence of each calibration test is to conduct such a laboratory experiment, the result of which depends to a greater extent on the desired interaction parameter and to a lesser extent on other variables. The conditions and result of the laboratory test are recorded. After that, the laboratory experiment is carried out numerically in the DEM model, while maintaining all the conditions of the experiment. By varying the sought-for interaction parameter, one finds such a value at which the result of numerical simulation corresponds to a laboratory experiment. After that, the found value of the required parameter is considered calibrated.
In this work, a representative sample of the investigated charge was prepared for calibration tests in accordance with the arithmetic mean proportion of the studied types of ores in it (mass fraction of ore M1-40%, M2-25%, T1-35%). Particles of the correct shape from the size class −12.5 +10.0 mm of a representative sample of the investigated charge were used. To determine the coefficients of friction, the surface of the particles used was sprayed with water from a spray bottle until the visual effect of wetting with water was achieved.

Restitution Coefficient (Ore-Steel)
This experiment consisted of alternately dropping 50 ore particles onto an inclined steel surface of a steel cube with sides of 145 mm, as shown in Figure 3. The output of the experiment was the average distance traveled by particles from the point of impact on a steel inclined surface to the point of impact on a horizontal surface.   The arithmetic mean rebound of the particles was 182 mm. Based on the results of varying the restitution coefficient in the model of this experiment, duplicated in Rocky DEM software, this rebound value corresponds to a recovery factor of 0.703.

Restitution Coefficient (Ore-Ore)
This experiment consisted of alternately dropping 50 ore particles onto an inclined surface of a flat ore sample. The general configuration of the experiment is shown in Figure 4. The output of the experiment was the arithmetic mean distance traveled by ore particles from the point of impact on the inclined surface of a flat ore sample to the point of impact on the horizontal surface [36].

Friction Coefficient (Ore-Ore)
The test was conducted to determine the angle of repose of the particles of the charge. The massif is formed after lifting from the h vertically located plastic tube filled with the test material. The mass o such that the height of the material column in the tube is at least four t diameter [37]. In the present test, the particles' mass was 6856.3 g. The i was 103 mm. The tube was lifted at a speed of approximately 2 cm/s, a of particles was then located under the tube. The arithmetic mean a 39.1º. According to the results of the numerical simulation in Rocky value of the angle of repose corresponds to a coefficient of static fric Figure 5 illustrates the numerical (a) and full-scale (b) of the experime The arithmetic mean rebound of the particles was 182 mm. Based on the results of varying the restitution coefficient in the model of this experiment, duplicated in Rocky DEM software, this rebound value corresponds to a recovery factor of 0.623.

Friction Coefficient (Ore-Ore)
The test was conducted to determine the angle of repose of the assemblage of the particles of the charge. The massif is formed after lifting from the horizontal surface a vertically located plastic tube filled with the test material. The mass of particles must be such that the height of the material column in the tube is at least four times the inner tube diameter [37]. In the present test, the particles' mass was 6856.3 g. The inner pipe diameter was 103 mm. The tube was lifted at a speed of approximately 2 cm/s, and the assemblage of particles was then located under the tube. The arithmetic mean angle of repose was 39.1º. According to the results of the numerical simulation in Rocky DEM software, this value of the angle of repose corresponds to a coefficient of static friction equal to 0.870. Figure 5 illustrates the numerical (a) and full-scale (b) of the experiment.
To estimate the effect of the dynamic friction coefficient of the ore on the dynamics of material movement in a rotating cylinder, a series of numerical experiments was carried out in Rocky DEM. It was found that the coefficient of dynamic friction did not significantly affect the trajectories of material. In this regard, the coefficient of dynamic friction was taken to be practically equal to the coefficient of static friction (0.855). Figure 6 shows 2 experiments from this series of experiments in which all the initial parameters were the same, except for the coefficient of dynamic friction between the ore particles, equal to 0.870 and 0.570, respectively.
such that the height of the material column in the tube is at least four times the inner tube diameter [37]. In the present test, the particles' mass was 6856.3 g. The inner pipe diameter was 103 mm. The tube was lifted at a speed of approximately 2 cm/s, and the assemblage of particles was then located under the tube. The arithmetic mean angle of repose was 39.1º. According to the results of the numerical simulation in Rocky DEM software, this value of the angle of repose corresponds to a coefficient of static friction equal to 0.870. Figure 5 illustrates the numerical (a) and full-scale (b) of the experiment. To estimate the effect of the dynamic friction coefficient of the ore on the dynamics of material movement in a rotating cylinder, a series of numerical experiments was carried out in Rocky DEM. It was found that the coefficient of dynamic friction did not significantly affect the trajectories of material. In this regard, the coefficient of dynamic friction was taken to be practically equal to the coefficient of static friction (0.855). Figure 6 shows 2 experiments from this series of experiments in which all the initial parameters were the same, except for the coefficient of dynamic friction between the ore particles, equal to 0.870 and 0.570, respectively. Figure 6. Two images from a series of numerical experiments to determine the influence of the dynamic friction coefficient between ore particles on the overall dynamics of the material in the rotating drum.

Friction Coefficient (Ore-Steel)
The test was carried out to determine the inclination angle of the steel surface at which the particle assemblage slid off it. The particles were in a thick paper ring 50 mm in height and 160 mm in diameter. The weight of the particles was 2454.0 g, and the steel surface and the particle assemblage were wetted with water.
The arithmetic mean sliding angle under the described conditions was 28.0°. According to the results of numerical simulation, this value of the sliding angle corresponds to the static friction coefficient at rest, equal to 0.553. Figure 7 shows the numerical (a) and laboratory (b) conditions of the experiment. By analogy with the "ore-ore" friction coefficients, the dynamic friction coefficient was taken as almost equal to the static friction coefficient (0.545). Figure 6. Two images from a series of numerical experiments to determine the influence of the dynamic friction coefficient between ore particles on the overall dynamics of the material in the rotating drum.

Friction Coefficient (Ore-Steel)
The test was carried out to determine the inclination angle of the steel surface at which the particle assemblage slid off it. The particles were in a thick paper ring 50 mm in height and 160 mm in diameter. The weight of the particles was 2454.0 g, and the steel surface and the particle assemblage were wetted with water.
The arithmetic mean sliding angle under the described conditions was 28.0 • . According to the results of numerical simulation, this value of the sliding angle corresponds to the static friction coefficient at rest, equal to 0.553. Figure 7 shows the numerical (a) and laboratory (b) conditions of the experiment. By analogy with the "ore-ore" friction coefficients, the dynamic friction coefficient was taken as almost equal to the static friction coefficient (0.545).

Friction Coefficient (Steel-Steel)
The test was performed to determine the angle of inclination of a steel surface, at which an array of three steel 80 mm balls, fastened to each other, slid off. Prior to the experiment, the balls and steel surface were wetted with water to achieve the friction conditions estimated for wet grinding. The average sliding angle obtained with 5 measurements was 22.4 • . According to the results of numerical simulation of this experiment, carried out in the Rocky DEM software, this value of the slip angle corresponds to a coefficient of static friction of 0.405. Figure 8 displays the numerical (a) and full-scale (b) of the experiment conditions. By analogy with the "ore-ore" friction coefficients, the sliding friction coefficient was taken as almost equal to the static friction coefficient (0.400). Figure 6. Two images from a series of numerical experiments to determine the influence of the dynamic friction coefficient between ore particles on the overall dynamics of the material in the rotating drum.

Friction Coefficient (Ore-Steel)
The test was carried out to determine the inclination angle of the steel surfac which the particle assemblage slid off it. The particles were in a thick paper ring 50 in height and 160 mm in diameter. The weight of the particles was 2454.0 g, and the s surface and the particle assemblage were wetted with water.
The arithmetic mean sliding angle under the described conditions was 28.0°. Acc ing to the results of numerical simulation, this value of the sliding angle correspond the static friction coefficient at rest, equal to 0.553. Figure 7 shows the numerical (a) laboratory (b) conditions of the experiment. By analogy with the "ore-ore" friction co cients, the dynamic friction coefficient was taken as almost equal to the static friction efficient (0.545).  The test was performed to determine the angle of inclination of a steel surfac which an array of three steel 80 mm balls, fastened to each other, slid off. Prior to experiment, the balls and steel surface were wetted with water to achieve the friction ditions estimated for wet grinding. The average sliding angle obtained with 5 meas ments was 22.4°. According to the results of numerical simulation of this experiment ried out in the Rocky DEM software, this value of the slip angle corresponds to a co cient of static friction of 0.405. Figure 8 displays the numerical (a) and full-scale (b) o experiment conditions. By analogy with the "ore-ore" friction coefficients, the sliding tion coefficient was taken as almost equal to the static friction coefficient (0.400). The interaction coefficients obtained in calibration tests are presented in Table 2

Description of the Model Conditions
The geometry of the mill was simplified to the surface of its lining. The calcu area with dimensions is shown in Figure 9. The interaction coefficients obtained in calibration tests are presented in Table 2.  1 The restitution coefficient of the pair "steel-steel" is taken according to reference data.

Description of the Model Conditions
The geometry of the mill was simplified to the surface of its lining. The calculated area with dimensions is shown in Figure 9. Table 3 shows the main characteristics of the semi-autogenous grinding process, adopted for modeling. The duration of the simulation time was 20 s. The number of particles in the model was 74 thousand.
Ore-steel 0.703 0.560 0.545 Steel-steel 0.900 1 0.405 0.400 1 The restitution coefficient of the pair "steel-steel" is taken according to reference da

Description of the Model Conditions
The geometry of the mill was simplified to the surface of its lining. The calcu area with dimensions is shown in Figure 9.

Results
The general dynamics of the material in the mill (a) and the distribution of the mass fraction of the liquid phase in the volume (b) are shown in Figure 10.  Table 3 shows the main characteristics of the semi-autogenous grinding process, adopted for modeling. The duration of the simulation time was 20 s. The number of particles in the model was 74 thousand.

Results
The general dynamics of the material in the mill (a) and the distribution of the mass fraction of the liquid phase in the volume (b) are shown in Figure 10. The steady-state mode of operation of the SAG mill was achieved when the weighted average coordinates of all particles and grinding media along the vertical and horizontal axes came to conventionally unchanged values (± 0.15 m) ( Figure 11). The analysis of the energy interactions of the material was carried out for the simulation time interval from 10 to 20 s. The steady-state mode of operation of the SAG mill was achieved when the weighted average coordinates of all particles and grinding media along the vertical and horizontal axes came to conventionally unchanged values (± 0.15 m) ( Figure 11). The analysis of the energy interactions of the material was carried out for the simulation time interval from 10 to 20 s. Minerals 2021, 11, x FOR PEER REVIEW 9 of 17 Figure 11. Average particle coordinates during simulation. Figure 12 shows the net power consumed by the mill, obtained from the simulation results after reaching steady-state mill operation (after 10 seconds). These values were recalculated in proportion to the difference in mill feed weight, considering the simulation of a 500 mm drum width from 5180 mm. The arithmetic mean value was 7017 kW. This value of energy consumption does not consider the imperfection in all units of a mill. To calculate the real energy consumption, it is necessary to consider the power losses in the engine (up to 8.5%), main bearings, girth gear, pinion shaft bearings (up to 2.5%) and gearbox (2%) [38]. Thus, it can be assumed that the total power consumed by the mill drive is 8026 kW, which is 7.96% different from the actual power consumption of the mill drive under the given grinding conditions (8720 kW).

SAG Energy Spectra
The main challenge in DEM modeling is interpreting the simulation results in terms of the efficiency of the grinding process.
Due to the fact that direct modeling of the ore particles' breakage with the formation of fragments is an extremely resource-intensive calculation, the process is evaluated by analyzing the statistics of collisions of particles with each other, balls and lining. The main tool in the analysis is a set of so-called "energy spectra" curves for the studied grades of ore size in the mill. The energy spectra is the distribution of the cumulative specific power (W/kg) attributable to impacts of various specific energy (J/kg or kWh/t) over narrow particle size classes.
The obtained energy spectra for the charges under study are shown in Figure 13. All obtained curves are divided into size fractions and types of ores. Net power, kW Time, s Figure 11. Average particle coordinates during simulation. Figure 12 shows the net power consumed by the mill, obtained from the simulation results after reaching steady-state mill operation (after 10 s). These values were recalculated in proportion to the difference in mill feed weight, considering the simulation of a 500 mm drum width from 5180 mm. The arithmetic mean value was 7017 kW. This value of energy consumption does not consider the imperfection in all units of a mill. To calculate the real energy consumption, it is necessary to consider the power losses in the engine (up to 8.5%), main bearings, girth gear, pinion shaft bearings (up to 2.5%) and gearbox (2%) [38]. Thus, it can be assumed that the total power consumed by the mill drive is 8026 kW, which is 7.96% different from the actual power consumption of the mill drive under the given grinding conditions (8720 kW).  Figure 12 shows the net power consumed by the mill, obtained from the simulation results after reaching steady-state mill operation (after 10 seconds). These values were recalculated in proportion to the difference in mill feed weight, considering the simulation of a 500 mm drum width from 5180 mm. The arithmetic mean value was 7017 kW. This value of energy consumption does not consider the imperfection in all units of a mill. To calculate the real energy consumption, it is necessary to consider the power losses in the engine (up to 8.5%), main bearings, girth gear, pinion shaft bearings (up to 2.5%) and gearbox (2%) [38]. Thus, it can be assumed that the total power consumed by the mill drive is 8026 kW, which is 7.96% different from the actual power consumption of the mill drive under the given grinding conditions (8720 kW).

SAG Energy Spectra
The main challenge in DEM modeling is interpreting the simulation results in terms of the efficiency of the grinding process.
Due to the fact that direct modeling of the ore particles' breakage with the formation of fragments is an extremely resource-intensive calculation, the process is evaluated by analyzing the statistics of collisions of particles with each other, balls and lining. The main tool in the analysis is a set of so-called "energy spectra" curves for the studied grades of ore size in the mill. The energy spectra is the distribution of the cumulative specific power (W/kg) attributable to impacts of various specific energy (J/kg or kWh/t) over narrow particle size classes.
The obtained energy spectra for the charges under study are shown in Figure 13. All obtained curves are divided into size fractions and types of ores.

SAG Energy Spectra
The main challenge in DEM modeling is interpreting the simulation results in terms of the efficiency of the grinding process.
Due to the fact that direct modeling of the ore particles' breakage with the formation of fragments is an extremely resource-intensive calculation, the process is evaluated by analyzing the statistics of collisions of particles with each other, balls and lining. The main tool in the analysis is a set of so-called "energy spectra" curves for the studied grades of ore size in the mill. The energy spectra is the distribution of the cumulative specific power (W/kg) attributable to impacts of various specific energy (J/kg or kWh/t) over narrow particle size classes.
The obtained energy spectra for the charges under study are shown in Figure 13. All obtained curves are divided into size fractions and types of ores. Minerals 2021, 11, x FOR PEER REVIEW 10 of 17 Figure 13. Energy spectra obtained by modeling.
Each narrow size fraction is characterized by its own average value of the specific minimum energy required to break a particle. The larger the particle, the lower the specific impact energy required to break it (hereinafter, referred to as the nominal breakage energy of the particles). This is due to the natural fracturing of mineral raw materials: statistically, large ore samples show more cracks of varying degrees of liberation than small ones.

Method for Determining the Nominal energy of Particles Breakage
To establish the facts of the destruction of particles of various sizes, a methodology is proposed for determining the nominal particle-specific breakage energy of the investigated samples in different size classes. This value allows us to determine the amount of energy accounted for by such collisions of particles that lead to their failure. This method is intended primarily for comparative analysis of SAG under various conditions.
The laboratory equipment for the study is proposed to use a drop weight tester used in the JK DWT [17] and SMC test [19] methodologies. The JK DWT test methodology includes a stage for determining the nature of damage in the low-energy region (parameter ta); however, in this test, several particles are used simultaneously without recording the number of destroyed ones. In the area of high specific energies in the JK DWT test (parameters A and b), an energy of ≥ 0.1 kWh/t always leads to particle breakage. In this regard, the results of the JK DWT test cannot be processed in such a way as to obtain from them the value of the nominal specific energy, which is the minimum necessary to destroy particles of a certain size. In this case, the drop weight tester is a ready-made tool for determining this value.
The study assumes the selection of particles of narrow size fractions, which can be obtained from intermediate size fractions when sampling for the JK DWT test (−75.0 + 63.0 mm, −53.0 + 45.0 mm, −37.5 + 31.5 mm and −26.5 + 22.4 mm). The average particle sizes in these sets are 69.0, 49.0, 34.5 and 24.5 mm, respectively. It has been calculated that the masses of particles in these size fractions will make it possible to achieve the minimum specific impact energy from 0.005 kWh/t for 24.5 mm and from 0.0005 kWh/t for 69.0 mm, provided that the light DW base (LB, 2.7 kg) is dropped at a height of 50 mm. Such a short drop height will not lead to the breakage of particles, i.e., will be lower than their minimum specific energy of breakage. By gradually increasing the drop height, it is necessary to achieve particle breakage. Figure 14 shows a conditional probability graph of particle breakage depending on the value of the specific impact energy. Cumulative specific power, kW/t Specific energy, kWh/t Each narrow size fraction is characterized by its own average value of the specific minimum energy required to break a particle. The larger the particle, the lower the specific impact energy required to break it (hereinafter, referred to as the nominal breakage energy of the particles). This is due to the natural fracturing of mineral raw materials: statistically, large ore samples show more cracks of varying degrees of liberation than small ones.

Method for Determining the Nominal Energy of Particles Breakage
To establish the facts of the destruction of particles of various sizes, a methodology is proposed for determining the nominal particle-specific breakage energy of the investigated samples in different size classes. This value allows us to determine the amount of energy accounted for by such collisions of particles that lead to their failure. This method is intended primarily for comparative analysis of SAG under various conditions.
The laboratory equipment for the study is proposed to use a drop weight tester used in the JK DWT [17] and SMC test [19] methodologies. The JK DWT test methodology includes a stage for determining the nature of damage in the low-energy region (parameter t a ); however, in this test, several particles are used simultaneously without recording the number of destroyed ones. In the area of high specific energies in the JK DWT test (parameters A and b), an energy of ≥ 0.1 kWh/t always leads to particle breakage. In this regard, the results of the JK DWT test cannot be processed in such a way as to obtain from them the value of the nominal specific energy, which is the minimum necessary to destroy particles of a certain size. In this case, the drop weight tester is a ready-made tool for determining this value.
The study assumes the selection of particles of narrow size fractions, which can be obtained from intermediate size fractions when sampling for the JK DWT test (−75.0 + 63.0 mm, −53.0 + 45.0 mm, −37.5 + 31.5 mm and −26.5 + 22.4 mm). The average particle sizes in these sets are 69.0, 49.0, 34.5 and 24.5 mm, respectively. It has been calculated that the masses of particles in these size fractions will make it possible to achieve the minimum specific impact energy from 0.005 kWh/t for 24.5 mm and from 0.0005 kWh/t for 69.0 mm, provided that the light DW base (LB, 2.7 kg) is dropped at a height of 50 mm. Such a short drop height will not lead to the breakage of particles, i.e., will be lower than their minimum specific energy of breakage. By gradually increasing the drop height, it is necessary to achieve particle breakage. Figure 14 shows a conditional probability graph of particle breakage depending on the value of the specific impact energy. To analyze the curves of energy spectra, it is necessary to determine th cific impact energy ( ) . , at which a particle of a certain size will be broke ability of 0.95.
The test methodology consists of the following algorithm of actions: 1) Preparation of sets of 20 particles of each size fractions: −75.0 + 63 45.0 mm, −37.5 + 31.5 mm and −26.5 + 22.4 mm. The selection is perform with the SMC test [20], but with a decrease in the tolerance of the mass of set from ± 30 to ± 5% of the average mass of a particle in the set. All parti the correct shape. At the end of the selection, four narrow size fractions shou within which the lightest particle differs in its mass from the heaviest by 10%; 2) Sequential breakage of all particles of each set in the JK DWT tester 2.1 The initial drop height of a light weight (LB, 2703.4 g) for all part from the top surface of the particle to the impact surface of the dropping w 2.2 The weight is dropped onto the particle, and after each drop the s creased by 10 mm until the particle is broken. The height is recorded in the which the particle is broken. A particle is considered broken when its la passes through a sieve with a cell size of the lower limit of the correspondin or has a mass of no more than 0.8 of the average particle masses in the set; 3) The average value of the installation heights is calculated, at whic of particles of each of the size fractions occurred. The resulting values a Table 4; To analyze the curves of energy spectra, it is necessary to determine the nominal specific impact energy E 0.95 cs(j) , at which a particle of a certain size will be broken with a probability of 0.95.
The test methodology consists of the following algorithm of actions:  [20], but with a decrease in the tolerance of the mass of particles in the set from ± 30 to ± 5% of the average mass of a particle in the set. All particles must have the correct shape. At the end of the selection, four narrow size fractions should be selected, within which the lightest particle differs in its mass from the heaviest by no more than 10%; (2) Sequential breakage of all particles of each set in the JK DWT tester; 2.1 The initial drop height of a light weight (LB, 2703.4 g) for all particles is 50 mm from the top surface of the particle to the impact surface of the dropping weight; 2.2 The weight is dropped onto the particle, and after each drop the set height is increased by 10 mm until the particle is broken. The height is recorded in the table as h i,j , at which the particle is broken. A particle is considered broken when its largest fragment passes through a sieve with a cell size of the lower limit of the corresponding size fraction or has a mass of no more than 0.8 of the average particle masses in the set; (3) The average value of the installation heights is calculated, at which the breakage of particles of each of the size fractions occurred. The resulting values are recorded in Table 4; (4) Calculation of the average drop height of the weight at which the particles breakage in the j-th set occurs: The average specific energy of particles' breakage of each set is determined by the equation: where: E 0.5 cs(j) -specific energy of particles breakage, at which with a probability of 0.5 a particle in the j-th set will be broken, kWh/t; Q w . -weight mass, kg; h j -average drop height of the weight at which the particles breakage in the j-th set occurs, m; Q pj -average particle mass in the j-th set, g; G-acceleration of gravity, m/s 2 .
The nominal specific energy of particles breakage of each set is determined by the equation: where E 0.95 cs(j) -specific impact energy, with a probability of 0.95 leading to the j-th set particle breakage, kWh/t; σ j -standard deviation of the achieved specific energy E cs of the j-th set of particles breakage. According to the method described above, for the tested samples of M1, M2 and T1 ores, the dependences of the nominal specific particle's breakage energy on their size were obtained. The results are shown in Table 5.
The results obtained from the analysis of the dependence of the nominal specific particle's breakage energy on the various size fractions of the three studied types of ores are presented in Figure 15.
The obtained value of nominal specific particle's breakage energy E 0.95 cs(j) cuts on the energy spectra of the j-th size fraction the useful collisions, occurring in the mill, from those which do not lead to breakage of the material, see Figure 16. Numerically, the total power transmitted to the j-th size fraction can be divided into two components: where W s(j) -is the cumulative specific power of the j-th size fraction, kW/t; W l(j) -specific cumulative power of the j-th size fraction, performed by impacts that do not lead to the particle's breakage, kW/t; and W u(j) -useful cumulative specific power of the j-th size fraction, performed by impacts leading to the particle's breakage, kW/t. This method assumes that the ratios of the infinitely narrow particle size fraction destruction rate for different operating modes must be equal to the corresponding ratios of useful energies W u(j) for these size fractions in the investigated SAG operating modes. Because the energy spectra curves are similar to each other (in specific particle size fractions), the difference in the selected nominal specific particle's breakage energy (within reasonable limits) will not have a very high impact on the final useful energies' ratio. This is the reason why the effect of particles weakening [39,40] during contacts with specific energies smaller than the nominal specific breakage energy may be not considered.   The results obtained from the analysis of the dependence of the nominal specific particle's breakage energy on the various size fractions of the three studied types of ores are presented in Figure 15. The overall efficiency of the grinding process can be determined as the weighted average over all size fractions, and the completeness of the transfer of the total consumed energy into the contact interaction of particles in the range of specific energies leading to their breakage: , % Thus, having determined the nominal specific particle's breakage energy E 0.95 cs(j) in the entire range of material sizes by calculation from the obtained dependencies, it is possible to process the data of energy spectra under specified conditions and determine the overall efficiency of the SAG process, see Table 6.
where ( ) -is the cumulative specific power of the j-th size fraction, kW/t; ( ) -specific cumulative power of the j-th size fraction, performed by impacts that do not lead to the particle's breakage, kW/t; and ( ) -useful cumulative specific power of the j-th size fraction, performed by impacts leading to the particle's breakage, kW/t. This method assumes that the ratios of the infinitely narrow particle size fraction destruction rate for different operating modes must be equal to the corresponding ratios of useful energies ( ) for these size fractions in the investigated SAG operating modes. Because the energy spectra curves are similar to each other (in specific particle size fractions), the difference in the selected nominal specific particle's breakage energy (within reasonable limits) will not have a very high impact on the final useful energies' ratio. This is the reason why the effect of particles weakening [39,40] during contacts with specific energies smaller than the nominal specific breakage energy may be not considered. The overall efficiency of the grinding process can be determined as the weighted average over all size fractions, and the completeness of the transfer of the total consumed energy into the contact interaction of particles in the range of specific energies leading to their breakage: Thus, having determined the nominal specific particle's breakage energy ( ) .
in the entire range of material sizes by calculation from the obtained dependencies, it is possible to process the data of energy spectra under specified conditions and determine the overall efficiency of the SAG process, see Table 6.  In Table 6, it is possible to track the share of each j-th size fraction in the value of the overall efficiency of the grinding process. In fact, the value E i,j shows the rate of breakage of the j-th size fraction relative to all others. For the simulation with the conditions in which the change in the particle size distribution of material is expected: the fine fractions' masses must decrease, and the coarse fractions' masses must increase, according to their E i,j /γ i,j ratios.
The obtained value of the efficiency of the semi-autogenous grinding process E en.gr can be used as an optimization criterion when comparing various technical and technological parameters of the SAG process in FEM-designed models. In comparative calculations with the presented method, it is expected that the resulting value of the grinding efficiency E en.gr will be inversely proportional to the specific power consumption of the SAG unit.

Discussion
With the development of computer hardware and software, fundamental discrete element modeling is increasingly used in the field of comminution processes modeling. This study presents the author's approach to the construction of a multiphase model of the SAG process and the processing of its results.
The developed integrated approach allows us to create dynamic models of AG and SAG mills. In the future, dynamic models can be used to increase the versatility of static models, as well as to obtain design points when creating digital twins of semi-autogenous grinding mills.
Further directions of research in this area: • standardization of laboratory methods for calibrating the parameters of material interactions-this kind of research should be easy to perform and aimed at finding specific parameters that govern the model of interaction of materials present in the model; • standardization and substantiation of laboratory methods for studying the strength of mineral raw materials, the results of which can be used to analyze models of the SAG process, designed with DEM; • development of methods for processing the results of DEM modeling of tumbling mills, aimed at predicting the productivity and grain size distribution of its product.

Conclusions
As a result of the work, a set of models for the interaction of ore and grinding media was determined, which makes it possible to design a stable multiphase DEM-CFD model of wet semi-autogenous grinding. A method for determining the coefficients required to create a DEM-CFD model of wet semi-autogenous grinding and a method for laboratory determination of the nominal specific energy mineral raw material particle's breakage in various sizes fractions using a JK DWT tester have been developed. This indicator can be used when processing the results of DEM modeling of the AG/SAG mill.
A method for calculating the efficiency of the semi-autogenous grinding process based on the results of numerical simulation by the discrete element method is proposed.