Identification of Optimal Mill Operating Parameters during Grinding of Quartz with the Use of Population Balance Modeling

It is known that ball milling is an energy intensive process and great efforts have been made over the years to improve energy efficiency. The use of population balance models (PBMs) can assist in the design of mineral processing circuits and the scale-up of laboratory mill results to full-scale. However, since each model has its own capabilities and limitations, it is believed that a combined use will provide more accurate information for the reliable description of the process. In this study, the simulation of grinding of quartz is investigated in order to identify the optimal mill operating parameters. With the use of population balance modeling the specific rate of breakage and the cumulative breakage parameters can be determined from mono-size, short grinding time batch tests. The determined breakage parameters were back calculated by minimizing the residual error between experimental and reproduced size distributions. By combining two software packages the back calculated breakage parameters were used for the prediction of the optimum ball filling volume. The proposed procedure can be also applied for the identification of optimal mill operating parameters for other minerals.


Introduction
Size reduction is one of the most fundamental operations in the mineral processing industry. The primary purpose of size reduction is to obtain an appropriate product size with the least possible energy consumption. However, it is known that a large share of the energy used in comminution is absorbed by the machine and only a small part is consumed for size reduction. This implies that the existing grinding machines are inefficient and engineers have focused their efforts not only on the design of new and more efficient machines, but also on the development of new approaches for data processing (Yahyaei M. et al., 2015;Umucu Y. et al., 2016).
The main theories that initially described, evaluated and improved the relationship between specific energy and size reduction in comminution are those of Rittinger P.R. (1867), Kick F. (1885) and Bond F.C. (1952). These theories are based on a specific size modulus of the mate-rial but they do not take into consideration the distribution modulus. Later, Charles R.J. (1957) and Stamboliadis E. (2002) used different breakage models by assuming that both the initial feed and the final product size distributions are expressed by the Gates-Gaudin-Schuhmann (GGS) equation.
During the last decades, considerable work has been done on the optimization of energy consumption in grinding mills using phenomenological grinding kinetics models based on population balance considerations. Population balance modeling is based on first order kinetics and uses two functions, namely the specific rate of breakage S i and the breakage function b i,j (Herbst J.A. and Fuerstenau D.W., 1980;Austin L.G. et al., 1984). These functions provide the fundamental size-mass balance equation for fully mixed batch grinding operations. Many researchers have underlined the advantages of these functions (Ipek H. and Göktepe F., 2011;Wang X. et al., 2011;Gupta V.K. and Sharma S., 2014), while the scale-up from laboratory data to full-scale mills has also been discussed in a number of recent studies (Fuerstenau D.W. et al., 2003;Deniz V., 2013).
Population balance modeling is the most widely used approach and assumes that the specific rate of breakage is constant and is not affected by the grinding time (Austin L.G. and Bhatia V.K., 1972;Austin L.G. and Bagga P., 1981;Nomura S. et al., 1991). Austin L.G. and Bagga P. (1981) also have shown that the breakage function parameters can be estimated from the size analysis of the products obtained after a short grinding time, using as initial feed a narrow particle size fraction, namely the one size fraction BII method.
The primary concept of simulation balance modeling, which suggests that the breakage function parameters must be estimated before any of the particles are reselected for further breakage, may result in calculation errors (Austin L.G. et al., 1984;Celik M.S., 1988;Bozkurt V. and Özgür I., 2007). The deviation from first-order kinetics is also well known in several cases. Celik M.S. (1988) stated that the breakage rate of anthracite in ball mills exhibited an acceleration, particularly noticeable for long grinding times, due to extensive agglomeration of fines. The acceleration or deceleration of breakage rates of all individual size classes can be determined during wet grinding using the G-H method (Rajamani R.K. and Guo D., 1992). Bilgili E. and Scarlett B. (2005) proposed a phenomenological theory which explains non-first order effects as a result of multi-particle interactions. Instead of the forward experimental estimation of breakage model parameters, the back-calculation method which calculates the model parameters by minimizing the error between experimental and predicted size distributions, is also widely employed. This method uses all available data at the same time and enables reduction of the calculation error. It can be applied in cases where limited data is available and thus the need for carrying out several experiments is minimized. In this regard, Austin L.G. et al. (2007) stated that it is always necessary to confirm that the experimental breakage model parameters reproduce the size distributions and the only way to accomplish this is to use the back-calculation method.
For the design of grinding circuits in mineral processing plants, many equations describing mill power and grinding energy are available. These equations can be categorized into two major groups. The first group describes the relationship between specific energy requirements and size reduction, while the second one relates power draft to mill size. The first group can be represented by the Bond F.C. (1960) test method. This method is widely used both for evaluating the grindability of a material and for determining the power required to drive the mill. However, it uses data from laboratory tests while many empirical correction factors were proposed for the scale-up to larger mills or other operating conditions. With respect to the second group, Herbst J.A. and Fuerstenau D.W. (1973) found out experimentally that a proportional relationship exists between the specific rate of breakage and the specific net mill power for a tumbling mill, which is not affected by the operating parameters. Nomura S. et al. (1994) also confirmed this relationship.
Two important factors should be determined before carrying out experiments in a mill, namely the ball filling J and the powder filling f c . J is the ball filling volume, which is the fraction of the mill filled by the grinding media bed (Eqn. (1)), while f c is the powder filling volume in the mill (Eqn. (2)) where ε is the bed porosity of balls and powder which is assumed to be 40 % for both (Austin L.G. et al., 1984). The fraction of space between the balls at rest that is filled with powder is called interstitial filling U and is calculated by Eqn. (3).
It is mentioned that when U approaches 1 the impacts are increasingly cushioned by excess particles between grinding media. Considering that a proportional relationship exists between the specific rate of breakage and the specific net mill power, the following Eqns. (4)-(8) can be used to determine the specific rate of breakage S i (1/min), for each size class i obtained after grinding, as a function of J. Bond F.C. (1960) proposed the following relationship for wet grinding in overflow ball mills to scale-up results from laboratory tests to larger mills Beeck R. (1970) proposed that for dry grinding of cement the following relationships, based on data from German and U.S plants respectively, could be used. Shoji K. et al. (1982) studied the specific rate of breakage in relation to the ball filling volume for quartz grinding using two tumbling ball mills with different capacity. They proposed the following relationship for mill diameter 0.6 m 1.20 2.3 1 e 1 6.6 Austin L. G. et al. (1984), taking into consideration the differences between various relationships, proposed that for batch dry grinding the following relationship should be used (mentioned here as 'Bond Modified') 5 (1 0.937 ) (1 5.95 ) Later, many other research efforts were made to define the relationship between S i and J under different operating conditions. These studies focused on the determination of S i using the grinding kinetic model. More specifically, Deniz V. and Onur T. (2002) investigated the breakage kinetics of pumice samples in relation to the powder filling volume in a ball mill, at J = 20 %, and found that the maximum normal breakage rate occurred at an interstitial filling U of approximately 0.4. A more recent study by Fortsch D.S. (2006) showed that the reduction of J leads to an increase in capital and installation costs of the milling equipment. He proposed J = 35 % as optimal ball filling volume for cement plants. On the other hand, Metzger M.J. et al. (2009) using silica sand as a test material showed that the lower the ball filling J the higher is the amount of the produced product. They also found that more efficient breakage occurs at ball filling volume J = 1.5 %.
The present study through batch kinetic experiments with the use of quartz as test material aims to assess the effect of ball filling volume J on the specific rate of breakage S i and identify the cumulative breakage function B i,j by keeping the powder filling volume f c constant. A new methodology is proposed, involving the use of population balance modeling, for the identification of the optimum ball filling during quartz grinding and the scale-up of laboratory mill results to full-scale.

