Review and Further Validation of a Practical Single-particle Breakage Model †

Particle breakage occurs in comminution machines and, inadvertently, in other process equipment during handling as well as in geotechnical applications. For nearly a century, researchers have developed mathematical expressions to describe single-particle breakage having different levels of complexity and abilities to represent it. The work presents and analyses critically a breakage model that has been found to be suitable to describe breakage of brittle materials in association to the discrete element method, either embedded in it as part of particle replacement schemes or coupled to it in microscale population balance models. The energy-based model accounts for variability and size-dependency of fracture energy of particles, weakening when particles are stressed below this value, as well as energy and size-dependent fragment size distributions when particles are stressed beyond it, discriminating between surface and body breakage. The work then further validates the model on the basis of extensive data from impact load cell and drop weight tests. Finally, a discussion of challenges associated to fitting its parameters and on applications is presented.

One application of single-particle breakage information that has grown significantly in recent decades has been to support mathematical models of comminution and mechanical degradation during handling. This information was initially used as part of phenomenological models of comminution devices, including ball mills, autogenous and semi-autogenous mills and cone crushers (Narayanan and Whiten, 1988;Napier-Munn et al., 1996), as well as of degradation during handling (Teo et al., 1990;Weedon and Wilson, 2000), although originally in a fairly simplistic way. Perhaps influenced by the conditions used in the standardized twin pendulum, and later also in the drop weight tester (JKDWT) by the Julius Kruttschnitt Mineral Research Centre (JKMRC) (Napier-Munn et al., 1996), attention was initially only given to the relationship between the energy applied and the resulting fragmentation at relatively high specific impact energies.
From the beginning of the 1990s, research on the discrete element method (DEM) applied to tumbling mills (Mishra and Rajamani, 1992;King and Bourgeois, 1993;Morrison and Cleary, 2006) has demonstrated that very few impacts occur at the high specific energies that are object of testing in the JKDWT, and that the majority of collisions occur at impact energies that are either able to only impose limited breakage to particles or break them only after repeated stressing. This led to recognition of the importance of deepening the investigation of conditions that result in primary breakage, including breakage probability, weakening by repeated impacts and progeny size distributions at very low specific impact energies (Tavares, 2007). This shed new light into the subject, which had been object of relatively scattered studies before then (Rumpf, 1973;Baumgardt et al., 1975;Krogh, 1980;Schönert, 1979;Yashima et al., 1987). Nevertheless, the importance of very low energy impacts was recognized over 70 years ago by Fred Bond, when developing the standard crushability test (Bond, 1947), in which particles are repeatedly stressed at increasing impact energies until primary breakage occurs (Tavares and Carvalho, 2007).
A definition that has a central role when dealing with particle breakage is that of particle fracture energy (Baumgardt et al., 1975;King and Bourgeois, 1993). It corresponds to the stressing energy below which the core of the particle is maintained unbroken, whereas potentially damaged, while the surface of the particle may lose mass due to abrasion-like breakage when it is subjected to a stressing event. This particle fracture energy varies from particle to particle as well as with particle size and defines the threshold between massive or body breakage at higher impact energies and surface breakage that occurs at lower stressing energies (Tavares and King, 1998).
As mathematical models of different mills, crushers and handling systems evolved, incorporating information from DEM, it became evident that more detailed information on particle breakage was necessary. It is worth mentioning that several approaches circumvented the need for this information by either using data from particle bed breakage tests (Datta and Rajamani, 2002) or from the process itself, by back-calculating the selection function from data (Capece et al., 2014;Wang et al., 2012). However, they are not as generally applicable given their more limited ability to decouple material from process (Tavares, 2017). As such, the most powerful approaches used in advanced mathematical models and simulations are particle-based, that is, require information on how individual particles respond to stresses, either when they are fully resolved in the DEM simulation or when they are not. It is important to recognize that in order to be useful in such models, particle response must be known under the variety of conditions to which it will be subjected to in a comminution device or a system that can cause mechanical degradation. These typically include widely varied particle sizes, stressing intensities, geometries and rates of stressing as well as modes of stress application.
A fairly large number of mathematical models describing particle breakage has been proposed in the literature, which may be potentially of use in association to advanced mathematical models and DEM simulations. In most of them the variable of interest is the stressing energy (Tavares and King, 1998;Vogel and Peukert, 2003;Morrison et al., 2007;Bonfils, 2017;Ballantyne et al., 2017), whereas in a few the variable of interest is impact velocity (Ghadiri and Zhang, 2002;Rozenblat et al., 2012) or force (Rodnianski et al., 2019). Also, several of these models have been proposed that relate the size distribution of fragments to the stressing energy applied, not discriminating between volume or surface breakage (Ballantyne et al., 2017;Ouchterlony and Sanchidrián, 2018). These models, however, are not meant to discriminate particles individually, being only suited to describe the response of the entire batch. Other models have been proposed to describe the so-called incremental breakage, lumping together surface breakage and weakening by repeated stressing (Morrison et al., 2007;Bonfils et al., 2016), being very limited when used in association to advanced mathematical models of comminution and degradation. Yet other models focused exclusively on surface breakage (Ghadiri and Zhang, 2002).
The present work analyzes critically and also presents further evidence of the validity of a comprehensive mathematical model proposed by the author and his co-workers over the last 25 years (Tavares, 1997;Tavares and King, 1998;Tavares and King, 2002;Tavares, 2009;Carvalho et al., 2015;Tavares and Chagas, 2021) and schematically shown in Fig. 1 to describe several aspects of particle breakage.

