Identification of the breakage rate and distribution parameters in a non-linear population balance model for batch milling

Powder Technology 208 (2011) 195–204

Contents lists available at ScienceDirect

Powder Technology j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / p ow t e c

Identification of the breakage rate and distribution parameters in a non-linear population balance model for batch milling Maxx Capece, Ecevit Bilgili ⁎, Rajesh Dave Department of Chemical, Biological, and Pharmaceutical Engineering, New Jersey Institute of Technology, 161 Warren St. 150 Tiernan Hall, Newark, NJ 07102, United States

a r t i c l e

i n f o

Article history: Received 9 September 2010 Received in revised form 13 December 2010 Accepted 15 December 2010 Available online 28 December 2010 Keywords: Non-Linear population balance model Inverse problem Non-linear optimization Milling Non-first order effects Back-calculation

a b s t r a c t Non-linear population balance models (PBMs), which have been recently introduced due to the limitations of the classical linear time-invariant (LTINV) model, account for multi-particle interactions and thus are capable of predicting many types of complex non-first order breakage kinetics during size reduction operations. No attempt has been made in the literature to estimate the non-linear model parameters by fitting the model to experimental data and to discriminate various models based on statistical analysis. In this study, a fully numerical back-calculation method was developed in the Matlab environment to determine the model parameters of the non-linear PBM. Not only does the back-calculation method identify the parameters of complicated non-linear PBMs, but also it gives the goodness of fit and certainty of the parameters. The performance of the back-calculation method was first assessed on computer-generated batch milling data with and without random error. The back-calculation method was then applied to experimental batch milling data exhibiting non-first order effects using both the LTINV model and two separate non-linear models. The back-calculation method was able to correctly determine the model parameters of relatively small sets of batch milling data with random errors. Applied to experimental batch milling data, the back-calculation method with a two-parameter non-linear model yielded parameters with reasonable certainty and accurately predicted the slowing-down phenomenon during dry batch milling. This study encourages experimenters to use advanced non-linear population balance models along with the back-calculation method toward estimating the breakage rate and distribution parameters from dense batch milling data sets. © 2010 Elsevier B.V. All rights reserved.

1. Introduction and mathematical formulation Population balance models (PBMs) as mathematical description of size reduction have been used extensively in literature [1–6]. Early PBMs for milling describe size reduction as a first-order rate process as developed by Sedlatschek and Bass [7]. Their first-order breakage hypothesis simply stated that the disappearance of particles of a given size due to breakage is proportional to the weight of particles of that size present. Later, concepts of breakage distribution and probability of selection-for-breakage [8,9] were included by Broadbent and Callcott [10–13], but they treated particle breakage as occurring in stages. In order to apply the PBM to time-continuous milling processes, Gaudin and Meloy [14] derived the time- and size-continuous mass density form of the equation for a well-mixed batch milling process. More commonly, the time-continuous size-discrete form seen in Eq. (1) is preferred as experimental data is inevitably in discrete form [2]. i−1 dMi ðt Þ = −Si Mi ðt Þ + ∑ bij Sj Mj ðt Þ dt j=1

N ≥i≥j ≥1 with Mi ð0Þ = Mini ⁎ Corresponding author. Tel.: + 1 973 596 2998; fax: + 1 973 596 8436. E-mail address: [email protected] (E. Bilgili). 0032-5910/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.powtec.2010.12.019

ð1Þ

In Eq. (1), i and j are the size-class indices and extend from sizeclass 1 containing the coarsest particles to size class N containing the finest particles usually in a geometric progression. Mi represents the mass fraction of particles in size-class i. Si is the specific breakage rate parameter and bij is the breakage distribution parameter, which describes the distribution of particles formed when a particle of size class j is broken. This equation is also known as the linear timeinvariant (LTINV) model because the specific breakage rate does not vary with time and its discretized value is only dependent on particle size of given size class. The following constraints also apply to Eq. (1) due to the conservation of mass: N

SN = 0;

∑

i=j + 1

bij = 1; bii = 0:

ð2Þ

Population balance models have the ability to simulate the evolution of the particle size distribution of a milling process, but also to elucidate the breakage mechanisms (e.g. fracture, cleavage, attrition) [6,15–17]. Numerous methods to determine Si and bij from experimental milling data for the LTINV model have been proposed using both direct measurements and back-calculation [18–24]. The so-called “direct measurement method” demands tedious breakage experiments on numerous mono-sized feeds to determine the parameters without

196

M. Capece et al. / Powder Technology 208 (2011) 195–204

resorting to complex non-linear optimization methods. Austin and Bhatia [18] outlined the experimental procedure for determining the breakage rate parameter and breakage distribution parameter assuming first-order breakage. Austin and Luckie [19] also detailed multiple methods, known as the BI, BII, and BIII methods, to determine the breakage distribution parameters. Aside from obvious complications of these methods such as time-consuming preparation of numerous mono-sized feeds, they also require certain assumptions such as negligible re-breakage of particles which may cause error in the calculation. Back-calculation, a technique which calculates the model parameters that best fit the model to the experimental data, is also widely employed and has significant advantages over the direct measurements. It allows the milling of a natural-sized feed instead of multiple mono-sized feeds and reduces the need for laborious preparation of material. Both linear [22] and non-linear [21,23] optimization techniques have been used with either Reid's analytical solution [25] to the batch milling equation or other approximate solutions [26] to back-calculate the model parameters. Despite the varying success of the above methods and application of the LTINV model, non-random deviation between experimental milling data and modeled predictions of Eq. (1) were noticed and became a source of some criticism [4,27–29]. Such deviations became pronounced at long milling times or when the presence of fines became significant. Austin and Bagga [4] investigated such deviation and observed that the specific rate of breakage decreased as fines accumulated and determined that the source of non-first order effects originated from a cushioning action provided by fine particles. To account for this, Austin and Bagga [4] and Austin et al. [27] introduced a time-dependence to the specific breakage rate parameter. They subsequently assumed that the specific rates of breakage for all particle sizes varied in the same way as the milling environment changed according to an acceleration– deceleration function κ(t). They solved the resulting linear time-variant (LTVAR) model by invoking the concept of a false time or equivalent first-order grind time, θ, as seen in Eq. (3). i−1 dMi ðθÞ = −Si ð0ÞMi ðθÞ + ∑ bij Sj ð0ÞMj ðθÞ dθ j=1

N ≥i≥j ≥1 with Mi ð0Þ = Mini

ð3Þ

