Ind. Eng. Chem. Res. 2010, 49, 5911–5924
5911
Comparison of High Resolution Schemes for Solving Population Balances Ala Eldin Bouaswaig* and Sebastian Engell Process Dynamics and Operations Group, Technische UniVersita¨t Dortmund, 44221 Dortmund, Germany
The physical properties of particulate products such as crystals and polymers are strongly dependent on their size distribution which is usually modeled by a population balance equation (PBE). The PBE is an integro-partial differential equation and due to its hyperbolic nature, steep fronts are encountered in the solution profile. When solved numerically, the discretization of the PBE has to be performed with care to avoid numerical diffusion, and dispersion and appropriate numerical methods have to be chosen for that purpose. In this work, an extensive investigation of the performance of several high resolution discretization methods when applied to the PBE is performed with the objective to determine the most promising among them in terms of accuracy and computation time. The comparison is performed on standard benchmark problems and on a practical PBE that arises from modeling the growth of a seed of particles in emulsion polymerization. 1. Introduction Particulate products such as emulsion polymers, crystals, biological cells, and aerosols are frequently encountered in daily life. Their physical properties are directly related to their molecular and particulate characteristics such as the size distribution which is usually modeled by a population balance equation (PBE).1-3 An analytical solution of the PBE only exists for some simplified cases,1 and when an analytical solution is not available, the PBE has to be solved numerically. Because of the hyperbolic nature of the PBE this task is challenging and has to be addressed with care. The numerical simulation of the PBE has attracted numerous researchers in the last half century and is still an active field of research. Among the numerical methods that were used for discretizing the PBE are finite element schemes,4-7 orthogonal collocation on finite elements schemes (OCFE),8-15 and fixed and moving pivot discretization schemes.16,17 More recently, high resolution schemes have gained interest.18-24 For a detailed list of the numerical methods that have been applied to discretizing the PBE the reader is referred to Ramkrishna,1 Vale and McKenna,21 and Kiparissides.25 Among the methods mentioned above, the high resolution finite difference/finite volume methods have the advantage that they were especially designed by researchers in the field of computational fluid dynamics (CFD) to handle hyperbolic conservation laws as encountered for example, in gas dynamics.20 They have proven to be very efficient for that purpose26-29 and also seem to be promising for solving PBEs. For example, Lim et al.18 applied the weighted essentially nonoscillatory (WENO) scheme to a crystallization process. The results of the scheme were compared to those of the modified method of characteristics (MOC), and it was found that the steep moving fronts were well captured by both methods. Motz et al.19 compared the performance of a standard finite volume scheme with piecewise constant profiles, a finite volume scheme extended by a flux limiter, and the space-time conservation element and solution element method (CE/SE) for growth problems. Among these methods the performance of the (CE/ SE) scheme was found to be superior for the investigated case studies, and it was also observed that the flux limited finite * To whom correspondence should be addressed: E-mail:
[email protected].
volume schemes perform better than the standard finite volume schemes. Gunawan et al.20 studied the application of flux limited finite volume algorithms for one-dimensional and multidimensional, batch and continuous population balance equations with nucleation and size-dependent growth. Time integration was performed explicitly, and the Courant-Friedrichs-Levy (CFL) condition was used to determine the time step to guarantee stability. From the results it became clear that the high-resolution methods significantly reduce the numerical diffusion associated with low order methods and avoid spurious oscillations that are observed when higher order methods are used. Furthermore, since a coarser grid can be used for the high-resolution methods to achieve the same accuracy as classical finite volume methods, the computational cost was reduced by orders of magnitude. The algorithms were implemented in the freely available software “ParticleSolver”. Vale and McKenna21 provided a review of the numerical methods that are used for discretizing the PBE. The authors recommended the use of the WENO scheme for discretizing the growth term in the PBE and the approach of Kumar and Ramkrishna30 for handling the aggregation terms. A brief comparison between the WENO scheme and the first order upwind scheme was provided. Qamar et al.22 investigated the use of the flux limited finite volume semidiscrete scheme of Koren31 for discretizing the PBE. Size-independent growth problems with and without stiff nucleation and sizedependent growth problems were studied for crystallization processes. Unlike Gunawan et al.,20 an ODE solver was used to simulate the semidiscrete equations. The performance of the flux limited methods was compared with that of the (CE/SE) method and of the commercial software for particulate processes PARSIVAL. It was found that the high resolution scheme performed very well and provided results comparable to those of the (CE/SE) method and of PARSIVAL. Recently, Qamar and Warnecke published two new numerical methods for the discretization of the PBE. The first method combines MOC for the growth terms with a finite volume (FV) scheme for the aggregation terms. The second method is a flux limited semidiscrete FV scheme.23 The methods were tested for pure growth, growth and aggregation, nucleation and growth, and simultaneous nucleation, growth, and aggregation. For handling the nucleation in the MOC/FV scheme an approach similar to that reported by Kumar and Ramkrishna30 was used. For the case studies investigated, the MOC/FV method
10.1021/ie9020057 2010 American Chemical Society Published on Web 05/10/2010
5912
Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010
generated results with significantly less numerical diffusion in comparison to the results of the FVS and the classical first order FV method. From the previous discussion it can be concluded that the use of high resolution schemes is increasingly attracting interest, but to the best of our knowledge, a comprehensive numerical analysis of the performance of these methods similar to the analysis performed on OCFE11,13 and on finite element methods5-7 has not yet been reported in the open literature. The objective of this paper is to perform a comparison among several high resolution schemes, including a recently published variation of the WENO scheme, known as WENO35Z.32 The comparison presented here is limited to numerical methods that are designed to handle hyperbolic conservation laws, that is, high resolution methods including flux limiter methods, adaptive stencil methods, and weighted stencil methods. The methods investigated in this work are stable, and the goal is to determine the most efficient ones among them in terms of accuracy and computation time. Inspired by the ideas of Vale and McKenna,21 the approach of Kumar and Ramkrishna30 for the aggregation terms is used for all numerical methods and the main difference is in the calculation of the growth term. As far as we are aware, this work is the first in which this combination between the pivot method and high resolution methods (suggested originally by Vale and McKenna21) is thoroughly investigated. Since it is common practice in engineering applications to use the method of lines for solving partial differential equation and off-shelf DAE solvers for solving the semidiscrete set of equations that result after discretization, the performance of all methods is analyzed when coupled with a standard DAE solver. The tests are performed on benchmark problems that comprise pure growth and growth and nucleation as well as growth and coagulation problems. The analytical solution for these problems are mainly taken from the literature; for two more growth and nucleation cases the analytical solutions are derived and presented in this work. In addition to the aforementioned benchmark problems, the numerical methods investigated in this paper are also tested on a practical case that arises from modeling the growth of a seed of polymer particles in emulsion polymerization. Since an analytical solution is not available for this example, the numerical results are compared to a semianalytical solution that is obtained by applying the MOC.33 This article is organized as follows: After an introduction of the PBE and its discretized form in section 2, the mathematical formulation of the numerical methods is presented in section 3. In section 4 the tested benchmark problems are summarized and the results are shown and discussed in section 5. Finally, section 6 gives a summary and conclusions. 2. Discretization of the PBE The PBE that describes the evolution of a particulate system undergoing particle coagulation, growth and nucleation when the volume is taken as the internal coordinate is given by the following hyperbolic PDE: ∂(Gv(V, t)nv(V, t)) ∂nv(V, t) + ) Rnucleation(V, t) + ∂t ∂V Rcoagulation(V, t)
(1)
Consider a grid defined by the points Vi ) i∆V, i ) 0, 1, ..., N, also called the cell centers. As shown in Figure 1, each cell Vi is surrounded by the cell boundaries Vi(1/2 ) Vi (
Figure 1. Stencils of a finite difference/finite volume method. Si(0), Si(1) are the k stencils (in this case k ) 2) that correspond to w ) 0, w ) 1, respectively.
∆V/2. On this grid, the discretized version of the PBE (eq 1) takes the following form: dnv,i ∂f + dt ∂V
|
V)Vi
) Γi
(2)
where f ) Gv(V,t)nv(V,t) and nv,i is the numerical approximation of the value of the density function nV at the point Vi. To maintain the conservative property of the PBE,32 the numerical flux ˆf(V, t) is implicitly defined by f(V, t) )
1 ∆V
∫
V+∆V/2
V-∆V/2
ˆf(V, t) dV
such that the derivative ∂f/∂V is calculated exactly by the following formula: ∂f ∂V
|
V)Vi
)
ˆfi+1/2 - ˆfi-1/2 ∆V
(3)
By substituting eq 3 in eq 2, the discretized form of the PBE (eq 1) is obtained: ˆfi+1/2 - ˆfi-1/2 dnv,i ) Γi + dt ∆V
(4)
Γi ) Rnucleation,i + Rcoagulation,i
(5)
Rcoagulation,i ) Rformation,i - Rdepletion,i
(6)
As described in sections 2.1 and 2.2, the method by which the right-hand side of eq 4 is approximated is identical for all the numerical methods that are compared in this work. However, the approximations of the fluxes ˆfi(1/2 differ. This is explained in detail in section 3. 2.1. Approximation of the Nucleation Term. The nucleation term in eq 5 is given by the following equation:
Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010
Rnucleation(V, t) ) δ(V - Vnuc)Nnuc(t)
(7)
That is, all the nucleated particles/crystals/elements are assumed to be of the size Vnuc. The discretization grid point Vnuc coincides with the first grid cell V0, and eq 7 is approximated as follows:
{
Nnuc(t) i)0 Rnucleation(Vi, t) ) ∆V 0 0
