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

ABSTRACT

In this work, numerical approximations for solving the one dimensional Smoluchowski coagulation equation on non-uniform meshes has been analyzed. Among the various available numerical methods, finite volume and sectional methods have explicit advantage such as mass conservation and an accurate prediction of different order moments. Here, a recently developed efficient finite volume scheme (Singh et al., 2015) and the cell average technique (Kumar et al., 2006) are compared. The numerical comparison is established for both analytically tractable as well as physically relevant kernels. It is concluded that the finite volume scheme predicts both number density as well as different order moments with higher accuracy than the cell average technique. Moreover, the finite volume scheme is computationally less expensive than the cell average technique.

Key words:
Aggregation; Particles; Population balance equation; Finite volume scheme; Cell average technique; Non-uniform grids

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. (2005, 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):

$$\frac{\partial n\left(t,u\right)}{\partial }=\frac{1}{2}{\int }_o^u\beta \left(v,u-v\right)n\left(t,u-v\right)n\left(t,v\right)dv-{\int }_o^{\infty }\beta \left(u,v\right)n\left(t,u\right)n\left(t,v\right)dv$$ (1)

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, 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., 2009, 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:

$$u_{\frac{1}{2}}=0,u_i=\frac{\left(u_{i+{\displaystyle \frac{1}{2}}}+u_{i-\frac{1}{2}}\right)}{2},∆u_i=u_{i+\frac{1}{2}}-u_{i-\frac{1}{2}}$$

Figure 1
One dimensional domain discretization.

We denote the number of particles in the i th cell by N_i,

$$N_i={\int }_{u_{i-\frac{{\displaystyle 1}}{{\displaystyle 2}}}}^{u_{i+\frac{1}{2}}}n\left(t,u\right)du$$ (2)

The distribution is then represented by a series of delta functions,

$$n\left(t,u\right)=\sum _{j=1}^I{N}_i\delta \left(u-u_j\right)$$ (3)

Substituting this form into (1) we obtain a series of ODE’s of the form:

$$\frac{d{N}_i}{dt}=B_i-D_i$$ (4)

with the birth term given by:

$$B_i=\sum _{u_{i-\frac{1}{2}}\le \left(u_j+u_k\right)<u_{i+\frac{1}{2}}}^{j\le k}\left(1-\frac{1}{2}{\delta }_{j,k}\right)\beta \left(u_j,u_k\right)N_j{N}_k$$ (5)

and the death term by:

$$D_i=\sum _{j=1}^I\beta \left(u_i,u_j\right)N_i{N}_j$$ (6)

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

$$a{u}_i+b{u}_{i+1}=B_i{\overline{u}}_{i}a+b=B_i$$ (7)

Here, $${\overline{u}}_i$$ is the mean size of new particles that form by aggregation in the ith cell,

$${\overline{u}}_i=\frac{V_i}{B_i}$$ (8)

and V_i is the flux of particle mass into the cell,

$$V_i=\sum _{u_{i-\frac{1}{2}}\le \left(u_j+u_k\right)<u_{i+\frac{1}{2}}}^{j\le k}\left(1-\frac{1}{2}{\delta }_{j,k}\right)\beta \left(u_j,u_k\right)N_j{N}_k\left(u_j+u_k\right)$$ (9)

Solving the system of equations (7), we get:

$$a=b_i{\lambda }_i^{+}\left({\overline{u}}_i\right),b=b_i{\lambda }_{i+1}^{-}\left({\overline{u}}_i\right)$$ (10)

where the functions $${\lambda }_i^{+}$$ and $${\lambda }_i^{-}$$ are given by:

$${\lambda }_i^{\pm }=\frac{\left(u-u_{i\pm 1}\right)}{\left(u_i-u_{i\pm 1}\right)}$$ (11)

In summary, after collecting all contributions to the birth at u_i, the modified birth takes the following form

$$B_i^{mod}=B_{i-1}{\lambda }_i^{-}\left({\overline{u}}_{i-1}\right)H\left({\overline{u}}_{i-1}-x_{i-1}\right)+B_i{\lambda }_i^{-}\left({\overline{u}}_i\right)H\left(x_i-{\overline{u}}_i\right)+B_i{\lambda }_i^{+}\left({\overline{u}}_i\right)H\left({\overline{u}}_i-x_i\right)+B_{i+1}{\lambda }_i^{+}\left({\overline{u}}_{i+1}\right)H\left(x_{i+1}-{\overline{u}}_{i+1}\right)$$ (12)