dθ = κðt Þdt with θð0Þ = 0 False time may be correlated to true grind time through a series of mono-sized feed milling experiments. The acceleration–deceleration function may be determined in the same manner. While this method can correctly predict a decrease (or increase) in specific breakage rate as milling progresses, it still lacked the ability to account for the source of the non-first order effects (i.e. cushioning of coarse particles by fine particles) in an explicit way. Similarly, other time-variant models to account for non-first order effects based on Kapur's method [22] do not either [28,30]. Bilgili and Scarlett [31] introduced a population balance framework to mathematically explain non-first order effects arising from multi-particle interactions for rate processes. In their model, the specific breakage rate is decomposed into an apparent breakage rate and a population dependent functional where the functional describes different types of non-first order breakage kinetics. Their non-linear population balance model in Eq. (4) is shown in size-discrete form. " # " # N i−1 N dMi ðt Þ = −ki Fi ∑ Piq Mq ðt Þ Mi ðt Þ + ∑ bij kj Fj ∑ Pjq Mq ðt Þ Mj ðt Þ dt q=1 q=1 j=1 N ≥i≥j ≥1 with Mi ð0Þ = Mini ð4Þ In Eq. (4), ki is the apparent specific breakage rate parameter, Fi is a functional of the weighted distribution of the mass fraction where Piq

expresses the contribution of the generic size q to the disappearance rate of particles of size i due to multi-particle interactions, and all other terms are identical to Eq. (1). The following constraints also apply to Eq. (4): ki ≥0; kN = 0;

N

∑

i=j + 1

bij = 1; bii = 0

ð5Þ

Fi ½ ≥0; Fi ½ →1 as ∀ Mq ðq≠i; t Þ→0: The choice of functional, as discussed by Bilgili and Scarlett [31] and Bilgili et al. [32], is partly empirical and depends on mill type, design variables, operation mode, operating variables, and material properties. However, it is well-established that short-time milling of a mono-sized feed is first-order [1]. It has also been shown in many milling studies [4,30,33,34] that non-first order kinetics is significantly contributed by multi-particle interactions. In other words, the breakage of a particle is affected by the surrounding population. Specifically, for dry milling in a ball mill, finer particles exert a “cushioning action” on the coarser particles, thus reducing the specific breakage rate of the coarser particles [4]. The functional chosen (e.g. F ≡ exp[],(1 + []− 1), etc.) and the weighting function P must accurately reflect these considerations and physical interactions. Finally, whatever functional chosen, it should explain Type I and Types II or III deviations from the LTINV model [31], while yielding parameters with higher statistical certainty than the LTINV model, as being demonstrated in this paper. Therefore, if a functional is chosen completely arbitrarily without taking into account any of the above considerations, the non-linear model will probably fail, resulting in inferior predictions and parameters as compared with those of the LTINV model. Numerical simulations using a simplified version of Eq. (4), known as Model B and shown in Eq. (6) in the sequel, were performed by Bilgili and Scarlett [31] and Bilgili et al. [32]. In their work, two specific forms of the functional (Eqs. (7) and (11) in this study) were considered and found to successfully describe different types of nonlinear kinetics. Furthermore, they successfully predicted the three types of milling behavior observed in literature including the decrease in specific breakage rate observed by Austin [27] that could not be explained by the LTINV model, i.e., Eq. (1). Further novelty of their framework extends from the fact that the non-linear model reduces to the linear model in the absence of multi-particle influence or in other words when F[ ] → 1. Due to the ability of this non-linear population balance model to encompass many different types of particle breakage kinetics phenomenologically, we argue that it can lead to better design, control, and optimization of size reduction processes. Because no analytical solution exists for the non-linear model, calculation of the breakage rate parameters, breakage distribution parameters, and the non-linear functional parameters remains a formidable challenge. Determining the model parameters is of upmost importance in order to use the PBM for process simulation, design, scaling, control, and optimization [1–3,15,35]. Having estimated the model parameters for a given set of operating conditions and material in a mill, one can predict the temporal evolution of the size distribution and final product size distribution for any given feed (initial) particle size distribution. This is especially important when non-first-order effects are significant, which result in a pronounced impact of the initial condition on the temporal evolution [31]. In fact, multiple feed size distributions can also be used to generate more dense data sets for determining PBM parameters with more statistical certainty. In addition, one can fit the same model to size distributions obtained from different operating conditions (e.g. volume percent media loading in ball milling) and assess the impact on the model parameters for a given material. This can allow a numerical optimization of the process parameters. Similarly, for the same operating conditions in the mill, one can determine the specific rate of breakage of different materials and assess their relative particle strengths. Finally, fitting the non-linear PBM

M. Capece et al. / Powder Technology 208 (2011) 195–204

to experimental data, also known as the inverse problem, may help researchers elucidate breakage mechanisms (e.g. fracture, cleavage, attrition) based on values of the breakage distribution parameters [6,15–17], especially when multi-particle interactions or non-linear effects are significant. While the non-linear model has exceptional modeling and predictive benefits, it appears that no attempt has been made to estimate the non-linear model parameters by fitting the model to experimental data. This study addresses the inverse problem of identifying the parameters of non-linear PBMs using batch milling data while elucidating non-first-order effects due to multi-particle interactions and discriminates various models based on statistical analysis. A fully numerical back-calculation method was developed that uses a nonlinear optimization technique in conjunction with a fast and accurate ordinary differential equation (ODE) solver for the solution of the non-linear PBM. Although simpler back-calculation methods have been applied to the LTINV model [21–23], this is the first time an advanced back-calculation method was applied to non-linear PBMs. This back-calculation method was performed on computer-generated batch milling data without random error to show the accuracy of the methods and techniques used. Similarly, the back-calculation method was applied to computer-generated data containing random errors to assess the method's general robustness. In addition, the method was applied to published experimental batch milling data [36] exhibiting slowing-down phenomenon during dry milling. We established that the fully-numerical back-calculation method is accurate and yields sufficient information to discriminate different models. The so-called two-parameter non-linear model was able to explain the slowingdown phenomenon with reasonable accuracy. Having many parameters in a more complex non-linear model, e.g. the three-parameter model, did not improve the goodness of the fit and the standard error of the fitted parameters were high. Hence, the two-parameter nonlinear model was discriminated as the best model to explain the slowing-down phenomenon during the batch dry milling. 2. Methods 2.1. Computer-generated batch milling data In order to test the accuracy and robustness of the back-calculation method developed in this study, the technique was first applied to computer-generated (artificial) batch milling data. Batch milling data exhibiting non-first-order effects can be generated artificially by solving the non-linear PBM with pre-chosen model parameter values. In this manner, the model parameters determined by the back-calculation method may be compared to the model parameters chosen to produce the computer-generated data. To further test the back-calculation method, random errors were computationally introduced into the computer-generated data to simulate errors that may occur during an experimental procedure. This will effectively allow a proper assessment of the numerical methods used to solve the PBM and the optimization techniques used to determine the model parameters. Finally, the backcalculation method will be performed on genuine experimental batch milling data which will be discussed in detail in the following section. Computer-generated data was produced by numerical solution of Eq. (6) presented below (known as Model B [32]) along with the following size-dependent functions for ki, bij, Fi and Piq in Eq. (7). A power-law function for ki and a normalized function for Bij were chosen due to their common use in literature [19,37] while Fi and Piq have previously been used in batch milling simulations [31,32,38]. "

