Characterization of Simultaneous Evolution of Size and Composition Distributions Using Generalized Aggregation Population Balance Equation

The application of multi-dimensional population balance equations (PBEs) for the simulation of granulation processes is recommended due to the multi-component system. Irrespective of the application area, numerical scheme selection for solving multi-dimensional PBEs is driven by the accuracy in (size) number density prediction alone. However, mixing the components, i.e., the particles (excipients and API) and the binding liquid, plays a crucial role in predicting the granule compositional distribution during the pharmaceutical granulation. A numerical scheme should, therefore, be able to predict this accurately. Here, we compare the cell average technique (CAT) and finite volume scheme (FVS) in terms of their accuracy and applicability in predicting the mixing state. To quantify the degree of mixing in the system, the sum-square χ2 parameter is studied to observe the deviation in the amount binder from its average. It has been illustrated that the accurate prediction of integral moments computed by the FVS leads to an inaccurate prediction of the χ2 parameter for a bicomponent population balance equation. Moreover, the cell average technique (CAT) predicts the moments with moderate accuracy; however, it computes the mixing of components χ2 parameter with higher precision than the finite volume scheme. The numerical testing is performed for some benchmarking kernels corresponding to which the analytical solutions are available in the literature. It will be also shown that both numerical methods equally well predict the average size of the particles formed in the system; however, the finite volume scheme takes less time to compute these results.


Introduction
Population balance equations (PBEs) describe the behavior of particle properties' changes due to phenomena such as nucleation, growth, aggregation and breakage [1]. Many applications of PBEs can be found in the area of chemical engineering [2], depolymerization [3,4], waste water treatment [5], bubble columns [6], physics [7] and pharmaceutical sciences [8][9][10][11][12], where these mechanisms have a significant impact on the particle properties in the system. In the aggregation process, two or more smaller particles merge at a specific rate to build larger size particles. Aggregation-driven wet granulation processes such as twin-screw granulation [11,13], high shear granulation [14] and fluidized-bed granulation [15] are extensively used in pharmaceutical manufacturing when characterized using PBEs, mostly assuming that only one property of the particles is changing. However, the aggregation mechanism is caused by the mixing-driven presence of liquid binder in the system. Therefore, for such processes, the ability to mechanistically characterize aggregation requires consideration of the degree of mixing between particles and the binder phase. Iveson [16] has shown that the univariate PBEs are inadequate to capture the actual particle behavior in granulation and therefore suggested the application of multivariate (higher dimensional) PBEs. Such PBEs can be used for the simultaneous prediction of several components, the compositional distribution and hence the accurate prediction of the mixing state in the system. Quantification of the degree of mixing in binary component aggregation has been discussed by Matsoukas et al. [17]. The study concluded that the scaling of the variance indicates that the mixing of components is not characterized by a time scale but by a size scale. This emphasizes the need for accurate prediction of particle size distribution evolution, which depends on the selection of a suitable numerical scheme to solve the PBEs and aggregation.
In this study, we focus on the relevance of different numerical schemes for multivariate PBE solution in terms of their capability to capture the degree of mixing during aggregation. In the aggregation mechanism, the total number of particles reduces with time, but the total mass of the system remains constant. A pure bivariate aggregation population balance equation is an integro-partial differential equation which can be written as follows [18].
Here, n(u, t) is the number distribution functions having size u ≥ 0 at time t (s) and u (µm) is a vector given in terms of the summation of properties like the mass or volume of the components. The bold notations are used to denote vector quantities, i.e., u = [u 1 , u 2 , · · · , u p ], where u p expresses the property of the particle in p th direction. Moreover, the birth term in Equation (1) represents the formation of new particles of properties u due to the aggregation of smaller particles of properties u − u and u . Similarly, the death term describes the loss of particles with properties u due to the collision of particles with properties u . The aggregation kernel a(u, u , t) describes the rate of merging of two particles of properties u and u . It can be noted that the aggregation kernel is non-negative and symmetric in nature with respect to size variable, i.e., a(u, u , t) = a(u , u, t). Moreover, for simplicity, we assume that the aggregation kernels chosen for this study are time-independent, i.e., only size dependency is taken into consideration.
Apart from finding the number distribution function n, different properties of the system, namely integral moments corresponding to the number distribution function, are also important for understanding the complete dynamics of the system [19]. For bivariate PBE, the {i 1 , i 2 , · · · , i p } th order moment corresponding to the number distribution function is defined as Here, µ 0,0,··· ,0 (t) denotes the total number of particles in the system, which is also known as zeroth-order moment, whereas µ 0,0,··· ,1,··· ,0 (t) gives the total mass of the pth property. Similarly, other moments can also be obtained.
Among all numerical schemes, the cell average technique (CAT) is well known for the accurate prediction of the number distribution function and their moments. However, the mathematical formulation for CAT is very complex; hence, it is computationally expensive [54,55]. The finite volume schemes (FVS) are well recognized in the literature as being able to accurately and efficiently obtain these results. The studies conducted in the last decade restricted their comparison only to the number distribution functions and their integral moments [18,35,46,47,49,56]. To examine the accuracy of the prediction of the component mixing for a higher-dimensional PBE, the cell average technique [46] and the finite volume scheme [57] are implemented to approximate a multi-dimensional pure aggregation PBE and compared. The verification of the numerical results is also conducted by comparing the average size of particles predicted in the system using the exact solution.
The paper is organized as follows: to start, a brief introduction of the existing CAT as well as the FVS for solving bivariate pure aggregation PBE on non-uniform meshes is provided in Section 2. In Section 3, the qualitative and quantitative numerical results, particularly various order moments and number density functions computed by both numerical schemes, are analyzed. Finally, Section 4 summarizes the conclusions and discussion of this study.

