An accurate and efficient discrete formulation of aggregation population balance equation

An efficient and accurate discretization method based on a finite volume approach is presented for solving aggregation population balance equation. The principle of the method lies in the introduction of an extra feature that is beyond the essential requirement of mass conservation. The extra feature controls more precisely the behaviour of a chosen integral property of the particle size distribution that does not remain constant like mass, but changes with time. The new method is compared to the finite volume scheme recently proposed by Forestier and Mancini (SIAM J. Sci. Comput., 34, B840 - B860). It retains all the advantages of this scheme, such as simplicity, generality to apply on uniform or nonuniform meshes and computational efficiency, and improves the prediction of the complete particle size distribution as well as of its moments. The numerical results of particle size distribution using the previous finite volume method are consistently overpredicting, which is reflected in the form of the diverging behaviour of second or higher moments for large extent of aggregation. However, the new method controls the growth of higher moments very well and predicts the zeroth moment with high accuracy. Consequently, the new method becomes a powerful tool for the computation of gelling problems. The scheme is validated and compared with the existing finite volume method against several aggregation problems for suitably selected aggregation kernels, including analytically tractable and physically relevant kernels.


(Communicated by Lorenzo Pareschi)
Abstract. An efficient and accurate discretization method based on a finite volume approach is presented for solving aggregation population balance equation. The principle of the method lies in the introduction of an extra feature that is beyond the essential requirement of mass conservation. The extra feature controls more precisely the behaviour of a chosen integral property of the particle size distribution that does not remain constant like mass, but changes with time. The new method is compared to the finite volume scheme recently proposed by Forestier and Mancini (SIAM J. Sci. Comput., 34, B840 -B860). It retains all the advantages of this scheme, such as simplicity, generality to apply on uniform or nonuniform meshes and computational efficiency, and improves the prediction of the complete particle size distribution as well as of its moments. The numerical results of particle size distribution using the previous finite volume method are consistently overpredicting, which is reflected in the form of the diverging behaviour of second or higher moments for large extent of aggregation. However, the new method controls the growth of higher moments very well and predicts the zeroth moment with high accuracy. Consequently, the new method becomes a powerful tool for the computation of gelling problems. The scheme is validated and compared with the existing finite volume method against several aggregation problems for suitably selected aggregation kernels, including analytically tractable and physically relevant kernels.
1. Introduction. Population balances appear in many branches of engineering and science to model the dynamics of particulate processes such as crystallization, aerosol formation, precipitation, polymerization, agglomeration etc. The one dimensional continuous population balance equation (PBE) for pure aggregation processes in a well mixed batch system is given by [10] as along with the given initial condition f (0, x) = f in (x), x ∈ (0, ∞). Here t ≥ 0 stands for time and f (t, x) represents the number density function. Without loss of generality, the above population balance equation is written in a nondimensional form and all variables are taken as dimensionless, [27]. The first right hand side term represents the birth (gain) of the particles of size x as a result of the aggregation of particles of sizes (x − x ) and x . In this paper, we shall refer to size as the particle volume or mass. The second term, known as the death (loss) term, describes the merging of particles of size x with other particles. The dynamics of the process are modelled by the coagulation kernel β. The kernel is non-negative and satisfies the symmetry condition β(t, x, x ) = β(t, x , x). Besides predicting the particle size distribution, estimation of moments of the distribution is also of great interest in several applications. Moments are integral properties of the distribution and some of them represent physical quantities that are easy to measure experimentally. Moreover, assessment of different numerical methods is usually based on their accuracy and efficiency to produce the particle size distribution as well as its moments. To this end, it is therefore important to formally define the jth moment of the particle size distribution as The zeroth and first moments are proportional to the total number and total mass of the particles in the system at time t, respectively.
For a class of aggregation kernels, a phase transition occurs at the so called gelation point where mass is lost from particles of finite size and appears in particles of infinite size, [6]. Such kernels are called gelling kernels and this phenomenon is known as gelation. From a macroscopic point of view the gelation effect is represented by a loss of mass in the solution. Mathematically, the gelation point can be described in terms of the first moment, which is proportional to mass, as [5] The presence of gelation has been investigated by [30] for both the continuous stirred tank (CST) and the batch mode of operation. According to their study, the sizeindependent kernel is a non-gelling kernel for both modes of operation, while the sum kernel is a gelling kernel for CST operation but a non-gelling kernel for batch operation. Further, the product kernel is a gelling kernel for both batch and CST mode of operation. The exact values of gelling time t gel for various cases can be found in [30].
It should be noted that the size variable x in the population balance equation (1) ranges from 0 to ∞. In order to apply a numerical method to solve the population balance equation, it is therefore necessary to fix the computational domain by changing ∞ to a finite number. In this work, we denote the maximum particle size by X and therefore ∞ appearing in (1) is replaced by X. Thus, the truncated population balance equation is given as along with the truncated initial condition f (0, x) = f in (x), x ∈ (0, X]. Without loss of generality, we continue with the same notation for the number density f of the truncated problem. The equation (2) is now suitable for numerical computation, however, it fails to hold the property of mass conservation.
Several numerical schemes based on finite volume methods [3,31,7,4,8,13], finite element methods [23,21,28], the method of moments [2,20,22,1], sectional approaches [9,11,19,32,17,18,15,12,26,16] have been developed for solving one dimensional aggregation population balance equation. The scheme given in [15] and other similar schemes referenced therein distribute the new-born particles in each cell to neighboring nodes so that the desired integral properties are preserved. In particular, the cell average technique computes the average position of all new-born particles in each cell and distributes them to neighboring nodes in one step, while the fixed pivot technique does the distribution for particles individually without computing the average position. As far as a typical finite volume approach for solving aggregation population balance equation is concerned, work in the literature based on finite volume methods is restricted to [7] and [8]. In our present work, the formulation of the proposed scheme preserves the desired integral properties by introducing some correction factors. We shall observe that the proposed scheme leads to a rather simple mathematical structure and consequently computationally less expensive with comparable accuracy.
This paper is based on the recently developed mass conserving finite volume scheme on nonuniform meshes by [8]. Before proposing the new discretization, we first briefly explain the finite volume scheme of [8]. For any numerical method, the first step is to divide the domain into small cells (intervals), Λ i := (x i−1/2 , x i+1/2 ], i = 1, . . . , I, with x −1/2 = 0 and x I+1/2 = X. The representative or the grid point of a cell, denoted by x i , is taken as the mean size of the cell, i.e., x i = (x i−1/2 + x i+1/2 )/2. Moreover, for a fully discrete formulation, the time domain needs to be discretized. Denoting the nth time step by ∆t n , the discrete time follows t n = n∆t n , n = 1, 2, . . .. Let f i (t) denote the average number density of the ith cell, i.e., The numerical values approximating f i (t) at time t n is denoted by f n i . At this juncture, we can give the definition of a mass conserving scheme. Definition 1.1. We call a numerical scheme mass conserving or first moment conserving if the following holds x i ∆x i f n i , for each n.
Next step for the derivation of the finite volume scheme is to integrate the equation (2) over the cell i as d dt Using the standard approach of a finite volume scheme, [8] arrive at the following fully discrete formulation where β n j,k = β(t n , x j , x k ), The above discretization, however, does not hold the mass conservation property. In order to ensure mass conservation, authors have modified the kernel of the birth term such that the total mass becomes conserved. The modified kernel of the birth term is taken as Thus, the resulting formulation takes the following form To this end, it is important to point out that the discrete formulation (5) shows the conservation of the first moment, however, no care is taken during the process of modification for other moments. Therefore, the question arises: How much does this modification affect the accuracy of, for example, zeroth, second or higher moments. It is shown later that this modification leads to the overprediction of the second as well as other higher moments and the underprediction of the zeroth moment. In this work we show that, if due care is taken during the modification process, a numerical scheme that provides better solution of several moments in addition to the mass conservation can be formulated.
The objective of this work is to present a new discrete formulation that predicts different moments of the particle size distribution more accurately with comparable accuracy of the complete particle size distribution than the existing formulation (5) without adding any further complexity to the formulation. Moreover, the method is proposed with the target of reducing the computational cost, in other words, a method that can produce marginally accurate results by using only a few number of cells. Furthermore, it is well known that the gelling problems are more difficult to solve and most of the numerical methods fail to conserve mass before reaching the gelling time. Therefore, one of the aims is to propose an efficient method that can be employed to identify the gelling time numerically.
2. New method. To get a discrete formulation of the integrated continuous population balance equation (3), there are two common ways to treat the underlying integrals. The one which is used in the work of [8] relies on approximating the integrands within the two representatives by interpolating polynomials. The other approach assumes the number density function as the point masses concentrated on representatives, with the measure of each representative being equal to the total number of particles in the interval containing the respective representative. The new method is based on the latter approach.
We start the formulation of the new scheme by defining the following sets of indices as As mentioned before, we deal integrals with the assumption of point masses concentrated on representatives, i.e., Substituting the above approximation of f (t, x) in equation (3) and pursuing the calculations from the work of [14], we end up with the following discrete equation Discretizing the time derivative using the first order Euler method, we get Similar to the discrete equation (4), the above discrete formulation (7) also suffers with the mass conservation, if applied on any meshes other than the uniform meshes of the type x i = ih, for some constant h. Before we propose some modifications to the above formulation, we first introduce here a rather general concept of preservation of moments. Note that, a conserved moment (first moment in our case) is a moment that remains constant in time. This is a property that is very easy to verify numerically, because the time derivative of a conserved moment is zero. The new concept of preservation allows us to compare schemes with respect to how well they reproduce the correct behavior of moments that are not conserved but change in time. We present the concept of preservation only for the zeroth moment, however, a similar study can be performed with respect to any other moment.
Integrating the population balance equation (1) from 0 to ∞ and further simplifying by suitable substitutions as well as change of order of integration, we get the following differential equation for the temporal change of the zeroth moment Using the same assumption of number density as the point masses concentrated on representatives, i.e., substituting f from (6) into the above equation, we obtain After discretizing the above equation in time variable, we find