Theoretical background
Let's consider a mass of material M in a ball mill that after breakage needs to be divided, by using x i screens, into i + 1 narrow size classes. Normally, for a size class i bounded between two successive screens x i and x i+1 containing a mass fraction m i (t) at time t, it is assumed that breakage follows a first-order law. That is, the rate of breakage of size i is proportional to the mass of size i in the mill hold up. According to this, the specific rate of breakage S i is independent of time and since the mill hold up M is constant the following Eqns. (9), (10) are obtained (Austin L.G. and Luckie P.T., 1972; Klimpel R.R. and Austin L.G., 1977): where m i (t) and m i (0), are the mass fractions for size class i, at times t and 0 respectively. Thus, a log-linear plot of m i (t) versus t should result in a straight line and S i can be determined from the slope of this line. The specific rate of breakage depends on the particle size of the mill hold up. Austin L.G. et al. (1984) used the following Eqn. (11) to describe this relationship: where x i (mm) is the upper size of size 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 specific rate of breakage at size x i = 1 mm and has the same units as S i (1/min). Q i is a correction factor, which is 1 for small particles (normal breakage) and less than 1 for too large particles to be nipped and fractured by the grinding media (abnormal breakage). In the abnormal breakage region each material size behaves as it is composed of a weak and a strong fraction (Ipek H. and Göktepe F., 2011). Q i is calculated by Eqn. (12): where μ is a parameter depending 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 rates of breakage decrease as size increases (Λ ≥ 0). The evolution of the specific rate of breakage vs. size is shown in Fig. 1. When a material of size j breaks once, the mass fraction of the broken products with size i can be represented by the breakage function b i,j , i > j. Thus, b 2,1 is the mass fraction that is reduced from size 1 into size 2, and so on. The breakage function is usually represented in cumulative form B i,j as, Austin L.G. and Bagga P. (1981) have shown that for short grinding times the values of B i,j can be estimated from the size analysis of the products, using as initial feed a narrow particle size fraction j (the one size fraction BII method): where P i (t) is the mass fraction less than size x i at time t. The above Eqn. (14) presupposes that a small mass fraction of particles with smaller sizes will re-break. Experience suggests that good results are obtained when the selected time of grinding results in a 20 % to 30 % material that is broken out from size j (Austin L.G. et al., 1984).
The cumulative breakage function B i,j can also be represented by the following empirical Eqn. (15) which is the sum of the two straight lines on log-log plot (Austin L.G. and Luckie P.T., 1972): where x j is the top size, B i,j is the cumulative breakage function, γ is the slope of the lower straight line part of curve, β is the slope of the upper part and Φ j is the intercept, as shown in Fig. 2. γ also, characterizes the relative mass of fines produced after breakage and therefore is related to the grinding efficiency. Generally, the breakage parameters Φ j , γ and β depend on the properties of the material. Thus, the parameters of the cumulative breakage function, which are determined by laboratory experiments, can be used directly for larger scale mills.