Numerical Methods and System Analysis
In this section, the mathematical formulations of the existing cell average technique [46] and the finite volume scheme [57] for solving a bivariate aggregation PBE (1) on non-uniform grids are outlined ( Figure 1). For developing the expressions of both numerical methods, it is assumed that particles within a grid cell are concentrated on its representative (or mean of the cell). Before giving the description of the numerical methods, it is necessary to fix the computational domain. For the numerical approximations, we define the size variables u in PBE (1) ranges from 0 to ∞; thus, a large sufficient vector such as u max = [u 1 max , u 2 max , . . . , u p max ] T is replaced with ∞ in the second integral of PBE (1). Thus, the original PBE (1) takes the following form: with the modified initial condition cell cell cell The above expression well suited to numerical simulations; however, it does not capture the property of mass conservation. During numerical computations, a sufficiently large u max is chosen to minimize loss of mass from the system. Further, we assume that the whole domain is divided into I = (I 1 , I 2 , I 3 , · · · , I p ), where I r is the number of grids in r direction for r = {1, 2, 3, · · · , p}. Now, for any r, the mesh points and the step size can be defined by

Cell Average Technique (CAT)
First, the mathematical explanation of the CAT developed by Kumar et al. [46] on the non-uniform grids is presented. Let us suppose that N i defines the number of particles in the ith cell, which can be expressed as where du = ∏ p r=1 du r and Now, let us express the number distribution function in terms of dirac delta functions, i.e., Substituting the above expression in the original PBE (1), we obtain the following set of ordinary differential equations: Here, the discrete forms of birth and death terms are expressed as and Here, j = [j 1 , j 2 , . . . , j p ] and k = [k 1 , k 2 , . . . , k p ] are p-dimensional vectors. The corresponding mass or flux of the particles, which takes birth in the ith cell, can be computed as follows: Further, once the aggregation events of the particles in the ith cell are completed, it is necessary to compute the average properties of all birth events in the cell using the following expression: It is assumed that the particle properties are concentrated on the representative of the cell; however, the possibility of the aggregating particles falling on the representative is low. Therefore, if the birth B i takes place in the ith cell and is not represented by the nodes, then the averaging properties, such as b 1 , b 2 , · · · , b p , are distributed to the neighboring nodes in such a fashion that the integral properties, particularly the zeroth and first order moments for our case, remain conserved. The assigned values b 1 , b 2 , · · · , b p can be calculated from the following relations: where u A i 1 , u A i 2 , · · · , u A ip are the coordinate vectors of the vertices of the ith cell. Hence, the final set of discrete equations can be written as It can be noted that the cell average technique is computationally very expensive as it is necessary to determine the average properties of the aggregating particles in each cell and then redistribution of particle properties is carried out to the neighboring representative of the cell. The distribution is carried out in such a way that the zeroth and first order moments of the system are conserved. A detailed description of the cell average technique can be found in Kumar et al. [46] and its schematic flowchart is provided in Figure 2a.
i=1 Compute from eq. (10), b 1 , b 2 and b 3 from eq. (11) and B i from eq. (7) Find vertices number u A , u B and u C of the triangle where falls Compute the birth (B i ) and death (D i ) terms from eq. (14) Add the weights from eqs. (16) and (17) to eq. (14) which are responsible for conservation of integral properties