Putting the above expression in equation (4), the final set of ODE’s formed is:

$$\frac{d{N}_i}{dt}=B_{i-1}{\lambda }_i^{-}\left({\overline{u}}_{i-1}\right)H\left({\overline{u}}_{i-1}-x_{i-1}\right)+B_i{\lambda }_i^{-}\left({\overline{u}}_i\right)H\left(x_i-{\overline{u}}_i\right)+B_i{\lambda }_i^{+}\left({\overline{u}}_i\right)H\left({\overline{u}}_i-x_i\right)+B_{i+1}{\lambda }_i^{+}\left({\overline{u}}_{i+1}\right)H\left(x_{i+1}-{\overline{u}}_{i+1}\right)-\sum _{j=1}^I\beta \left(u_i,u_j\right)N_i{N}_j$$ (13)

where the Heaviside step function is defined by:

$$H\left(u-u_i\right)=\left\{\begin{array}{c}1,u-u_i>0\\ \frac{1}{2},u-u_i=0\\ 0,u-u_i<0\end{array}\right.$$ (14)

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:

$$\frac{d{N}_i}{dt}=\sum _{\left(j,k\right)\in y^i}\beta \left(u_j,u_k\right)N_j{N}_k\frac{∆u_j∆u_k}{∆u_i}{\phi }_{i,j,k}-\sum _{j=1}^I\beta \left(u_i,u_j\right)N_i{N}_j∆u_j$$ (15)

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:

$${𝛾}^i=\left\{\left(j,k\right)N\times N:u_{i-\frac{1}{2}}\le \left(u_j+u_k\right)<u_{i+\frac{1}{2}}\right\}$$

The factor φ_{i, j, k} is defined as:

$${\phi }_{i,j,k}=\frac{u_j+u_k}{u_i}$$

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,

$$I_{agg}\left(t\right)=1-\frac{{\mu }_0\left(t\right)}{{\mu }_0\left(t=0\right)}$$ (16)

to track the extent of aggregation. Since mass is constant, the index I_agg has a simple relationship to the mean size

$$\overline{u}\left(t\right)=1-\frac{{\mu }_1\left(t\right)}{{\mu }_0\left(t\right)}$$

$$\frac{\overline{u}\left(t\right)}{u_0}=\frac{1}{1-I_{agg}}$$ (17)

where u₀ 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:

$${∆}_j\left(t\right)=\sum _{j=1}^I\left|N_j^{ana}\left(t\right)-N_j^{num}\left(t\right)\right|u_i^j/\sum _{j=1}^I\left|N_j^{ana}\right|u_i^j$$ (18)

where the superscripts ana and num denotes the analytical and numerical values for the number distribution function. Moreover, Δ₀ and Δ₁ shows the errors in the zeroth and first order moments, respectively.

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, β(u, v) = 1, the sum kernel, β(u, v) = u + v, and the product kernel, β(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 size-independent kernel, i.e., β(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).

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

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., β(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
Different order moments and number density for the sum kernel, i.e., β(u, v) = u + v.

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.

Moreover, the comparison based on computational time shows that FVS is less expensive than CAT, as the times taken by FVS and CAT are 0.2701 s and 0.4746 s, 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.

Table 1
Weighted errors in various order moments using the sum kernel.

Δ	CAT (with 40 grids)	FVS (with 40 grids)	CAT (with 50 grids)	FVS (with 50 grids)
Δ0	0.02801	0.02669	0.01716	0.01520
Δ1	0.02801	0.02669	0.01716	0.01520
Δ2	0.02801	0.02669	0.01716	0.01520

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., β(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 (μ₂ (t)) indicates that the second order moment is more accurately approximated by the FVS, whereas the CAT shows more under-prediction.

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

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. 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.

Table 2
Weighted errors in various order moments using the product kernel.

Δ	CAT (with 40 grids)	FVS (with 40 grids)	CAT (with 50 grids)	FVS (with 50 grids)
Δ0	0.03498	0.03130	0.02187	0.02123
Δ1	0.03498	0.03130	0.02187	0.02123
Δ2	0.03498	0.03130	0.02187	0.02123

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

$$n^{in}\left(u\right)=\frac{1}{\sqrt{2\pi }}\frac{N_0}{\sigma u}{e}^{-{\ln }^2\left(\frac{u/u_0}{2{\sigma }^2}\right)}$$

with N₀ = 1, $$u_0=\raisebox{1ex}{$\sqrt{3}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$$ and $$\sigma =\sqrt{\ln \left(\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.\right)}$$. Piskunov et al. (2002) found the normalized zeroth order and second order moments for different number of cells.

Brownian kernel, β_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:

$${\beta }_b\left(u,v\right)=\left(u^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}+v^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\right)\left(u^{\raisebox{1ex}{$-1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}+v^{\raisebox{1ex}{$-1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\right)$$

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³ 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.

Table 3
Comparison of the zeroth and second order moments for the Brownian kernel.

T	CAT		FVS		GA
	μ0 (t)	μ2 (t)	μ0 (t)	μ2 (t)	μ0 (t)	μ2 (t)
0	1.00	1.00	1.00	1.00	1.00	1.00
1	0.359	2.74	0.343	3.27	0.326	4.09
5	9.26 × 10-2	11.14	8.70 × 10-2	14.20	8.68 × 10-2	16.54
10	4.55 × 10-2	24.12	4.49 × 10-2	31.24	4.50 × 10-2	32.11
50	7.40 × 10-3	194.58	7.8 × 10-3	153.86	9.26 × 10-3	157.14

Coalescence kernel, β₊ - Here we choose the Coalescence Kernel, i.e., β₊(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⁵, 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.

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

T	CAT		FVS		GA
	μ0 (t)	μ2 (t)	μ0 (t)	μ2 (t)	μ0 (t)	μ2 (t)
0	1.00	1.00	1.00	1.00	1.00	1.00
1	0.461	3.35	0.454	3.36	0.443	3.88
5	7.75 × 10-2	39.81	7.73 × 10-2	44.95	7.05 × 10-2	54.83
10	2.32 × 10-2	1.94 × 102	2.35 × 10-2	2.18 × 102	2.07 × 10-2	2.75 × 102
50	7.82 × 10-4	1.51 × 104	7.99 × 10-4	1.68 × 104	6.84 × 10-4	2.27 × 104

Gravitational kernel, β_g - Similar to the case of the Coalescence kernel, the same comparison is also done for the Gravitational kernel, i.e., β_g(u, v) = (u^1/3 + v^1/3)² |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⁶.

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.

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

T	CAT		FVS		GA
	μ0 (t)	μ2 (t)	μ0 (t)	μ2 (t)	μ0 (t)	μ2 (t)
0	1.00	1.00	1.00	1.00	1.00	1.00
1	0.721	1.58	0.722	1.57	0.832	1.47
5	0.207	10.02	0.205	9.02	0.314	9.02
10	4.8 × 10-2	68.88	4.9 × 10-2	66.65	7.90 × 10-2	67.07
50	7.86 × 10-5	7.12 × 104	9.20 × 10-5	2.88 × 104	1.34 × 10-4	1.83 × 104

Gelling kernel, β_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:

$${\beta }_{gc}\left(u,v\right)={\left(u^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}+v^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\right)}^2\left|u^{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}-v^{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}\right|$$

The degree of aggregation for the comparison is taken as 0.50 and the domain with u_max = 6 × 10⁵ 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.

Table 6
Comparison of the zeroth and second order moments for the Gelling kernel.

T	CAT		FVS		GA
	μ0 (t)	μ2 (t)	μ0 (t)	μ2 (t)	μ0 (t)	μ2 (t)
0.0	1.00	1.00	1.00	1.00	1.00	1.00
0.1	0.864	1.35	0.871	1.29	0.926	1.20
0.2	0.730	2.13	0.745	1.73	0.847	1.54
0.3	0.596	10.05	0.626	2.42	0.761	2.21
0.4	0.448	1.40 × 103	0.512	3.56	0.668	3.71

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 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³ 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.

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

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² corresponding to the binary breakage function. For the comparison, the size domain range from u_min = 10^-10 to u_max = 10⁵ 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.

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

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 kernels. The moment errors indicate that the results 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.

REFERENCES

• Aldous, D. J. Deterministic and stochastic models for coalescence (aggregation and coagulation): A review of the mean-field theory for probabilists. Bernoulli 5(1), 3-48 (1999).
• Alexopoulos, A. H. and Kiparissides, C. A., Part II: Dynamic evolution of the particle size distribution in particulate processes undergoing simultaneous particle nucleation, growth and aggregation. Chemical Engineering Science 60(15), 4157-4169 (2005).
• Attarakih, M., Jaradat, M., Drumm, C., Bart, H., Tiwari, S., Sharma, V., Kuhnert, J. and Klar, A., A multivariate sectional quadrature method of moments for the solution of the population balance equation. Computer Aided Chemical Engineering 28, 1551-1556 (2010).
• Attarakih, M. M., Drumm, C. and Bart, H. J., Solution of the population balance equation using the sectional quadrature method of moments (SQMOM). Chemical Engineering Science 64(4), 742-752 (2009).
• Bennett, M. and Rohani, S., Solution of population balance equations with a new combined lax-wendro_/crank-nicholson method. Chemical Engineering Science 56(23), 6623- 6633 (2001).
• Bove, S., Solberg, T. and Hjertager, B. H., A novel algorithm for solving population balance equations: the parallel parent and daughter classes. Derivation, analysis and testing. Chemical Engineering Science 60(5), 1449-1464 (2005).
• Filbet, F. and Laurenҫot, P., Numerical simulation of the smoluchowski coagulation equation. SIAM Journal on Scientific Computing 25(6), 2004-2028 (2004).
• Forestier-Coste, L. and Mancini, S., A finite volume preserving scheme on nonuniform meshes and for multidimensional coalescence. SIAM Journal on Scientific Computing 34(6), B840-B860 (2012).
• Goldschmidt, M. J. V., Hydrodynamic modelling of fluidised bed spray granulation. Twente University Press (2001).
• Hasseine, A. and Bart, H. J., Adomian decomposition method solution of population balance equations for aggregation, nucleation, growth and breakup processes. Applied Mathematical Modelling 39(7), 1975-1984 (2015).
• Hounslow, M., Ryall, R. and Marshall, V., A discretized population balance for nucleation, growth, and aggregation. AIChE Journal 34(11), 1821-1832 (1988).
• Hulburt, H. and Katz, S., Some problems in particle technology: A statistical mechanical formulation. Chemical Engineering Science 19(8), 555-574 (1964).
• Kettner, C., Peglow, M., Metzger, T. and Tsotsas, E., Distributed product quality in fluidized bed drying. In Proceedings of the 15th International Drying Symposium, Budapest, Hungary (2006).
• Kostoglou, M. and Konstandopoulos, A. G., Evolution of aggregate size and fractal dimension during brownian coagulation. Journal of Aerosol Science 32(12), 1399-1420 (2001).
• Kumar, J., Kaur, G. and Tsotsas, E., An accurate and efficient discrete formulation of aggregation population balance equation. Kinetic and Related Models 9(2), 373-391 (2016).
• Kumar, J., Peglow, M., Warnecke, G., Heinrich, S. and Mörl, L., Improved accuracy and convergence of discretized population balance for aggregation: The Cell Average Technique. Chemical Engineering Science 61(10), 3327-3342 (2006).
• Kumar, J., Warnecke, G., Peglow, M. and Heinrich, S., Comparison of numerical methods for solving population balance equations incorporating aggregation and breakage. Powder Technology 189(2), 218-229 (2009).
• Kumar, J., Warnecke, G., Peglow, M. and Tsotsas, E., A note on sectional and finite volume methods for solving population balance equations. In Micro-Macro-interaction, pp 285-297. Springer Science and Business Media (2008).
• Kumar, S. and Ramkrishna, D., On the solution of population balance equations by discretization I. a fixed pivot technique. Chemical Engineering Science 51(8), 1311-1332 (1996).
• Lee, K. and Matsoukas, T., Simultaneous coagulation and break-up using constant- N Monte Carlo. Powder Technology 110(1), 82-89 (2000).
• Leyvraz, F., Scaling theory and exactly solved models in the kinetics of irreversible aggregation. Physics Reports 383(2-3), 95-212 (2003).
• Lin, Y., Lee, K. and Matsoukas, T., Solution of the population balance equation using constant-number monte carlo. Chemical Engineering Science 57(12), 2241-2252 (2002).
• Litster, J., Smit, D. and Hounslow, M., Adjustable discretized population balance for growth and aggregation. AIChE Journal 41(3), 591-603 (1995).
• Madras, G. and McCoy, B., Reversible crystal growth-dissolution and aggregation- breakage: Numerical and moment solutions for population balance equations. Powder Technology 143, 297-307 (2004).
• Mandarins, N. V., Daoutidis, P. and Srienc, F., Numerical solution of multi-variable cell population balance models: I. finite difference methods. Computers and Chemical Engineering 25(11), 1411-1440 (2001).
• Mantzaris, N. V., Liou, J. J., Daoutidis, P. and Srienc, F., Numerical solution of a mass structured cell population balance model in an environment of changing substrate concentration. Journal of Biotechnology 71(1), 157-174 (1999).
• Marchal, P., David, R., Klein, J. and Villermaux, J., Crystallization and precipitation engineering- I. an efficient method for solving population balance in crystallization with agglomeration. Chemical Engineering Science 43(1), 59-67 (1988).
• Marchisio, D. L. and Fox, R. O., Solution of population balance equations using the direct quadrature method of moments. Journal of Aerosol Science 36(1), 43-73 (2005).
• Mountain, R. D., Mulholland, G. W. and Baum, H., Simulation of aerosol agglomeration in the free molecular and continuum ow regimes. Journal of Colloid and Interface Science 114(1), 67-81 (1986).
• Peglow, M., Kumar, J. and Mörl, L., Investigation of coalescence kinetics of microcrystalline cellulose in fluidised bed spray agglomeration: experimental studies and modelling approach. Brazilian Journal of Chemical Engineering 22(2), 165-172 (2005).
• Peglow, M., Kumar, J., Warnecke, G., Heinrich, S. and Mörl, L., A new technique to determine rate constants for growth and agglomeration with size-and time-dependent nuclei formation. Chemical Engineering Science 61(1), 282-292 (2006).
• Piskunov, V. and Golubev, A., The generalized approximation method for modeling coagulation kinetics-part 1: Justification and implementation of the method. Journal of Aerosol Science 33, 51-63 (2002).
• Piskunov, V., Golubev, A., Barrett, J. and Ismailova, N., The generalized approximation method for modeling coagulation kinetics-part 2: comparison with other methods. Journal of Aerosol Science 33, 65-75 (2002).
• Ramkrishna, D., Population balances: Theory and applications to particulate systems in engineering. Academic Press (2000).
• Ruckenstein, E. and Pulvermacher, B., Growth kinetics and the size distributions of supported metal crystallites. Journal of Catalysis 29(2), 224-245 (1973).
• Scott, W. T., Analytic studies of cloud droplet coalescence I. Journal of the Atmospheric Sciences 25(1), 54-65 (1968).
• Singh, M., Kumar, J., Bück, A. and Tsotsas, E., A volume-consistent discrete formulation of aggregation population balance equations. Mathematical Methods in the Applied Sciences 39(9), 2275-2286 (2015).
• Smit, D., Hounslow, M. and Paterson, W., Aggregation and gelation I. analytical solutions for cst and batch operation. Chemical Engineering Science 49(7), 1025-1035 (1994).
• Smith, M. and Matsoukas, T., Constant-number monte carlo simulation of population balances. Chemical Engineering Science 53(9), 1777-1786 (1998).
• Vemury, S., Kusters, K. A. and Pratsinis, S. E., Time-lag for attainment of the self-preserving particle size distribution by coagulation. Journal of Colloid and Interface Science, 165(1), 53-59 (1994).
• Ziff, Robert M. and McGrady, E.D., The kinetics of cluster depolymerisation. Journal of Physics A: Mathematical and General 18(15), 3027-3037 (1985).

Publication Dates

• Publication in this collection
  Dec 2018

History

• Received
  26 Jan 2017
• Reviewed
  15 May 2017
• Accepted
  03 July 2017