Breakage probability 2.1 Body breakage probability distributions
Arguably, the central part of the model is the breakage probability ( Fig. 1). It represents the relationship between the applied stressing energy and the resulting proportion of particles that have undergone body or massive breakage when particles contained in a batch are stressed under identical conditions. Breakage in this context is defined arbitrarily as loss of at least 10 % of the particle mass (Dan and Schubert, 1990;Tavares, 2007). Estimates of breakage probability may be obtained experimentally from tests in which lots of particles contained in narrow size ranges are subjected to impacts by propelling them, one-by-one, against a target at a controlled velocity or dropping them under gravity. In this case the mass-specific impact energy is given by is the average velocity at the instant of collision against the target, which can be obtained, for instance, from measurements using a high-speed camera (Cavalcanti et al., 2021).
Alternatively, an estimate of breakage probability may be obtained by impacting particles using a drop weight tester. Assuming no loss of momentum during drop of the impactor (free-fall), the specific impact energy may be given by where m is the mass of the drop weight, p m is the average weight of particles in the lot tested, g the acceleration due to gravity and h the net drop height, that is, the distance between the bottom of the striker and the top of the particle.
In either case, the estimate of breakage probability for a given mass-specific impact energy E is given by the proportion broken, represented by the ratio where N is the number of stressed particles and N b is the number of particles that underwent body breakage as a result of application of the specific impact energy E k . Although Eqn.
(3) relates numbers of particles, it is assumed to be equivalent to the mass-based probability, given their narrow size range. In order to remove ambiguity in identifying which particles underwent breakage during execution of the test, each should be weighed individually prior to the test and also after each impact, to verify if loss of at least 10 % of their mass occurred. The measurement after the impact becomes obviously unnecessary if the particle suffered evident disintegration. Either when particles are projected against a target or when a drop weight impacts individual particles, the specific impact energy E k is actually a random variable, subjected to variability owing to the velocity of impact and, in the case of the later, also due to variability in particle weight. In addition to that, the proportion broken given by Eqn. (3) is only an estimate of the breakage probability, with the uncertainty reducing as the number of testparticles increase, provided no bias existed in selecting particles for the test. The magnitude of the sampling error has been estimated (Tavares, 2009) on the basis of the Binomial distribution, so that the confidence interval for the proportion broken is given by where α is the significance and z is the reduced Gaussian distribution.
An example of data collected following this approach is presented in Fig. 2. Obtained from impact of iron ore pellets against a target, it exhibits the typical S shape and asymmetry, evident from the skewness of the distribution to the right.
A more direct and less time-consuming approach is to measure directly the energy required to break individual particles, using either slow compression devices (Sikong et al., 1990;Tavares, 2007;Rozenblat et al., 2011;Ribas et al., 2014;Campos et al., 2021) or instrumented drop weight test devices, such as the impact load cell (King and Bourgeois, 1993;Tavares and King, 1998;Bonfils, 2017;Lois-Morales et al., 2020), analyzing the distribution of data. While requiring more specialized equipment, this approach has the advantage of demanding significantly smaller sample volumes. In this approach, the force- Fig. 2 Proportion broken as a function of specific impact energy from propelling fired iron ore pellets contained in the size range 12.5-9.0 mm against a steel target (data from Cavalcanti et al., 2021). Vertical error bars are given by Eqn. (4) and horizontal error bars were estimated from measurements of impact velocity using a high-speed camera. deformation profile during loading of each particle is recorded and the mass-specific particle fracture energy E is given by numerical integration of the curve up to the deformation or force responsible for fracture where P(E) represents the breakage probability or the cumulative distribution, E is the mass-specific fracture energy of the particle, which corresponds to the maximum stressing energy that it can sustain in a collision and not break, E max is the upper truncation value of the distribution, E 50 and σ 2 are the median and the variance of the distribution, respectively. The upper truncation value is often represented by the ratio E max /E 50 . In the case when this ratio is equal to infinity, E* = E and Eqn. (7) becomes the lognormal distribution.
Besides stressing energy, several other variables also influence breakage probability (Fig. 1). The stiffness of the surface in contact is such a variable, which is analyzed in greater detail in 2.2, whereas the important effect of particle size is analyzed in section 2.3. The effect of impact angle on the breakage probability has been the object of some studies (Dan and Schubert, 1990;Salman et al., 2003;Tavares et al., 2018;Cavalcanti et al., 2021). In a comprehensive campaign of experiments in which iron pellets were propelled against a target Cavalcanti et al., 2021) it was found that breakage probability, and therefore body breakage, may be described considering only the normal component of the impact velocity, at least for impacts at angles that range from 90° (normal) to 30°. The effects of other variables, such as particle shape, stressing rate and mode of stressing (single or double) ( Fig. 1) have been studied in the literature (Tavares, 2007;Bonfils, 2017;Saeidi et al., 2017;Lois-Morales et al., 2020) but the lack of a proper understanding on how they affect the several aspects of breakage for brittle materials in general made it impossible to properly account for them in the model. As such, it is advised that the user fits the model parameters using data that most closely match conditions that dominate in the device that is meant to be simulated.
Adapting Eqns. (7) and (8) to be used to describe discrete particles, such as those embedded in a DEM simulation (Tavares et al., 2020b; is quite straightforward, since it only requires randomly sampling the distribution every time a particle is created.

Effect of stiffness on body breakage probability
Analyzes in section 2.1 assumed that the stiffness of a particle is significantly smaller than that of the surface of the testing device or machine which is in contact with the particle. In several instances, however, that is not the case, so that a correction must be used. Given the local nature of the contacts involving application of stresses to particles and assuming validity of the Hertz contact theory (Tavares and King, 1998;, the proportion e of the stressing energy involved in any given event that is used to deform the particle is given by (Becker et al., 2001;Tavares and Carvalho, 2011; where k p is the particle stiffness (Tavares and King, 1998), k s is the Hertzian stiffness of the surface in contact with the particle. For instance, Eqn. (9) may be used to estimate the proportion of the energy recorded using Eqn. (5) in the test that was actually dissipated in breakage of the particle. Since steel surfaces are typically used in the testing devices, with Young's modulus (Y) of 209 GPa and Poisson's ratio of 0.3, their stiffness [k s = Y/(1 -ν 2 )] is 230 GPa. The apparent particle stiffness may be estimated using a convenient procedure in an earlier work (Tavares and King, 1998), whereas a more accurate, but more laborious, measure may be obtained using a procedure that was recently proposed (Angulo et al., 2020). If the particle stiffness is significantly lower than that of the surface, then e ≅ 1. However, in a collision involving two particles of the same material, then e = 1/2, with the energy being equally shared between the colliding particles.
In order to assess the validity of Eqn. (9), tests were conducted in an impact load cell device on limestone particles using different types of impactors. In this modified impact load cell device a rod rather than a drop ball is used as impactor (Tavares and King, 2004), so that tests were carried out comparing the specific fracture energies measured us-ing steel or hard polymer (Polyether ether ketone or PEEK) impacting rods. The values of Hertzian stiffness of the impactors were estimated as 230 GPa for steel and 4.8 GPa for PEEK, this later calculated on the basis of data from manufacturer's catalogs (Young's modulus of 4.0 GPa with Poisson's ratio of 0.4). The median apparent stiffness for the limestone particles, estimated using the procedure proposed by Tavares and King (1998) was 6.4 GPa. Fig. 4 shows that the energies required to break particles when the softer impactor (PEEK) is used are significantly higher than when steel surfaces are used. These higher values for the polymeric striker are consistent with the larger dissipation of the energy which results from using a striker with lower stiffness (Eqn. 9). Considering now the application of Eqn. (9) considering the stiffness values, Fig. 4 shows the dashed line, which nearly superimposes to data for the steel striker, showing that the actual fracture energies of the limestone particles were independent of the striker used.
It is important to bear in mind that Eqn. (9) only accounts for the elastic part of the deformation . Whenever stressing events occur with significant contribution from plastic deformation, Eqn. (9) provides only a conservative estimate. An example of this was the impact of iron ore pellets against an aluminum alloy plate, in which Eqn. (9) underestimated the energy dissipated on the impact surface, leading to incorrect predictions (Cavalcanti et al., 2021).
This effect of stiffness of the parts in contact can be dealt with readily using Eqn. (9) either in the case of applying the model to continuum distributions, such as in microscale population balance model formulations (Tavares, 2017), and in discrete forms imbedded in DEM (Tavares et al., 2020a;.

Size effect on body breakage probability
Besides stressing energy, the most critical variable influencing breakage probability is particle size. When dealing with brittle materials, it is expected that the stress or energy require for their fracture increases as particle sizes reduce (Schönert, 1979;Yashima et al., 1987;Tavares, 2007). However, such variation depends on material, with Fig. 5 showing the distributions of particle fracture energies for a sample of gneiss.
Several expressions have been proposed in the literature to describe the size effect on the median value of the distribution (Vogel and Peukert, 2004;Rozenblat et al., 2012). In the present model the effect of particle size on E 50 is described using an expression inspired in reliability theory (Tavares and King, 1998;Tavares and Neves, 2008), where E ∞ , d o and φ are model parameters that must be fitted to experimental data and d i is the representative size of particles contained in size class i. k p and k s are the values of Hertzian stiffness of the particle and of the surface of the device used in measuring the breakage characteristics, respectively. As such, Eqn. (10) already incorporates the effect of stiffness of the tools, given by Eqn. (9). Data from a variety of materials are presented in Fig. 6, which also shows the model fit. It is evident that the model is able to represent variations that range from power-law relationships of fracture energies with particle size up to cases in which a constant fracture energy is achieved at coarse sizes. Parameters of Eqn. (10) for these and other materials are presented in Table 1. Values of parameter φ vary from as low as about 0.7 to as high as 2. However, since it describes primarily the power-law increase of the median fracture energies with the reduction in particle size, this parameter is more prone to uncertainty, due to the challenge associated with estimations of fracture energies at fine sizes, as discussed in section 5.1. Values of parameter d o , on the other hand, varied significantly, from values as low as 0.5 mm to values as high as 200 mm. When this value lies in the range of particle sizes tested (Table 1), it may be interpreted as a characteristic dimension of the material microstructure (Tavares and King, 1998;Tavares and Neves, 2008). On the other hand, values of d o above these make it simply a fitting parameter with no particular physical meaning, implying that no clear plateau would have been reached in the measurements of fracture energies at coarse sizes.
Fitting all three parameters in Eqn. (10) typically requires data from at least five widely different particle size classes, ideally spanning at least two size decades. This is not typically a straight-forward task, since it would almost certainly require using multiple impact load cell or compression testing devices, given the different resolutions required.
Besides the effect of particle size on the median particle fracture energy, it is also worth examining its effect on the standard deviation of the distribution (Eqn. 7). For several materials, the variation of the standard deviation of the distribution may be considered small or negligible. An example of this is presented in Fig. 7, which shows data from Fig. 5 plotted as a function of ratio between the specific particle fracture energy of individual particles and the median value E 50 for each size. Since data approximately superimpose, forming a mastercurve, it may be assumed that the standard deviation of the distribution is invariable with particle size in the range studied. Values of the variance σ 2 have been found to vary, exceptionally, from as low as 0.15 Fig. 5 Distribution of specific particle fracture energies for different particle sizes of gneiss (Santa Luzia quarry). Symbols represent experimental data and lines the fit to the lognormal distribution (Eqn. 7).

Fig. 6
Variation of median particle fracture energies of selected materials as a function of representative particle size.
to as high as 1.2, with typical values of about 0.7. For some materials, however, it has been found that the variability of the distribution increases for finer particle sizes (Barrios et al., 2011) and expressions that are analogous to Eqn. (10) have been successfully used to describe the data (Carvalho and Tavares, 2013). The ratio E max /E 50 also shows the same general trend as the variance with size. These variations are typically associated to rocks/ores, which behave as composites at coarser sizes, responding as individual minerals at finer sizes, with the corresponding wide range of fracture energies of the components.

Weakening by repeated stressing events
Breakage models in the literature have been proposed to describe the so-called incremental damage (Morrison et al., Table 1 Values of parameters in the size model for selected materials (Eqn. 10).

Material
Source  2007; Bonfils et al., 2016), which is given by the proportion of fines generated as a function of repeated stressing. In these models, however, weakening by repeated stressing is confounded with the fines associated to surface breakage. Thus, a more rigorous approach is to describe the increase in breakage probability with the repeated impacts separately. Indeed, when the specific stressing energy involved in a collision is smaller than the fracture energy of the particle, the particle will not undergo breakage (Fig. 1), but will sustain internal crack-like damage, which will make it more amenable to break in a future stressing event. This weakening has been described on the basis of a model based on continuum damage mechanics (Kachanov, 1958), through which the specific fracture energy of the particle will be reduced to (Tavares and King, 2002;Tavares, 2009) where E′ is the fracture energy of the particle after the stressing event and D is the damage variable, which is given by (Tavares, 2009) where γ, the only fitting parameter in the model, is the damage accumulation coefficient, which characterizes the amenability of a material to sustain damage prior to catastrophically breaking, and eE k is the specific energy in the collision available to the particle. Eqn. (12) is implicit but can be easily solved iteratively with the initial guess D = 0. The constant γ is typically estimated through tests in which lots of particles are repeatedly impacted at a given specific energy and the proportion broken recorded (Eqn. 3). From the distribution of specific particle fracture energies, given by Eqns. (7) and (8) and the model Eqns. (11) and (12) the optimal value of the constant may be estimated. This is illustrated in Fig. 8, which compares the model fit to experiments. Tavares (2009) carried out a comprehensive study, which demonstrated that the constant γ is typically independent of particle size, with more limited data showing independence of particle shape. Also, it was observed that values of γ varied from about 1.5 for materials that are highly amenable to accumulate damage prior to breakage, up to 8.1 for materials with exceptionally little amenability to weaken by repeated stressing, with average values being in the range of 3 and 4. Finally, it has then been found that the value of γ is not directly influenced by the mechanical strength of the particle, but by its microstructure: it is often lower for materials with complex microstructures, which are more capable of sustaining damage prior to collapsing and higher for more brittle materials (Tavares, 2000).
Eqns. (11) and (12) are useful when the model is applied to individual particles, such as when it is embedded in a DEM simulation (Tavares et al., 2020a;). An alternative formulation of the model, useful when it is applied in a continuum, such as associated to microscale population balance formulations (Tavares and Carvalho, 2009;Carvalho and Tavares, 2013;Oliveira et al., 2020), is given by which should be solved in association to where E in Eqns. (13) and (14) represent the fracture energy of the particle after the impact. The model does not have a minimum threshold, that is, as E k reduces, D vanishes in Eqns. (12) and (14). In practice, however, a minimum value of E k is often selected, for computation efficiency, when using the model, so as to prevent computing damage in collisions with negligibly small magnitudes.

Surface breakage mass loss
In the breakage model the description of surface breakage is certainly in an earlier stage of development than the models described in sections 2.1 to 2.4. Surface breakage has been widely investigated as part of attrition studies and an elegant model has been proposed by Ghadiri and Zhang (2002). An adaptation of the model for the case when the controlled variable is the specific impact energy gives where κ is the attrition parameter, fitted from data, and d i is the representative particle size. This model has been suc- cessfully applied to both natural and recycled aggregates (Cunha et al., 2014;Moreno-Juez et al., 2021) and iron ore pellets (Boechat et al., 2018;. Unfortunately, it has only been validated for particles contained in a relatively limited (22-9 mm) range of sizes. Selected results are presented in Fig. 9 for tests in which particles are dropped one-by-one by gravity and their mass loss recorded after each impact, showing the validity of the model. Evidently, this model does not account for the fact that surface breakage tends to deviate from linearity as the number of impacts increases and neither for the dependency on particle shape, so it may only be used as a first approximation of a rather complex problem. In an attempt to account for the effect of impact angle on surface breakage,  used the total energy dissipated in the collision, given by the sum of the normal and the shear energy loss estimated in DEM in place of E k in Eqn. (15), obtaining good agreement to data from impact of iron ore pellets.

Breakage distribution function
As mentioned in the introductory section, several models have been proposed in the literature to describe the progeny size distribution from single-particle breakage. Several of them benefit from the fact that when cumulative progeny size distributions are plotted as a function of the relative fragment size, that is, the ratio between the sieve size and the representative size of the parent particles, they superimpose, independently of parent particle size. This representative parent particle size is given by the geometric mean sizes of the initial sieve sizes containing the tested lot.
An illustration of this is presented in Fig. 10, which shows data for different initial sizes for syenite. In spite of the scatter, data for the different initial sizes superimpose reasonably well, demonstrating that such normalization of data in respect to the initial parent size is valid. This behavior has been explained on the basis of the fractal nature of breakage of brittle particulate materials (Tavares, 2007). However, for some materials such normalization is only valid above a minimum fragment size. This is illustrated in Fig. 11 for a sample of gneiss. These so-called nonnormalizable fragment size distributions have been identified to be associated to some brittle materials in which discontinuities concentrate in well-defined size scales, Fig. 9 Average percentage of mass loss per impact of selected materials dropped by gravity onto a thick steel plate. Gneiss and granulite: 22.4-19 mm; recycled aggregate: 14-10 mm; iron ore pellet: 12.5-9 mm.

Fig. 11
Progeny size distributions from testing gneiss particles from Santa Luzia quarry contained in different size ranges in drop weight testers (modified from Tavares and Neves, 2008). creating an inflection point in the cumulative fragment size distribution (Tavares and Neves, 2008;Powell et al., 2014). As such, data contained in sizes below this value should be treated separately, as addressed later in the section.
One convenient method to deal with data such as those in Fig. 10 is to select particular markers on each curve. This is the so-called t 10 procedure originally proposed by Narayanan and Whiten (1988). As part of this procedure the percentage of material passing different 1/nth fractions of the initial particle size, called t n s, are estimated from the cumulative size distribution, typically interpolated using cubic splines. Commonly used values of n are 2, 4, 10, 25, 50 and 75. With the aim of representing the progeny size distribution in greater detail, additional markers may also be used, including 1.2 and 1.5 (Carvalho et al., 2015;Ballantyne et al., 2017). In order to illustrate this, breakage data in Fig. 10 are presented in normalized form in Fig. 12, which shows that data for each n value follows well-defined trends as a function of t 10 . The principle of the method is that, if the user has an estimate of the t 10 value of the size distribution, likely from an expression that relates the initial particle size and the specific impact energy, the entire size distribution can be reconstructed (Napier-Munn et al., 1996;Tavares, 2017).
Several functional forms have been proposed over the years to describe cumulative progeny size distributions from drop weight tests, including the cubic splines (Narayanan and Whiten, 1988;Napier-Munn et al., 1996), the upper truncated Rosin-Rammler (King, 2001;Tavares, 2004) and combined exponential functions (Ballantyne et al., 2017). Given its ability to describe data with great flexibility, the incomplete beta function is used as part of the present model (Carvalho et al., 2015), given by x where the cumulative mass of the particles passing a screen (t n ) with relative size x is calculated by the given t 10 value of the distribution. In the model two parameters (α and β) must be fitted to each value of n selected. As such, with the seven n values (1.2, 1.5, 2, 4, 25, 50 and 75) a total of 14 parameters must be fitted from data, which is substantial. Fig. 12 shows that data can be well described using Eqn. (16). The equation has also demonstrated to successfully fit the data for several other materials, as was also evident in previous publications (Barrios et al., 2011;Carvalho and Tavares, 2013;Carvalho et al., 2015). Unfortunately, given the lack of a proper physical description of fragmentation that would allow extrapolation of data, characterizing the response of a given material requires extensive testing over a range of stressing energies and sizes, which is a time-consuming task. Although it may be argued that each material has its well-defined breakage pattern, so that they will present different t 10 versus t n relationships such as the one presented in Fig. 12 (Tavares and Neves, 2008), it is worthwhile to search for a mean distribution that would be able to describe data for a variety of materials with reasonable confidence.
As such, data from drop weight tests of a total of 40 materials have been gathered (Tavares and Chagas, 2021) and plotted in normalized axes as those in Fig. 12, with results being presented in Fig. 13 for selected values of n. A large scatter in the data appears, but it is also evident that data can be reasonably well represented using a single set of curves. The scatter in the graph in a logarithmic scale suggests the structure of the errors and the validity of the logarithmic objective function used in fitting the model parameters (Carvalho et al., 2015).
A more detailed examination of the data (Tavares and Chagas, 2021) allowed to conclude that the different types of materials (rocks, ores, cement clinkers and coals) did not appear particularly segregated in each of the plots in Fig. 13, which suggests the validity of the generalization proposed. The only exceptions were iron ores (Tavares and Chagas, 2021), which were not included in model fitting. Also, deviations were typically limited to the extremes of the size distributions, that is, either at the coarse (t 1.2 and t 1.5 ) or the fine (t 50 or t 75 ) end.
The list of parameters that resulted from fitting the data in Fig. 13 is presented in Table 2, covering the 40 materials (Tavares and Chagas, 2021). These parameters may be used as default in case of the absence of enough fragmentation data that would allow fitting of the various parameters in Eqn. (16) for a material of interest. Fig. 14 is obtained by plotting the various t n values as a function of normalized fragment size for selected t 10 values using Eqn. (16) with parameters given in Table 2.
Although not included in the present work (Cavalcanti Fig. 12 Progeny size distributions presented in the normalized t 10 versus t n format for syenite from Vigné quarry (data from Fig. 10).
et al., 2021), data from projecting particles against a target could also be well described using the incomplete beta function and parameters from Table 2.
The normalizable breakage function can then be obtained from Eqn. (16) using the ratio between the fragment size D i and the parent particle size d j , so that As mentioned earlier in the section, there are instances in which the normalization approach proposed by the t 10 procedure is not valid below a given fragment size D*. This is the case of the so-called non-normalizable progeny size distributions, which appear typically when materials have a well-defined dimension in which discontinuities concentrate (Tavares, 2000). Indeed, non-normalizable breakage functions have occasionally been used with the population balance model equations for materials such as iron ores (Faria et al., 2019;Campos et al., 2019). Table 3 shows examples of materials whose breakage description required the non-normalizable part of the breakage distribution. This behavior can be described in the form of a powerlaw relationship below a minimum normalization size D*, Fig. 13 Relationship between t 10 and selected t n s for a variety of materials, with the exception of iron ores. Solid lines represent the fit to the incomplete beta function (data from Tavares and Chagas, 2021).
which shows that the non-normalizable part of the distribution is attached to the fine end of the normalizable progeny size distribution obtained using the incomplete beta function (Eqn. 16). Values of the exponent η are shown in Table 3 for selected materials. It is worth comparing these to the average slope at fine sizes from the model for the 40 materials (Fig. 14), which varies from 0.53 at low t 10 values to 0.63 (Tavares and Chagas, 2021). Adapting the present model to describe breakage of discrete particles in DEM has been carried out by different methods, depending on the shapes of particles. In the case of polyhedral particles (Tavares et al., 2020a) these have been used as part of the Voronoi tessellation scheme to create the appropriate fragment size distribution and are implemented in the commercial software Rocky DEM. On the other hand, in order to represent breakage of spherical particles in a DEM simulation, families of spherical progeny have been created for a variety of values of t 10 using a procedure described elsewhere (Tavares and Chagas, 2021). This later has been implemented in the commercial software Altair EDEM .
The described normalization procedure only accounts partially for the effect of particle size. The additional effect of particle size associated to the increase in resistance to breakage as particle sizes reduce, as well as the effect of stressing energy, are described using expressions presented in detail in section 3.2. In addition to those, studies have shown that other variables, such as impact geometry (Tavares, 2007;Saeidi et al., 2016), stressing rate (Saeidi et al., 2017) and mode of stressing, if single or double impact (Tavares, 2007), are known to have an influence on breakage distribution. However, the lack of a proper understanding of their effect made it unfeasible to take them into account explicitly in the model.
At stressing energies that are insufficient to cause body breakage, particles will undergo surface breakage. Carvalho and Tavares (2011) proposed to describe it using the Gaudin-Schuhmann equation where d A is the top size of surface breakage products, which together with λ, must be fit from data. The variation of these as a function of parent particle size and stressing energy has not been studied in as much detail in the literature. From tumbling tests d A has been found to vary typically from about 0.20 to 0.45 mm (Tavares and Neves, 2008;Carvalho and Tavares, 2011). This parameter will likely be strongly influenced by parent particle size, but the data are not yet comprehensive enough to characterize such dependence (Cavalcanti et al., 2021).

Effect of particle size and stressing energy on body breakage distribution
In the original t 10 method (Narayanan and Whiten, 1988;Napier-Munn et al., 1996), size-independent relationships between the stressing energy and the t 10 value were proposed. In the present model the relationship accounts simultaneously for the effects of stressing energy and particle size, this later incorporated in the median fracture energy of particles (Tavares, 2009), where A and b′ are model parameters fitted to experimental data, in which A corresponds to the maximum value of t 10 that can be achieved when breaking a material in a single stressing event, whereas E 50b is the median specific fracture energy of particles that underwent breakage. The higher the specific impact energy E k in comparison to the specific fracture energy of the particles, the higher the value of t 10 and the finer the progeny size distribution. By using the median specific fracture energy of the particles that underwent breakage, Eqn. (20) is able to account for the particle size effect in breakage. Data from a variety of materials have been fitted to Eqn. (20) and results are listed in Table 4. Values of parameter A, representing the maximum t 10 value achieved in a single stressing event, were found to vary between about 40 to 75, with typical values around 50 %. This demonstrates that a limiting maximum fineness can be achieved in a single stressing event on individual particles, regardless of the intensity of the stressing energy. The product of the parameters (A*b′), which corresponds to the derivative of Eqn.
(20) at the origin, represents the amenability of the material to fragmentation. Since the toughness of the material is already captured in the value of E 50b , this product actually characterizes how effectively the stressing energy is used to create new fracture surfaces and fragments. Its values range from as low as 0.5 % to as high as 10 % ( Table 4). As such, some grouping in respect to type of material has been observed (Tavares, 2000), in which materials with simple microstructures, such as individual minerals or crystals, present values below about 1 %, whereas most rocks and ores present values in the range from about 1 and 2 %. Values higher than 2 % have been observed for ores and rocks with poor intergranular bonding (Tavares, 2000), coals, as well as some porous materials, such as cement clinkers.
Typical results from fitting the model to data are presented in Fig. 15, which shows the good agreement between the two. In some cases, however, the scatter is larger when particles of several sizes are included in the fit, potentially a consequence of effects that are not accounted for in the model, such as variation in standard deviation of the distribution with particle size, impact velocity and particle shape. The direct connection in the model between the breakage probability and the fragment distribution, offered by the normalizing value of E 50b in Eqn. (20), is considered an important feature, which renders internal consistency in the model.
Adapting the model given by Eqn. (20) to describe breakage of discrete particles results in the modified equation where k is a variable introduced because the mean and the median values of the distribution do not coincide (Tavares and Chagas, 2021). This expression has been used in as-  sociation to discrete implementations of the model for the case of either spherical (Tavares and Chagas, 2020; or polyhedral particles (Tavares et al., 2020a) in DEM.

Model implications 4.1 Primary body breakage function
In several mills and crushers, the breakage distribution function is widely used to describe the average size distribution from each stressing event inside the devices, playing a central role in traditional population balance models (Austin et al., 1984;King, 2001). In practice, rather than a true material property, it represents the result of model fitting, at best estimated on the basis of the size distribution from short grinding times (Austin et al., 1984).
For the case in which the particle is stressed with the minimum energy to produce body breakage, the resulting fragment distribution is called primary single-particle breakage function (Saeidi et al., 2016). In the few studies that attempted to measure directly this primary breakage distribution function of individual particles, Bourgeois (1993) estimated in about 2 % the value of t 10 for quartz particles under double impact (Saeidi et al., 2016), whereas the corresponding value was about 1 % for iron ore pellets subjected to single impact (Cavalcanti et al., 2021). Tavares (2009) then demonstrated, for a sample of limestone, that it did not vary after each repeated impact. The importance of this primary single-particle breakage function has been demonstrated by Saeidi et al. (2016), who showed that it was possible to predict the response of a material in a drop weight tester to higher impact energies on the basis of this primary breakage function, besides Eqn. (10) and empirical functions representing probability of selection of fragments to secondary breakage. This has also been demonstrated to be valid in DEM environment, where the size distribution from higher impact energies could be predicted using the breakage model in the case of polyhedral particles in Rocky DEM (Tavares et al., 2020a) and Altair EDEM .
The present model allows estimating this primary breakage function by considering the case in which eE k ≅ E in Eqn. (21), giving Given the values of the constants A and b′ (Table 4) and since at such low impact energies the constant k in Eqn. (22) may be considered equal to about 0.95 (Tavares and Chagas, 2021), then the resulting t 10 values corresponding to primary breakage range from as low as 0.5 to, exceptionally, higher than 5 %. As such, from comparison with the experimental measurements by Bourgeois (1993) and Cavalcanti et al. (2021), values estimated using Eqn. (22) are likely an underestimation of the actual value. Nevertheless, they may be considered generally valid, given the challenge in estimating directly the primary body breakage function.
With the aim of illustrating the variations encountered for different materials, Fig. 16 compares primary body breakage distribution functions for selected materials estimated using this approach. It shows the variations in shape and position of the primary fragment size distributions, representing differences in fragmentation pattern, with the single-phase materials (minerals) generating less fines than rocks/ores (Tavares, 2000).

Predicting the outcome of a single stressing event
Arguably, the simplest demonstration of the application of a breakage model is to predict the outcome of a single impact on individual particles, such as in drop weight tester (Napier-Munn et al., 1996;Tavares, 2007). The outcome of this single impact event, considering both body and surface breakage, may be described by the mass balance given by and w i are the weight fractions of the material contained in size class i after and before the impact event, respectively. ξ i is the fractional generation of surface breakage products, a ij is the surface breakage distribution function and b ij is the energy-specific body breakage distribution function. The term [1 -P i (eE k )][1 -ξ i (eE k )] on the right-hand side of Eqn. (23) represents the mass fraction of material that did not suffer body breakage during the stressing event and remained in the original size, whereas the term represents the proportion of material from size classes i and coarser that suffered body breakage and generated fragments reported to size class i. The term gives the cumulative proportion of surface breakage products that reports to size class i. P i (eE k ) is the probability that a particle contained in size class i will suffer body breakage when it captures energy eE k from a collision event.
When varying the specific impact energy applied to a batch of particles, breakage can transition from no body breakage at very low impact energies, where all particles undergo only surface breakage, to body breakage of all particles at high impact energies. In between these are impact energies in which a fraction of the particles underwent body breakage, whereas the remainder only suffered surface breakage. This is illustrated in Fig. 17 for particles of Pedra Sul granulite. In the figure, stressing energies above 360 J/kg were able to produce body breakage in the first impact on all particles, whereas those at 9 J/kg produced only surface breakage. For stressing energies between these, only a fraction of the particles underwent body breakage, resulting in size distributions that represent the transition between surface and body breakage. Unfortunately, validation data are available only for the higher impact energies, since at the lower impact energies the small amount of surface breakage products makes collection of accurate data a challenging task.

Breakage testing and model fitting
The breakage model contains a variety of parameters that must be fitted on the basis of experiments. Although some parameters, in particular those related to the fragment size distribution, may be assumed as default, since their variation for different materials may be relatively limited (section 3.1), experiments are required to estimate the several remaining ones. The originally proposed and most direct approach to estimate parameters from the model was based on the use of the impact load cell device (Tavares, 2007). It is an instrumented drop weight tester in which the minimum energy required for primary fracture (Tavares and King, 1998) as well as the net energy absorbed (Tavares, 1999) may be estimated for each individual particle. Besides that, repeated impact tests have then been used to estimate the damage accumulation parameter (section 2.3), whereas the fragment size distributions from different impact energies served as the basis to fitting the parameters in the model described in 3.2. Slow compression tests can also be used to estimate the fracture energy distribution of particles. Data from non-instrumented devices can also be used to estimate the parameters in several of the model equations. For instance, the breakage probability distributions may be estimated from either standard drop weight tests (Baumgardt et al., 1975) or tests in which particles are projected against a target (Dan and Schubert, 1990;Cavalcanti et al., 2021). These later are particularly relevant when modeling with the aim of simulating material degradation during handling.
Although the detailed analysis of the experimental procedures and devices used to fit parameters in the breakage model is beyond the scope of the present work and is well covered in other publications (Tavares and King, 1998;Tavares, 2007;Rozenblat et al., 2012;Ribas et al., 2014;Bonfils, 2017), it is worthwhile to highlight some important issues related to the tests. Indeed, one key challenge when estimating the fracture energy of individual particles using Eqn. (5) is the identification of primary fracture. This primary fracture or primary body breakage point is given by the condition in which the force suddenly drops, given the inability of the particle to continue transmitting stresses. Such condition is relatively straightforward to identify in some cases, such as in Fig. 18a, which is often found from testing regularly-shaped brittle particles. In these cases separation of newly created fracture surfaces is rapid, resulting from propagation of a major crack through the entire particle, facilitating the identification of the fracture point, given the sudden drop in load. In other cases, however, identifying the instant of fracture can be far more challenging and amenable to subjective judgement. A few examples of such cases are also illustrated in Fig. 18.  Fig. 18b shows the case in which the particle suffered chipping at approximately 100 ms, but that primary breakage is still relatively straightforward to identify. This response is often found while testing highly irregularly-shaped particles, in which a small fraction of their volume becomes subjected to concentration of tensile stresses early in the test, causing the appearance of such superficial cracks. Yet in other cases either multiple peaks appear (Fig. 18c) or no clear fracture point may be identified (Fig. 18d), requiring some subjectivity in identifying the fracture load. Challenges such as these appear in association to testing more deformable and/or less brittle materials. This can be associated to cases in which either a single crack or several cracks grow sub-critically, being responsible for not such a clear reduction in load at primary fracture. However, since no relative motion of the fragments occur, due to the limited brittleness of the particle, the fragments continue to be able to sustain load, preventing a drop in load-bearing ability of the particle and, therefore, identification of the instant of primary breakage. Techniques that can be used to assist in interpreting these results will be analyzed in greater detail in a future publication.
A few useful guidelines when conducting experiments with the aim of estimating the distribution of particle fracture energies are: • Particles should be contained in very narrow size ranges, ideally resulting from sieving following a 4 2 (37) series; • Particle selection for the test should be conducted randomly, so as to guarantee that each lot tested has approximately the same average weight; • Sample preparation should be carried out carefully, so as to avoid the use of crushing processes that may be responsible for generating flaky material, such as roll crushing (Gupta, 2016). Ideally either no crushing or the adoption of the size reduction process that is meant to be used in the targeted industrial operation should be preferred. Also, size reduction processes which are known to induce cracks in particles, sub as high-pressure roll grinding (Tavares, 2005), should also be avoided, whenever possible, in this task; • Differences in particle shapes introduce variability in the distribution (Tavares, 2007;Lois-Morales et al., 2020).
One possibility to deal with this is to segregate particles contained within the 4 2 (37) series according to mass and test separately. If that approach is not used, particles should be randomly selected and the size of the tested lot increased accordingly with the variability; • The number of particles tested in each size range must be at least 30, but should increase as the variability of particle composition and shape in the sample also increases, so that 50 or 100 are typically optimal; • In the case of tests in impact load cell devices, it is advisable to select an average specific impact energy for the test (Eqn. 2) that is just high enough to break the toughest particle in the lot. The use of specific impact energies that are much higher than those must be avoided, in particular in the case of testing less brittle materials, since they can lead to greater difficulty in identifying properly the load responsible for primary fracture (Fig. 18d).
While particles in the size range from 100 mm (Fandrich et al., 1998;Tavares, 2007) down to 0.25 mm (Tavares and Cerqueira, 2007) have been tested in impact load cell devices, the greatest challenge lies in measuring fracture energies of particles contained in the lower range and even finer sizes, given the tedious nature of testing single particles and the limitation in the resolution of the devices. Some researchers have demonstrated the potential of measuring the strength and fracture energies of fine particles in special devices (Steier and Schönert, 1972;Yashima et al., 1979;Sikong et al., 1990;Antonyuk et al., 2005;Nguyen et al., 2009;Yan et al., 2009). Ribas et al. (2014) showed the potential of using a commercially-available micro compression tester to measure the fracture energy distributions of single particles in the 45-37 μm size range, whereas Campos et al. (2021) recently applied it to test iron ore particles contained in the range of 45 to 150 μm, with good results. Nevertheless, due to the challenge in measuring breakage characteristics of single particles down to such fine sizes, some indirect estimation approaches have been proposed. Barrios et al. (2011) proposed testing particles as fine as 150 μm contained in monolayer particle beds, combining them to data from single-particle breakage tests at coarser sizes, and estimating the breakage parameters from back-calculation, with reasonable success. More recently, Tavares et al. (2020b) proposed combining measurements of particle fracture energies of particles at relatively coarse sizes to data from impact tests using a rotary breakage tester. This allowed estimating selected material breakage parameters of a copper-gold ore down to 350 μm (Tavares et al., 2020b).
Yet another alternative to estimate breakage parameters in Eqn. (10) is to back-calculate them from batch grinding tests, assuming that predictions using the UFRJ mechanistic mill model (Tavares and Carvalho, 2009;Carvalho and Tavares, 2013) are strictly accurate. This approach has been used to estimate breakage parameters of an iron ore from batch milling data  as well as of limestone to stirred milling data (Oliveira et al., 2020). In these cases, the parameters from Eqns. (7), (12), (16) and (20) were taken from similar ores.
Evidently, errors can appear in estimates from backcalculation, since they rely on the validity of the various assumptions. In the particular case when only batch grinding results are used, errors may be associated to the fact that any limitation in the mill model in describing the phenomenon would inevitably contaminate the estimates of material breakage parameters, partially contradicting the very approach of the advanced models in fully decoupling material from process. In addition to that, this approach also requires that some parameters be taken from similar materials, which also involves a potential source of bias, since each material is intrinsically unique in regard to breakage response.
As discussed in section 3.3, more effort should be dedicated to modeling the primary breakage function from single-particle breakage, given its central role in the application of the model in particular to discrete systems in DEM. One device that could be useful in this task is the rigidly-mounted roller mill, which crushes particles nearly individually between rigidly-mounted rollers with carefully controlled gaps . In spite of challenges in measuring the energy required to break particles contained in the finest sizes studied (45 μm), it provided good estimates of the size distribution of particles under stressing conditions that approach the minimum required for fracture .

Model applications
The breakage model described has been used either in complete form or in parts in simulating different systems of interest using a variety of approaches. In its discrete form it has been embedded in DEM as part of different particle replacement schemes. Initially, a cutdown version of the model, which did not incorporate description of weakening by repeated impacts (section 2.4), was implemented in EDEM (Jiménez-Herrera et al., 2018;Barrios et al., 2020), demonstrating some ability to describe breakage in both confined and unconfined particle beds. A more complete version of the model, incorporating a new stochastic description of particle replacement (Tavares and Chagas, 2021) was implemented in an updated version of EDEM , demonstrating very good ability to predict unconfined breakage of particle beds. Both approaches were based on spherical particle replacement, which has the advantage of computational efficiency. This later implementation also incorporated description of the unresolved fines, which were added upon completion of the simulation to the size distribution given by the resolved spheres. The more complete version of the model has also been implemented in the commercial software package Rocky DEM (Tavares et al., 2020a). Such implementation, based on breakage of polyhedral particles, provided very good predictions of breakage in unconfined particle beds (Tavares et al., 2020a) and cone crushing (André and Tavares, 2020).
The breakage model has also been widely used as part of advanced models of crushers and mills. For instance, Austin (2004a;2004b) used the damage model from section 2.4 in a phenomenological model of high-speed hammer milling. The full version of the breakage model has been used to predict degradation during handling of materials . It has also been used, coupled to DEM, to model crushing in a vertical shaft impact crusher (Cunha et al., 2014), as well as of degradation of recycled aggregate in a high-shear mixer (Moreno-Juez et al., 2021). The modeling approach has also been used to predict degradation of iron ore pellets during handling using slightly customized mathematical expressions Tavares et al., 2018). The widest use of the model has been in association to a microscale population balance formulation, called UFRJ model, to describe grinding in media mills, including ball mills (Tavares and Carvalho, 2009;Carvalho and Tavares, 2013), semi-autogenous mills  and stirred mills (Oliveira et al., 2020;. A version of the ball mill model of continuous milling with some simplifications is available in the state-of-the-art mining and mineral processing plant simulator IES (Integrated Extraction Simulator).

Summary and conclusions
A practical model describing breakage of single particles has been proposed over the last 25 years. It discriminates body and surface breakage, accounting for the body breakage probability through the upper-truncated lognormal distribution. It is able to predict breakage by repeated impacts owing to a model describing weakening based on continuum damage mechanics. Body-breakage progeny size distribution is described using the incomplete beta function, which also considers the exceptional cases in which the material exhibits non-normalizable breakage response. With an empirical expression that describes the stressing energy and size-dependent breakage intensity, it can predict the size distribution of fragments from any stressing event.
This breakage model has been implemented and is now an integral part of particle replacement schemes embedded in commercial DEM software (Rocky DEM and Altair EDEM), through which it may be used to simulate a number of industrial processes of interest. It has also been incorporated in microscale population balance models describing several mill types, namely ball mills, stirred mills, autogenous and semi-autogenous, being even available (ball mill model) in the plant simulator IES.
While being able to account for several important variables that influence particle breakage, challenges are identified in application and fitting of the model parameters for particular applications. Being a single-particle-based model, challenges exist in estimating breakage parameters for fine particles. Establishing the validity of novel testing devices as well as of back-calculation approaches, is an important task. Further, modifications in the model in order to make it applicable to describe breakage of non-brittle materials (Tavares and Almeida, 2020) is another worthwhile task for advancing model applicability. Finally, models describing surface breakage, which are particularly relevant in predicting mechanical degradation during handling, require more attention, since they are still at a comparatively earlier stage of development and validation. e proportion of the impact energy dissipated in the particle (-) E mass-specific fracture energy of the particle (J/kg) E k specific impact energy (J/kg) E max upper truncation value of the distribution (J/kg) E 50 median specific fracture energy (J/kg) E 50b median specific fracture energy of particle that broke in an event (J/kg) E* variable used to incorporate the upper truncation in the lognormal distribution of particle fracture energies (J/kg) E′ fracture energy of the particle after the stressing event (J/kg) E ∞ parameter describing the variation of the particle fracture energy with particle size (J/kg) F force (N) g acceleration due to gravity (m/s 2 ) h net drop height (m) k variable introduced to correct the bias of using the median instead of the mean value of the distribution (-) k p apparent particle stiffness (GPa) k s apparent surface stiffness (GPa) m mass of the drop weight (kg) m p weight of each individual particle (kg) proportion broken in impact of magnitude E k (-) t n proportion passing 1/n th of the parent particle size in the fragments of a stressing event (-) t 10 proportion passing 1/10 th of the parent particle size in the fragments of a stressing event (-) v (52) average velocity at the instant of collision against the target (m/s) w i weight fraction of material contained in size class i (-) x relative fragment size (-) Y Young's modulus (GPa) z reduced Gaussian distribution (-) α significance level of the confidence interval α n parameter of the incomplete beta function for t n (-) β n parameter of the incomplete beta function for t n (-) γ damage accumulation coefficient (-) Δ deformation (m) Δ c deformation at primary fracture (m) η parameter describing the shape of the nonnormalizable part of the breakage distribution (-) κ attrition parameter (kg/mmJ) λ parameter in the non-normalizable body breakage function (-) ν Poisson's ratio (-) ξ i proportion of attrition products (-) σ standard deviation of the lognormal distribution of fracture energies (-) φ parameter describing the variation of the particle fracture energy with particle size (-)