A COMPARATIVE STUDY OF NUMERICAL APPROXIMATIONS FOR SOLVING THE SMOLUCHOWSKI COAGULATION EQUATION

1 BIOMATH, Department of Mathematical Modelling, Statistics and Bioinformatics, Faculty of Bioscience Engineering, Ghent University, 9000 Ghent, Belgium. 2 Department of Mathematics, Indian Institute of Technology Kharagpur, Kharagpur-721302, India. 3 Pharmaceutical Process Analytical Technology, Department of Pharmaceutical Analysis, Faculty of Pharmaceutical Sciences, Ghent University, 9000 Ghent, Belgium.


INTRODUCTION
Population balance equations (PBEs) generally describe the changes taking place in the particles due to certain mechanisms like coagulation (aggregation), breakage, growth and nucleation.Many other applications which involve these mechanisms are crystallization, granulation, fluidized bed etc., as illustrated in Marchal et al. (1988), Ramkrishna (2000), Peglow et al. (2005Peglow et al. ( , 2006) ) and Kettner et al. (2006).Mathematically, population balance equations are classified in the class of integro-partial differential equations of the hyperbolic or parabolic type.In the present work, our main concern is to solve the one dimensional Smoluchowski coagulation population balance equation, which in the continuous form can be written as (Hulburt and Katz, 1964): Taking ''size'' to refer to the mass v of a particle, n(t, u) is the size distribution at time t such that n(t, u)dt is the number concentration of particles with mass in (v, v + dv).The kernel n(v, u) is the rate constant for coagulation between sizes v and u and is determined by the physics of the specific process.The coagulation equation describes a wide range of physical phenomena such as polymerization, colloidal aggregation or agglomeration of granular materials and its solution has been the subject of a substantial body of literature (Aldous, 1999;Leyvraz, 2003) , , Except for a small class of kernels that lead to analytic solutions, the coagulation equation must be solved numerically.Several numerical techniques have been presented, including the Monte Carlo method (Smith and Matsoukas, 1998;Lee and Matsoukas, 2000;Lin et al., 2002), finite difference method (Mantzaris et al., 1999(Mantzaris et al., , 2001)), method of successive approximations (Ramkrishna, 2000), the method of moments (Madras and McCoy, 2004;Alexopoulos and Kiparissides, 2005;Marchisio and Fox, 2005;Attarakih et al., 2009Attarakih et al., , 2010)), finite volume schemes (Bennett and Rohani, 2001;Filbet and Laurenҫot, 2004;Bove et al., 2005;Kumar et al., 2008;Forestier-Coste and Mancini, 2012;Kumar et al., 2016), and Adomain Decomposition Method (Hasseine and Bart, 2015).These numerical methods merely focused on predicting the number distribution function accurately; however, the consistency of these solutions with experiments may be checked via the zeroth and first moments.Therefore, computing different order moments is also of some interest.
There are other numerical methods developed by Hounslow et al. (1988), Litster et al. (1995), Kumar and Ramkrishna (1996) and Kumar et al. (2006) which account for the accurate prediction of the number distribution function along with various order moments.Among these methods the 'Cell Average Technique (CAT)' developed by Kumar et al. (2006) is the most accurate and efficient method for solving the aggregation PBE.This technique relies on distributing the average properties of particles entering into a cell due to the aggregation to the neighboring nodes.Due to averaging of particles properties, the results computed for particle property distribution and different order moments were improved.Previously, a study was conducted by Kumar et al. (2009) for comparing the cell average technique (Kumar et al., 2006) and finite volume scheme (Filbet and Laurenҫot, 2004).The results revealed that both methods predict the number distribution functions with comparable accuracy, but the various order moments are more accurately predicted by the cell average technique as compared to the finite volume scheme.The major drawback of both numerical approximations is that the mathematical formulations of both numerical methods were highly complex.In contrast to this study, our aim is to compare a newly developed finite volume scheme (Singh et al., 2015) with the cell average technique (Kumar et al., 2006) for obtaining the solution of an aggregation PBE on non-uniform grids.
The paper is structured as follows: our exercise begins with a brief introduction of the existing cell average technique as well as the finite volume scheme on non-uniform meshes.We then discuss the numerical results to prove the efficiency as well as accuracy of both the numerical methods.Finally, we summarize the conclusions of this study.