#

"

#

N i−1 N dMi ðt Þ = −ki Fi ∑ Piq Mq ðt Þ Mi ðt Þ + ∑ bij kj Fj ∑ Pjq Mq ðt Þ Mj ðt Þ dt q=i j=1 q=j

N ≥i≥j ≥1 with Mi ð0Þ = Mini ð6Þ

 m x ki = A i x0 !μ !v x x Bij = φ i−1 + ð1−φÞ i−1 ; where bij = Bij −Bi + xj xj " #−1   N xq α Fi = 1 + λ ∑ 1− Mq xi q=i

197

1j

and bN j = BN j

ð7Þ

Eq. (6) is similar to the non-linear PBM of Eq. (4) except that in the latter case, it is assumed that the specific breakage rate of particles of size i is affected only by the mass fraction of particles of size classes j≥i. Additional forms of the non-linear PBM are discussed in detail by Bilgili et al. [32]. The constraints of Eq. (5) apply likewise to Eq. (6). The solution of Eq. (6) yields the mass fraction distribution from which the cumulative mass fraction undersize distribution was calculated so as to compare the model predictions and the experimental data. The cumulative mass fraction undersize distribution can easily be calculated by use of Eq. (8) where Qi represents the mass fraction of material less than a given size class. N

Qi = ∑ Mi i

ð8Þ

First, batch milling data without random error was generated. Each size class was represented by its upper edge size, xi, that progresses from x1 = 1000 μm downward with a geometric progression ratio s of 21/16. The number of size classes N was selected to be 160. A Gaussian feed (initial particle size distribution) with a mean size of 600 μm and a standard deviation of 100 μm was generated in Matlab using the function “normpdf”. The resultant mass density distribution was converted to a mass frequency distribution. Cumulative mass fraction undersize data were generated for milling times of 10, 20, 30, 40, 60, and 80 min. The parameters of the size-dependant functions used in the solution of the model (A = 0.5 min− 1, m = 1.0, x0 = 1000 μm, φ = 0.3, μ = 1.3, ν = 5.8, α = 3.0, λ = 16) were chosen to display a deceleration in breakage rate (slowing-down) as milling progressed. The complete specified set of coupled, non-linear ODEs that consists of Eqs. (6) and (7) with the chosen parameter set was integrated using a stiff ODE solver “ode15s” [39] in Matlab v7.9. The relative error tolerance and absolute error tolerance for the integrator were set to 10− 4 and 10− 6 respectively. Second, batch milling data with random error was generated. The error-free mass fraction distribution data generated by the procedure described above was used as the data set to which errors were applied. A Gaussian distribution with a mean of 0.05 and a standard deviation of 0.01 was generated in Microsoft Excel. A value from the distribution was randomly chosen using Excel's random selection function, “rand”, and was then randomly assigned as a positive or negative number using Excel's random function, “randbetween”. The value chosen from the Gaussian distribution, now being either positive or negative, was multiplied by the mass fraction of a single size class and added to the original mass fraction of that same size class as shown in Eq. (9). 

Mi = Mi + Mi × Errori

ð9Þ

In this equation, Mi* is the mass fraction now with error, Mi represents the mass fraction of the error-free computer-generated data, and Errori is the value randomly selected from the Gaussian distribution. This procedure was repeated individually for all 160 size classes in which a value from the Gaussian distribution was randomly chosen and randomly assigned as a positive or negative number and incorporated in the error-free mass fraction distribution according to Eq. (9). Given that

198

M. Capece et al. / Powder Technology 208 (2011) 195–204

some error values of the Gaussian distribution may increase or decrease a given error-free mass fraction value more or less than 5%, it was assured that the average relative error applied for all 160 size classes was 5%. Thus, this set of data is referred to as computer-generated data with 5% error. After introduction of the random error, the mass fraction distributions were normalized to 1 and the cumulative mass fraction undersize distributions were generated. This procedure was repeated identically for 10% relative error except that a Gaussian distribution was generated with a mean of 0.10 and a standard deviation of 0.01. This data set is referred to as computer generated data with 10% error. It is important to reiterate that all the data were generated with the same model parameter set specified above; however, the actual mass fraction distribution was different due to the level of random error introduced (0%, 5%, and 10%). 2.2. Experimental batch milling data Experimental batch milling data were taken from a study on the kinetics of fine dry milling of quartz performed by Austin et al. [36]. In the referenced study, quartz, 600 μm by 425 μm, was milled to a product containing 80% by weight less than 10 μm in a tumbling ball mill for a maximum of 512 min. Particle size distributions were determined by dry sieving and laser diffraction which were combined to give a size distribution including both the coarsest and finest particles down to 1 μm. In their study, milling was not restricted to short milling times which would otherwise minimize the effects of multi-particle interactions. As a result, the specific breakage rate of quartz decreased significantly as milling progressed due to the generation of finer particles in the mill. Due to the observed non-first order kinetics, a LTVAR model introduced by Austin et al. [4,27], which uses the concept of false time as was previously discussed (refer to Eq. (3)), was fit to the quartz data [36]. Although this data set is not very dense, meaning 9 particle size distributions and about 15 data points within each size distribution, it still offers an opportunity to apply a non-linear PBM and demonstrate the back-calculation method developed in this study. The cumulative mass fraction undersize distribution at 1 min milling time was used as the initial feed distribution. Accordingly, all subsequent milling times were reduced by 1 min. 2.3. Back-calculation method While direct measurement method may be used to identify the model parameters as was previously discussed, it requires numerous breakage tests performed on mono-, binary-, and natural-sized feeds. This approach entails tedious experimentation along with timeconsuming material preparation. The alternative method to identify the model parameters is through back-calculation method utilizing non-linear optimization. The sum of squared residuals between the model predictions and the experimental data is minimized by a nonlinear optimizer and the best set of parameters is chosen. Various models can be fit in this way and can be discriminated based on the goodness of fit to the experimental data and the degree of certainty of the fitted parameters. In this study, the non-linear optimizer “fmincon”, part of the Matlab v7.9 optimization toolbox, was used to minimize the calculated error between the model and the experimental data. “fmincon” uses the algorithm, “interior-point”, to solve the non-linear optimization problem at hand and is suitable for general non-linear optimization. Readers are referred to Byrd et al. [40] for details regarding this particular algorithm. “fmincon” requires the minimum of an initial parameter set and a stopping criterion (e.g. maximum number of iterations). A simple sum of squared residuals SSR, seen in Eq. (10), was used as the optimizer's objective function. Z

