Reliability of the Non-linear Modeling in Predicting the Size Distribution of the Grinding Products Under Different Operating Conditions

During the modeling of grinding systems, population balance modeling (PBM) which considers a constant breakage rate has been widely used over the past years. However, in some cases, PBM exhibited some limitations, and time-dependent approaches have been developed. Recently, a non-linear framework which considers the traditional linear theory of the PBM as a partial case was introduced, thus allowing the estimation of product particle size distribution in relation to grinding time or the specific energy input to the mill. In the proposed model the simplified form of the fundamental batch grinding equation was transformed into the well-known Rosin–Rammler (RR) distribution. Besides, the adaptability and reliability of the prediction model are among others dependent upon the operating conditions of the mill and the adjustment of the RR distribution to the experimental data. In this study, a series of grinding tests were performed using marble as test material, and the adaptability of the non-linear model was investigated using three loads of single size media, i.e., 40, 25.4, and 12.7 mm. The results indicate that the proposed model enables a more accurate analysis of grinding, compared to PBM, for different operating conditions.


Introduction
The main objective of grinding in mineral processing is to produce a desired product size and liberate the mineral(s) of interest from the gangue, so that they can be separated from each other using several physical methods. Grinding is most frequently performed in rotating cylindrical steel vessels known as tumbling mills and can be achieved by several mechanisms including impact or compression and attrition. These mechanisms deform particles beyond certain limits determined by their degree of elasticity and cause them to break. In these systems, the combined action of repeated impact and attrition over time causes size reduction. It has been found, however, that impact forces mostly reduce coarser particle sizes, whereas finer particle sizes are primarily reduced by attrition [1]. According to estimates, 3-4% of the global electrical energy and almost 50% of the mining energy consumption are used for mineral comminution [2][3][4][5][6]. Also, grinding is an inefficient process and is characterized by high CO 2 emissions and increased processing cost [7]. The European Union has set targets to reduce the greenhouse gas (GHG) emissions by 40% and 80-90% in 2030 and 2050 respectively, compared to the 1990 levels. Especially for the mining sector, this target can be achieved through the implementation of innovative and efficient ore processing technologies and the adoption of circular economy principles [8].
Over the past decades, several attempts were made to improve grinding efficiency and reduce energy requirements during size reduction. In this respect, Walker et al. [9] proposed a general relationship (Eq. (1)) where each of the existing comminution theories, namely, Rittinger [10], Kick [11], and Bond [12], can be represented as a partial case.
where dε is the infinitesimal specific energy (energy per unit mass) required to reduce by dx the size of a particle with size x, while C is a constant related to the material type, and m is a constant indicating the order of the process. If exponent m is replaced by 2, 1, or 1.5, then Eq. (1) leads to the relationships of Rittinger, Kick, and Bond, respectively. First, Charles [13] and later Stamboliadis [14] extended the existing theories of comminution and proposed Eq. (2) to calculate the specific energy ε required to obtain a particulate material described by the Gates-Gaudin-Schuhmann (GGS) distribution.
where k and x G are the distribution and size modulus (100% passing-screen size), respectively, as derived by the GGS distribution (Eq. (3)).
where P x is the mass (in %) finer than size x.
It is noted that if Rosin-Rammler's (RR) Eq. (4) is used to describe the particle size distribution (PSD) of the grinding products, and following Eq. (5) can be used instead of Eq. (2).
where x R and n are the particle size modulus (63.2% passingscreen size) and the distribution modulus (index of uniformity), respectively. Also, it has been shown that the exponents (1−m and 1−m´) of Eqs. (2) and (5) can be considered almost equal [15].
Other equations proposed in the literature to describe the energy-size relationship include Hukki's [16] equation, where the exponent (1−m) is dependent on the size of the comminuted material and the degree of fineness, and Morrell's [17] equation which is an alternative form of Bond's equation. Furthermore, Kapur [18] verified with experimental data that the particle size distributions of comminuted materials acquire self-similar character, indicating that distribution curves corresponding to different grinding periods collapse into a single spectrum when the cumulative finer mass is plotted as a function of dimensionless size. Gupta [19] validated the self-similarity approach based on batch grinding distribution data; however, Bilgili [20] pointed out that the self-similarity does not exist if grinding exhibits nonfirst-order behavior.
In recent years, kinetic models deriving from population balance considerations have been developed in order to optimize energy consumption in grinding circuits. The population balance model (PBM) provides the fundamental size-mass balance equation for fully mixed batch grinding operations, and several studies highlighted its advantages for the design, optimization, and control of grinding circuits [21][22][23][24]. The PBM, also referred to as the linear timeinvariant (LTINV) model, assumes that the breakage rate (or selection function), which is given by Εq. (6), does not vary with grinding time [25,26], where m i (t) is the mass fraction for size class i at time t, and S i is the breakage rate function of size class i.
The integration of Eq. (6) results in Eq. (7) which shows that the log-linear plot of m i (t) versus t gives a straight line, the slope of which determines the value of S i .
Austin et al. [15] proposed the following empirical relationship (Eq. (8)) of the breakage rate with respect to particle size x of size class i, which has been adopted by many researchers [23,27,28], where α T and α are model parameters depending on the milling conditions and the material type, respectively, while q x is a correction factor which defines the region of breakage; q x equals 1 in the normal breakage region, whereas q x <1 in the abnormal region. It is noted that abnormal breakage in ball mills is defined as a deviation from the first-order kinetics and occurs mainly for particles present in the feed that are too large to be properly nipped by the media. In this case, the breakage rate drops steadily and tends to zero.
Even though the linear theory of PBM has been widely used to analyze, control, and optimize grinding processes, several researchers critically addressed its limitations. In light of this, the non-linear population balance modeling (also referred to as the linear time-varying (LTVAR) model) was introduced to explain the deviations from the linear kinetic approach. Acceleration or deceleration of the breakage rate has been experimentally observed in both dry and wet grinding systems, the degree of which depends on the material type and operating conditions [29][30][31][32].
Recently, Petrakis and Komnitsas [33] developed a nonlinear framework for the prediction of the particle size distribution of the grinding products, where the linear theory of PBM was considered as a partial case. Based on this approach, the fundamental batch grinding mass balance equation was transformed into the well-known Rosin-Rammler (RR) distribution (Eq. (4)), thus allowing the determination of the breakage rate parameter α Τ as a function of the specific energy ε consumed for size reduction, as shown in Eq. (9), where M p and M are the power of the mill and the total mass of the feed material, respectively, and b (equals 1−m´) is the exponent of Eq. (5). In Eq. (9), n/b is the acceleration-deceleration parameter; for n/b>1, acceleration of the breakage occurs, and for n/b<1, the breakage rate decelerates, while for n/b=1, the breakage rate remains constant during grinding and corresponds to the linear timeinvariant (LTINV) model. Also, the magnitude of n/b defines the degree of acceleration or deceleration of the breakage rate. By substituting α T from Eq. (9) in Eq. (8) and setting C´=M p /(M·C R n/b ), the breakage rate S x is obtained as a function of specific energy input ε in the mill, as seen in Eq. (10), Finally, the following Eq. (11) is proposed to determine the particle size distribution (in the form of cumulative percent undersize), P x (ε), as a function of specific energy input for small particles in the normal breakage region.
The objective of the present study is to validate the non-linear framework as its accuracy depends on the mill operating conditions, the modeling assumptions, and/or the regression method of the RR distribution that fits the experimental data. For this, several grinding tests were carried out using marble as test material, and the effect of the media (balls) size on the grinding efficiency was investigated through the use of the non-linear model. In each series, the accuracy of the RR distribution model was determined by employing well-known metrics, such as the correlation coefficient R 2 values [34,35], the root-mean-square error (RMSE), and the modified index of agreement (IoA′) [36] using different regression methods.
The root-mean-square error (RMSE) (Eq. (12)) is a frequently used measure of the difference between values predicted by a model and the values actually observed. These individual differences are also called residuals, and the RMSE serves to aggregate them into a single measure of predictive power.
The modified index of agreement (IoA′) (Eq. (13)) ranges between 0 and 1 and is a standardized measure of the mean error and expresses the agreement directly; the optimal value is 1.
The assessed measures contain variables that denote the following features: Ο i and P i are the observed data and predictions respectively, _ O i and _ P i are the mean of the observed data and predictions respectively, while N is the number of the observed-predicted pairs.
Their effect on the adaptability and reliability of the nonlinear model to predict the product PSD was also investigated. The limitations of the traditional linear model for batch grinding are critically discussed, and suggestions are made for a more accurate analysis of the grinding operation. It is also mentioned that the present study has a certain degree of novelty since it presents, critically analyses, and validates a non-linear grinding model that can be used as a simple and reliable tool to predict the size distribution of the grinding products.