NUMERICAL METHODS
In this section, the mathematical formulations of the existing cell average technique (Kumar et al., 2006), as well as the finite volume scheme (Singh et al., 2015) for solving the aggregation PBE (1) on non-uniform meshes are provided.The idea of both the numerical methods is based on the assumption that particles within a grid cell are concentrated on its representative.

Cell Average Technique (CAT)
A finite one dimensional computational domain with upper limit u max < ∞ is taken and divided into I number of smaller cells having u i as representative volume, for i = 1, 2, 3, ..., I (see Figure 1).The whole idea of the CAT is to convert the continuous eq. ( 1) into a set of I ordinary differential equations (ODEs).The mesh points and the step are: We denote the number of particles in the i th cell by N i , ( ) The distribution is then represented by a series of delta functions, ( ) ( ) ∑ Substituting this form into (1) we obtain a series of ODE's of the form: with the birth term given by: (2) The CAT is based upon the concept of averaging all properties of the particles in a cell and apportioning the property between the two nodes that bracket the cell.Accordingly, the mass of new particles that are formed inside the ith cell is re-distributed between nodes i and i + 1 by the respective factors a and b whose values are determined such that the zeroth and first moments of the population are strictly conserved.These factors are obtained by solving where the Heaviside step function is defined by: ( ) ( ) and the death term by: ( ) Here, is the mean size of new particles that form by aggregation in the ith cell, and V i is the flux of particle mass into the cell, Solving the system of equations ( 7), we get: where the functions and are given by: ( ) ( ) In summary, after collecting all contributions to the birth at u_i, the modified birth takes the following form ( ) ( ) Putting the above expression in equation ( 4), the final set of ODE's formed is: The detailed formulation of the cell average technique can be found in Kumar et al. (2006).

Finite Volume Scheme (FVS)
CAT is quite expensive because of the added cost of averaging the particle mass and reassigning its value to neighboring nodes.Recently we have formulated a finite volume scheme (FVS) whose mathematical formulation is much simpler than the CAT as no averaging of properties and distribution to the neighboring nodes is required.The discretization of the size domain on FVS is the same as in CAT, with each size interval (u i-1/2 , u i+1/2 ) represented by the midpoint u i , and with the number of particles N i in the interval represented by the integral in equation ( 2).In the FVS method the coagulation equation in equation (1) becomes a set of ordinary equations for N i of the form: The first summation on the right-hand side is taken over all pairs (j, k) such that the sum of midpoints u j + u k falls in interval i.This condition may be expressed formally as: ) The factor j i, j, k is defined as: , , and ensures that mass is conserved.As shown in Singh et al. (2015), this scheme conserves the total mass, as long as no mass leaves the upper boundary of the size domain.In reality, mass may leak through the upper boundary if the size domain does not extend to sufficiently large values.
For the cell average technique, there is always a possibility of numerical diffusion when the particle properties are distributed to the neighboring nodes after the averaging of the particleproperties (Kumar et al., 2006).But, for the finite volume scheme, no such distribution is required as only a correction function is added to the formulation for ensuring the mass conservative property.

RESULTS AND DISCUSSION
To establish the accuracy of the FVS relative to the CAT, we ran various test cases, including kernels with known analytical solutions, as well as kernels that appear frequently in the literature of particle technology, including the Brownian kernel.We use the aggregation index I agg ,