N

SSR = ∑ ∑

z=1 i=1

  Mod Exp 2 Qzi −Qzi

ð10Þ

Q refers to the cumulative mass fraction undersize distribution, the superscripts Exp and Mod refer to the experimental data and model predicted data respectively, and the subscripts z and i refer to the product size distribution for a particular milling time and the size class respectively. This optimization technique also requires, during each iteration, the discretized PBM of Eq. (6) to be solved numerically using the ODE integrator, “ode15s” with its relative and absolute tolerances set to 10− 4 and 10− 6 respectively. The optimizer stopping criteria include the termination tolerance on the constraint violation, termination tolerance on the function value, and termination tolerance on the parameters, which were all set to 10− 9. Maximum iterations and function evaluations were both set to 106. “fmincon” also allows constraints (boundaries) to be placed on the parameters. All parameters except for two, φ and λ, were constrained between 0.01 and 15.0. The parameter, φ, was constrained between 0.01 and 1.00 due to the form of the breakage distribution function chosen. The parameter λ, which is a scaling factor for the non-first order effects, was constrained between − 25.0 and 25.0 to allow for the prediction of either deceleration or acceleration of the breakage rate. Due to non-linear structure of the objective function, comparatively large number of parameters to fit with the non-linear model, and possibly poor initial parameter guess, convergence may be excessively slow or may not be achieved. In order to circumvent this problem, the LTINV model, which has 5 parameters due to the forms of ki and bij chosen, was fit to the experimental data first. The parameter set obtained from back-calculation using the LTINV model was then used as the initial guess in the subsequent fitting of the non-linear model. The LTINV model of Eq. (1) and the non-linear model of Eq. (6) with the size-dependent functions for ki, bij, Fi and Piq in Eq. (7) were used to fit the computer-generated data as well as the experimental milling data. The above back-calculation procedure was performed on both the computer-generated data and experimental milling data similarly except for one difference in discretizing the size domain. For backcalculation of the experimental milling data, the top most size class, x1, was specified as 850 μm and the geometric progression ratio s was specified as 21/12.2. The number of size classes N remained 160. In addition, a second form of Fi and Piq seen in Eq. (11), which had been used previously in milling simulations [31], was only used in the fitting of the experimental data along with Eq. (6). 2

3 N

Fi = exp4λ ∑

q=i

1+

1 

 Mq 5

xq = xi α β

ð11Þ

This non-linear model is shortly referred to as the three-parameter non-linear model as the functional Fi contains three parameters while the remaining five parameters of this model are common to the LTINV model. Likewise, the non-linear model using the form of Fi and Piq of Eq. (7) is referred to as the two-parameter non-linear model. In order to compare models based on the goodness of fit, the sum of squared residuals SSR was calculated using Eq. (10). To account for the degrees of freedom for the differing models, the standard error of the residuals SER was also calculated as follows: rffiffiffiffiffiffiffiffiffiffiffiffi SSR SER = D−K

ð12Þ

where D is the number of comparisons made between the model predictions and the experimental data (i.e. number of experimental or computer-generated data points) and K is the total number of parameters of the model. SSR and SER were used to quantify the accuracy or goodness of the model fit for the genuine experimental data, while %error of back-calculated parameter values, in addition to SSR and SER, were used to quantify the same for the computergenerated data.

M. Capece et al. / Powder Technology 208 (2011) 195–204

199

Table 1 Back-calculation results for computer-generated batch milling data without error. Parameter

A (min− 1) m φ μ ν α λ SSR SER a b

Parameter valuea

Initial parameter setb

LTINV model Back-calculated value

% |Error| of back-calculated value

Back-calculated value

% |Error| of back-calculated value

0.500 1.000 0.300 1.300 5.800 3.000 16.000

5.000 0.100 3.000 13.000 0.580 0.300 0.000

0.245 2.176 0.292 5.079 1.063 – – 2.335×10− 1 1.564×10− 2

51.049 117.649 2.710 290.682 81.675 – –

0.500 1.000 0.300 1.300 5.800 3.000 16.000 8.804×10− 14 9.611×10− 9

2.854×10− 4 1.429×10− 5 2.409×10− 5 5.036×10− 5 1.735×10− 4 4.425×10− 4 1.027×10− 4

Non-linear model

Parameter values used to generate the computer simulated batch milling data. Initial parameter guess for the optimizer. Initial guess for α and λ were only used for back-calculation of the non-linear model.

To be able to discriminate different models, the standard error of the parameter SEP was calculated and evaluated in addition to SSR and SER. While SSR and SER give an indication of how well the model fits the experimental data, SEP is an important statistic in determining the degree of certainty in the estimated parameters. A high SEP in comparison with the estimated parameter value indicates low certainty of the parameter estimate (or a wide confidence interval), which may be due to a lack of sufficiently accurate dense data set and/ or the use of a poor model. SEP was calculated by taking the square root of the diagonals of the parameter covariance matrix C determined at the back-calculated parameter set m′ as detailed by Aster et al. [41] and shown in Eq. (13). h    i−1  0 2 0 T 0 C m = SER J m J m ¼

¼

ð13Þ

where T and − 1 stand for matrix transpose and inverse operators respectively. The Jacobian matrix J was also computed at the parameter set determined by the back-calculation method. Each row of the Jacobian matrix is comprised of the partial derivatives of the residual for a single experimental observation and model prediction with respect to each model parameter. The partial derivatives were estimated by using a first-order finite difference method also detailed by Aster et al. [41]. A parameter step size of 10− 4 was used in the finite difference method for the Jacobian matrix calculation. The SEP values obtained from reducing the step size from 10− 3 to 10− 4 differed by less than 0.1%; hence, no further reduction in step size was needed. 3. Results and discussion 3.1. Back-calculation using computer-generated batch milling data The non-linear PBM parameters were back-calculated first from error-free computer-generated data in order to determine the accuracy of the numerical methods used. This is by no means a trivial exercise as Berthiaux and Dodds [24] showed that some non-linear optimization techniques considered in their study may give incorrect results. A successful result is only attained if the back-calculation method returns the same parameters used to generate the data. As can be seen from Table 1, the optimizer converged to the correct set of parameters for the two-parameter non-linear model, which had been used to generate the data, with excellent accuracy. The small differences are likely due to rounding errors and other numerical errors associated with the optimization and ODE integration. On the other hand, the fitted parameters of the LTINV model had significant errors. Comparison of the SSR and SER values for the LTINV model and the two-parameter non-linear model shows a superior fit for the latter. These findings imply that the LTINV model may not accurately