Materials and Methods
The material used in the present study is marble with a density of 2.7 g/cm 3 and low porosity (~0.3%) obtained from a quarry in west Crete, Chania region, Greece. The mineralogical and chemical analyses were carried out using X-ray diffraction and X-ray fluorescence techniques, respectively. The results indicated that marble consists of calcite (95%) and dolomite (4%), while its main chemical composition (in the form of oxides) is (wt.%) CaO 53.6, SiO 2 1.2, and Al 2 O 3 1.4. The calculated loss on ignition (LOI) is 43.5% after heating the material at 1050 °C for 4 h.
Three series of grinding tests were performed using a laboratory-scale ball mill (Sepor, Los Angeles, CA, USA) with a volume of 5423 cm 3 operating at 66 rpm, which corresponds to 70% of its critical speed ( Table 1). The ball charge consisted of balls with various sizes and a density of 7.85 g/cm 3 . Parameters J and f c , expressed as ball and material filling volume, were kept constant at 20 and 4%, respectively. As a result, the fraction of the space between the balls at rest that is filled with material (interstitial filling) U was 50%. The power of the mill was calculated using the formula presented in a previous paper [37].
In this test series, three loads with the same mass of single size balls, i.e., 40, 25.4, and 12.7 mm in diameter, were used. Four mono-size fractions of marble, i.e., −3.35 + 2.36 mm, −1.70 + 1.18 mm, −0.850 + 0.600 mm, and −0.425 + 0.300 mm, were prepared for the tests, and each one was ground in the mill for various grinding times t (0.5, 1, 2, 4 min). The products obtained after each grinding test were wet-sieved using a series of screens for the determination of particle size distribution. It is underlined that before preparing the feed material for the mill, the crushed product (less than 4 mm) was first pre-ground in the ball mill for 2 min, as recommended by Gupta and Sharma [22], in order to avoid abnormal behavior during the initial grinding period. All grinding tests were carried out in duplicate, and the values given in the paper for all parameters are average. It is mentioned that the variation of measurements was in all cases around ±2.3%.
The analysis of the grinding data obtained for different operating conditions is carried out using the linear theory of PBM and the one-size fraction BII method proposed by Austin et al. [15], along with the non-linear model proposed by Petrakis and Komnitsas [33].
The cumulative undersize distribution of the grinding products is simulated by the RR distribution. In light of this, two regression methods are used in order to define which one best describes the experimental data and to assess its influence on the accuracy of the predicted non-linear model. These methods involve (i) linear regression that provides the best-fit straight line when particle size distribution data are plotted on loglog(100/(100−P x )) versus logx (linear form of Eq. (4)) and (ii) non-linear regression using the Solver tool of Microsoft Excel that minimizes the sum of squared error between experimental and estimated particle size distributions.