Analytically Tractable Kernels for Pure Aggregation PBE
Here we test the numerical methodologies against three kernels whose analytical solutions are known in the literature.These are the constant kernel, b(u, v) = 1, the sum kernel, b(u, v) = u + v, and the product kernel, b(u, v) = uv.For the initial condition n(0, u) = e -u , the exact results of number density as well as various order moments were found by Scott (1968).
Constant kernel -In this section, the accuracy of the numerical methods is investigated for the sizeindependent kernel, i.e., b(u, v) = 1, corresponding to the initial condition n(0, u) = e -u .The forward simulation is run for 0.90 degree of aggregation on a non-uniform grid.It consists of 30 cells, ranging from u min = 2 × 10 -6 to u max = 600 whereas time ranges from 0 to 18.
The comparison of different order moments along with the number density obtained by both numerical methods with the exact results are shown in Figure (2).It is revealed that the zeroth and first order moments predicted by both the methods are in very good agreement with the exact results.Moreover, the second order moment is slightly better predicted by the FVS as compared to the CAT.In addition, the number density function computed by FVS and CAT coincides with the exact number density.However, the results obtained for distinct order moments by numerical methods are comparable, but the computational time taken by FVS is 0.2374s whereas CAT took 0.5746s.This shows that the finite volume scheme is significantly more efficient and accurate than the CAT.The difference is not that large for a single simulation, but becomes important in cases where many model simulations are required (e.g., Monte Carlo based optimization).
Sum kernel -Similar to the constant kernel, the moments and number density obtained using FVS and CAT on non-uniform grids are compared for the sum kernel (i.e., b(u, v) = u + v) in Figure (3).Here, the simulations are run up to 0.80 degree of aggregation.The domain considered for the comparison ranges from u min = 2 × 10 -6 to u max = 2600, divided into 30 parts, whereas time ranges from 0 and 1.2.
Figure (3) indicates that the results for the zeroth and first order moments predicted by both numerical methods coincide with the exact moments.Although both methods under-predicted the second order moments, the moment obtained using CAT show more deviation than FVS.Further, the comparison of the number density computed using the numerical methods with the exact number density is shown in Figure 3(d).The comparison affirms that the number density approximated by FVS and CAT coincides with the exact number density.

( ) ( ) (
) to track the extent of aggregation.Since mass is constant, the index I agg has a simple relationship to the mean size ( ) ( ) ( ) where u 0 is the mean size at time 0. All computations were carried out using MATLAB on an i5 CPU with 2.67 GHz and 8 GB RAM.The same grid is used in both the FVS and CAT methods.To enhance the comparison, we also calculated the relative weighted errors in various order moments using the expression: where the superscripts ana and num denotes the analytical and numerical values for the number distribution function.Moreover, D 0 and D 1 shows the errors in the zeroth and first order moments, respectively.Furthermore, the quantitative weighted errors in the different order moments computed using the expression ( 18) are shown in Table 1.It reveals that the weighted errors in the various order moments decrease as the number of grids increases.Hence, the accuracy of the moments can be improved to the desired levels by choosing a denser grid.
Product kernel -Having compared the different order moments and number density obtained using different numerical schemes with the exact results for constant and sum kernels, a similar comparison was made for the product kernel, i.e., b(u, v) = uv.The computational domain consists of 30 size bins, ranging from u min = 10 -6 to u max = 1300.Moreover, the product kernel is a gelling kernel for any arbitrary particle size distribution (Smit et al., 1994).Therefore, the simulation in this case was performed until the product kernel showed the gelling behavior.The kernel exhibits gelling behavior when the new forming particles aggregate at a greater frequency than their parents.
The comparison of different order moments predicted by both numerical schemes is shown in Figure 4.One can see that the zeroth order moments predicted by FVS and CAT show good agreement with the exact moment.Furthermore, it is also observed that the first order moments predicted by FVS show more accurate prediction as compared to CAT.Similarly, the comparison of the second order moment (m 2 (t)) indicates that the second order moment is more accurately approximated by the FVS, whereas the CAT shows more under-prediction.
Analogous to the constant and sum kernels, the behavior of number density obtained by both numerical schemes against the exact number density is also tested graphically.Figure 4(d) shows that the number density found by FVS and CAT coincides with the exact result.
Table 1.Weighted errors in various order moments using the sum kernel.Table 2. Weighted errors in various order moments using the product kernel.
Table 3.Comparison of the zeroth and second order moments for the Brownian kernel.
However, on the computational basis, FVS computes all numerical results with less computational time, i.e., the FVS and CAT took 0.2890 s and 0.4831 s, respectively.It can be concluded that the FVS is more suitable for solving the pure aggregation PBE with the kernel which exhibits gelling behavior.In addition to the above, it can be observed that the number density functions obtained by both numerical methods are equal in accuracy to all analytically tractable kernels, but the higher order moments, in particular the second order moments predicted by the CAT, deviate more from the exact results.This is because the CAT only accounts for the conservation of the zeroth and first order moments.However, in contrast to the CAT, the FVS is even consistent with the higher order moments.
Similar to the sum kernel, the quantitative weighted errors shown in Table 2 are also calculated for the product kernel.It can be seen that the errors in the moments reduce to a certain level by choosing a more refined grid.Hence, it can be concluded that the errors in the various order moments certainly depend on how dense a grid is chosen for a discretization.