fit experimental data that exhibit significant multi-particle interactions such as cushioning action of the finer particles on the coarser particles [4]. This conclusion is further supported by examining the SEP and relative SEP values listed in Table 2. The high SEP and relative SEP values for the LTINV model as compared to the two-parameter non-linear model show more uncertainty in the estimated parameters of the LTINV model. A more rigorous assessment of the back-calculation method is achieved by performing the procedure on computer-generated data with error. Random errors are inevitably incurred in any experiment. To simulate this and to test the robustness of the back-calculation method, 5% and 10% relative error was introduced into the error-free milling data. Table 3 shows results for the back-calculation method performed on computer-generated data with 5% error for the LTINV model and the two-parameter non-linear model. It can be deduced from Table 3 that the back-calculation method with the twoparameter non-linear model is robust even with data containing random errors. In addition, the SSR and SER values indicate that the non-linear model predicts the computer-generated batch milling data with much higher accuracy than the LTINV model. The SEP and relative SEP values in Table 4 suggest an increase in uncertainty for the parameters of the two-parameter non-linear model when random errors were introduced. The SEP and relative SEP values for the LTINV model still remained large. Similar conclusions can be drawn for the back-calculation results for the computer-generated data with 10% error, as shown in Tables 5 and 6. The parameters of the twoparameter non-linear model were calculated with high accuracy even when 10% relative error was introduced. A comparison of the relative SEP values in Tables 2, 4, and 6 leads to the conclusion that the certainty of the parameters decreases when random errors increase, as intuitively expected. The predictions of the LTINV model and the two-parameter nonlinear model were further investigated with their respective graphical fittings. Fig. 1 shows the LTINV model fit to the computer-generated data with 10% error. Data presentation in cumulative mass fraction undersize, as opposed to mass fraction distribution, tends to attenuate Table 2 SEP for computer-generated batch milling data without error. Parameter

A (min− 1) m φ μ ν α λ a

LTINV model

Non-linear model

SEP

Relative SEPa (%)

SEP

Relative SEPa (%)

0.038 0.016 0.078 0.193 0.015 – –

15.681 0.726 26.646 3.798 1.366 – –

1.124×10− 7 1.867×10− 8 2.243×10− 8 6.411×10− 8 2.351×10− 8 1.710×10− 7 5.126×10− 8

2.248×10− 5 1.867×10− 6 7.478×10− 6 4.932×10− 6 4.054×10− 7 5.699×10− 6 3.204×10− 7

Relative to the parameters determined by back-calculation.

200

M. Capece et al. / Powder Technology 208 (2011) 195–204

Table 3 Back-calculation results for computer-generated batch milling data with 5% error. Parameter

A (min− 1) m φ μ ν α λ SSR SER a b

Parameter valuea

Initial parameter setb

LTINV model Back-calculated value

% |Error| of back-calculated value

Back-calculated value

% |Error| of back-calculated value

0.500 1.000 0.300 1.300 5.800 3.000 16.000

5.000 0.100 3.000 13.000 0.580 0.300 0.000

0.244 2.164 0.294 5.102 1.062 – – 2.396×10− 1 1.584×10− 2

51.116 116.426 2.131 292.441 81.688 – –

0.518 0.994 0.298 1.311 5.904 2.869 15.975 1.190×10− 3 1.117×10− 3

3.644 0.617 0.683 0.881 1.801 4.350 0.157

Parameter values used to generate the computer simulated batch milling data. Initial parameter guess for the optimizer. Initial guess for α and λ were only used for back-calculation of the non-linear model.

any fluctuations in particle size distribution. As a consequence, visually small deviations between the model predicted and experimentally determined cumulative mass fraction distributions can actually be significant. In light of this, Fig. 1 shows that the LTINV model prediction deviates from the computer-generated data exhibiting non-first order effects. The LTINV model does however predict the general trend of temporal variation still offering some predictive capability. Fig. 2, in contrast, shows a perfect fit of the two-parameter non-linear model to the computer-generated data with 10% error for all milling times. The findings in this section are significant for several reasons. The numerical methods and non-linear optimization technique used in this study are both highly accurate and robust. This back-calculation method can be used as a reliable technique to determine the nonlinear model parameters even for a relatively small set of experimental data containing significant random error. Considering the advancements in particle size measurement techniques including in-line techniques [42], the data fit in this study, six time sets, represent a relatively small data set. Furthermore, the SEP values were shown to be small even for data with 10% random error instilling high confidence in results of the back-calculation method. Additional findings in this section showed that the LTINV model still retains some predictive capability in fitting batch milling data exhibiting non-first order effects. However, the high relative SEP values also show that the parameters of the LTINV model were not calculated with great certainty possibly due to the inability of this model to account for nonfirst order effects. 3.2. Back-calculation using experimental batch milling data The back-calculation method outlined in Section 2.3 was applied to the experimental dry milling data on quartz, which was previously investigated by Austin et al. [36]. The LTINV model of Eq. (1), the nonlinear model of Eq. (6) with the size dependent functions of Eq. (7) (two-parameter non-linear model), and a second non-linear model

Table 4 SEP for computer-generated batch milling data with 5% error. Parameter

A (min− 1) m φ μ ν α λ a

Non-linear model

LTINV model

Non-linear model

SEP

Relative SEPa (%)

SEP

Relative SEPa (%)

0.039 0.016 0.078 0.195 0.015 – –

15.953 0.737 26.631 3.816 1.378 – –

1.346×10− 2 2.163×10− 3 6.165×10− 3 2.771×10− 3 7.543×10− 3 1.910×10− 2 5.445×10− 3

2.598 0.218 2.069 0.211 0.128 0.666 0.034

Relative to the parameters determined by back-calculation.