Simulation of the Product Size Distribution
As mentioned, the more reliable RR distribution is selected for the description of the particle size distribution of the mill products obtained after different grinding times (or specific energy inputs) [7,[38][39][40]. Besides, the proposed non-linear model is used with the parameters acquired by fitting the RR distribution to the experimental data, while the goodness of fit defines the accuracy of the model. Figure 1a-c, which is used as an example, shows the RR plots of particle size distributions for the −3.35+2.36 mm feed fraction at different grinding times using linear regression analysis, when three loads of single size balls are used, i.e., 40 mm, 25.4 mm, and 12.7 mm. The results indicate that even though in some cases the PSD data points deviate from the straight line, especially for the coarse particles, the RR distribution can describe with high reliability the PSD across the entire range of particle sizes. Regarding the effect of media size on the grinding of the same coarse fraction, Fig. 1d presents the RR plots after 4 min of grinding for the three different ball sizes used. It is seen that for the coarse feed fraction, −3.35 + 2.36 mm, larger balls are required for efficient grinding, while the use of smaller balls, i.e., 12.7 mm, results in a much coarser product. This is in line with the general principle that larger grinding media are required for the efficient breakage of coarser particles, while smaller media having larger surface area break more efficient finer particles [41]. The RR plots (not shown) for the smaller feed fractions of marble confirm this behavior.
The correlation coefficients (R 2 ) were also calculated in order to assess the accuracy of the RR distribution when linear regression analysis was used ( Table 2). The results confirm that the straight lines of RR distribution describe the particle size distributions of the available data set well, and the R 2 values range from 0.962 to 0.998. Table 2 also presents the variation of R 2 values with grinding time for each feed fraction tested, when three loads of single size balls were used, using non-linear regression analysis. It is revealed that in this case, higher R 2 values (ranging from 0.968 to 0.999) are obtained. Since the predicted non-linear model uses the RR distribution to simulate the particle size distribution of the grinding products, it is not clear at this stage which is the most appropriate regression method that should be used to improve the accuracy of  The application of the non-linear regression to the RR distribution data enables the calculation of the distribution parameters, i.e., particle size modulus x R (63.2% passing-screen size) and distribution modulus n (index of uniformity). Table 3 presents the variation of the distribution parameters with grinding time (or specific energy) for the three loads of single size balls used. The results show that the operating conditions, e.g., grinding media size, the amount of energy required per mass of material (specific energy, ε), and the feed size affect the particle size distribution. In general, for a given load of single size balls, the n values increase with decreasing feed fraction indicating that the product size distribution becomes narrower when finer particle sizes are fed to the mill. The influence of ball size on the n values depends on the feed fraction tested. In particular, for coarse fractions, there is a tendency for n values to increase with increasing ball size (e.g., the average n values increase from 0.68 to 1.15 for the −1.70 + 0.850 mm feed fraction), while this trend becomes less evident for the −0.425 + 0.212 mm fraction (a slight increase from 0.96 to 1.07 was observed). It is noted that n values are incorporated in the predicted non-linear model (as shown in Eqs. (9)-(11)); thus, if the variation is known, it is possible to foresee how the model evolves over time under different milling conditions. Figure 2a-d shows the specific energy required versus the size modulus x R (values are shown in Table 3) on loglog scale for the four feed fractions tested when three loads of single size balls are used. It is observed from the trendlines that the energy required to reduce each feed fraction and produce a particulate material with a size modulus x R is described by the Charles relationship (Eq. (5)); thus, the exponent 1−m´ is determined from the slope of the straight line. It is reminded that the exponent 1−m´ was set equal to the parameter b used to determine (i) the breakage rate parameter α Τ (Eq. (9)), (ii) the breakage rate S x (Eq. (10)), and (iii) the particle size distribution P x (Eq. (11)), as a function of specific energy input. Figure 2a clearly shows that from the energy-saving point of view, the 25.4-mm diameter balls improve energy efficiency for the coarse feed fraction −3.35 + 2.36 mm, whereas the process becomes inefficient when smaller or larger balls are used. This finding is consistent with the results of previous studies which indicate that there is an optimum grinding media size that maximizes the ball milling efficiency [37,42,43]. The optimum ball diameter depends on, among others, the feed-to-product size ratio, the mill dimensions, and the breakage parameters. Also, the results show that (Fig. 2b) as the feed size decreases, the largest balls tested (40 mm) exhibit the lowest energy efficiency, while for even finer feed fractions (Fig. 2c, d), smaller balls (25.4 or 12.7 mm) are required for efficient grinding.  7) to the experimental data and according to the linear theory of PBM. The breakage rate S i can be determined from the slope of the straight line under the assumption that grinding follows first-order kinetics and the rate remains constant during grinding. In addition, by accepting the previous assumption, these results enable the elucidation of the effect of grinding media size on marble grinding kinetics. It is observed that the coarse feed fraction −3.35 + 2.36 mm exhibits a higher breakage rate with the use of 40-mm or 25.4-mm balls, while its breakage rate is significantly lower when smaller balls (12.7 mm) are used (see also Table 4). As seen in Fig. 2a, the energy efficiency for this size fraction is improved with the use of 25.4-mm diameter balls, indicating that the specific energy rather than the grinding time should be used in breakage kinetics approaches for the design, control, and optimization of grinding systems. Similar results, pertinent to the effect of ball diameter on the breakage rate of finer feed fractions are obtained.

