Breakage, the process of fragmenting large particles into smaller ones, is fundamental across multiple industries such as pharmaceuticals, mining, and chemical manufacturing. [1][2][3][4][5] In these sectors, controlling particle size is critical for optimizing product properties like dissolution rates, bioavailability, and material handling. Breakage is often implemented via milling, grinding, or other comminution methods to achieve desired particle size distributions. Efficient modeling of breakage behavior is crucial for predicting the evolution of particle size distribution (PSD), enabling industries to design more effective and scalable processes. Modeling breakage processes has traditionally been approached through empirical and mechanistic methods. [6][7][8] Empirical models often rely on fitting experimental data to predefined functional forms, which may lack the flexibility to generalize across different systems. [8] Population balance equations (PBEs) provide a general framework for modeling particle interactions, including breakage rates and distribution functions. [9] While the overall structure of PBEs is mechanistic in that it enforces mass conservation, the breakage kernels themselves can either be derived from first principles (making the model fully mechanistic) or be empirical/phenomenological (making the model data-driven). The evolution of particle size distribution is described by a population balance equation that includes birth and death rates, which determine the particle sizes. [10][11][12] While PBEs have been widely used to model particle size reduction, solving these equations efficiently remains computationally challenging, particularly for kernels that exhibit strong nonlinearity, multi-modal distributions, or time-dependent behavior. Our approach provides a generalized framework that remains computationally efficient across a broad range of kernel types without the need for iterative numerical solvers.
Extensive research on population balance equations (PBEs) modeling aggregation, breakage, and related particulate phenomena has spurred the development of diverse numerical and semi-analytical techniques, each with distinct advantages and limitations. Various techniques have been used to solve the population balance equations described in the past: method of moments [13], analytical solutions [14], numerical solutions [15], method of discretization using fixed pivot and moving pivot technique [16][17][18], and Monte Carlo [7]. Though there are many studies on the modeling of fragmentation, the numerical methods are often subject to instabilities and loss of accuracy and are computationally expensive. [19,20] Meanwhile, the analytical solutions are constrained to specific kernels, making it difficult to expand them to other fragmentation kernels. Early discretization-based methods for volume-conserving schemes on non-uniform meshes and for 2D aggregation equations achieve conservation of key quantities but often suffer from grid sensitivity and challenges in accurately capturing higher-order moments. [21,22] Semi-analytical approaches, including the homotopy perturbation method and homotopy analysis method, yield elegant series solutions yet are typically confined to short-time domains unless rigorous convergence analyses are performed. [23,24] Recent advancements have further extended these methods: efficient analytical strategies for nonlinear particle aggregation and coupled aggregation-breakage models enhance long-term accuracy, while novel formulations for hyperbolic fragmentation and modified variational iteration techniques improve computational precision. [25][26][27][28] Comparative studies also reveal trade-offs between variational iteration and Adomian decomposition methods, as highlighted in recent reviews of semi-analytical techniques. [29,30] Moreover, Laplace transform-based approximations offer promising accuracy for pure aggregation and breakage problems, while new homotopy-based approaches and multi-dimensional numerical comparisons between fixed pivot and finite volume schemes underscore the ongoing challenge of balancing computational efficiency, physical fidelity, and long-term accuracy in modeling complex particulate systems. [31][32][33] More recently, a fixed pivot discretization was used to formulate the PBEs as a set of coupled ODEs with subsequent eigendecomposition. [34] In contrast, our approach uses a direct graph-theoretic framework and Sylvester's expansion to deconstruct the milling matrix while providing a generalized eigen-decomposition that yields closed-form analytical solutions. These generalized eigenformulations ensure the applicability to a wide range of kernels. Although recursive relationships for numerically generating eigenvalues and eigenvectors have been documented in the literature, an explicit analytical solution has remained elusive. [35,36] In this study, an explicit analytical form is derived and employed to facilitate the identification of breakage kernels, a task that has received limited attention. Notably, while previous works obtained eigenvalues and eigenvectors from experimental data, these efforts did not extend to kernel identification. [10] However, our work, which uses closed-form analytical expressions for eigenvalues and eigenvectors, significantly simplifies the process of kernel identification.
Experimental validation of breakage models has often been performed in tandem with PBE-based simulations. [9] By integrating experimental data from processes such as high-shear wet milling, researchers have been able to tune and validate PBEs to accurately predict PSDs across different milling conditions. [9] They have demonstrated the utility of PBE-based modeling in scaling up high-shear rotor-stator wet milling, using experimental data to optimize parameters such as breakage rate and distribution functions. Such studies underline the importance of combining mechanistic models with experimental validation to ensure that models reflect real-world conditions. However, these approaches often rely on specific kernel forms or simplifying assumptions that limit their applicability across different types of breakage phenomena or the use of toolboxes to solve the PBEs with specific breakage kernels. Present work not only streamlines the predictive process but also enhances reproducibility and applicability across various domains, thereby addressing some of the inherent limitations of current solutions.
In this work, we propose an innovative analytical framework capable of accurately predicting particle breakage for both simple and complex kernels while ensuring numerical stability over prolonged time scales. We begin by outlining a conceptual approach for solving the population balance equation with a focus on fragmentation dynamics. Thereafter, we derive an analytical solution for pure breakage using Sylvester's formula and benchmark our results against the established solution by Ziff and McGrady. [37] Notably, our method is versatile, accommodating both linear and non-linear kernels. The proposed framework is further corroborated by experimental data obtained from wet milling experiments conducted on two high-shear mills, namely the IKA Magic Lab and Quadro Ytron. Overall, the manuscript details the development of the model, the derivation of the analytical solution with generalized eigen formulations, its validation through comparison with literature, and kernel identification for experimental data.
Lab scale wet mills: IKA Magic Lab and Quadro Ytron are used for milling compound A. IKA Magic Lab fitted with three (1 Medium and 2 Fine) heads was used, and the milling experiments were performed with a single rotation rate of 10,000 rpm (~15.6 m/s). For the Quadro Ytron, the measurements were provided in the tip speed of the blade, where 15.7 m/s was used. A 150 ml of the slurry with 30 Volumes concentration was used in IKA, whereas 1500 ml of the slurry was used with the same concentration in the Quadro wet mill. Since the IKA wet mill is relatively smaller, a flow rate of 1000 ml/min was used, whereas in the Quadro wet mill, 10,000 ml/min was used as the flow rate. Samples were collected after each time required, filtered, and stored aside for further characterization. Particle Size distributions of the wet samples were determined by a laser diffraction technique using the Malvern Mastersizers MS3000. This equipment provides a volume-based particle size distribution of equivalent spheres in the range of 0.02 -2000 µm. The dispersant used was n-heptane.
Consider a pure breakage process or a milling operation in which particles of size x are fragmented into smaller particles while simultaneously being generated from the disintegration of larger particles (y > x). The evolution of number density n(x,t) of a particle of size x at time t is provided by the population balance equation as shown in the Eq. ( 1):
where V(y) denotes the average number of fragments produced from a particle of size y, s(y) is breakage frequency, and P(x|y) represents the probability density function for producing a fragment of size x from a particle of size y. The probability density function must satisfy the normalization condition such that
and the conservation of mass requires
The product V(y) P(x|y) dy is known as the breakage distribution function B, which characterizes the size distribution of fragments emerging from a parent particle of size y. Given an initial number density, n o (x), the population balance Eq. ( 1) can be solved using an appropriate discretization method. For discretization, consider a finite discretization of x into m bins of average sizes x 1 , x 2 , …, x m arranged in ascending order. The edges of the bin of average size x i are v i and v i+1, as shown in Fig. 1. Although the primary objective is to develop a generalized analytical solution for the population balance equation, the discretization process is employed solely to facilitate the eigendecomposition and Sylvester's expansion that underpins the analytical framework. Importantly, the analytical solution itself does not depend on the particular discretization scheme used. In this work, a straightforward discretization method was adopted to construct the milling matrix. However, this framework is sufficiently flexible to accommodate alternative strategies, such as Fixed Pivot techniques, that more effectively handle discontinuities by adapting the mesh resolution as needed. [16][17][18] Thus, if a discontinuous kernel is encountered, one can refine the discretization approach without altering the core analytical insights and predictions of our model. Such works that focus purely on the analytical solution exist without employing a rigorous discretization technique. [10] The discretized form of the population balance Eq. ( 1) becomes
where,
The normalization condition and mass conservation require that
The discretized Eq. ( 4) can be expressed in matrix form as
Eq. ( 7) can be written as:
Where M is the milling matrix. The solution to the Eq. ( 8) for the initial number distribution, N ⇀ 0 is given by
Direct evaluation of the matrix exponential at each time point in Eq. ( 9) is computationally expensive, with the cost generally scaling as O(n 3 ) for a square matrix of size n. [38] In contrast, Sylvester's expansion expresses the exponential as a sum over eigen values and eigen vectors as shown in Eq. (10), which reduces the computational cost significantly. Using Sylvester's expansion, the number distribution function represented in the Eq. ( 9) can be decomposed as
where λ i is the eigenvalue, X ⇀ i is the right eigenvector and Y ⇀ i is the left eigenvector of the milling matrix M. Since smaller particles do not generate larger ones, the milling matrix M resulted in the upper triangular matrix. Hence, the eigenvalues are the diagonal elements of the milling matrix M.
However, the direct application of our method to discontinuous breakage processes, for example, when a Dirac Delta daughter distribution is involved, the milling matrix may become degenerate due to repeated eigenvalues. In such cases, the solution can be formulated using the Jordan canonical form. Specifically, if the milling matrix is expressed as M = PJP -1 , the analytical solution becomes exp(Mt) = Pe Jt P -1 , where e Jt is evaluated based on the structure of Jordan blocks. This extension seamlessly incorporates discontinuous breakage into our generalized analytical framework. The analytical expression for the components of the right and left eigenvectors of the milling matrix are derived by the method of successive approximation, (detailed in the supplementary information) and can be validated numerically using singular value decomposition. The elements of the normalized right eigenvector are defined as follows:
The modal matrix of normalized right eigenvectors is given as:
where T j,i k represents the relative rate of k th order transition from bin i to j. Here, the k th order transition corresponds to k discrete jumps from bin i to j. The total number of possible k th order transitions from i to j is given by the combinations L= i-j-1 C k-1 . The relative rate T j,i k obtained with the method of successive approximation is: Fig. 1. One dimensional domain discretization with representative size x i and size group boundary v i .
where W k l is an l th set of pairs (u,v) involved in k th order transition, for example, in the case where a particle of size 6 produces a fragment of size 2, k may range from 1 to 4. For 1st order transition, there should be a single jump from a particle of size 6 to a particle of size 2. Hence, there is only one possibility (L = 6-2-1C1-1 = 1), which is shown in Fig. 2a for k = 1. For 2nd order transition, there must be two steps from a particle of size 6 to reach a particle of size 2. This is possible in three (L = 6-2-1 C 2-1 = 3) ways, i.e., 6→5→2, 6→4→2 and 6→3→2, which are represented in Fig. 2b for k = 2. Similarly, for the 3rd order transition, there are three possible steps, as shown in Fig. 2c (for k = 3), and for the 4th order transition, there is only one possible step, as shown in Fig. 2d (for k = 4). The maximum order of transition is the difference between bigger and smaller particle sizes ( in this case, 6-2=4).
To evaluate T j,i k , let us consider the same example. For k = 1, we get L= 1,
Similarly for k = 2, we get L= 3,
For k = 3, we get L= 3,
We consider microscopic reversibility that occurs in the system, and eventually the system attains equilibrium. This means that there is a reversible process where particles will combine at the micro level, eventually forming bigger particles. This comes into account by considering the left eigenvectors of the milling matrix. The components of normalized left eigenvector for the milling matrix and the modal matrix of normalized left eigenvectors are defined as follows:
where τ j,i k denotes the relative rate of k th order transition from bin i to j. The total number of possible k th order transitions from i to j is given by
The relative rate τ j,i k obtained from the method of successive approximation as
We confirmed that our approach for evaluating the model matrices produces results identical to those obtained by MATLAB's eig function for the normalized right and left eigenvectors. Subsequently, the general solution of Eq. ( 10) can be compressed into matrix form, such as
The solution in the Eq. ( 22) can be expressed as Fig. 2. Graphical representation of transition of a particle of size 6 to a particle of size 2 for a first-order transition (a), second-order transition (b), third-order transition (c), and fourth-order transition (d).
where
To assess the validity of our analytical framework, a specific case of binary breakage problem was considered, and the findings were compared with an existing analytic solution of Ziff and McGrady [37], which applies to a specific kernel. Ziff and McGrady [37] reported that for pure breakage with kernels s i = s(x i ) = x 2 i and B i,j = 2 Δx xj the analytical solution (refer case F(x,y) = x + y in Ziff and McGrady [37]) is
A discrete counterpart of this solution can be written as
In this study, we utilize a widely employed power-law formulation for the breakage rate and the distribution kernel. [2] This choice is driven by the need to capture the underlying physics of the breakage process while providing insights into system behavior under varied conditions. Notably, for certain breakage parameters, the chosen breakage rate function and kernel align with kernels used by Ziff and McGrady [37], allowing for direct comparison and validation of our generalized approach. Specifically, the breakage rate and the associated probability distribution function are defined as s(x) = x a (26)
where a and q are breakage parameters and Γ denotes gamma distribution. Note that P(x|y) is of the form P(x|y) = 1 y f
)
, where f (z) is a symmetric beta distribution function with parameter q. Here q governs the shape of the distribution, indicating whether breakage predominantly produces smaller fragments (q < 0), larger fragments (q > 0), or uniformly distributed fragments (q = 0). For a= 2 and q= 0, Eqs. ( 26) and ( 27) reduce to the specific kernels utilized by Ziff and McGrady [37]. Under these conditions, our generalized analytical solution was directly compared to the Ziff and McGrady model to validate its accuracy and consistency. All simulations were performed on a laptop featuring an AMD Ryzen 7 5800 H processor (3.2 GHz, 16 GB RAM) using MATLAB (version-R2023a). Our method exhibits remarkable computational efficiency, solving the PBE in just 5 seconds of CPU time per simulation shown in Fig. 3.
To evaluate the analytical solution, the initial particle size distribution was assumed to follow a Gaussian distribution with size ranging from 0 to 1 and a standard deviation equal to one-tenth of the mean. However, the proposed model is applicable to any particle size distribution, as demonstrated in the experimental validation, where the size distribution extends up to 2000 µm. The results demonstrated excellent agreement between the generalized analytical solution and the analytical solution in the literature [37], as shown in Fig. 3. Furthermore, the left and right eigenvector matrices predicted by our model matched the normalized eigenvectors matrices obtained using MATLAB's built-in function, eig, confirming the reliability of our eigenvector calculations. Fig. 3(a-e) illustrates the number density function at different time points, comparing the developed analytical solution with Ziff's solution. The two solutions are virtually indistinguishable, validating the accuracy of our generalized analytical approach. This close agreement emphasizes the robustness of the model, showing that it can replicate established analytical results while maintaining flexibility for more complex scenarios.
Further validation was done with an initial particle size distribution where the number density exponentially decays with size (N 0 = e -x ).
The number density at a time t = 4, obtained from our approach is compared with the Ziff's solution and a Padé approximant of the series solution derived from projected differential transform method with Laplace transformation (LPDTM-PA) from the literature. [28] Both techniques yield almost identical number density versus particle size plots (Fig. 3f), demonstrating that both can reach Ziff's solutions under the test conditions. However, our approach, which discretizes the particle size domain and employs an eigen-decomposition of the milling matrix, offers a simpler implementation and greater generalizability, allowing it to be easily adapted to a broader range of breakage kernels and process conditions.
In Fig. 3 and Fig. 4(a), the initial particle size distribution at t = 0 follows a Gaussian distribution. As the time progresses, the number density function steadily shifts toward smaller particle sizes, indicating that the system continuously evolves due to ongoing breakage rather than reaching a steady state. This trend is evident as the number of smaller particles increases, ultimately resulting in the highest density at the smallest particle size (x = 0.002), which can be attributed to the kernel's independence from the size of the daughter fragments when q = 0. Consequently, a decrease in the Sauter mean diameter (ratio of third to the second moment), which is the average measure of particle size, can be seen in Fig. 4(b) which decreases with time. Additionally, analysis of the integral moments confirms that while the zeroth moment (representing the total number of particles) increases with time, the first moment (indicative of the total mass) remains nearly conserved, reinforcing the physical consistency of the model. Equation for the r th moment of the number density function n(x, t) is defined as:
Although our analytical solution model is discretized, the total number (using the first moment) is comparable to Ziff's analytical solution, as shown in Fig. 4 (c). The negligible variation in the first moment demonstrates the generality and robustness of our approach. Notably, for a discretization comprising 100 points, the maximum observed error relative to Ziff's solution was approximately 3 % at extended time scales, with additional error percentages illustrated in Fig. 4(d). Although the discretization error may depend on the chosen mesh resolution, by increasing the number of bins or employing adaptive discretization schemes, the error can be minimized without compromising the internal consistency of the model. While a large number of grid points was used for error analysis and the final solution, linear breakage PBEs may be accurately resolved with as few as 20 -30 grid points. However, the analytical framework is flexible enough to incorporate alternative discretization strategies for enhanced resolution when necessary. These considerations, while important, do not alter the core analytical insights presented herein and will be further explored in future work. Furthermore, a comparison of both the zeroth and first moments for an even longer time horizon (t = 100) is provided in Figure S1 of the supplementary information.
After binary breakage was validated and number conservation ensured, the influence of the parameters q and a on the particle size distribution (PSD) over time was examined. The relationship between these breakage parameters and the resulting PSD for any given scenario was explored using the analytical solution. As illustrated in Fig. 5a, an increase in a, which defines the size dependence of the breakage rate (s(x) = x a ), was found to cause the number density function to become narrower, more symmetric, and essentially time-invariant because of the sharper breakage mechanism. This indicates that, with higher a larger particles are more likely to break compared to smaller ones. Physically, this reflects a system dominated by size-selective breakage, where higher a leads to rapid stabilization of PSD due to the accumulation of smaller particles. Lower a, on the other hand, results in a slower, less selective breakage process, allowing larger particles to persist longer. On the other hand, parameter q, governs the shape of the daughter particle size distribution, influencing whether breakage produces smaller or larger fragments preferentially. When q > 0, the kernel skews towards an increased number of fragments of very small particles and a narrowing of the secondary peak in the PSD with an increase in q, as shown in Fig. 5b. This behavior signifies a selective process where specific fragment sizes dominate, reducing the diversity of intermediate sizes. Together, a and q provide insights into the mechanistic aspects of breakage processes, with a controlling the size-dependence of breakage rates and q shaping the relative size distribution of the resulting fragments. Tuning these parameters allows engineers to optimize PSD for desired applications, such as rapid size reduction or the generation of controlled fragment sizes for specific material properties. Overall, it can be understood from Fig. 5 that our analytical solution can be extensively used to study the breakage model with various kernels and parameters robustly. For cases where the rate kernel is size-independent, implying that all particles (including smaller ones) break at the same rate, the analytical solution can be expressed in terms of the Jordan canonical form. In scenarios where the kernel is not size-independent and the milling matrix is readily diagonalizable, our approach remains accurate.
The analytical model was further validated against experimental data obtained from wet milling processes using two distinct high-shear mills: the IKA Magic Lab and Quadro Ytron. The PSD quantiles (D90, D50, and D10), representing the 90th, 50th, and 10th percentiles of cumulative particle size distribution, were tracked over time and compared to theoretical sizes, as shown in Figs. 6a and 6b. The close agreement between the experimental data and the model demonstrates the robustness of the analytical solution in capturing the physical processes underlying particle breakage.
To enable this comparison, the breakage rate and breakage distribution functions described by Austin et al. were adapted and simplified, as outlined in Eqs. 29-31.[39] These functions were then incorporated into the analytical model, allowing for a more detailed exploration of the breakage dynamics and alignment with experimental observations.
) α (29)
B i,j = f i,j -f i+1,j (31) where s i is the specific breakage rate of particles of size i, s N is the specific breakage rate of the coarsest particle, x i is the diameter of particles in bin i, f i,j is the cumulative mass-based breakage distribution function, B i,j is the mass-based breakage distribution function that specifies the masses of child particles of size i formed when particles of size j are broken. α, ϵ are the material and mill-dependent parameters whereas S N depends on mill conditions and type. Our model seamlessly integrates such complex kernels into the analytical framework by parameterizing them through the milling ma-Fig. 5. Effect of breakage parameters a and q for t = 1.
trix and eigenvector analysis. This flexibility allows our model to adapt to a wide variety of breakage kernels. The kernel identification with Austin et al.'s model demonstrates its effectiveness in capturing realistic milling dynamics where breakage is influenced by both material properties and operating conditions. For the IKA Magic Lab, optimized parameter values obtained are s N = 8272.7, α= 4.076, and ϵ= 1.25 and for Quadro, values are s N = 1124.2, α= 4.11, and ϵ= 1.31. The breakage frequency values for both the wet mills show that IKA has a more aggressive breakage mechanism compared to the Quadro, indicating a stronger size dependence on the coarsest particle. It can be hypothesized that the geometrical differences between the two mills are a major contributing factor to such differences in the breakage mechanisms. For example, the number of generators for the IKA magic lab is three, whereas for the Quadro wet mill, it's just one. Moreover, this can be observed with the experimental data as well, where IKA has smaller particle sizes at all times during milling compared to Quadro.
In this work, we developed a generalized analytical solution for population balance equation (PBE) describing pure breakage processes. By discretizing the PBE and utilizing Sylvester's expansion, we successfully constructed a robust framework capable of handling different breakage kernels. Validation of the analytical solution against Ziff and McGrady's established model along with the kernel identification using experimental data from two different types of wet mills (IKA Magic Lab and Quadro Ytron) and Austin et al.'s breakage kernels demonstrated the model's accuracy, applicability and reliability.
In conclusion, the generalized nature of this solution, combined with simplicity, ease of applicability, and minimal computational requirement, provides a powerful tool for researchers and process development engineers, offering predictive capabilities for particle size evolution in systems such as milling and grinding. Its ability to adapt to complex kernels ensures broad applicability across industries such as pharmaceuticals, mining, and materials processing. Furthermore, integrating real-time experimental feedback with the model could enable adaptive process control, further enhancing its utility for optimizing particle processing technologies. Future work could explore extending this approach to more complex systems involving coupled phenomena, such as aggregation, dissolution, and/or reactive processes. Overall, this work represents a significant advancement in analytical solutions for PBE, paving the way for more precise and scalable approaches to modeling particle breakage systems. original draft, review & editing, Visualization, Validation, Methodology, Investigation, Formal analysis, Data curation.