On the Similarity of Austin Model and Kotake–Kanda Model and Implications for Tumbling Ball Mill Scale-up †

The aim of this theoretical investigation is to seek any similarities between the Austin model and the Kotake–Kanda (KK) model for the specific breakage rate function in the population balance model (PBM) used for tumbling ball milling and assess feasibility of the KK model for scale-up. For both models, the limiting behavior for small particle size-to-ball size ratio and the extremum behavior for a given ball size are described by “power-law.” Motivated by this similarity, specific breakage rate data were generated using the Austin model parameters obtained from the lab-scale ball milling of coal and fitted by the KK model successfully. Then, using the Austin’s scale-up methodology, the specific breakage rate was scaled-up numerically for various mill diameter scale-up ratios and ball sizes of 30–49 mm and coal particle sizes of 0.0106–30 mm. PBM simulations suggest that the KK model predicts identical evolution of the particle size distribution to that by the Austin model prior to scale-up. Upon scale-up, the differences are relatively small. Hence, modification of the exponents in the Austin’s scale-up methodology is not warranted for scale-up with the KK model. Overall, this study has established the similarity of both models for simulation and scale-up.


Introduction
Tumbling ball mills have been commonly used in a multitude of industries, especially the minerals industry . Population balance modeling (PBM) enables engineers to simulate, design, control, and scale-up various particulate processes (Randolph and Larson, 1988;King, 2001), including the tumbling ball milling Prasher, 1987). The size-discrete, time-continuous form of the PBM (Sedlatschek and Bass, 1953) for a wellmixed ball mill operating in batch mode is expressed as where i and j are the size-class indices running up to N s (sink size class). Here, t, M i , S i , and b ij represent time, mass fraction, specific breakage rate parameter, and breakage distribution parameter, respectively. Eqn. (1) also assumes first-order (linear) breakage kinetics and has been successfully used for simulating the milling of a wide variety of materials in different types of mills King, 2001). Readers are referred to Bilgili and Scarlett (2005) and Capece et al. (2011) for an extensive review and a PBM framework for treatment of nonlinear breakage kinetics. It has been experimentally well-established that S i varies significantly with design and operating conditions for a given material to be ground in ball mills . The celebrated Austin model describes how S i varies as a function of particle size x i : Here, a T and μ T are parameters that strongly depend on the milling conditions, the material, and the ball size, whereas α and Λ are positive constants that largely depend on the properties of the material. Austin et al. (2007) assert that Λ can be satisfactorily taken as 3 for most materials. Some researchers assumed Λ = 3 (e.g., Bwalya et al., 2014) or fitted Λ to find it within ±10 % of this value (Rogers et al., 1986;Petrakis et al., 2017); yet others found notably different values: 4.74 (Chimwani et al., 2013) and 2.4 (Shahcheraghi et al., 2019). Traditionally, these parameters have been estimated via either direct calculation of S i by one-size-fraction method (Austin and Bhatia, 1971-72) or the optimization-based back-calculation method (Klimpel and Austin, 1977) using experimental data from small or pilot scale ball mills. As balls are an integral element of the ball mill design, the impact of ball size d B on a T and μ T has been extensively studied (Austin et al., 1976a;Katubilwa and Moys, 2009). The Austin model can be recast into the following form with explicit ball size dependence (e.g., Katubilwa and Moys, 2009, with slight differences in formalism): (3) where a 0 and μ 0 are the values of a T and μ T for specific ball size d B,0 in the small/pilot scale batch milling experiments, while η and ξ are constant ball-size exponents. For a mixture of balls, the specific breakage rate is calculated from the mass fraction of each ball size M B, p and its corresponding S i, p as follows: Scale-up of tumbling ball mills from small-scale batch mills to industrial scale mills has been successfully implemented within the context of PBM Yildirim et al., 1999;Chimwani et al., 2014;De Oliveira and Tavares, 2018). In this approach, empirical equations are used for determining the dependence of the breakage parameters on ball/mill dimensions and operating variables at the lab-scale and large-scale industrial mills. It is important to emphasize that scale-up from a batch mill to continuous large-scale mills entails modification of the breakage parameters (S i and b ij ) as well as consideration of the residence time distribution and internal-external classification in the continuous mills (Herbst and Fuerstenau, 1980;Austin et al., 1984;Rogers and Austin, 1984;De Oliveira and Tavares 2018). Although each one of these aspects is important, this study focuses on the scale-up of the specific breakage rate S i .
Upon scale-up, S i is known to change significantly, whereas b ij or its cumulative counterpart B ij is generally assumed to remain invariant unless ball size/type changes significantly (Herbst and Fuerstenau, 1980;Austin et al., 1984). In accounting for scale-up related changes, the Austin's methodology scales-up S i in Eqn. (2) via correction factors K 1 -K 5 that depend on the design parameters of the dry ball mill and operating conditions as follows: where D, J, U, and ϕ C denote the mill diameter, the fractional mill filling by the balls, the void filling fraction, and the fraction of actual rotation speed compared to the critical speed N C . The subscript T refers to the parameters used in the small-scale batch test, which are used to fit the model parameters in Eqns. (2) and (3). Typical scale-up exponents are N 0 = 1, N 1 = 0.5, N 2 = 0.2, N 3 = 2, and N 4 = 0.2 . However, other values of N 3 were also used such as 1 (Yildirim, 1999) and 1.2 (Austin et al., 2007). One should note that N 0 and N 3 should be taken as ξ and η of Eqn.
(3), respectively, from a study on the ball size impact on S i at the test (T) scale, which is a preferred approach (Mulenga, 2017). While the Austin model in Eqn.
(2) or (3) has been the most widely used kinetic model for ball milling, Kotake et al. (2002Kotake et al. ( , 2004 developed the following kinetic model for the breakage of narrowly sized feed with size x f : which we briefly refer to as the KK (Kotake-Kanda) model for the sake of simplicity. The parameters C 1 and C 2 strongly depend on the milling conditions and material properties, whereas α is a positive constant that largely depends on the properties of the material; m and n are material-dependent, ball size exponents. Kotake et al. (2002Kotake et al. ( , 2004 and Deniz (2003) demonstrated the fitting capability of this model for ball milling of a multitude of minerals. The KK model, Eqn. (10), is appealing as it has 5 parameters as opposed to the Austin model with 6 parameters, when ball sizes are explicitly considered. The milling studies, where the back-calculation approach involved more than a few breakage parameters to be estimated, have shown that the accuracy decreases as the number of parameters to be estimated increases (Klimpel and Austin, 1977;Kwon and Cho, 2021). However, to the best knowledge of the author, the KK model has not been used for the PBM simulation of ball milling of a natural feed with a wide size particle size distribution (PSD), and no scale-up has been performed using it. More importantly, the KK model has never been compared to the Austin model. The present theoretical investigation seeks answers to the following questions: (i) Are there any similarities between the Austin model and the KK model despite their different mathematical forms? (ii) Is the Austin's scale-up methodology (scale-up correction factors) applicable to the KK model? (iii) Do these models predict similar time-wise evolution of PSD in ball mills before and after scale-up? To be able to answer these questions, we first generated synthetic specific breakage rate data using the Austin model, whose parameters were obtained from the lab-scale ball milling of coal. Next, this data was fitted by the KK model to determine its parameters. Then, using the Austin's scale-up methodology, the specific breakage rate was scaled-up numerically (virtual scale-up) using both models for various mill diameter scale-up ratios D/D T from 4 to 10, ball sizes d B of 30-49 mm, and coal particle sizes x i of 0.0106-30 mm. PBM simulations assumed plug flow and invariance of b ij in the large-scale mills and incorporated the scaled-up S i values. The PBM simulations were performed for the small-scale mill and the large-scale mills. Upon scale-up, we find that the differences in both models' predictions are relatively small. Instead of just using N 1 = 0.5 and N 2 = 0.2 for the virtual scale-up with the KK model, which is based on the Austin's scale-up methodology, we also fitted these exponents to the S i values generated by the Austin kinetic model, keeping the other scale-up exponents the same. Comparison of the PSD evolution prediction by the KK model (with/without modified N 1 and N 2 ) and the Austin model will allow us to find the answers to the above-mentioned questions.

Experimental data and model calibration
The KK model can be cast into the following form, which is valid for any arbitrary particle size x i Note that both Eqns. (3) and (11) are explicit in d B and x i , which enables us to discern any similarities or differences between the two models by analyzing the limiting behavior for small particle size-to-ball size ratio and the extremum behavior for a given ball size. To demonstrate that the KK model can predict the same S i as that predicted by the Austin model, we need to have S i parameters from an experimental data set on batch ball milling. Hence, we considered the experimental work on the batch ball milling of a South African coal (Katubilwa and Moys, 2009). They fitted and/or set the Austin model parameters as follows: a 0 = 0.48 mm -0.81 •min -1 , μ 0 = 19.27 mm, d B,0 = 38.8 mm, α = 0.81, Λ = 3, ξ = 1, and η = 1.96. Using this set of parameters, a synthetic data set for S i was generated by the Austin model using Eqn. (3) for ball sizes of 30-49 mm with increments of 1 mm and coal particle sizes of 0.0106-30 mm with a geometric progression ratio of 2 1/13 . As the purpose of this study is to examine theoretical similarity of the KK model to the Austin model, we did not fit the KK model to the experimental data directly, but to the synthetic S i data generated by the Austin model to minimize any bias. Both models have been shown to fit various ball milling data successfully in the milling literature; the question we seek to answer here is not how well they can fit specific experimental data. Therefore, the parameters of the KK model, Eqn. (11), were estimated by minimizing the sumof-squared residuals between the KK model prediction and the synthetic data generated by the Austin model for the South African coal. In comparison to no transform and logarithmic transform, a square-root transform of the response (S i ) led to the best overall fitting and simulation results when the KK model with the estimated parameters were compared with the Austin model. Hence, only the KK model parameters obtained by fitting after the square-root transform of the response (S i ) are presented. The optimization was performed via the Marquardt-Levenberg method in Minitab version 21.1.

Virtual scale-up
A virtual scale-up of the specific breakage rate S i was performed from the small-scale batch mill to the largescale continuous ball mills. Using the Austin's scale-up methodology, the specific breakage rate S i obtained experimentally for the Austin model (Katubilwa and Moys, 2009) and the KK model with the estimated parameters (refer to Section 2.1) were scaled-up numerically for various mill diameter scale-up ratios D/D T from 4 to 10, ball sizes d B of 30-49 mm, and coal particle sizes x i of 0.0106-30 mm. For simplicity, K 4 and K 5 are set to 1 in this study for the scale-up with both models regardless of any particular values of J, U, and ϕ c values because K 4 and K 5 are independent of d B /d T ≡ d B /d B0 and D/D T , which were varied in this study. Hence, the scaled-up S i , i.e., S i *, was calculated by The above equations incorporate K 1 , K 2 , and K 3 directly, which depend on d B /d B0 and D/D T . First, Eqn. (12) was used to generate the S i * by the Austin model for all scale-up scenarios with various d B /d B0 and D/D T . These scaled-up data were taken as the "synthetic scale-up data" against which the scale-up with the KK model was compared.
Then, two different virtual scale-up approaches for the KK model were used. In the first approach, the standard scale-up exponents of the Austin's scale-up methodology, i.e., N 1 = 0.5 and N 2 = 0.2, along with N 4 = 0.2 were used in Eqn. (13) to predict the S i *. However, as there is no prior scale-up study with the KK model, one cannot assume that the Austin's scale-up exponents N 1 and N 2 are directly applicable for use in Eqn. (13). Hence, in the second approach, we fitted S i * values generated by Eqn. (12) for each D/D T from 4 to 10 and ball sizes d B of 30-49 mm using Eqn. (13) and minimizing the following standard error SE: which allowed us to estimate N 1 and N 2 . An additional fit was performed to minimize SE for all D/D T cases simultaneously. Here, N and N p refer to the total number of data points and parameters estimated, respectively. Then, these simultaneously estimated values of N 1 and N 2 were used to predict S i * for each scale-up scenario separately.

PBM simulations
The scale-up from a small batch mill or a pilot-scale mill to a large-scale industrial continuous mill is a complex endeavor, entailing considerations of the residence time distribution (RTD), internal classification, and external classification in the continuous mill circuits (Herbst and Fuerstenau, 1980;Austin et al., 1984;Rogers and Austin, 1984;De Oliveira and Tavares, 2018) besides scale-up of the specific breakage rate. However, any consideration of these factors in the scale-up process would complicate the analysis to discern the differences between the Austin model and KK model in the scale-up of S i . Hence, we followed the simplest approach for the modeling of an open-circuit continuous mill by assuming plug flow behavior and disregarding any internal classification in the mill. This simple approach purposefully isolates the problem of "breakage kinetics" in scale-up and allows us to focus on the discrimination of the KK model from the Austin model within the context of simulation and scale-up. The plug-flow is a viable assumption for some industrial scale ball mills (Austin, 1973). For example, a 10 tanks-in-series model represented the RTD data for a large-scale cement ball mill (Austin et al., 1975). The plug-flow assumption was also adopted by other researchers in the ball mill scale-up (Chimwani et al., 2014). Developing a sophisticated PBM for a continuous ball mill circuit with internalexternal classification is not within the scope of this study.
Note that Eqn.
(1) is applicable to mills operating in batch mode as well as plug-flow continuous mode. One can replace time t with contact time t* = y/u for any axial position y in the mill with u denoting the average axial velocity of the particles: u = L/τ = LF/M h , where τ, L, F, and M h stand for the space-time (the average residence time), length of the mill, inlet mass flow rate, and mass hold-up at steady state, respectively. So, after replacing t with t*, one can solve the set of ordinary differential equations (ODEs) in Eqn. (1) with the initial condition M i (0) = M i,ini , with M i,ini representing the feed PSD, and predict the steadystate PSD at t* = τ. Henceforth, although we use the above formalism, for the sake of simplicity, we will still denote time as t in the results, but keeping in mind, t represents retention time in the batch mill and contact (residence) time in the plug-flow continuous mill at the steady state.
In the PBM simulations, we used the S i values in Section 2.1 for the small-scale batch mill simulations and S i * values, as described in Section 2.2, for the continuous milling. The following cumulative breakage distribution parameters B ij were taken from Katubilwa et al. (2011) about the same South African coal for which S i parameters were reported in Section 2.1: Here, 0 < ϕ < 1 is a breakage constant, and γ and β are breakage exponents. For the South African coal, ϕ, γ, and β were reported to be 0.51, 0.53, and 3.2, respectively. It is well-known that these parameters are largely material dependent, and are relatively insensitive to the milling conditions, except the ball size Prasher, 1987;King, 2001). In a more elaborate model, one could incorporate the impact of ball size on B ij by making γ dependent on the ball size (Austin et al., 2007). In the present simulations and scale-up, we assumed that B ij remain invariant to ball size and scale-up, similar to what Chimwani et al. (2014) implemented. Using the B ij values from Eqn.
The initial PSD for the small batch mill simulation and the feed PSD for the continuous mill simulation were taken the same: a Gaussian PSD with a mean size of 20 mm and a standard deviation of 2 mm. The Gaussian PSD was generated using the function "normpdf" in MATLAB version 9.11. Considering the time scale of the batch milling experiments (Katubilwa and Moys, 2009), we simulated 0.5, 1, 2, 4, and 8 min of milling time (or residence time). Another rationale for the selection of 8 min is the following: powerstation boilers usually demand 60-70 % of coal particles below 75 μm (Prasher, 1987). Our exploratory simulations suggest 8 min residence time in the large-scale mills could meet this demand approximately. In all PBM simulations, the number of size classes N s and the geometric progression ratio were set as 320 and 2 1/13 , respectively, which yields grid-independent simulation results. The set of ODEs in Eqn.

Austin model vs. Kotake-Kanda model
Despite the differences in their mathematical forms, both the Austin model and the KK (Kotake-Kanda) model have similar limiting behavior and extremum behavior. Let us first consider the limiting case of x i << d B , for which Eqns.
(3) and (11) reduce to respectively, because the dominator of Eqn. (3) and exponential term in Eqn. (11) approach 1. Clearly, both models provide power-law dependence of S i on x i and d B in this limit. Note that x i << d B has not been strictly defined and x i /d B < O(10 -2 ) is reasonable for this limiting behavior. To prevent a possible confusion, we assert that this limiting behavior is not a practically irrelevant asymptotic limit. As demonstrated by Austin et al. (1976b) and De Oliveira and Tavares (2018), the whole milling kinetics could be governed by this limiting behavior if the feed or initial size is much smaller than the ball sizes used. The KK model can be made equivalent to the Austin model by imposing as well as the sameness of α in both equations in terms of the limiting behavior.
Let us now consider the extremum behavior of S i by setting dS/dx = 0 after taking the equivalent particle size-continuous forms of Eqns. (3) and (11). This analysis yields the particle size x m at which S i is at a maximum and has the following S m : for the Austin model and the KK model, respectively. Similar to the limiting behavior, both models predict a power-law dependence of x m and S m on d B . However, the exponents for the limiting behavior and the extremum behavior are different. One can make both models equivalent in terms of the extremum behavior alone, i.e., identical x m and S m , by setting the KK parameters as n = η, m = -ξ, and same α in both equations. An interesting theoretical question arises as to what happens if the equivalence conditions, i.e., either Eqn. (17) for the limiting behavior or Eqns. (22)-(23) for the extremum behavior are imposed on the KK model. We note that the conditions imposed on C 1 are incompatible because e α (1-α/Λ) = 1 is overly restrictive and cannot be satisfied, except for α = 0, which is untenable. If one imposes the equivalence in terms of the extremum behavior, then C 1 provided by Eqn. (22) will be similar to C 1 provided by Eqn. (17) in terms of the power-law dependence with the same exponent ξ. However, Eqn. (22) will overpredict S i by a factor of e α (1-α/Λ) in the limit x i << d B . This factor equals 1.641 (64.1 %) for the coal considered in Section 2.1, which is not acceptable.
The upshot of the above theoretical analysis is that the KK model is similar, but not equivalent or identical, to the Austin model in terms of the limiting behavior and the extremum behavior. It is impossible to ensure equivalency of these models in terms of both limiting behavior in the size domain x i << d B and the extremum behavior in the neighborhood of x i = x m . Owing to the power-law dependence of S i on x i and d B in the size domain x i << d B and that of x m and S m on d B , we hypothesize that both models could predict similar S i while satisfying the limiting behavior and the extremum in some approximate statistical sense. Note that they may differ significantly in the size domain x i << d B and near the extremum. This hypothesis can be tested by first fitting the KK model parameters to the S i generated by the Austin model and then analyzing the resulting S i profiles as well as these models' prediction of the temporal evolution of the PSD. To this end, the synthetic S i data generated by the Austin model, whose parameters were obtained by fitting to the ball milling data on a South African coal (refer to Section 2.1), were fitted by the KK model. The parameters of the KK model are presented in Table 1. The sum-of-squared errors was 4.08 min -1 , and the standard error SE of the fit was 3.69 × 10 -2 min -0.5 . The parameters had low standard errors: all below 8 %. Clearly, the KK parameters in Table 1 will not satisfy the requirements of Eqn. (17) and Eqns. (22)- (23), which, respectively, ensure equivalency in terms of the limiting behavior and the extremum behavior. These strict conditions are not required for the models to be similar both qualitatively and quantitatively when experimental errors in determining S i are considered (see later discussion). Fig. 1 presents the specific breakage rate S i data generated by the Austin model for the South African coal and its fitting by the KK model, which led to the estimated parameters in Table 1. While there are minor deviations, the KK model estimated the S i profiles and the extremum behavior reasonably well. The KK model exhibited a slightly depressed, shifted peak as compared with the Austin model. Both models predict a cliff or falling off S i for x i > x m , which can be explained by the inability of the balls to nip or capture the particles effectively. In view of Eqn. (16), the values of ξ = 1 and m = -1.444 suggest that the smaller ball sizes are effective for x i << d B because for the smaller particles, the number frequency of the collisions with the beads govern the specific breakage rate, which is suggested to be scaled by the inverse of d B , as signified by ξ = 1 . Overall, both models' phenomenological representation of the breakage kinetics appears to be qualitatively and quantitatively similar. The location of the extremum points in Fig. 1 are given by x m -S m , which depend on the ball size d B and are obtained via Eqns. (18) (24) and (25). According to the Austin model, each ball size most effectively breaks a certain particle size. Large balls are the best for breaking coarse particles; small balls are the best for breaking small particles. For example, 46 mm balls are most effective for breaking ~20 mm particles; yet, they are the worst for breaking <~8 mm particles; 30 mm balls are most effective for the latter. Practically, for a feed with wide PSD, the common practice is to use a mixture of different ball sizes Prasher, 1987;Katubilwa and Moys, 2009). As can be seen from Fig. 2, the KK model provides a similar ball size dependence of S m and x m to that by the Austin model although there is some deviation. The mean relative deviation for S m and x m were 8.18 % and 7.31 %, respectively, while the maximum relative deviation for S m and x m were 13.2 % and 16.1 %, respectively. One interesting trend is that as d B increased, the S m deviation of the KK model from the Austin model decreased, while the x m deviation increased, which counter each other. This observation along with <10 % mean relative deviation overall imply that the deviations may have small impact on the PBM simulation and scale-up results. This is indeed the case, as will be demonstrated next through PBM simulations before and after scale-up.
To develop a deeper understanding of the ball size effects, one must resort to Discrete Element Modeling (DEM) of tumbling ball mills. DEM provides significant fundamental insights and microdynamic information such as collision energies among the balls-particles-mill wall/liners and their frequencies (Tavares, 2017;Rodriguez et al., 2018). In an excellent critical analysis, Rodriguez et al. (2018) pointed out some fundamental learnings from the DEM: (i) Collision energies involving events with only particles (particle-particle and particle-liner) are not sufficiently high to cause breakage of particles and (ii) a fraction of the collisions inside the ball mill do not involve particles; hence, they do not directly contribute to particle breakage. In agreement with the above, the multi-scale PBM-DEM modeling of the experimental tumbling ball mill (Kotake et al., 2002) by Capece et al. (2014) predicted the S i profile similar to those in Fig. 1. Moreover, their simulations suggest that for a given narrow feed size, an effective ball size that maximizes the specific breakage rate exists. This effective ball size appears to maximize collision frequency with the particles (ball-particle collisions) with sufficiently high impact energy above a minimum threshold impact energy that is particle size-dependent. Such minimum impact energies required to break individual particles under impact loading should be determined experimentally (Tavares and King, 1998;Vogel and Peukert, 2003;Tavares, 2007).
While the analysis of the data in Figs. 1 and 2 suggests that the KK model is similar to the Austin model and they both provide similar S i , x m , and S m , one would argue about the deviations. At this juncture, it is imperative to mention the errors involved in the determination of S m and     x m . For example, Kotake et al. (2002;2004) described an averaging procedure to determine S m and x m because of the experimental scatter around the maxima. The experimental deviations were certainly greater than the <~10 % S m and x m deviations of the KK model from the Austin model observed here. Moreover, the 85 % of the fitted S 1 values by the KK model were within a band of ±20 %, with 15 % of the values outside this band (Kotake et al., 2004); similarly notable deviations can be seen in Deniz (2003). Most importantly, a long-due study (Shahcheraghi et al., 2019) clearly demonstrated that at 95 % confidence level, the experimentally determined S i parameters of the Austin model had 21.7-34.1 % uncertainty, whereas the B ij parameters had 7.8-17.8 % uncertainty. Hence, the deviations shown in Figs. 1 and 2 between the Austin model and the KK model would most likely be within the range of experimental uncertainty when the two models are fitted to the same experimental data with repeats. This has never been done in the experimental ball milling literature. In this study, we will confirm the similarity of the two models via PBM simulations next. A batch ball milling process was simulated by PBM, as described in Section 2.3 to assess if the deviations between the two models' S i values could cause significant difference in the simulated time-wise evolution of the PSD. To the best knowledge of the author, this is the first study in which a batch milling process has been simulated by the KK model. Fig. 3 clearly illustrates that despite the S i deviations illustrated in Fig. 1, the two models simulated the temporal evolution of the PSD almost identically. The maximum deviations in the PSDs estimated by the KK model and the Austin model are within a few percentages. In view of the above discussion and the simulation results, it is clear that the PBM with either the KK model or the Austin model yields a similar temporal evolution of the PSD in a batch mill prior to scale-up. The KK model has one less parameter than the Austin model (5 vs. 6); its use may be advantageous for parameter estimation.

Virtual scale-up with the KK model
Historically, the Austin kinetic model has been used along with Austin's methodology for ball mill scale-up . As the KK model has never been used for ball mill scale-up, it is doubtful if it can be integrated with the Austin's scale-up methodology. In Section 2.2, we developed the scaled-up specific breakage rate parameter S i * using this methodology. To the best knowledge of the author, this is the first attempt to use the KK model for scale-up, albeit virtual. This virtual scale-up will allow us to assess if the KK model could provide (or fit to) S i * similarly to the Austin model upon scale-up. Table 2 presents the virtual scale-up scenarios and the standard error SE, as per Eqn. (14). Four diameter scale-up ratios D/D T with ball sizes of 30-49 mm and coal particle sizes of 0.0106-30 mm were considered. The Austin model,via Eqn. (12), was used to produce the synthetic scale-up S i * data to which the S i * generated by the KK model, via Eqn. (13), will be compared (SU1-4 in Table 2). In SU1-SU4, standard values of the scale-up correction exponents, i.e., N 1 = 0.5 and N 2 = 0.2 were used. In SU5-SU8, N 1 and N 2 were estimated by fitting Eqn. (13) to the synthetic scale-up Scale-up no. The SE of simultaneous fitting with SU9-12 was 5.66 × 10 -2 min -0.5 . Initial PSD data for each D/D T separately. In SU9-SU12, they were fitted to the synthetic scale-up data for all D/D T data together.
Fitting of the scaled-up specific breakage rate parameter S i * by the KK model to the synthetic scale-up S i * data (Austin model) either one-at-a-time for each scale-up scenario (SU5-8) or simultaneously for all scenarios (SU9-12) led to almost identical exponents of the scale-up correction factor: N 1 = ~0.46 and N 2 = ~0.24, which deviated only by 8 % and 20 % from the standard values of the Austin's scale-up methodology, i.e., N 1 = 0.5 and N 2 = 0.2 (SU1-4). Hence, one set of modified N 1 and N 2 can be used for the scale-up of S i of the KK model. When the SE of SU9-12 was compared with that of SU1-SU4, the modified N 1 and N 2 reduced SE by 12 %, 19 %, 24 %, and 28 %, respectively, for D/D T values of 4, 6, 8, and 10.
At the small-scale, the fitting of the KK model to the synthetic S i data for its parameter estimation has a SE of 3.69 × 10 -2 min -0.5 . Upon scale-up, the SE values in Table 2 for SU1-4 are 57 %-130 % higher than that associated with the parameter estimation. The situation is somewhat better when the newly estimated, modified N 1 = 0.455 and N 2 = 0.241 are used. The SE values for SU9-12 (upon scale-up) are 37 %-65 % higher than that associated with the parameter estimation. Clearly, the S i differences between the KK model and the Austin model become more pronounced upon scale-up whether the standard Austin scale-up correction exponents or the modified exponents are used. While the SE increased upon scale-up, as intuitively expected, in absolute terms, the SE values are still low which will become apparent once the deviations are presented in visual form next.
Figs. 4 and 5 present the synthetic S i * data by the Austin model and S i * provided by the KK model (SU1 and SU3) and fitted by the KK model (SU9 and SU11). The data from the other scale-up scenarios were not presented in graphical form for the sake of brevity. It suffices to state that Figs. 4 and 5 are representative of the general scale-up trends for all scenarios.
A cursory look at Figs. 4 and 5 vs. Fig. 1 reveals the remarkable impact of the scale-up on the specific breakage rate. Upon scale-up and an increase in D/D T , the specific breakage rate increases, which accords well with the experimental ball milling literature (Austin, 1973;. N 1 captures the increase in mill power as mill diameter increases (Austin, 1973). S m and x m also increase, which suggests the maximum point moves to the right and its peak is heightened. With an increase in D/D T , the cliff after x i > x m becomes less pronounced and the abnormal breakage region shrinks. In fact, for d B = 46 mm, the maximum point almost disappears. Obviously, for a higher D/D T value, the balls are dropped from a greater height, and they hit the particles at a much higher impact force. According to a semi-theoretical analysis (Hayashi and Tanaka, 1972), the impact force for free-falling balls scales with D 0.5 . A recent advanced DEM simulation study (De Carvalho et al., 2021) provides a more fundamental explanation. Their DEM simulations point out a significant expansion of the collision frequency-collision energy envelope as the ball milling is carried out at larger scales: from batch to pilot and then to industrial mills. They determined the specific collision frequency to be 2404 collisions/kg•s, 2746 collisions/kg•s, and 3795 collisions/kg•s at the batch, pilot, and industrial scales, respectively. Overall, higher impact forces, collision energies, and specific collision frequency explain the higher specific breakage rates upon scale-up fundamentally.
The differences in S i * between the Austin model and the KK model, with/without the modified N 1 and N 2 exponents, illustrated in Figs. 4 and 5 seem to be unremarkable and most likely statistically insignificant when actual ball milling data are considered. The readers are referred to the discussion and the references in Section 3.1. However, proving this point with comparative experimental study is beyond the scope of this study. We pose a different, yet quite relevant question: can the S i * differences depicted in Figs. 4 and 5 translate into significant differences in terms of the PSD evolution with residence time in a continuous mill? Fig. 6 illustrates that the PSD shifts from the Gaussian feed PSD to the left monotonically along the axial direction of a continuous mill with D/D T = 4, which is scaled from the small batch mill whose PSD evolution is shown in Fig. 3. The time refers to retention time in the batch mill (Fig. 3) and the residence time in the continuous mill ( Fig. 6), which corresponds to different axial locations in the mill. A cursory look at Fig. 6 vs. Fig. 3 reveals that the product becomes much finer upon scale-up, which is in agreement with the higher specific breakage rate parameter upon scale-up (S i * in Fig. 4 vs. S i in Fig. 1). A larger diameter mill with D/D T = 8 (Fig. 7) exhibits similar behavior; the PSDs are finer as compared with those in Fig. 6, which again can be attributed to the respective S i * values in Figs. 4 and 5.
Let us now answer the question posed earlier: the Austin model and the KK model, with either the standard scale-up exponents N 1 = 0.5 and N 2 = 0.2 or the modified scale-up exponents estimated in this study, i.e., N 1 = 0.455 and N 2 = 0.241, predict very similar PSDs in a plug flow continuous mill. Although the differences increased slightly upon scale-up to a larger mill with D/D T = 8 (Fig. 7), both the Austin model and the KK model estimated similar PSDs, with a max. deviation of 4 % at 10 μm in the cumulative mass fraction.
An interesting finding from the continuous mill simulations presented in Figs. 6 and 7 is that the modification of the Austin's scale-up exponents N 1 and N 2 was not warranted for the integration of the KK model with the Austin's scale-up methodology. The KK model with the modified scale-up exponents showed less deviation from the Austin model than that with the standard exponents at the earlier residence time values. However, this trend was reserved at the later residence time points and at the exit of the mill (t = τ = 8 min). Despite the lower SE values of the S i * associated with the modified scale-up exponents in Table 2, these differences in the SE values did not seem to make a significant difference in terms of the PSD change along the mill axis and the product PSD, as proven by the simulation results in Figs. 6 and 7. Longer residence times were not considered for two reasons. First, the product PSD (t = 8 min) in the large-scale continuous mills reached the desired target approximately for coal applications: 60-70 % below 75 μm. Second, the Austin model parameters were estimated using the one-size-fraction data from a ball milling study on South African coal (Katubilwa and Moys, 2009;Katubilwa et al., 2011). In these studies, the smallest feed size fraction studied contained 300-425 μm coal particles. Hence, both the Austin and the KK models extrapolate S i of the finer particle range (< 300 μm) based on the coarser particles' breakage kinetics. While this is a common practice, we suspect that the α parameter in both models could be sensitive the particle size domain in the model calibration. Hence, the slightly higher deviation of the KK model from the Austin model for particles smaller than 100 μm (refer to Fig. 7) could be indirectly linked to the model calibration and the slight differences in the α value (α = 0.81 for the Austin model and α = 0.95 for the KK model).

Concluding remarks and outlook
This theoretical study presents a comparative analysis of the Austin kinetic model and the Kotake-Kanda (KK) kinetic model in terms of their prediction of the specific breakage rate parameter as a function of particle size and ball size. Both models have been shown to be similar in terms of their limiting behavior and the extremum behavior because both behaviors are governed by powerlaw relationships. Mathematical requirements have been formulated to ensure equivalency of both models under these conditions; however, ensuring equivalency for both behaviors through these derived conditions is shown to be severely restrictive or incompatible. Hence, to investigate their similarities further, synthetic specific breakage rate S i data were generated by the Austin model and fitted by the KK model reasonably well. The PBM simulations of a batch (reference) mill have demonstrated almost identical PSD evolution, confirming the similarity of the two models. Then, a virtual scale-up from the batch reference mill to a continuous plug-flow mill with different D/D T has been performed. Using the Austin's scale-up methodology along with the Austin model, the scaled-up specific breakage rate parameter S i * has been determined. The KK model has also been scaled-up using the Austin's scale-up methodology with and without modified scale-up correction exponents. Also, the large-scale continuous mill operation has been simulated via PBM with the S i *. This study has overall demonstrated that (i) the KK model is similar, but not equivalent, to the Austin model, (ii) the KK model and the Austin model provide similar PSD evolution after scale-up, (iii) the modification of the Austin's scale-up exponents N 1 and N 2 is not warranted, and (iv) the Austin's scale-up methodology is applicable to the KK model as well as the Austin model.
This theoretical study has generated additional questions and ideas for further research. A well-designed experimental ball milling study, using both one-size fraction method and the optimization-based back-calculation method on several natural feed PSDs, should be conducted. The models should be discriminated based on their goodnessof-fit, parameter uncertainties and confidence intervals, and statistical significance. Such studies could reveal if the KK model has some advantages in terms of parameter uncertainties over the Austin model as the former has one less parameter than the Austin model. Such a future study is especially important to the calibration of the KK model (parameter estimation) as the currently used method with one-size-fraction milling experiments is conducive to accumulation of experimental errors. Finally, it is hoped that the Kotake-Kanda model will be used within the Austin's scale-up methodology for ball mill scale-up in the future.

Acknowledgements
The author would like to thank his PhD student Mr. Nontawat Muanpaopong for preparing the figures and proof-reading the manuscript. This theoretical research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Data Availability Statement
The data analysis file and the data supporting the findings of this study are available in J-STAGE Data.