Finite Volume Scheme (FVS)
Now, we provide the mathematical formulation of the finite volume scheme [57] for solving a generalized aggregation PBE. The finite volume scheme is based on the idea of preserving the total number of particles and conserving the total mass in the system by just adding two correction factors in the formulation. For developing the mathematical expression of the FVS in a similar domain, it is necessary to define the following set of indices: where u i−1/2 and u i+1/2 denote the boundaries of the ith cell, respectively, and u i represents the mean value of the ith cell. The representation of the ≡ i defines all those pairs of cell indices j and k with particle properties u j and u k , such that the addition of their particle properties u j + u k will fall in any cell having representative u i after the aggregation of particles. Similar to the CAT, for the FVS, it is presumed that the point masses are concentrated on the representatives. Therefore, by proceeding in a similar way as in the CAT, the expression for the FVS is obtained, which is given by We know that the finite volume schemes are well known for the conservation of the various properties. The above formulation (14) takes into account the preservation of the zeroth order moment but not the conservation of the total mass (first order moment) of the system. Nevertheless, this can be resolved by introducing only two weights in the above equation, leading to the following expression: where the correction factors ϕ j,k and φ i,j are defined as and where θ(u i ) denotes the sum of the components of the vector u i , i.e., θ(u i ) = θ(u i 1 ) + θ(u i 2 ) + · · · + θ(u i p ) and l i,j is the index of the cell where (u i + u j ) falls. The theoretical proof of the preservation of the total number of particles as well as total mass of the particles in the system along with the CFL conditon can be found in Kaur et al. [57] and its flowchart is depicted in Figure 2b.

Kernel Selection
The efficiency and accuracy of the numerical methods is tested against the analytically tractable kernels corresponding to 2D PBE. In particular, two kernels, namely constant (size-independent) and sum (size-dependent) kernels, are chosen for the comparison.

Size-Independent Kernel
For a size-independent kernel a(u, v, u , v ) = 1, the analytical solution corresponding to the size-independent kernel was formulated by Gelbard and Seinfeld [22] and is listed in Table 1.

Size-Dependent Kernel
Mathematically, the sum kernel can be expressed as a(u, v, u , v ) = u + v + u + v and is heavily dependent on the size of the particles. The analytical solution in this case is formulated by Fernández-Díaz and Gómez-García [21] and is listed in Table 2. Table 2. Analytical solutions of number density functions for size-dependent kernel.

Parameter
Value φ total mass of the particles in the system

Model Initialization and Post-Processing
For the initial condition n(u, v, 0) = 16uv exp(−2u − 2v), the exact results of number density as well as various order moments for constant and sum kernels are provided in the literature by Gelbard and Seinfeld [22] and Fernández-Díaz and Gómez-García [21], respectively. Before comparing the numerical results, it is important to define the degree of aggregation I agg : which describes the decrease in the number of initial primary particles due to the aggregation process. At time t = 0, I agg = 0, and as it approaches large values, I agg → 1, with all primary particles forming one large particle. Furthermore, the weighted relative error is also calculated to test the accuracy of the number distribution quantitatively [49]: The superscripts ana and num represent the exact and numerical solutions, respectively. Here, ∆ 0,0 represents the relative error in the distribution of the number of particles over the whole size domain. The two solutions corresponding to the two numerical schemes may identify the same prediction for the total number of particles, even though the distribution of particle populations may disagree considerably. This is well captured by ∆ 0,0 . Similarly, the ∆ 1,0 expresses the relative error in the distribution of the total mass of the system. All the simulations and computations for both schemes were carried out using "MATLAB" on a i5 CPU with 2.67 GHz and 16 GB RAM.