using the form of Fi of Eq. (11) (three-parameter non-linear model) were used in the back-calculation method. The two non-linear models used in this section are very similar in several respects. They both assume that the specific breakage rate of particles of size i is affected only by the mass fraction of particles of size classes j≥i. Both non-linear models also use the same size dependant functions for ki and bij. The distinguishing factors between the two non-linear models are the form of Fi used and the number of parameters associated with Fi. In this study, a two-parameter form of Fi seen in Eq. (7) and a three-parameter form seen in Eq. (11) were chosen. To study the predictive capability of the three above mentioned models, they were fit to the data of Austin et al. [36] by the backcalculation method developed in this study. Based on the SSR and SER values presented in Table 7, both non-linear models fit the experimental data better than the LTINV model. The SER value for the LTINV model was more than 1.5 times those for both non-linear models. The two-parameter non-linear model using the size dependent functional of Eq. (7) was chosen to be the best fit due to the lowest SSR and SER, although the 3-parameter non-linear model fit the data comparatively well. Furthermore, both versions of non-linear model predicted the slowing-down phenomenon (deceleration in breakage rate) as shown by the acceleration-deceleration scaling factor λ. Being a model dependent parameter, λ reveals the severity of the multi-particle interactions and deviations from the LTINV model based on its magnitude and the type of deviation based on its sign. The functional of Eq. (7), due to its form, predicts slowing breakage (deceleration) for λ N 0 while the functional of Eq. (11) predicts slowing breakage for λ b 0. Both models will predict the absence of non-first order effects when λ → 0 (i.e. F[ ] → 1), but both showed a significant effect with values of 3.382 and −1.870 respectively. These results confirm that the non-linear models have better fitting capability over the LTINV model and that they are capable of predicting the correct non-linear breakage behavior (slowing-down). The relative SEP values for each model in Table 7 are also important to discriminate the different models. Overall, the two-parameter nonlinear model has greater certainty of the estimated parameters as compared with the LTINV model and the three-parameter non-linear model. However, unlike the high SEP values associated with the LTINV model, which is partly due to the inability of the model to account for the non-first order effects, the high SEP values of the two-parameter non-linear model may be due to the lack of sufficient data. SEP can be reduced by using more dense data sets in order to gain more certainty in the parameters. The three-parameter non-linear model shows a significantly higher SEP than the two-parameter non-linear model, which may be due to the use of an overly complex model. The parameters of the three-parameter model cannot be determined with certainty given the limited amount of experimental data. The parameter β cannot be predicted with any certainty. Although both

M. Capece et al. / Powder Technology 208 (2011) 195–204

201

Table 5 Back-calculation results for computer-generated batch milling data with 10% error. Parameter

A (min− 1) m φ μ ν α λ SSR SER a b

Parameter valuea

Initial parameter setb

LTINV model Back-calculated value

% |Error| of back-calculated value

Back-calculated value

% |Error| of back-calculated value

0.500 1.000 0.300 1.300 5.800 3.000 16.000

5.000 0.100 3.000 13.000 0.580 0.300 0.000

0.244 2.184 0.280 5.214 1.062 – – 2.387×10− 1 1.581×10− 2

51.118 118.424 6.620 301.111 81.696 – –

0.539 1.006 0.303 1.315 5.908 2.726 16.049 3.400×10− 3 1.889×10− 3

7.841 0.608 0.970 1.155 1.866 9.130 0.309

Non-linear model

Parameter values used to generate the computer simulated batch milling data. Initial parameter guess for the optimizer. Initial guess for α and λ were only used for back-calculation of the non-linear model.

non-linear models fit the experimental data almost equally well as seen by the SER, the two-parameter non-linear model proves to be a superior model due to the higher certainty in the estimated parameters. Further comparisons can be made between the LTINV model and the two-parameter non-linear model in graphical form. Fig. 3 illustrates that the LTINV model fits experimental data well at short milling times, however, is unable to do so at times greater than 15 min. This is not surprising since it is well known that a majority of short-time milling experiments validate first-order breakage kinetics [1,2,31,32]. At longer milling times, the LTINV model deviates significantly from what was experimentally determined. This finding is in agreement with the study of Austin et al. [36] in which a slowingdown factor was required for times greater than 15 min in order to accurately fit the experimental data. In contrast to the LTINV model, the two-parameter non-linear model fits data reasonably well for most milling times without requiring a slowing-down factor or substituting true milling time with a false milling time, as used in [36] (see Fig. 4). The back-calculation of F[ ] within the context of the non-linear PBM allows a quantitative and explicit assessment of the multi-particle interactions during a milling process, which is another novel aspect of this paper. A more detailed examination of the non-first order effects originating from particle cushioning as milling progresses for the twoparameter non-linear model is presented in Fig. 5. The value of the functional F[ ], which expresses the effect of multi-particle interactions, is plotted versus particle size for various milling times. As milling progresses and fine particles accumulate, the functional value decreases translating to a decreased specific breakage rate and more multi-particle interaction. This effect is most pronounced at the longest milling times in accordance with the lowest values for F[ ] and the most fines present in the mill. Furthermore, coarse particles are more affected, whereas the finest particles experience the least nonfirst order effects as F[ ] → 1. As shown, multi-particle interactions play an important role and its effects can be quantified through use of the non-linear PBM. The methodology adopted here did not require

Fig. 1. Comparison of computer-generated batch milling data with 10% error and the LTINV model prediction. Back-calculated model parameters: A = 0.244 min− 1, m = 2.184, φ = 0.280, μ = 5.214, ν = 1.062.

Table 6 SEP for computer-generated batch milling data with 10% error. Parameter

A (min− 1) m φ μ ν α λ a

LTINV model

Non-linear model

SEP

Relative SEPa (%)

SEP

Relative SEPa (%)

0.040 0.016 0.085 0.207 0.014 – –

16.336 0.732 30.183 3.962 1.334 – –

2.381×10− 2 3.643×10− 3 1.057×10− 2 4.732×10− 3 1.297×10− 2 3.130×10− 2 8.764×10− 3

4.416 0.362 3.488 0.360 0.219 1.148 0.055

Relative to the parameters determined by back-calculation.

Fig. 2. Comparison of computer-generated batch milling data with 10% error and the two-parameter non-linear model prediction. Back-calculated model parameters: A = 0.539 min− 1, m = 1.006, φ = 0.303, μ = 1.315, ν = 5.908, α = 2.726, λ = 16.049.

202

M. Capece et al. / Powder Technology 208 (2011) 195–204

Table 7 Back-calculation results and SEP for experimental batch milling data. Parameter

A (min m φ μ ν α β λ SSR SER a

−1

)

LTINV model

Two-parameter non-linear model

Three-parameter non-linear model

Back-calculated value

SEP

Relative SEPa(%)

Back-calculated value

SEP

Relative SEPa(%)

Back-calculated value

SEP

Relative SEPa (%)

1.082 1.136 0.737 1.053 5.122 – – – 1.405×10− 1 3.180×10− 2

0.342 0.012 0.267 0.094 2.348 – – –

31.582 1.070 36.187 8.916 45.839 – – –

1.369 1.091 0.522 0.816 3.195 7.327 – 3.382 5.768×10− 2 2.052×10− 2

0.277 0.015 0.187 0.085 0.500 2.763 – 1.307

20.239 1.351 35.872 10.362 15.663 37.704 – 38.634

1.660 1.075 0.468 0.832 3.443 1.525 0.181 − 1.870 5.831×10− 2 2.071×10− 2