JITENDRA KUMAR, GURMEET KAUR AND EVANGELOS TSOTSAS
The truncated form of above equation can be written as Defining the numerical zeroth moment asμ n 0 := To this end, we give the following definition: We call a numerical scheme number preserving or zeroth moment preserving if the numerical solution follows the equation (10) at each time step t n on the kernel.
Let us now try to understand the above preservation property with the help of a simple example of constant aggregation kernel, β(t, x, x ) = β 0 . The equation (8), in this case, reduces to Discretizing the time variable using the first order Euler method, we obtain where µ n 0 approximates the zeroth moment of the truncated problem with µ 0 0 =μ 0 0 . Notice that µ n 0 is obtained without discretizing the space (size) variable x, and, in general, it is different than the numerical first momentμ n 0 . The only error in µ n 0 accumulates due to the time discretization, whereasμ n 0 suffers from both errors in space and time. Also, note that the error due to the time discretization can easily be controlled either by taking small time step size or integrating the model equation over time by a higher order solver.
For β n i,j = β 0 , the equation (11) reduces tô Comparing equations (13) and (14), we find thatμ n 0 = µ n 0 . This indicates that a number preserving numerical scheme must produce very accurate results for the zeroth moment in the case of a size-independent kernel. The condition (10) respects the underlying assumptions on number density function (6) and signifies the consistency of the discrete formulation with respect to the zeroth moment.
To conclude, the preservation property (10) controls the growth of the zeroth moment more accurately. Note that the condition (10) and the discrete formulation (7) are derived under the same assumption on particle size distribution. Indeed, one can easily show that the discrete formulation (7) satisfies the condition (10). Hence, any numerical method which is based on the assumption that particles within the cell are concentrated on its representative must be consistent with the condition (10). It should be emphasized here that any modification to the original discrete equation (7) should not violate the condition (10), otherwise the growth of zeroth moment may differ significantly from its true value. We will observe later in numerical comparisons for various kernels that a zeroth moment preservative scheme predicts the right behavior of the zeroth moment. The discrete formulation (7) is number preserving, however, as pointed out earlier, mass conservation is an essential characteristic of a numerical scheme that has to be fulfilled. Thus, the aim is to propose a method that has both the properties of mass conservation and number preservation. This can be achieved by introducing some correction factors in the formulation (7).
In what follows, we propose the following new finite volume scheme (NFVS) as where w b j,k and w d i,j are the weights responsible for the number preservation and mass conservation. These weights are defined as Here l ij is the index of the cell where (x i + x j ) falls. Notice that the index l ij is symmetric with respect to its subindices, i.e., l ij = l ji . Next, we verify the earlier stated properties for the proposed method. Proof. The above result can be proved as follows. Multiplying the discrete formulation (15) by x i ∆x i and summing over all i, we obtain where T is given as In order to prove the mass conservation we need to show that T = 0. With the definition of l ij and I i , the representative x i in the first term may be replaced by x l jk , and, therefore, the above expression can be rewritten as By the definition of I i , (j,k)∈I i represents the combination of the cells j and k having properties x j and x k such that their sum x j + x k fall in some cell i. Further, it can be noted that the sum x j + x k may fall outside the computational domain.
Then the corresponding indices belong to the set I * , i.e., (j, k) ∈ I * for which w b j,k = 0 by the definition. So, So, the equation (16) modified to Further, I j=1 represents all possible combinations of the cells i and j. Here two cases will arise: either the sum of properties of the cells i and j, i.e., x i + x j will fall either inside or outside of the computational domain. So, I j=1 = (j,k)∈I i + (j,k)∈I * , and therefore, equation (17) becomes Note that, due to the symmetry of β, and the relation l ij = l ji , we have .
Substituting the weights w b j,k and w d i,j into the preceding expression of T and using the above identity, we find that T = 0. This proves the mass conservation property of the new method.
Proposition 2. The proposed finite volume scheme (15) is zeroth moment preserving.
Proof. Multiplying the discrete formulation by ∆x i , summing over all i, and proceeding similarly to the earlier calculations, we arrive at the following expression where T = 1 2 Combining the two sums in the above expression we obtain It is easy to show that w b i,j − 2w d i,j = −1. This implies Substituting T in the equation (18), we get the desired result.
Proposition 3. The weights w b j,k and w d j,k are equal to 1 for uniform grids. Proof. Consider the uniform grids of the type x i = ih where h is a constant. This implies: Thus, with the definition of the index l jk , we have This readily implies that both the weights w b j,k and w d j,k are equal to 1. It it worth to mention that for uniform grids of the type x i = ih the discrete formulation takes a simpler form The above equation is also known as the standard discrete form of the coagulation population balance equation. Also, notice that the discretized formulation by [8] also reduces to the same form for such uniform grids. It should also be noted that the uniform grids and so the above simplified formulation are not suitable for aggregation problems. This is due to the requirement of a large number of grid points to discretize a domain of moderate size in the case of uniform grids. Therefore, to reduce the CPU time, some other special nonuniform grids are required to discretize large domains. In practice, uniform meshes on logarithmic scale or geometric meshes are often used to deal with such problems.
Since we are working with an explicit scheme, one has to impose some restrictions on time step to ensure positivity of the solution (often called the CFL condition). Similar to the condition proposed by [8], the following restriction on the time step ensures the positivity of the solution: Here B n i and D n i are discrete birth and death terms, respectively. They are given as Various other issues regarding, stability, order of convergence, and other mathematical analysis of the scheme are not the primary concern of this work and will be discussed elsewhere . The rest of the paper is devoted to the discussion of numerical results.
Remark 1. In this work, a fully discrete scheme is formulated and the main focus of discussion is centered on the discretization of the size variable. The proposed numerical scheme may be written in a semidiscrete form where the time variable remains continuous. Then, the resulting system of ordinary differential equations may be integrated in time using a higher order solver.
3. Numerical results. In order to validate the new finite volume scheme (NFVS) and to illustrate the improvements over the standard finite volume scheme (SFVS) by [8], numerical results of both schemes are compared with known analytical solutions and also for some physically relevant problems that do not have analytical solutions. In particular, some physically relevant kernels are considered and results of zeroth and second moments are compared with the generalized approximation (GA) method by [24]. [25] compared several methods and concluded that the GA method is more accurate for the computation of moments. This method relies on approximating the size distribution by a sum of I delta functions. There are two unknowns, namely the particle size distribution and the floating nodal point coordinates, need to be evaluated in this method. The mathematical formulation involves a set of 2I differential equations, I equations for each unknown.
As pointed out earlier, both the methods turn into the same formulation for uniform grids and hence we carry out all comparisons on nonuniform grids. The cell boundaries are generated with the rule x i+1/2 = rx i−1/2 for some real number r > 1, and their midpoints are taken as representatives of the cells. Since both the schemes produce the same results for fine grids, most of the numerical comparisons are performed for coarse grids to compare the methods and discuss their deviation from the true values. To ensure the positivity of the solution, time step is calculated using the CFL condition (19) at each time step. In order to know the extent of aggregation in numerical computations, it is convenient to express the zeroth moment in a dimensionless form. The index (degree) of aggregation, a dimensionless form of the zeroth moment, is defined as I agg = 1 − µ 0 (t)/µ 0 (0), where µ 0 (t) is the zeroth moment of the distribution at time t. The value of the index of aggregation mentioned in this work is reached at final time of computation.
The comparisons are based on a number of features including different normalized moments as well as particle size distribution. Normalization of moments is done by dividing the values at different times by the initial value of the moment. In particular, for the case of analytically tractable kernels, comparisons are performed to assess the accuracy of numerical results to predict the first three moments and the particle size distribution. Whereas, only zeroth and second moments are compared for application oriented kernels. Furthermore, it can be perceived from the formulation of the two methods that they both have similar complexity, which is also observed in CPU time required by the computer to run the simulations. So, both the methods take comparable CPU time to solve a particular problem and therefore we do not consider this attribute for comparison.
From a computational point of view, it is more difficult to handle the gelling problems. Therefore, we consider both gelling and non-gelling problems to show the effectiveness of the new technique developed here. In the gelling case, we emphasize more the prediction of moments to show that the new technique is able to predict results very near to the gelling point. For the non-gelling case, the evaluation of the complete particle size distribution along with its moments is analyzed.
3.1. Analytically tractable kernels. In this section, numerical results of the new scheme are compared with the existing finite volume scheme for test problems with known analytical results. We consider three types of kernels: size-independent, sum and product kernels. The analytical results can be found for various initial conditions and different kernels in [29]. Unless stated otherwise, all numerical test cases of this section consider the following exponential initial condition The minimum and maximum sizes of the computational domain are taken to be 10 −6 and 10 3 , respectively. The number of cells for these test problems is 30. The degree of aggregation varies for each problem and therefore it is mentioned individually in each case.
3.1.1. Size-independent kernel. Figure 1 shows the numerical results for the simplest case of size-independent kernel, β(t, x, y) = 1. The computation is carried out for a very large degree of aggregation, I agg = 0.95. The prediction of the zeroth moment is shown in Figure 1a. As expected, the numerical values of the zeroth moment using the NFVS overlap with its true values. Whereas, a consistent underprediction is observed in the numerical results of the SFVS. The tendency of loosing the first moment by the existing finite volume scheme is shown in Figure 1b. However, the loss in this case is not significant because this kernel does not rapidly produce particles of big sizes. Nonetheless, the prediction of the first moment by the new method is very accurate as it remains constant throughout the simulation time. The existing method, even at short times, overpredicts the second moment at great extent. On the other hand, the new method predicts extremely accurate results. This shows that prediction by the new technique is considerably better than by  Figure 1d, the particle population of each cell is plotted against the representative size of the respective cell. The numerical results show the improvement of the overprediction of particles at large particle sizes by the new technique over the existing technique. Clearly, the numerical results of particle distribution by the new method are highly accurate. Thus, it can be concluded that the new technique assigns the particles within the intervals more accurately without increasing the complexity in numerics.
3.1.2. Sum kernel. Figure 2 compares results of the two techniques for the sum kernel, β(t, x, y) = x + y. The degree of aggregation is taken to be 0.80. The deviation of the zeroth moment from its true values using the SFVS is more prominent than in the earlier case. The reason is the dependency of the aggregation kinetics on sizes of the particles. Since the whole population is captured by a small number of representatives and the aggregation kinetics depends on these representatives, the error propagation due to the discrete nature of the problem is faster than in the size-independent case. Nevertheless, it is interesting to observe that the decay of the zeroth moment is well captured by the new method and the prediction follows very closely the true values. The preservation property of the zeroth moment is well reflected in these results. The loss of the first moment in this case is significant by the SFVS as shown in Figure 2b. Interestingly, the new method captures the first moment very accurately. Figure 2c includes a comparison of the exact and numerical solutions for the variation of the second moment. Once again, the SFVS overpredicts the results and a diverging behavior is observed at high degree of aggregation. In contrast to the results of the SFVS, the new method slightly underpredicts the results of the second moment. Comparison of the particle distribution is made in Figure 2d. Contrary to the case of size-independent aggregation, the prediction is poor at large particle sizes by both the methods. The SFVS overpredicts the particle population whereas the new method underpredicts the results at large particle sizes. Moreover, one can also observe an underprediction of the numerical results by the SFVS in some part of the distribution. Overall, it seems that the results of the new method are more accurate than those of the SFVS.
3.1.3. Product kernel. In Figure 3, the numerical results computed with both the schemes are compared with the exact solution for the product kernel. The computation is made only up to I agg = 0.20. A very clear deviation of numerical results of the zeroth moment by the SFVS from the true values can be seen in Figure 3a. It is more prominent than in the earlier two cases because of stronger size dependency in this case. Once again, even for such a gelling kernel, the figure shows that the zeroth moment predicted by the new technique is in excellent agreement with the analytical results in the entire time range. Figure 3b shows that the SFVS technique starts loosing mass after a certain time while the new technique conserves mass up to the end of the simulation. This shows that the new method is more suitable for computing solutions of gelling problems, whereas the SFVS fails after a certain extent of aggregation. The comparison of the second moment is plotted in Figure  3c. The values of the moment diverge by the SFVS whereas a good prediction is observed by the new method. In Figure 3d, numerical and exact results of the particle size distribution are depicted. It is difficult to compare these results because of significant mass loss by the SFVS. However, a similar trend as in case of the earlier two kernels is observed. The results of the new method are underpredicting, whereas those of the SFVS are overpredicting.
Note that the product kernel is a gelling kernel for any arbitrary charge particle size distribution, see [30], however, the degree of aggregation at gelation depends on the shape of the charge distribution. In order to find gelation time numerically, we perform one more test case by taking mono-disperse particles of size unity as the initial condition. The exact value of the degree of aggregation at gelation is known to be 0.5 in this case, see [30]. The computational domain with the maximum particle size 3×10 3 is divided only into 30 nonuniform cells. Since our main focus in this case is to find the gelation time, only the first moment is compared. Figure 4 compares numerical solutions of the first moment with both the schemes. Since the degree of aggregation is more important in order to analyse the gelation phenomenon, the first moment is plotted with respect to the degree of aggregation instead of time. Figure 4 shows that after a certain degree of aggregation (I agg ≈ 0.45) the SFVS starts loosing mass while the new method conserves mass very closely to the gelling point (I agg ≈ 0.50). This, once again, shows that the new method is a powerful tool to identify the gelation time numerically.
In summary, the new technique, in each case, provides an excellent prediction of the zeroth and second moments as well as conservation of the first moment. In general, the SFVS tends to overpredict the number density at large particle sizes and, consequently, some mass loss and a diverging behavior of the second moment are observed. The results are obtained with a coarse grid. A finer grid can obviously be used to improve the accuracy of the numerical results to any desired level. The new method presented in this paper seems to be more adequate to compute moments very near to the gelation point at a low computational cost.   [25]. They considered the following initial condition with N 0 = 1, x 0 = √ 3/2 and σ = ln(4/3). As done in the previous section, the volume domain is discretized by the rule x i+1/2 = rx i−1/2 . These comparisons are based on the results of normalized zeroth and second moments. It should be mentioned here that [25] calculated zeroth and second moments for different number of nodal points. We have chosen their values corresponding to the parameter I = 4. The first three test cases are run for a very high degree of aggregation (≈ 0.99) that can be estimated from the value of the zeroth moment at final time. Whereas, the Numerical results obtained using SFVS, NFVS, and generalized approximation (GA) method are compared. Table 1 gives the zeroth and second moments obtained by the three methods at different dimensionless times. The maximum value of the size which ensures mass conservation is taken as 4 × 10 3 . The computational domain is discretized into 30 nonuniform cells. As can be seen from the table, the values of the zeroth moment using the new method are very close to the values obtained by the GA method, whereas, similar to the earlier observations, the SFVS underpredicts these values. Furthermore, prediction of the second moment by the NFVS is highly accurate. On the other hand, the SFVS once again overpredicts these results. However, as evident from the table, deviation of numerical values using the SFVS is not large because size dependency in Brownian kernel is not very strong. As a result, in this case the particle spectrum remains narrow and accurate results may be expected by numerical methods.