Physically Relevant Kernels for the Pure Aggregation PBE
In this section, comparisons of normalized zeroth order and second order moments obtained using the generalized approximation (GA) method (Piskunov and Golubev, 2002) with FVS and CAT are made.For the comparison, various practical oriented kernels like Coalescence, Gravitational and Gelling kernels have been considered.Piskunov et al. (2002) chose the following initial condition Brownian kernel, b b -First, the Brownian kernel is taken into the consideration for comparing the normalized zeroth and second order moments predicted by both numerical approximations.Mathematically, the Brownian kernel is of the following form: with N 0 = 1, and .Piskunov et al. (2002) found the normalized zeroth order and second order moments for different number of cells.

( ) (
)( ) It can be seen that this kernel is heavily dependent on the size of the particles and many of its applications related to aerosol (including colloidal sol, liquid droplets, air bubbles etc.) can be found in Mountain et al. (1986) and Kostoglou and Konstandopoulos (2001).For the comparison, the computational domain ranges from u min = 10 -6 to u max = 5 × 10 3 and is divided into 30 cells.
The values of the normalized zeroth and second order moments approximated by both numerical methods and the GA method are provided in Table 3.It reveals that, as we move from time t = 0 to 50, FVS is very consistent in predicting the zeroth order moment, whereas the CAT shows more under prediction from the GA method.Furthermore, the values of the second order moment computed by the FVS are closer to the GA method as compared to the CAT.
Coalescence kernel, b + -Here we choose the Coalescence Kernel, i.e., b + (u, v) = u 2/3 + v 2/3 for the comparisons of quantitative normalized zeroth order and second order moments obtained by all numerical schemes with the GA method.The degree of aggregation for the comparison is taken as 0.99.For the testing, the computational domain is u max = 8 × 10 5 , whereas the number of grid points is taken as 60.The application of this kernel can be found in Ruckenstein and Pulvermacher (1973).
Table 4 shows that the results of the zeroth order moments found by both numerical methods are in good agreement with the values obtained by the GA method.However, the value of the second order moment accessed by FVS is closer to the value found by the GA method, whereas the table clearly reveals that the second order moment found by the CAT slightly deviated from the moment obtained by using GA method.This shows that the FVS is a better approximation for predicting zeroth and second order moments than the CAT.
Gravitational kernel, b g -Similar to the case of the Coalescence kernel, the same comparison is also done for the Gravitational kernel, i.e., b g (u, v) = (u 1/3 + v 1/3 ) 2 |u 1/6 -v 1/6 |.The Gravitational kernel is more complex than the Coagulation kernel and highly dependent on size.For the comparison, the same degree of aggregation is taken as the Coalescence kernel case, whereas the number of grid points is taken as 100.For the comparison, the computational domain is considered as u max = 2 × 10 6 .
Table 5 shows the comparison of the numerical values of the zeroth and second order moments found by FVS and CAT with values obtained the GA method.The results reveal that the FVS predicts both order moments with more accuracy than the CAT.Thus it can be concluded that FVS provides a better approximation of the different order moments than the CAT.
Gelling kernel, b c -In this subsection, a very complex kernel, namely the Gelling kernel, is considered for conducting the comparison.The mathematical representation of the Gelling kernel is given as: The degree of aggregation for the comparison is taken as 0.50 and the domain with u max = 6 × 10 5 is divided into 50 non-uniform cells.It can be seen from Table 6 that the zeroth order moment predicted by the FVS is in good agreement with the GA method, whereas CAT deviates more from the exact moment.One can also observe from Table 6 that the second order moment obtained by the FVS nearly coincides with the GA method as compared to the CAT.Furthermore, at time t = 0.4s, Table 6 reveals the divergence of the second order moment obtained by the CAT due to the gelation effect, whereas the same moment obtained by the FVS shows good agreement with the moment found by the GA method.