0.511 0.011 0.247 0.079 0.460 1.275 3.387 1.994

30.776 0.984 52.658 9.541 13.358 83.607 1870.427 − 106.650

Relative to the parameters determined by back-calculation.

In this study, a fully numerical back-calculation method was developed and applied to batch milling data to determine the model parameters using a LTINV model and two separate non-linear models. Unlike other studies of similar scope in the literature, a more detailed statistical analysis of the parameters has been made to discriminate the models. The back-calculation method was first applied to errorfree computer-generated batch milling data in order to determine the method's numerical accuracy. The back-calculation method was then applied to computer-generated batch milling data with 5% and 10% random error to test the robustness of the method. Benchmarking the back-calculation method in this manner has revealed it to be highly accurate, robust, and capable of determining the model parameters with reasonable certainty for the relatively small data sets under consideration. The use of the back-calculation method was then extended to genuine experimental batch milling data from literature. The data taken from a study by Austin et al. [36] exhibited non-first order

breakage kinetics resulting from cushioning of coarse particles by fine particles. The earlier study fit the experimental data by incorporating a slowing-down factor to account for a deceleration in breakage rate. Besides being laborious, their methodology is of limited predictive capability because the specific breakage rate is assumed independent of the particle size distribution contrary to other published experimental data [4,28–30,33,34]. The back-calculation method developed in this study is the first attempt to apply a non-linear PBM [31,32] to some experimental data. While allowing a detailed and explicit analysis of multi-particle interactions and having a predictive capability for various initial size distributions [31,32,43], the twoparameter non-linear PBM appears to be better than both the traditional LTINV model and the LTVAR model. It led to a more accurate fit with higher certainty of the estimated parameters, while accounting for multi-particle interactions explicitly and without requiring time-consuming and laborious mono-sized breakage tests. Despite the promising results obtained with the application of the fully-numerical back-calculation method to the two-parameter nonlinear model, it is still left to future research to fully investigate and expand upon the methods used in this study. Other non-linear models with different forms of the functional can be considered for different milling systems. Potential applications include dry and wet batch milling performed in different equipment as well as continuous operation mode (see e.g. [44]) and milling operations such as roll

Fig. 3. Comparison of the experimental batch milling data (from [36]) and the LTINV model prediction. Back-calculated model parameters: A = 1.082 min− 1, m = 1.136, φ = 0.737, μ = 1.053, ν = 5.122.

Fig. 4. Comparison of the experimental batch milling data (from [36]) and the twoparameter non-linear model prediction. Back-calculated model parameters: A = 1.369 min− 1, m = 1.091, φ = 0.522, μ = 0.816, ν = 3.195, α = 7.327, λ = 3.382

milling experiments on multiple mono-dispersed samples, unlike the LTVAR model employed by Austin et al. [36]. Readers are referred to [31,32] for other advantages of the non-linear PBM over the LTVAR model.

4. Summary, conclusions, and outlook

M. Capece et al. / Powder Technology 208 (2011) 195–204

t x x0 Z

203

time, min particle size, μm normalizing reference particle size, μm total number of milling time samples, #

Greek letters α exponent of the weighting function, dimensionless β inflection parameter of the weighting function, dimensionless θ false time, min κ acceleration-deceleration function, dimensionless λ acceleration-deceleration parameter, dimensionless μ breakage distribution exponent, dimensionless ν breakage distribution exponent, dimensionless φ breakage distribution constant, dimensionless

Fig. 5. Population dependent functional value for various milling times calculated using the two-parameter non-linear model, which was fit to the experimental batch milling data (from [36]). Dotted line represents particles which are no longer present in the mill due to breakage.

milling and particle bed compression tests, which can be mathematically treated as a series of elementary breakage events (see e.g. [45]). In addition, significant improvement can be made to reduce the uncertainty of the parameters determined by the back-calculation method. Denser particle size distribution data sets and more milling experiments with different initial size distributions (excluding breakage tests on mono-dispersed samples) can help decrease the uncertainty of the non-linear model parameters. It is expected that the current study will stimulate extensive use of the non-linear PBMs in fitting batch milling systems that exhibit significant non-first-order effects. Notation A bij Bij C D Errori F[ ] J K ki LTINV LTVAR m m′ Mi N ODE PBM Piq Qi s SEP SER Si SSR

breakage rate constant, min− 1 breakage distribution parameter, dimensionless cumulative breakage distribution parameter, dimensionless parameter covariance matrix number of experimental or computer-generated data points, # random error value chosen from Gaussian distribution, dimensionless a non-linear function, dimensionless Jacobian matrix number of parameters, # apparent specific breakage rate parameter, min− 1 linear time-invariant linear time-variant breakage rate exponent, dimensionless vector of the parameter set determined by back-calculation mass fraction in size class i, dimensionless total number of size classes, # ordinary differential equation population balance model weighting parameter, dimensionless cumulative mass fraction undersize, dimensionless ratio of the upper edge to the lower edge of a size class, dimensionless standard error of the parameter standard error of residuals, dimensionless specific breakage rate parameter, min− 1 sum of squared residuals, dimensionless

Subscripts i, j, q size class indices, # ini initial N index for the finest size class z index for milling time sample

Superscripts Exp experimental data Mod model predicted data T matrix transpose operator * indicates data with introduced error −1 matrix inverse operator

Acknowledgments The authors gratefully acknowledge financial support from the National Science Foundation Engineering Research Center for Structured Organic Particulate Systems (NSF ERC for SOPS) through the Grant EEC0540855. References [1] C.L. Prasher, Crushing and Grinding Process Handbook, Wiley, Chichester, 1987. [2] L.G. Austin, A review: introduction to the mathematical description of grinding as a rate process, Powder Technol. 5 (1971) 1–17. [3] L.G. Austin, J. Shah, J. Wang, E. Gallagher, P.T. Luckie, An analysis of ball-and-race milling. Part I. The Hardgrove mill, Powder Technol. 29 (1981) 263–275. [4] L.G. Austin, P. Bagga, An analysis of fine dry grinding in ball mills, Powder Technol. 28 (1981) 83–90. [5] A. Heim, R. Leszczyniecki, Determination of parameters for wet-grinding model in perl mills, Powder Technol. 41 (1985) 173–179. [6] S.L.A. Hennart, W.J. Wildeboer, P. van Hee, G.M.H. Meesters, Identification of the grinding mechanisms and their origin in a stirred ball mill using population balances, Chem. Eng. Sci. 64 (2009) 4123–4130. [7] K. Sedlatschek, L. Bass, Contribution to the theory of milling processes, Powder Metall. Bull. 6 (1953) 148–153. [8] B. Epstein, The mathematical description of certain breakage mechanisms leading to the logarithmico-normal distribution, J. Franklin Inst. 244 (1947) 471–477. [9] B. Epstein, Logarithmico-normal distribution in breakage of solids, Ind. Eng. Chem. 40 (1948) 2289–2291. [10] S.R. Broadbent, T.G. Callcott, A matrix analysis of processes involving particle assemblies, Phil. Trans. R. Soc. (London) A249 (1956) 99–123. [11] S.R. Broadbent, T.G. Callcott, Coal breakage processes: I. A new analysis of coal breakage processes, J. Inst. Fuel. 29 (1956) 524–528. [12] S.R. Broadbent, T.G. Callcott, Coal breakage processes: II. A matrix representation of breakage, J. Inst. Fuel. 29 (1956) 528–539. [13] S.R. Broadbent, T.G. Callcott, Coal breakage processes III. The analysis of a coal transport system, J. Inst. Fuel (London) 30 (1957) 13–17. [14] A.M. Gaudin, T.P. Meloy, Model and a comminution distribution equation for single fracture, Trans. AIME. 223 (1962) 40–43. [15] M.J. Hounslow, The population balance as a tool for understanding particle rate processes, Kona 16 (1998) 179–193.

