Improved Modeling of the Grinding Process through the Combined Use of Matrix and Population Balance Models

The mechanistic approach has proven so far to be flexible and successful for simulation of the grinding process. The basic idea underlying mechanistic models, namely the matrix and population balance models, is based on the identification of natural events during grinding. Since each model has its own capabilities and limitations, their combined use may offer additional advantages on this aspect. In this study, the matrix model and the selection function, namely the probability of breakage of the population balance model, were combined through a MATLAB code to predict the size distribution of the grinding products of quartz, marble, quartzite and metasandstone. The modeling results were in very good agreement with the particle size distributions obtained after grinding the feeds in a ball mill.


Introduction
Comminution is one of the most important processes of the mining industry and, more specifically, of the materials processing industry.However, it is an energy intensive process, especially at the stage of grinding, and consumes 3-4% of the electricity generated worldwide [1].In addition, a large share of this energy is absorbed by the device and only a small fraction, given that the efficiency of Semi-Autogenous Grinding (SAG) mill is about 15%, is used for comminution [2].Considering these factors, it is deduced that a slight increase in process efficiency may have a large impact on the operating cost of the mill, the environment and the conservation of resources.
Scientific research, which resulted in an increase of the grinding efficiency, started with the empirical relationships of Rittinger [3], Kick [4], Bond [5] and Walker and Show [6], which described the energy-particle size relationship.Then, theories which incorporated the concept of particle size distribution were developed [7,8], while in recent years research focused on the use of the mechanistic approach that was based on the recognition of physical events taking place during grinding [9].The mechanistic approach includes the matrix model, which considers grinding as a series of breakage events and the kinetic model, namely the population balance model, which considers grinding as a continuous process.Other researchers focused their efforts on the advantages of the Attainable Region (AR) approach [10][11][12], while recently Shi and Xie used the t 10 procedure to simulate batch [13] and continuous [14] grinding in a ball mill.
Mechanistic models consider that breakage of particles in a machine is a complex process, which cannot be analyzed by a single breakage event.Thus, it is necessary to split the process of comminution into individual sections and define the operations involved in each one.First, Epstein [15] considered that the process involves repetitive steps and each step can be described by a probability function, namely a selection or breakage rate function and a distribution function.Based on the concept of Epstein, first Broadbent and Callcott [16] considered the process of comminution in matrix form.In the matrix model, comminution involves a series of single breakage events, where the feed to each event is the product of the previous one [17].The longer the period of grinding, the greater the number of breakage events.The particle size distributions of the feed and products of such a process can be determined by dividing the material into predetermined n size classes.The particles in all size classes have a probability of being broken and this probability may change as the size of the particle changes.With this model, after a number of u breakage events, the comminution process can be described by the following Equation ( 1), p = (B where p and f are the product and the feed size distributions (1 × N vectors), B is the breakage function (N × N lower triangular matrix), S is the selection function (N × N diagonal matrix) and I is the identity matrix.The notation of the breakage function, the selection function and the identity matrix is shown in Figure 1, where b i,j is the mass fraction of the breakage product of size j that reports to a finer size i [18,19].
on the concept of Epstein, first Broadbent and Callcott [16] considered the process of comminution in matrix form.In the matrix model, comminution involves a series of single breakage events, where the feed to each event is the product of the previous one [17].The longer the period of grinding, the greater the number of breakage events.The particle size distributions of the feed and products of such a process can be determined by dividing the material into predetermined n size classes.The particles in all size classes have a probability of being broken and this probability may change as the size of the particle changes.With this model, after a number of u breakage events, the comminution process can be described by the following Equation ( 1), where p and f are the product and the feed size distributions (1 × N vectors), B is the breakage function (N × N lower triangular matrix), S is the selection function (N × N diagonal matrix) and I is the identity matrix.The notation of the breakage function, the selection function and the identity matrix is shown in Figure 1, where bi, j is the mass fraction of the breakage product of size j that reports to a finer size i [18,19].
where Bi, j is the mass fraction of the breakage product of size j which falls below size i.The perfect mixing model, which can be considered as a special case of the population balance model, has also been used to optimize and scale up several different grinding circuits [21,22].The basic concept of this model was developed by Whiten [23] and is the simplest case of multi-segment models, which consider only one segment [17].The Julius Kruttschnitt Mineral Center (JKMRC) used this model to simulate industrial autogenous (AG) and semi-autogenous (SAG) mills for over 30 years and attempts were made to empirically correlate the breakage rate with the operating conditions and design characteristics of the mills [24].
In recent decades, emphasis was given to the optimization of energy consumption in grinding mills using phenomenological kinetic models based on population balance considerations.The population balance model is a discrete-size, continuous-time model which is based on first order kinetics and uses two functions, namely the breakage rate (or selection function) Si and the breakage function bi, j [18,25].These functions define the fundamental size-mass balance equation for fully mixed batch grinding operations.Many researchers have underlined the advantages of these functions [26][27][28], while the scale-up from laboratory to industrial mills has also been discussed in a number of recent studies [1,[29][30][31][32].Furthermore, several studies investigated the effect of the variation of kinetic model parameters, when different mill operating conditions were used.These conditions include ball size [33][34][35][36], media shape [37,38], mill speed [28,39] and powder loading [40].
As far as the breakage rate Si (min −1 ) is concerned, Austin et al. [13] have proposed an empirical relationship (Equation (3)), which has been widely accepted [41][42][43], Broadbent and Callcott proposed a modified form of the Rosin-Rammler (RR) distribution [20] which can be used to describe the particle size distribution of the comminution products (Equation (2)), where B i,j is the mass fraction of the breakage product of size j which falls below size i.The perfect mixing model, which can be considered as a special case of the population balance model, has also been used to optimize and scale up several different grinding circuits [21,22].The basic concept of this model was developed by Whiten [23] and is the simplest case of multi-segment models, which consider only one segment [17].The Julius Kruttschnitt Mineral Center (JKMRC) used this model to simulate industrial autogenous (AG) and semi-autogenous (SAG) mills for over 30 years and attempts were made to empirically correlate the breakage rate with the operating conditions and design characteristics of the mills [24].
In recent decades, emphasis was given to the optimization of energy consumption in grinding mills using phenomenological kinetic models based on population balance considerations.The population balance model is a discrete-size, continuous-time model which is based on first order kinetics and uses two functions, namely the breakage rate (or selection function) S i and the breakage function b i,j [18,25].These functions define the fundamental size-mass balance equation for fully mixed batch grinding operations.Many researchers have underlined the advantages of these functions [26][27][28], while the scale-up from laboratory to industrial mills has also been discussed in a number of recent studies [1,[29][30][31][32].Furthermore, several studies investigated the effect of the variation of kinetic model parameters, when different mill operating conditions were used.These conditions include ball size [33][34][35][36], media shape [37,38], mill speed [28,39] and powder loading [40].
As far as the breakage rate S i (min −1 ) is concerned, Austin et al. [13] have proposed an empirical relationship (Equation (3)), which has been widely accepted [41][42][43], where x i (mm) is the upper size of class i, x 0 is the standard size (1 mm), α T is a parameter depending on milling conditions and α is a characteristic parameter depending on the material properties.More specifically, α T is the breakage rate for size x i = 1 mm and has the same units as S i (min −1 ).Q i is a correction factor, which is 1 for small particles (normal breakage) and less than 1 for large particles that need to be nipped and fractured by the grinding media (abnormal breakage).Q i is calculated by Equation ( 4), where µ is a parameter that depends on milling conditions and denotes the particle size when the correction factor is 0.5.Λ is a positive number that depends on the material type and shows how rapidly the breakage rate decreases as size increases (Λ ≥ 0).In addition, the particle size x m for which the breakage rate takes the maximum value is related with the parameter µ as follows [34,44], Equation (5) states that x m is directly proportional to µ considering that, for the same material, Λ and α are constants.
Kinetic and perfect mixing models received much more attention than the matrix model for simulating purposes.A noteworthy exception is the study of Deniz [45], in which a computer simulation of a matrix notation of the Broadbent and Callcott model was carried out and the obtained particle size distributions were compared with the experimental ones.Unlike the matrix model, the population balance model takes into consideration the abnormal breakage of large particles.In a kinetic model of this type there is an optimum size for which the breakage rate takes the maximum value, while above this size it decreases sharply (Equations ( 3) and ( 4)).Furthermore, it is mentioned that the matrix model is not directly related to the material type and the operating conditions (e.g., mill diameter, ball filling volume and mill speed).Thus, it is believed that the combined use of the matrix and population balance models will overcome a number of limitations of the traditional matrix model and result in a more realistic simulation of the grinding process.
The present study, through batch grinding experiments using quartz, marble, quartzite and metasandstone as test materials, aims to predict the particle size distributions of the products under different operating conditions.A MATLAB code was used to determine the kinetic parameters that minimize the error between the experimental Rosin-Rammler (RR) distributions and the model, for three different grinding times.Then, the same parameters were used to estimate the size distributions of different feed sizes.

Materials and Methods
The materials used in the present experimental study obtained from three different Greek sites, namely quartz from Assiros, near Thessaloniki, marble from west Crete, Chania area, and quartzite and metasandstone from west Crete, Kissamos area.The mineralogical analysis shows that marble mainly consists of calcite, while quartzite and metasandstone of quartz and some mica.Therefore, the chemical analysis shows that marble mainly consists of CaO present in calcite and dolomite, while quartzite and metasandstone of SiO 2 present in quartz.The quartz samples used were of high purity white quartz.Porosity measurements of the raw materials, obtained with the use of Archimedes method [46], indicated that quartz has near zero porosity, marble and quartzite less than 1%, and metasandstone almost 8%.Also, the density of metasandstone is significantly lower than that of the other materials (Table 1).The experiments were performed in a laboratory ball mill with a constant speed of n = 66 rpm (1.1 Hz) which is 70% of its critical speed.The mill charge consisted of stainless steel balls with density ρ b = 7.85 g/cm 3 .The parameters that describe the mass of grinding media J and material f c , respectively, are calculated from Equations ( 6) and (7).J is the ball filling ratio, which is the fraction of the mill filled by the media bed, while f c is the filling ratio of the material fed in the mill.

J =
Volume o f solid balls mill volume where Φ is the bed porosity of balls and material (assumed to be 40%) [18].The fraction of space between the balls at rest, U, that is filled with material bed (interstitial filling) can be calculated from Equation ( 8), The effect of two operating conditions, namely the interstitial filling U and the ball diameter d, were investigated by carrying out two series of tests (Table 2).In the first series three different U values (50%, 100% or 150%) were investigated, with f c constant at 4%, corresponding to 345 g of quartz or 351 g of marble.In the second series, balls with constant weight of almost 5.1 kg and three different diameters (40 mm, 25.4 mm or 12.7 mm) were used independently.The samples were homogenized by the cone and quarter method and 6 kg of each material type was used for the tests.Each material was crushed with the use of a jaw crusher to a size of −4 mm.The product of the crusher was wet screened at 150 µm, while the +150 µm fraction was dried and screened using a series of screens, with an aperture ratio √ 2 to obtain the feed fractions.As a result, five mono-sized fractions for each material (−3.35 + 2.36 mm, −1.7 + 1.18 mm, −0.850 + 0.600 mm, −0.425 + 0.300 mm and −0.212 + 0.150 mm) were obtained, while the natural products obtained after 0.5 min of grinding were considered as initial feed.Then, grinding tests were performed at various grinding times t (0.5, 1.5, and 3.5 min) corresponding to three different grinding events (1, 3 and 7), with a step of 0.5 min.The obtained products were wet sieved using a series of screens with an aperture ratio of √ 2 for the determination of particle size distributions.For the determination of the kinetic parameters and for modeling the process, a MATLAB code was developed using the following steps:  3) and ( 4)).Given that the selection function elements of the matrix model are pure numbers, the determined S i (min −1 ) values multiplied by the time interval dt of each individual event are assumed to be 0.5 min.

•
Determination of the product p size distribution after u breakage events, with the use of Equation ( 1).

•
Non-linear optimization of the objective function F through the "fmincon" function of the MATLAB software, and determination of the optimum α T , α, µ and Λ values that minimize the error between experimental size distributions RR and the model, for three different grinding times.The objective function F is defined with the use of Equation ( 9), where Q = loglog(100/(100 − P)) and P is the cumulative mass % finer than size x, as defined by the RR distribution [15].The exponents Mod and Exp refer to the experimental data and the estimated results of the model, respectively, while u and i are indicators that refer to the grinding steps and sizes classes, respectively.A flowchart of the MATLAB code is shown in Figure 2.

•
After the estimation of the optimal kinetic parameters for a given feed size of a raw material, the next steps are (i) to verify whether these parameters can be used for different feeds of the same material; and (ii) to assess the fits between the experimental and estimated size distributions.The same MATLAB code was used to evaluate the effects of the material type and different milling conditions on the estimated parameters.Screen views of the code and optimization results (experimental and estimated RR size distributions) are presented in Figures 3 and 4.

First Series of Tests -Effect of Different Interstitial Filling Values
The investigation of the effect of different interstitial fillings U on the simulation of the grinding process focused on the study of marble and quartz as test materials under the conditions given in Table 2.According to the experimental order, the −0.850 + 0.600 mm feed was used for the optimization of the objective function F and the estimation of the kinetic parameters α T , α, µ and Λ of the selection function.Then, the values of parameters were used for the estimation of the particle size distributions of the −3.35 + 2.36 mm, −1.7 + 1.18 mm and −0.425 + 0.300 mm feed fractions.Figures 5-7 show, as an example, the experimental and estimated size distributions of marble for each initial feed size and different interstitial filling values U (50%, 100% and 150%).From the results obtained it is seen that a very good agreement between experimental and estimated values, within the limits of experimental error, can be observed.The accuracy of the estimated distributions was determined by using the correlation coefficient R, as seen in Tables 3-6 for marble for each initial feed and U = 50%.The same conclusions can be drawn for quartz using different U values.

First Series of Tests -Effect of Different Interstitial Filling Values
The investigation of the effect of different interstitial fillings U on the simulation of the grinding process focused on the study of marble and quartz as test materials under the conditions given in Table 2.According to the experimental order, the −0.850 + 0.600 mm feed was used for the optimization of the objective function F and the estimation of the kinetic parameters αΤ, α, μ and Λ of the selection function.Then, the values of parameters were used for the estimation of the particle size distributions of the −3.35 + 2.36 mm, −1.7 + 1.18 mm and −0.425 + 0.300 mm feed fractions.Figures 5-7 show, as an example, the experimental and estimated size distributions of marble for each initial feed size and different interstitial filling values U (50%, 100% and 150%).From the results obtained it is seen that a very good agreement between experimental and estimated values, within the limits of experimental error, can be observed.The accuracy of the estimated distributions was determined by using the correlation coefficient R, as seen in Tables 3-6 for marble for each initial feed and U = 50%.The same conclusions can be drawn for quartz using different U values.(c) (d) The estimated parameters αΤ, α, μ and Λ of the selection function for marble and quartz for different U are shown in Table 7.These parameters enable the assessment of the relation of the matrix model with the material type and the change of interstitial filling U value.As far as the parameter αT (min −1 ) is concerned, it is seen that the corresponding values for marble are much higher than those for quartz indicating that marble is ground with a higher rate.Also, from the same table, it is seen that for the same material the value of αΤ increases when U decreases and this suggests that grinding is more effective at lower U values.However, a recent study has shown that there is an optimum value of U for which the breakage rate obtains the maximum value and grinding becomes more efficient [30].7.These parameters enable the assessment of the relation of the matrix model with the material type and the change of interstitial filling U value.As far as the parameter α T (min −1 ) is concerned, it is seen that the corresponding values for marble are much higher than those for quartz indicating that marble is ground with a higher rate.Also, from the same table, it is seen that for the same material the value of α T increases when U decreases and this suggests that grinding is more effective at lower U values.However, a recent study has shown that there is an optimum value of U for which the breakage rate obtains the maximum value and grinding becomes more efficient [30].

Second Series of Tests-Effect of Different Ball Diameter
The investigation of the effect of ball diameter on the simulation of the grinding process focused on quartzite and metasandstone as test materials under the conditions given in Table 2. Following the procedure described previously, the optimum values α T , α, µ and Λ of the −0.850 + 0.600 mm initial feed fraction were used for the prediction of the particle size distributions of the −3.35 + 2.36 mm, −1.7 + 1.18 mm and −0.425 + 0.300 mm initial feeds.Figures 8-10 present as an example, the experimental and estimated size distributions of quartzite for each initial feed size and different ball diameter d (12.7 mm, 25.4 mm and 40 mm).From the results obtained it is seen that in general a very good agreement between experimental and estimated values, within the limits of experimental error, can be observed.However, poorer fit of the experimental data can be seen when the coarse size, −3.35 + 2.36 mm, was used as feed.This may be due to the abnormal grinding behavior of large particles that need to be nipped and fractured by the grinding media, as also reported by other researchers [18,26,47].The experimental and estimated size distributions for quartzite at each step are shown in Tables 8-11.The accuracy of the estimated distributions was determined by using the correlation coefficient R. The same conclusions can be drawn for metasandstone when balls of different diameter d are used.

Second Series of Tests -Effect of Different Ball Diameter
The investigation of the effect of ball diameter on the simulation of the grinding process focused on quartzite and metasandstone as test materials under the conditions given in Table 2. Following the procedure described previously, the optimum values αΤ, α, μ and Λ of the −0.850 + 0.600 mm initial feed fraction were used for the prediction of the particle size distributions of the −3.35 + 2.36 mm, −1.7 + 1.18 mm and −0.425 + 0.300 mm initial feeds.Figures 8-10 present as an example, the experimental and estimated size distributions of quartzite for each initial feed size and different ball diameter d (12.7 mm, 25.4 mm and 40 mm).From the results obtained it is seen that in general a very good agreement between experimental and estimated values, within the limits of experimental error, can be observed.However, poorer fit of the experimental data can be seen when the coarse size, −3.35 + 2.36 mm, was used as feed.This may be due to the abnormal grinding behavior of large particles that need to be nipped and fractured by the grinding media, as also reported by other researchers [18,26,47].The experimental and estimated size distributions for quartzite at each step are shown in Tables 8-11.The accuracy of the estimated distributions was determined by using the correlation coefficient R. The same conclusions can be drawn for metasandstone when balls of different diameter d are used.The estimated parameters αΤ, α, μ and Λ of the selection function are shown in Table 12, for quartzite and metasandstone and for different balls diameter d.As far as the parameter αT (min −1 ) is concerned, it is seen that the corresponding values for metasandstone are much higher than the ones for quartzite indicating that metasandstone is ground with a higher rate.Also, from Table 12 it is seen that for the same material the value of αΤ increases when d decreases and this suggests that grinding becomes more effective when smaller balls are used.Regarding the effect of ball diameter d on μ parameter it is revealed that the value of μ increases when d increases as well.Since μ is directly proportional to the optimum size xm (Equation ( 4)), these findings indicate that small balls are more effective for grinding fine particles, whereas large balls are required to break the coarser ones [18,35].12, for quartzite and metasandstone and for different balls diameter d.As far as the parameter α T (min −1 ) is concerned, it is seen that the corresponding values for metasandstone are much higher than the ones for quartzite indicating that metasandstone is ground with a higher rate.Also, from Table 12 it is seen that for the same material the value of α T increases when d decreases and this suggests that grinding becomes more effective when smaller balls are used.Regarding the effect of ball diameter d on µ parameter it is revealed that the value of µ increases when d increases as well.Since µ is directly proportional to the optimum size x m (Equation ( 4)), these findings indicate that small balls are more effective for grinding fine particles, whereas large balls are required to break the coarser ones [18,35].

Conclusions
In this study, the grinding of quartz, marble, quartzite and metasandstone was simulated under different operating conditions, through the combined use of Broadbent and Callcott matrix and population balance models.Also, with the use of a MATLAB code, the breakage rate parameters (α T , α, µ, Λ) of an initial feed that minimize the error between the experimental Rosin-Rammler (RR) and the model distributions were determined for three different grinding times.Then, the same parameters were used to estimate the size distribution of other feed sizes under the same operating conditions.The results showed a very good agreement between experimental and estimated particle size distributions.
In light of the results obtained from the study of the breakage rate parameters for different interstitial filling U values, it can be concluded that for different α T values marble is ground with a higher rate when compared with quartz.On the other hand, the value of µ parameter seems to remain constant when U changes and this suggests that the optimum size x m at which the breakage rate obtains the maximum value is independent of U. Regarding the study of the effect of different ball diameter d on the grinding process simulation, it is shown that the value of α T increases with decreasing d and this suggests that grinding is more effective when smaller balls are used.Also, the findings of this study support the general concept that small balls are more efficient for grinding finer particles, whereas large balls are required to break coarser ones.Under the same grinding conditions (U = 50% and d = 25.4 mm) the breakage rate of the examined raw minerals follows the order: metasandstone > marble > quartzite > quartz.
Finally, it is deduced that the proposed approach offers the following advantages: • Reduces the number of tests required and therefore reduces the cost of grinding.

•
Predicts more accurately the distribution of grinding products by taking into account the effect of the type of the feed and the grinding conditions used.

•
It can be also applied to feeds with broader size distribution, in contrast with the population balance model, which is mainly applied for mono-sized fractions.

•
Can be easily applied to full-scale operations.

Figure 1 .
Figure 1.Notation of the breakage function, the selection function and the identity matrix.

Figure 1 .
Figure 1.Notation of the breakage function, the selection function and the identity matrix.

Figure 2 .
Figure 2. Flowchart of the MATLAB code.

Figure 3 .
Figure 3. Screen view of the MATLAB code.

Figure 4 .
Figure 4. Screen view of experimental and estimated RR particle size distributions.

Figure 2 .
Figure 2. Flowchart of the MATLAB code.

Figure 2 .
Figure 2. Flowchart of the MATLAB code.

Figure 3 .
Figure 3. Screen view of the MATLAB code.

Figure 4 .
Figure 4. Screen view of experimental and estimated RR particle size distributions.

Figure 3 . 17 Figure 2 .
Figure 3. Screen view of the MATLAB code.

Figure 3 .
Figure 3. Screen view of the MATLAB code.

Figure 4 .
Figure 4. Screen view of experimental and estimated RR particle size distributions.

Figure 4 .
Figure 4. Screen view of experimental and estimated RR particle size distributions.

Table 1 .
Porosity and density of the raw materials used.

Table 2 .
Grinding conditions used in the first and second series of tests.

Table 7 .
Selection function parameter values for different U.

Table 12 .
Selection function parameter values for different ball diameter.