Analytically Tractable Kernels for the Pure Breakage PBE
In this section, the numerical results obtained using FVS and CAT for the pure breakage PBE (Ziff and McGrady, 1985) are compared with the exact results.Binary breakage functions corresponding to linear and quadratic selection functions are chosen to demonstrate the comparison.For the binary breakage function, only two objects are produced per fragmentation event.The analytical results for these kernels are provided by Ziff and McGrady (1985) for exponential    ( ) Comparison of the zeroth and second order moments for the Gelling kernel.
and monodisperse initial conditions.For our study, similar to the aggregation problem, the exponential initial condition, i.e., n(0, u) = e -u is considered.For the simulations, the time limits vary from 0 to 30 and the number of cells was taken to be 40.
Linear Selection Function -Here, the comparison of the numerical results for the pure breakage PBE corresponding to the linear selection function S(u) = u and binary breakage function b(u, v) = 2/v are conducted.For the comparison, the size domain range from u min = 10 -10 to u max = 10 3 is considered.
Figure 5 demonstrates the comparison of the numerical results with the exact results.It reveals that the zeroth order moment is more accurately predicted by the FVS as compared to the CAT.Moreover, both numerical schemes show the conservation of the total mass of the system at various times and matches well with the exact result.Additionally, the second order moment computed by the FVS shows better prediction than the CAT.Similar to the aggregation problem, the number density functions are obtained for the breakage problem as shown in Figure 5(d).The numerical results reveal that both schemes capture the number density functions with comparable accuracy.
Quadratic Selection Function -Similar to the previous case, the accuracy of the numerical schemes is checked by performing a comparison for a quadratic function, i.e., S(u) = u 2 corresponding to the binary breakage function.For the comparison, the size domain range from u min = 10 -10 to u max = 10 5 is considered.
The comparison of the numerical approximations in terms of various order moments, as well as number density functions, is shown in Figure 6.It reveals that the zeroth order moment is far better predicted by the FVS, whereas the CAT shows more deviation from the exact moment.However, both schemes computed the first order moment with equal accuracy, coinciding with the exact result (see Figure 6(b)).Furthermore, the second order moment calculated by the FVS shows better accuracy than the CAT.Additionally, the number density functions are equally well approximated by both numerical schemes.

CONCLUSIONS
This work presents a performance study of the cell average technique and recently developed finite volume scheme for solving the pure aggregation population balance equation on non-uniform meshes.The qualitative and quantitative comparisons of the cell average technique and finite volume scheme in terms of different order moments are shown for both analytically tractable as well as physically relevant obtained by the finite volume scheme are in better agreement with the exact results in comparison to the cell average technique.Moreover, the finite volume scheme not only conserves the total mass of particles in the system (first moment) and is consistent in the total number of particles (zeroth moment), but it also predicted higher order moments along with number density distributions more accurately than the cell average technique.Furthermore, the finite volume scheme is more efficient from a computational point of view.Therefore, the finite volume scheme is an attractive approximation for the one-dimensional aggregation population balance equation.

Figure 2 .
Figure 2. Different order moments and number density for the constant kernel, i.e., b(u, v) = 1.

Figure 3 .
Figure 3. Different order moments and number density for the sum kernel, i.e., b(u, v) = u + v.

Figure 4 .
Figure 4. Different order moments and number density for the product kernel, i.e., b(u, v) = uv.

Figure 5 .
Figure 5. Different order moments and number density for the linear selection function.

Figure 6 .
Figure 6.Different order moments and number density for the quadratic selection function.

Table 4 .
Comparison of the zeroth and second order moments for the Coagulation kernel.

Table 5 .
Comparison of the zeroth and second order moments for the Gravitational kernel.