Average Size Particles
The average size of particles formed in the system along u and v components is determined using the following expression:ū Moreover, the total average size of particles formed in the system is given bȳ

Quantification of Mixing
For identifying the mixing of components in a bivariate PBE, Matsoukas et al. [17] provided a theory to show that the mixing of components is calculated using the χ parameter corresponding to the variance of excess binder. The underlying principle relies on the assumption that if the components are perfectly mixed among all aggregates, the amount of binder in an aggregate of size v would be φv. If the actual amount of binder in the aggregate is u i , the difference m − φv defines the excess binder χ.
For kernels independent of both composition and initial conditions, the χ 2 parameter should be constant at all times for constant and sum kernels. The expression for a χ 2 parameter is given below: where η = 0.5. A detailed description of the theory of mixing of components can be found in Matsoukas et al. [17] and Matsoukas and Marshall Jr [19].

Results and Discussion
The efficiency and accuracy of the numerical schemes for 2D PBE solutions are tested using exact solutions applying analytically tractable kernels. For qualitative and quantitative testing of numerical schemes, CAT and FVS are compared in terms of various order moments and number density functions approximated by both numerical schemes. The results are presented separately for two selected kernels, namely constant (size-independent) and sum (size-dependent) kernels.

Comparison of Moments and Number Density Prediction
For a constant (or size-independent) kernel, the computational domain is taken from u min , v min = 6 × 10 −5 to u max , v max = 21. This domain is partitioned into 20 × 20 non-uniform cells and time ranging from 0 to 60. Numerical simulations achieved a degree of aggregation I agg (t) = 0.97 during this time on a non-uniform grid. Figure 3 shows various order moments predicted by both numerical methods compared to the exact moments. This reveals that the zeroth and first-order (Figure 3a) moments are equally well predicted by both methods and match well with the exact moments. In other words, both numerical approximations ensure the preservation of the total number of particles and the conservation of the total mass of the particles in the system. Additionally, the second-order moments (µ 20 & µ 11 ) computed by both numerical methods are in good agreement with the exact moments (Figure 3b). However, the third-order moment, namely µ 30 , is more accurately computed by the FVS than by the CAT (Figure 3e). Similar behavior is also shown by the other third-order moment (µ 21 ) as FVS shows a more precise result than the CAT, overlapping with the exact result.
In addition to this, using the flat representation concept [58], the number of particles in each cell of the computational domain obtained numerically is compared with the exact results qualitatively. The idea of the flat representation is to sort the cells' indices from 1 to I 1 × I 2 in ascending order, find the particle population in each cell by the relation N k = n k 1 ,k 2 ∆u k 1 ∆u k 2 and plot against its pivot index k. The plots given in Figure 3e,f show that the CAT more accurately predicts the particle population in each cell, whereas the FVS shows more deviation from the exact solution.

Comparison of Average Particle Size and Mixing State Prediction
The average size of particles formed and mixing of components (variance of excess binder) is calculated using Equations (21) and (22). In order to test the accuracy of prediction by the two methods, the numerical results for the system are compared with the exact results in Figure 4. It can be observed that the results for the average size of particles are equally well computed by both numerical methods (see Figure 4a). This is due to the fact that the accuracy of the average size particles depends on the accuracy of the zeroth and first-order moments and both methods predicted these results with equal precision. Moreover, the variance of excess binder plotted in Figure 4b shows that the CAT is in better agreement with the analytical solution than the finite volume scheme. However, the accuracy of these results can be improved to any level by considering a well-refined grid (as shown in Figure 4c,d).  As a further comparison, the quantitative weighted errors (19) existing in the number distribution functions are also found using the numerical methods (see Table 3). The table below reveals that the weighted relative errors in the distribution of the various moments found using the CAT are larger than the FVS. One can also observe that the weighted errors in various order moments can be enhanced by opting for a more refined grid. In terms of CPU time, the FVS took 30.35 s to obtain all numerical results, whereas the CAT took 35.22 s for a size-independent kernel. Therefore, from the above discussion, it can be concluded that the accurate prediction of the various order normalized moment depends enormously on the behavior of the weighted sectional moments errors but not always on the behavior of the number distribution function.