3.2.2.
Coagulation kernel β + . We extend the assessment by including the following kernel This kernel is computationally more expensive than the Brownian kernel as the kernel depends on its argument more strongly. To ensure mass conservation, the maximum size is taken as 8 × 10 5 . The computational domain is divided into 60 nonuniform cells. The final distribution in this case is much wider than that for Brownian kernel as can be seen from the numerical results for the second moment, summarized in Table 2. Once again, more or less the same observations as before are found in this case. The results of the zeroth moment by both the methods are in good agreement with the values obtained using the GA method. The values of the second moment using the NFVS are quite close to those obtained by the GA method. Whereas, the table clearly shows the over-prediction of the second moment by the SFVS.   Table 3. Comparison of the first and second moments for the gravitational kernel
Computation with this kernel is more difficult than in the earlier two cases because of the stronger size dependency. As mentioned by [25], the population balance equation becomes very stiff and difficult to solve and consequently computation becomes very expensive in this case. The maximum value of the size is taken as 2×10 6 , and the domain is divided into 100 nonuniform cells. Table 3 shows the results of zeroth and second moments using all three techniques. Even for such a strongly size dependent case, the results of the zeroth moment using the NFVS are in good agreement with the GA method, while GA results are underpredicted using the SFVS. The results of second moment using the SFVS technique are overpredicting at large times. However, the values using the NFVS seem to be much more accurate than the SFVS ones.
Like for the product kernel, second and higher moments of the distribution diverge at a finite time t gel for this kernel. The computational domain is divided into 50 cells and the maximum particle size goes up to 6 × 10 5 . Table 4 shows the values of the zeroth and second moments obtained using the SFVS and NFVS. The Table 4. Comparison of the first and second moments for the gelling kernel SFVS NFVS GA t µ 0 (t) µ 2 (t) µ 0 (t) µ 2 (t) µ 0 (t) µ 2 (t) 0 After that a sharp decay of the zeroth moment is observed. However, the prediction of the zeroth moment using the NFVS is matching very well with the results of the GA method. The values of the second moment using the SFVS increase rapidly after a very short time t = 0.1 and clearly diverge at larger times. The values of the second moment using the NFVS are slightly smaller at larger times than those obtained from the GA method. Overall, the results of the new finite volume scheme are comparable in quality to those obtained from the GA method.

4.
Conclusions. In this study, a new approach based on a finite volume method is presented for solving the aggregation population balance equation. The proposed method is computationally very efficient and easy to code. A new feature of the numerical method is presented by introducing preservation property of different moments which is beyond the requirement of merely conservation of the first moment. In particular, the concept of preservation of the zeroth moment is introduced and its importance is reflected through various numerical results. The scheme is validated and compared with the existing finite volume scheme against several test problems that include analytically tractable as well as physically based kernels. The new method produces acceptable results using only a few grid points. Moreover, it is found that the new method is capable of solving gelling problems whereas the standard finite volume scheme fails because of mass loss and consequently blow up of second moment for large times. The numerical results show the ability of the new method to predict very well the time evolution of the second moment as well as the complete particle size distribution.