Grinding Kinetics
With the assumption that marble grinding exhibits a first-order behavior under different operating conditions, the model parameters α T and α, along with the breakage rate S x with respect to particle size x of size class i, can be determined using Eq. (8). From Table 4, which shows the S x values and model parameters when three loads of single size balls are used, it is observed that the highest breakage rate S x for the feed size fractions tested, i.e., −3.35 + 2.36 mm, −1.70 + 1.18 mm, −0.850 + 0.600 mm, and −0.425 + 0.300 mm, is 2.04, 1.56, 0.95, and 0.49 min −1 . The results indicate that the coarse, intermediate, and fine fractions exhibit higher breakage rates when 40-mm, 25.4-mm, and 12.7-mm balls are used, respectively. The characteristic parameter α, which depends on the material properties, remains constant (0.90), while the parameter α T , which is the (normal) rate of breakage at feed size 1 mm, varies with ball diameter (ranging from 0.73 to 1.20).
The major limitation of the linear theory of PBM is the use of the slope of the best-fit straight line of the first-order plots to obtain the breakage rate S x from experimental data. Several breakage rate values have been experimentally observed during grinding; thus, the values that can be used in large-scale mills cannot be accurately estimated [44]. In light of this, Table 5 shows the S x values of feed sizes tested for various grinding time intervals in the three series of tests, assuming that at each interval, grinding exhibits first-order kinetics. It is obvious that S x values vary with grinding, and therefore, the values determined by fitting a straight line to the data points ( are questionable. In light of this, Petrakis and Komnitsas [33] reported that the most appropriate way to obtain reliable results is the use of a relationship involving the breakage rate S x and the specific energy ε consumed for size reduction, such as Eq. (10). This is discussed in the following section.

Evolution of the Breakage Rate with Specific Energy Input
Having discussed in detail the drawbacks associated with the use of first-order kinetics to accurately determine the breakage rate parameters, this section indicates a reliable approach to obtain meaningful S x values. In this context, Fig. 4a-c shows the evolution of breakage rate S x as a function of specific energy at each feed fraction tested with the use of the three loads of single size balls, as derived from Eq. (10) of the proposed non-linear model. It is mentioned again that the incorporation of specific energy rather than grinding time has the advantage of extrapolating the obtained data to larger-scale mills; the S x values are independent of mill dimensions and operating variables in the normal grinding region [22,45]. The results indicate that, contrary to the linear theory of PBM, S x varies with specific energy input, and the trend of variation (acceleration or deceleration) depends for each feed fraction on the diameter of the balls used. In particular, for the coarse feed fraction tested, −3.35 + 2.36 mm, with the use of 40-mm balls, a remarkable acceleration of S x during the initial grinding period (up to 1.16 kWh/t corresponding to 2-min grinding time) is observed, which then remains constant at almost 2.5 min −1 ; a maximum S x value of 2.47 min −1 was observed. When 25.4-mm balls are used, S x accelerates gradually until the end of the grinding period selected (4-min grinding time) and reaches 2.34 min −1 . On the other hand, the use of smaller balls (diameter 12.7 mm) causes gradual deceleration of S x for the coarse fraction tested. A significant acceleration of the S x for the −1.70 + 1.18 mm feed fraction is observed during the grinding period when 40-mm balls are used, and the S x reaches the value of 1.9 min −1 after 4 min of grinding. However, the use of 25.4-mm balls results in a small decrease in the breakage rate during the grinding period, ranging from 1.9 to 1.8 min −1 . In other combinations of the feed size fraction and ball diameter, only deceleration is observed,  which is clearly seen in Fig. 4a-c. It is underlined that the deviations from the linear theory depend on factors such as grinding time (or specific energy), operating conditions, and the heterogeneity of the feed material. In the case of marble, which is a homogeneous material, the breakage rate of each particle size is affected by the accumulation of fine particles as grinding proceeds due to multi-particle interactions. For example, when large balls (40 mm) are used, the increase of fines in the mill results in an increase in the breakage rate of the coarse fraction −3.35 + 2.36 mm, as seen in Fig. 4a. Fine particles cover the surface of the feed material, and the ball impact forces are transmitted through the fines to the coarse particles. On the other hand, the generation of many more fines is not beneficial for the breakage rate, and this is attributed to the cushion effect that takes place; S x remains almost constant as grinding continues. Similar observations can be made for the feed fraction 1.70 + 1.18 mm (Fig. 4a); however, the difference in the magnitude of S x for each specific energy is much smaller. In general, the results shown in Fig. 4a-c indicate that in terms of breakage rate, grinding of coarse and intermediate fractions, namely, −3.35 + 2.36 mm and −1.70 + 1.18 mm, can be carried out efficiently with the use of 40-mm or 25.4-mm balls, while for any other different combination of feed size fraction and ball diameter, the grinding efficiency decreases. The evolution of the breakage rate with a specific energy input of both homogeneous and heterogeneous feeds using the proposed non-linear model was also investigated in a previous recent study [33]. The main advantages of the above kinetic analysis are that (i) it takes into account the specific energy input instead of the grinding time and (ii) it determines with high accuracy the evolution of breakage rate S x as a function of specific energy. The data obtained from the laboratory experiments indicate that the optimum conditions which enable scaling up to larger mills can be reliably determined.

Use of the Non-linear Model
The evolution of the particle size distribution with specific energy (or grinding time) for marble can be determined with the use of the non-linear model (Eq. (11)) proposed by Petrakis and Komnitsas [33]. Parameter n incorporated in this equation is the distribution modulus of the RR distribution, and b (equals 1−m´) is the exponent of the Charles relationship (Eq. (5)), while the n/b ratio defines the degree of acceleration or deceleration of the breakage. As mentioned earlier in this paper, the n values can be determined using either linear or non-linear regression analysis, and the method used defines the accuracy of the RR distribution model. The accuracy of this model was determined using the correlation coefficient (R 2 ) [35,46] which showed that general higher R 2 values are obtained when non-linear regression analysis was used. In addition, it was revealed that the estimated n values were significantly different under the same operating conditions after the two methods were applied. Thus, the investigation of the degree in which these methods affect the adaptability and reliability of the prediction of non-linear models in determining the PSD of the grinding products is of great interest. Another issue to consider is the replacement of n values with the breakage parameter α (with a constant value of 0.90, Table 4) of the first-order kinetic approach, as indicated by the predicted non-linear model. Figure 5 shows the goodness of fit for the predicted nonlinear model which is assessed with the use of the correlation coefficient R 2 at various grinding times, when linear or nonlinear regression analysis is applied to the RR model data.
The results show that in most cases, the R 2 values are higher than 0.95, indicating quite high reliability of the non-linear model prediction in determining the particle size distribution of the grinding products. This is also observed by the high average R 2 values shown in Fig. 6. This figure also shows that the non-linear regression results have higher average R 2 values, except for the grinding time of 2 min. This can be attributed to the relatively low R 2 value (0.90) obtained for this grinding time, as seen in Fig. 5. This latter issue was further explored in order to justify this discrepancy. Other outliers shown in Fig. 5 were also investigated. The non-linear model prediction advantage over the linear one is further supported by the lower RMSE of 5.24 wt.% for the former and 6.54 wt.% for the latter. On the other hand, the IoA′ for both approaches was calculated equal to 0.89. The range of RMSE for the different grinding timesteps using the non-linear model is between 4.63 and 5.47 wt.% compared to 4.84 to 7.85 wt.% for the linear model. This denotes a more accurate performance for the non-linear model. In addition, this finding is further supported by the IoA′ values for the same timesteps, which show a range of 0.89-0.92 for the non-linear model and 0.86-0.92 for the linear one. An interesting finding is that, for the grinding time of 2 min, the RMSE of the non-linear model provides accuracy that is somewhat lower than the other steps; an RMSE decrease of 0.3% is observed when the specific grinding time is excluded from the measurements, while the IoA′ is maintained at the same levels (0.89). Figure 7 shows on log-log scale the evolution of S i values, obtained from the first-order plots of the linear theory of PBM, in relation to the upper feed particle size, when three loads of single size balls are used. It is shown that in all cases, S i increases with size up to a certain point (referred to as optimum feed size) and then decreases sharply, indicating that abnormal breakage takes place. The optimum feed size was found to be 6.0 mm, 3.0 mm, and 1.2 mm when The accuracy of the proposed non-linear model increases when normal breakage occurs; thus, feed size fractions in the abnormal breakage region could be excluded from the evaluation. The latter does not reduce the model's applicability since grinding in the abnormal breakage region operates inefficiently and should be avoided. Therefore, Fig. 8a shows the R 2 values at different grinding times for the predicted non-linear model when linear or non-linear regression analysis is applied to the RR equation data, after the exclusion of outliers, i.e., R 2 values corresponding to the abnormal breakage region. The results confirm that higher R 2 values are obtained when non-linear regression analysis was applied to the RR data. The average R 2 values for the different grinding timesteps are in the range of 0.982 to 0.989 and 0.984 to 0.993 for the linear and non-linear regression analysis, respectively. Similarly, the average RMSE and IoA′ of the latter is 4.97 wt.% and 0.90 correspondingly compared to 6.75 wt.% and 0.88 for the first, supporting further the advantage of the non-linear model. Also, this figure indicates that more reliable results are obtained when increased grinding times are used. This is consistent with the results shown in Table 2 and indicates that the adaptability and reliability of the predicted non-linear model are improved. Figure 8b shows the correlation coefficient R 2 values at different grinding times for the predicted non-linear model, when the n values of the RR equation are used or replaced by the breakage rate parameter α (with a constant value of 0.90) of the first-order kinetic approach. The results confirm that the use of parameter α instead of n values results also in high R 2 values; thus, the particle size distribution of the grinding products can be predicted with high accuracy (average R 2 values are in the range 0.983-0.991). The corresponding average RMSE and IoA′ are 5.67 wt.% and 0.89 wt.%, having similar values with those when n values are used with the RR equation.
It is known that during grinding, the mill and the feed material interact in a complex way which is a function of the operating conditions. The quality of the products depends on the type of mill, the mode of operation (dry or wet), and the properties of the material. Therefore, a deep understanding of the whole process requires the identification of each individual component that contributes to the grinding effect. In light of this, Pérez-Alonso and Delgadillo [47] studied the grinding process in ball mills using the discrete element method (DEM). The DEM approach was found to be a reliable tool to predict the PSD of the products using three components, namely, PBM, impact energy distribution of the mill, and breakage characteristics of particles. In another study, batch grinding experiments using various materials were carried out, in order to simulate the grinding process and predict the product PSDs under different operating conditions [48]. In that study, the breakage parameters of the linear theory of PBM were estimated using a MATLAB code, while the PSDs of the grinding products were described by the RR distribution. Hlabangana et al. [49] investigated the effect of operating conditions, namely, media/material filling and mill rotational speed on the milling efficiency of a laboratory-scale ball mill using the attainable region (AR) method. The findings showed that a correlation between ball size and feed size distribution is necessary to achieve optimum milling efficiency. Petrakis and Komnitsas [50] established potential correlations between the properties of various materials and the breakage parameters of the linear theory of PBM. In particular, it was found that the breakage rate parameter α, which is constant for the same material and independent of grinding conditions, is well correlated with P-wave velocity, Schmidt rebound value, and tangent modulus of elasticity when inverse exponential functions are used. From the energy-saving point of view, the dimensional properties of grinding products, namely, mass, surface area, length, and particle distribution in relation to the energy input in a ball mill were investigated [51]. In addition, the effect of material type and grinding conditions on the relationship between the dimensional properties and the specific energy input was also investigated, and valuable results were obtained regarding the energy requirements both during the initial grinding stages and at higher energy levels. Liao et al. [52] investigated the effect of grinding media type on the fine-grinding performance in wet milling. Comparative experiments between cylpebs and ceramic balls were performed under various conditions, and suggestions were made to improve grinding efficiency. From the comprehensive analysis of the results of the present study and the performance metrics, it is deduced that the proposed non-linear model demonstrates very good ability to predict the particle size distribution of products during marble grinding. Overall, the non-linear model results slightly overestimate the predictions, a finding which can be confirmed in similar research works on grinding tests in mineral processing [53,54]. Since its reliability is mainly affected by the adjustment of the RR distribution to the experimental data, it is suggested that for any material to be tested, including minerals or ores, the most suitable regression method should be used for the implementation of this model. A previous recent study showed that the non-linear model can also reliably predict the product size distributions of heterogeneous materials, such as limonitic laterites [33].

Conclusions
In the present study, the adaptability and reliability of the non-linear model when marble was ground in a laboratory ball mill under different operating conditions were investigated. This model, which considers the linear theory of the PBM as a special case, enables the calculation of the breakage rate parameters as a function of specific energy, as well as the determination of the particle size distribution of the grinding products. By taking into account that these parameters of the non-linear model are affected among others by the accuracy of the RR distribution model to simulate the grinding products, one of the objectives of this study was to clarify the most reliable regression method. In addition, further investigation involved the determination of the effect of the regression method on the non-linear model prediction of the grinding products.
Average R 2 values in the range of 0.962-0.998 and 0.968-0.999 were obtained when respectively linear or nonlinear regression analysis was applied to the RR distribution data. The RMSE and IoA′ metrics confirmed that the nonlinear regression analysis can more accurately describe the experimental data.
Experimental results obtained after grinding the test material for different loads of single size balls showed that the breakage parameters are not constant, as indicated by the linear theory of PBM, but vary with specific energy (or grinding time). The use of the proposed non-linear model provides meaningful and reliable values of the breakage rates for different grinding times. In this context, the acceleration or deceleration of the breakage rate during the grinding period for any combination of the feed size fraction and ball size can be determined.
The results also showed that the reliability of the nonlinear model to predict the particle size distribution of the grinding products depends on the regression (linear or non-linear) analysis applied to the RR distribution data. In general, the accuracy of the non-linear model is improved when non-linear regression analysis is applied. In addition, it was found that the presence of some low R 2 values was attributed to particle sizes exhibiting abnormal breakage behavior. After the exclusion of outliers, the average R 2 values for the different grinding timesteps were in the range of 0.982 to 0.989 and 0.984 to 0.993 for the linear and nonlinear regression analysis, respectively. The advantage of the non-linear model prediction was also confirmed by the reduced RMS (5.24 wt.%) when the non-linear regression analysis was used instead of the linear one (6.54 wt.%). Finally, the goodness of fit was found to be quite high when the n values of the RR model were replaced by the breakage rate α of the linear theory in the proposed equation of the non-linear model. This is considered significant since α is a characteristic parameter of the feed material and remains constant during grinding; thus, its use results in a reduction of the number of experiments that need to be performed for the application of the non-linear model as well as for the prediction of its accuracy.
Funding Open access funding provided by HEAL-Link Greece.

Conflict of Interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.