204

M. Capece et al. / Powder Technology 208 (2011) 195–204

[16] C. Varinot, S. Hiltgun, M.-N. Pons, J. Dodds, Identification of the fragmentation mechanisms in wet-phase fine grinding in a stirred bead mill, Chem. Eng. Sci. 52 (1997) 3605–3612. [17] E. Bilgili, R. Hamey, B. Scarlett, Nano-milling of pigment agglomerates using a wet stirred media mill: elucidation of the kinetics and breakage mechanisms, Chem. Eng. Sci. 61 (2006) 149–157. [18] L.G. Austin, V.K. Bhatia, Experimental methods for grinding studies in laboratory mills, Powder Technol. 5 (1972) 261–266. [19] L.G. Austin, P.T. Luckie, Methods for determination of breakage distribution parameters, Powder Technol. 5 (1972) 215–222. [20] H. Berthiaux, C. Varinot, J. Dodds, Approximate calculation of breakage parameters from batch grinding tests, Chem. Eng. Sci. 51 (1996) 4509–4516. [21] A. Devaswithin, B. Pitchumani, S.R. de Silva, Modified back-calculation method to predict particle size distributions for batch grinding in a ball mill, Ind. Eng. Chem. Res. 27 (1988) 723–726. [22] P. Purker, R. Agrawal, P.C. Kapur, A G-H scheme for back-calculation of breakage rate functions from batch grinding data, Powder Technol. 45 (1986) 281–286. [23] R.R. Klimpel, L.G. Austin, Determination of selection-for-breakage functions in the batch grinding equation by nonlinear optimization, Ind. Eng. Chem. Fundam. 9 (1970) 230–237. [24] H. Berthiaux, J. Dodds, A new estimation technique for the determination of breakage and selection parameters in batch grinding, Powder Technol. 94 (1997) 173–179. [25] K.J. Reid, A solution to the batch grinding equation, Chem. Eng. Sci. 20 (1965) 953–963. [26] P.C. Kapur, P.K. Agarwal, Approximate solutions to the discretized batch grinding equation, Chem. Eng. Sci. 25 (1970) 1111–1113. [27] L.G. Austin, J. Shah, J. Wang, E. Gallagher, P.T. Luckie, An analysis of ball-and-race milling. Part I. The Hardgrove mill, Powder Technol. 29 (1981) 263–275. [28] R.K. Rajamani, D. Guo, Acceleration and deceleration of breakage rates in wet ball mills, Int. J. Miner. Process. 34 (1992) 103–118. [29] T.P. Meloy, M.C. Williams, Problems in population balance modeling of wet grinding, Powder Technol. 71 (1992) 273–279. [30] R. Verma, R.K. Rajamani, Environment-dependent breakage rates in ball milling, Powder Technol. 84 (1995) 127–137.

[31] E. Bilgili, B. Scarlett, Population balance modeling of non-linear effects in milling processes, Powder Technol. 153 (2005) 59–71. [32] E. Bilgili, J. Yepes, B. Scarlett, Formulation of a non-linear framework for population balance modeling of batch grinding: beyond first-order kinetics, Chem. Eng. Sci. 61 (2006) 33–44. [33] C. Tangsathitkulchai, Acceleration of particle breakage rates in wet batch ball milling, Powder Technol. 124 (2002) 67–75. [34] D.W. Fuerstenau, A.-Z.M. Abouzeid, Effect of fine particles on the kinetics and energetics of grinding coarse particles, Int. J. Miner. Process. 31 (1991) 151–162. [35] L.G. Austin, K. Shoji, P.T. Luckie, The effect of ball size on mill performance, Powder Technol. 14 (1976) 71–79. [36] L.G. Austin, M. Yekeler, T. Dumm, R. Hogg, Kinetics and shape factors of ultrafine dry grinding in a laboratory tumbling ball mill, Part. Part. Syst. Charact. 7 (1990) 242–247. [37] R.R. Klimpel, L.G. Austin, The back-calculation of specific rates of breakage and non-normalized breakage distribution parameters from batch grinding data, Int. J. Miner. Process. 4 (1977) 7–32. [38] E. Bilgili, On the consequences of non-first-order breakage kinetics in comminution processes: absence of self-similar size spectra, Part. Part. Syst. Charact. 24 (2007) 12–17. [39] L.F. Shampine, M.W. Reichelt, The Matlab ODE suite, SIAM J. Sci. Comput. 18 (1997) 1–22. [40] R.H. Byrd, M.E. Hribar, J. Nocedal, An interior point algorithm for large scale nonlinear programming, SIAM J. Optim. 9 (1999) 877–900. [41] R. Aster, B. Borchers, C. Thurber, Parameter Estimation and Inverse Problems, Elsevier Academic Press, 2005. [42] D.M. Scott, A. Boxman, C. Jochen, In-line particle characterization, Part. Part. Syst. Charact. 15 (1998) 47–50. [43] E. Bilgili, B. Scarlett, Nonlinear effects in particulate processes, Nonlinear Anal. 63 (2005) e1131–e1141. [44] E. Bilgili, B. Scarlett, Numerical simulation of open-circuit continuous mills using a non-linear population balance framework: incorporation of non-first-order effects, Chem. Eng. Technol. 28 (2005) 153–159. [45] E. Bilgili, A unified approach to mathematical treatment of multi-particle interactions during size reduction: extending the nonlinear population balance model, AIChE Annu. Meet., Paper No: 712a, Philadelphia, PA, Nov., 2008.