Material and methods
The material used in the present experimental study is quartz obtained from the quarry of Assiros, in the vicinity of Thessaloniki, N. Greece, with S i O 2 content 99.1 % and density 2.65 g/cm 3 . 100 kg of the material was sampled from several locations to obtain representative samples. The microscopic examination using polished and thin sections showed that the material is compact, homogeneous, and does not show the presence of open cracks, discontinuities or alteration phenomena.
The experiments were performed in a laboratory ball mill that operated at a constant speed of Ν = 66 rpm (1.1 Hz) corresponding to 70 % of its critical speed ( Table  1). The mill charge consisted of 25.4 mm (1 inch) stainless steel balls with density ρ b = 7.85 g/cm 3 .
Three series of tests were performed at J 6.7 %, 10 % and 20 %, while f c that remains constant at 4 %, corresponds to 345 g of quartz.
The sample was homogenized by the cone and quarter method and 6 kg of the material was used for the tests. This quantity was crushed at a size of -4 mm using a jaw crusher. 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 monosized fractions of quartz (-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 prepared. The grinding tests were performed at various grinding times t (0.5, 1, 2 and 4 or 8 min). The products obtained from each test were wet sieved using a series of screens with a ratio of 2 for the determination of particle size distribution. All tests were carried out in duplicate and average values are given in the results section.
The mill power is calculated by the following formula (Stamboliadis E. et al., 2011):  where Μ is the total mass (kg) of the feed material and balls, Ν is the rotational speed (Hz) and D is the internal diameter of the mill (m). It should be also noted that the following issues were taken into consideration: i. A fresh feed charge was used in each test.
ii. Following the recommendation of Gupta V.K. and Sharma S. (2014), the crushed material was preground in the mill for 2 min in order to avoid abnormal breakage behaviour. iii. The feed charge was sampled and screen analyzed in the same way as the mill product. iv. The mill was layer loaded with balls and feed charge, following the standard practice (Gupta V.K. et al., 1985). The energy consumed by the mill is proportional to the time t and is given by Eqn. (17): where P is the power of the mill (W), t is the grinding time (s) and Ε is the net energy consumed (J). If m is the mass of the feed material (kg) then the specific energy for grinding ε' (J/kg) is given by Eqn. (18): The specific energy in kWh/t is calculated by dividing Eqn. (18) by 3,600.
Two software packages were used, namely Moly-Cop Tools TM version 1.0 and MODSIM TM academic version 3.6.25 (King, 2001), to identify breakage parameters and model the process.
Moly-Cop Tools TM includes a set of easy-to-use excel spreadsheets, that can be used to characterize and evaluate the operating efficiency of any given mineral ore grinding circuit, following standardized methodologies and widely accepted evaluation criteria. More specifically, the spreadsheet (BALLPARAM_BACH) uses all the differential equations of a population balance model for the calculation of S i and B i,j parameters. The software identifies the best combination of parameters that minimize the residual error between experimental and reproduced size distributions, while the objective function SSE is expressed by the following formula: where R i(exp) represents the particle size distribution of the mill product (as % retained on screen x i ), R i(mod) represents the response of each corresponding R i(exp) for a given set of model parameters, n is the number of screen intervals used, w i is the weight factor for each screen that quantifies the relative quality and reliability of each particular screen size value with respect to the other values (high values indicate more reliable measurements) and W is the total sum of the weight factors. The software input requires a minimum of four size distributions obtained at four grinding times t 1 , t 2 , t 3 and t 4 . For the back calculation of the breakage parameters, the experimentally calculated parameters (α T , α, μ, Λ, Φ j , γ and β) are used as initial input parameters in the Moly-Cop Tools TM software. As a rule of thumb, the parameters which are generally regarded as characteristics of the material should be kept constant and only those associated with mill conditions should be changed. In this study, the input parameters are α, μ, Λ, Φ j , γ and β and the initial value of α Τ is the one calculated to give the best fit between experimental and reproduced size distributions. After this, the fit is improved by adjusting the parameters μ, Λ, Φ j , and γ, while the parameters that depend on the material, i.e., α and β are kept constant.
MODSIM TM is a software that uses simulation models for the design of mineral processing circuits. The main unit operations include size reduction, crushing and grinding, classification and concentration based on particle mineralogical composition and solid-liquid separations. Furthermore, MODSIM TM which can up-scale laboratory mill results to full-scale mill, provides standard grinding mill units, e.g. MILL and GMSU based on population balance modeling, and thus the calculated experimental breakage parameters can be used as input.
The back calculated breakage parameters from the Moly-Cop Tools TM software are used as input in the unit model MILL of the MODSIM TM software and the distributions of the products are compared. MILL is the ball milling model that is based on standard Austin's models (Eqns. (14) and (15)) and takes into consideration the breakage parameters as does the Moly-Cop. It is mentioned that the MILL model does not take into consideration the operating conditions (e.g. mill diameter, ball filling volume and mill speed) and a different unit model should be used to elucidate the ball filling effect on the grinding process. For this purpose, the Austin's scaled-up unit model GMSU of MODSIM TM was used in this study. The back calculated breakage and other parameters, e.g. D = 0.204 m, N = 66 rpm, d = 25.4 mm were used as input in the GMSU model. Size distributions varied in relation with J. A total of six J values were examined and each of the obtained distributions was used as input in the Moly-Cop Tools TM software. Then, the optimum value of α Τ (1/ min) that minimizes the least-squares objective function (Eqn. (19)) was determined by keeping all other parameters constant. Fig. 3 shows the relationship between the remaining mass fraction of feed size and grinding time for J = 20 % using a normal-log plot. The feed size is the upper size of the tested size class. The results indicate that breakage follows a first order law (i.e. the specific rate of breakage S i is independent of time) and the specific rate of breakage can be determined from the slope of the straight lines. Each line corresponds to a different feed size tested. This figure also shows an increase of S i values when the feed size fraction is increased, except for the largest tested size (-3.35+2.36 mm), which exhibits lower breakage rate compared to the smaller fraction (-1.7+1.18 mm). This reveals that there is an increase in the rate of breakage up to a certain particle size and that grinding media cannot break coarser material efficiently. The same trend is obtained for J = 10 % and J = 6.7 %.

Determination of the specific rate of breakage
Regarding the effect of J on breakage rate S i , Fig. 4 shows that the value of S i increases with increasing J suggesting that breakage is more effective at higher J values.
Variations of specific rates of breakage S i for different feed size fractions are shown in Fig. 5 on a log-log scale.
It is shown that S i increases up to 1.7 mm (upper size class), but above this size breakage rates decrease sharply for all J values, since the particles are too large to be nipped and fractured by the grinding media used. It appears that the increase of S i is almost linear up to the optimum size of 1.7 mm and the slope of this linear part remains constant. In addition, it is also seen that the optimum feed size to obtain a maximum S i is independent of the ball filling volume J. It is believed that this size depends on the material type, the size of the balls and the type of the mill used.
The determined S i values for each feed size were used in Eqns. (11) and (12) and the parameters α T , α, μ, and Λ were calculated by using a non-linear regression technique. The calculated values are shown in Table 2.     finer than the optimum size (1.7 mm), while as described previously breakage takes place at the normal region. By definition, the values of B i,j are determined from the size distributions obtained at short grinding times, by fitting the experimental data to Eqn. (14). The Φ j , γ and β parameters were determined using Eqn. (15) and are shown in Table 2. Since β is generally accepted as characteristic constant of the material, it was decided to keep it constant at the value of 5.80, as indicated by Austin L.G. et al. (1984). This was considered necessary in order to reduce the number of parameters used and thus to investigate the sensitivity of only two parameters (Φ j , and γ). From the results obtained it is observed that changes in ball filling volume do not entail any significant changes in the values of parameters Φ j , and γ. This observation is in accordance with earlier studies (Fuerstenau D.W. et al., 2004;Ozkan A. et al., 2009;Chimwani N. et al., 2013), which also state that the B i,j parameters are independent of the mill operating conditions.

Back-calculation of the breakage parameters
From the calculated breakage parameters S i and B i,j ( Table 2) the size distribution of the grinding products can be reproduced using the spreadsheet BALLPARAM_ BACH of Moly-Cop Tools TM . Fig. 7 presents the reproduced values of cumulative mass passing for the (-0.850+0.600 mm) feed size after 0.5 min grinding time. The reproduction provides a coarser size distribution compared to the experimental values and this is also observed at longer grinding times. It is believed that the deviation between experimental and reproduced values in Fig. 7 is due to the value of α, which is a characteristic parameter of the material. In order to overcome this shortcoming, a readjustment of the grinding parameters is required. For quartz, the value of α proposed in literature is 0.80 (e.g. Austin L.G. et al., 1984), which is quite different from the calculated values (about 1.15) shown in Table 2. On the contrary, α T is a parameter depending on mill operating conditions and characterizes the grinding process. Thus, the variability of this parameter and the magnitude of its effect on the balling filling should be investigated.
In this study, the input parameters are α = 0.80, μ = 2.11, Λ = 2.90, Φ j = 0.76, γ = 1.02 and β = 5.80 and the initial value of α Τ is the one calculated to enable the best fit between the experimental and the reproduced size distributions. The fit can be improved by readjusting the values of the parameters μ, Λ, Φ j , and γ, while the parameters that depend on the material, i.e., α and β are kept constant. Table 3 summarizes the best combination of the back calculated parameters. These values are closer to those established in literature by other researchers (Yekeler M. et al. 2001) who obtained an average value of γ = 1.20 during wet grinding of quartz.

Simulation of the grinding process and validation of results
It was shown earlier that the value of S i increases with increasing J and this suggests that breakage is more effective as the ball mill filling volume increases. This conclusion is supported by three series of tests carried out with different J values 6.7 %, 10 % and 20 %.
In order to predict the effect of ball filling volume both the MODSIM TM and the Moly-Cop Tools TM software packages were used. According to this procedure, the back calculated breakage parameters given in Table 3 and other parameters shown in Table 1, i.e., D = 0.204 m, N = 66 rpm, d = 25.4 mm were used as initial input in the unit model GMSU of MODSIM TM . Fig. 9 shows the particle size distributions obtained from the model GMSU.
The results indicate an agreement between the experimental and the simulated particle size distributions for the tested ball filling volume J = 20 %, since both software packages use the standard Austin's kinetic models.
In order to simulate the grinding process, six J values (5 %, 10 %, 20 %, 30 %, 40 % and 50 %) were used and each of the obtained, from the unit model GMSU, size distributions was used as input in Moly-Cop Tools TM . Then, the optimum value of α Τ (1/min) that minimizes the least-squares objective function (Eqn. (19)) was obtained, by keeping all other parameters constant. Fig. 10 shows very good agreement between the experimental and the simulated α Τ values in relation to the ball filling volume J. However, the simulation indicated that there is an optimum ball filling volume (J = 30 %), for which the specific breakage rate obtains a maximum value. The predicted α Τ values were validated by additional experiments for J = 40 % and the results are given in Fig. 11. It can be seen from this figure that even above the optimum ball filling volume there is perfect match between experimental and simulated results, indicating that the approach used was able to predict accurately the variation between the α Τ values and the ball filling J.
Having successfully validated the simulated results, we attempted to compare them with those obtained from other earlier studies (Eqns. (4)-(8)). Fig. 12 presents the   predicted specific rates of breakage in relation to the ball filling volume J. In order to effectively compare the specific rates of breakage, as predicted from different models, the results are expressed as the ratio between S i value and a reference value S ref , taken as 50 % of the ball filling volume J. The approach proposed in the present study is indicated with the dashed line. It is clearly shown that Eqns. (4), (5), (6) and (8) are inversely proportional to the ball filling volume J and this is valid also for J values higher than the optimum value of about J = 30 %, as predicted in the present study. Eqn. (7) has also an optimum value J = 40 %, although the specific values of predicted breakage rates are lower than the ones proposed in this study. However, it can be said that Eqns. (4), (5), (6) and (8) are applicable for higher values of J (in the range of 30 % to 40 %), while the results derived by Shoji K. et al. (1982) (Eqn. (7)) and the present study are applicable over a wider range of ball filling volume J. The minor difference between the latter two studies may be due to different milling conditions and the material type.

Conclusions
In this study the grinding of quartz was simulated with the use of population balance modeling. The specific rate of breakage and the cumulative breakage parameters were determined from mono-size, short grinding time batch tests using a non-linear regression approach. The experimentally determined parameters (α T , α, μ, Λ, Φ j , γ and β) were compared for three J values (6.7 %, 10 % and 20 %). It is observed that the breakage rate increases with increasing J, suggesting that grinding in this case becomes more efficient. Grinding of different feed sizes shows that the specific rate of breakage increases up to the size of 1.7 mm, while for coarser sizes it decreases sharply. In addition, the change of J has minor effect on the parameters of the breakage function B i,j . The Moly-Cop Tools TM software was used to identify the best combination of the breakage parameters (α T , α, μ, Λ, Φ j , γ and β), for which the reproduced size distributions compare well with the experimentally determined ones. It was revealed that the experimental breakage parameters did not predict very well the size distributions and in order to improve the results a readjustment of the value of α, which is a characteristic parameter of the material, was required. The results of this study propose that for the more accurate prediction of the size distributions during grinding, the parameters which characterize the material (i.e. α and β) should be kept constant and only those associated with mill conditions should be modified accordingly.
The simulation of grinding of quartz was enabled by combining MODSIM TM and Moly-Cop Tools TM software packages. Using the back-calculated breakage parameters, a good match was observed between the experimental and the simulated α Τ values, in relation to J. The simulated results indicated that efficient grinding occurs at J = 30 %, for which the specific rate S i is maximized. Validation of the results was also done for a higher J value, 40 %.
The approach proposed provides more accurate information for the reliable description of the grinding process and enables the scale-up of laboratory results to larger scale mills. Future work will focus on the evaluation and exploration of other important mill operating parameters, e.g. ball size and mill speed.