Comparison of Moments and Number Density Prediction
The same numerical methods that were applied using the constant kernel are repeated using the size-dependent kernel (sum kernel) and the results are compared. The computational domain taken from u min , v min = 6 × 10 −5 to u max , v max = 30 is divided into 20 × 20 non-uniform cells with time varying from 0 to 20. The simulations attained a degree of aggregation of 0.90 at the end time. Figure 5 represents the comparison of results obtained using both numerical approximations with the exact results. It can be seen that the zeroth order moments obtained using the FVS and the CAT match well with each other, though both show under-prediction of the exact result ( Figure 5a). However, these numerical results can be improved by choosing a more refined grid. The first-order moment predicted by both numerical methods shows good agreement with the exact result, i.e., the total mass of the particles in the system is conserved. Furthermore, Figure 5b reveals that the second-order moments (µ 20 & µ 11 ) computed by the FVS show less deviation from the exact moments as compared to the CAT. Similarly, it can be seen from Figure 5c,d that the third-order moments, namely µ 21 and µ 30 , predicted by the FVS, overlap with the exact results, whereas the CAT shows more deviation from the exact moments. The deviations in the moments may be due to the numerical diffusion taking place in the CAT while distributing the averages of the particles to the neighboring nodes, but no such distribution of particles is required in the case of FVS.
Similar to the constant kernel, the particle population in each cell obtained numerically is compared with the exact result using a flat representation for the size-dependent kernel (see Figure 5e,f). These figures reveal that the CAT shows more accurate results than the FVS. In addition to this, the quantitative weighted relative errors in the different order moments computed using both numerical methods are shown in Table 4. The results reveal that the weighted errors existing in the various order moments are larger for the CAT, whereas FVS shows fewer errors. It can also be observed that the errors in the higher-order moments increase substantially more than the lower-order moments for the CAT, whereas the FVS shows stable results. The weighted relative errors in the various moments decrease when a more refined grid is taken into consideration. Consequently, from the above discussion, it can be again concluded that the accurate prediction of the various normalized moments certainly depends on the behavior of the sectional weighted relative errors in the moments.

Comparison of Average Particle Size and Mixing State Prediction
The average size particles and variance of the excess binder are compared in Figure 6 for the sum kernel. One can easily observe that the results are similar to the previous case. The average size particles are equally well obtained by both numerical methods, as illustrated in Figure 6a.
However, the excess binder variance is more accurately obtained by the cell average technique, even though the second-order moments are more accurately computed by the finite volume scheme (see Figure 5b). In contrast, the FVS estimates the second-order moments more accurately than the CAT. However, it shows a large deviation from the exact results, as shown in Figure 6b. Using a refined grid, the accuracy of the results can be improved to a large extent, as shown in Figure 6c,d. In a computational sense, the CAT took 20.09 CPU time to approximate the numerical results. However, the FVS required less time to predict these results (approximately 17.91 s).

Conclusions
This study demonstrates that the accurate prediction of second-order moments does not always predict the variance of the excess binder, which is very important for computing the extent of component mixing in aggregation-driven processes such as pharmaceutical granulation. Two numerical methods, namely the cell average technique and finite volume scheme, have been implemented for solving a multi-dimensional aggregation population balance equation to show these results. Many applications, including granulation and crystallization, merely focus on the accuracy of the various order moments and number distribution function. For this purpose, the finite volume scheme can be used for various practical applications such as granulation, crystallization and bubble columns due to its simpler mathematical formulation and higher accuracy than the cell average technique. However, to compute the extent of the mixing of binder in the system, the accurate prediction of excess binder variance (χ 2 parameter) is very important. Even though the cell average technique is less accurate than the finite volume scheme in second-and third-order moment prediction, it is recommended for computing the χ 2 parameter as it requires less computational CPU time. The accuracy of the χ 2 parameter can be improved to the desired level; however, this is at a computation expense as more grid points would need to be considered in the computational domain.