Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 125128, 11 pages

http://dx.doi.org/10.1155/2013/125128

Research Article

Weighted Essentially Nonoscillatory Method for Two-Dimensional Population Balance Equations in Crystallization

Chunlei Ruan,1,2 Kunfeng Liang,³ Xianjie Chang,² and Ling Zhang²

¹School of Materials Science and Engineering, Zhengzhou University, Zhengzhou 450002, China

²School of Mathematics and Statistics, Henan University of Science and Technology, Luoyang 471023, China

³Vehicle & Motive Power Engineering School, Henan University of Science and Technology, Luoyang 471023, China

Received 8 April 2013; Revised 15 July 2013; Accepted 10 August 2013

Academic Editor: Chuangxia Huang

Copyright © 2013 Chunlei Ruan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Population balance equations (PBEs) are the main governing equations to model the processes of crystallization. Two-dimensional PBEs refer to the crystals that grow anisotropically with the change of two internal coordinates. Since the PBEs are hyperbolic equations, it is necessary to build up high resolution schemes to avoid numerical diffusion and numerical dispersion in order to obtain the accurate crystal size distribution (CSD). In this work, a 5th order weighted essentially nonoscillatory (WENO) method is introduced to compute the two-dimensional PBEs. Several numerical benchmark examples from literatures are carried out; it is found that WENO method has higher resolution than HR method which is well established. Therefore, WENO method is recommended in crystallization simulation when the crystal size distributions are sharp and higher accuracy is needed.

1. Introduction

Crystallization operation has a wide range of applications in the chemical and pharmaceutical industrial production. A lot of solid products exist in the form of crystals or derived by crystallization separation [1]. Crystallization kinetics is the theory foundation of the development of crystallization technology, which provides the main basis for the selection and design of crystallizer as well as the quality control of crystalline products. Population balance equations (PBEs) [1] are the widely accepted equations to model the dynamic crystallization behavior. Therefore, it is of great importance to solve the PBEs accurately and to predict the crystal size distribution (CSD) precisely since it may provide further understanding of crystal growth as well as the development of optimal control strategies for these processes.

Commonly, population balance modeling for crystallization processes only considers one inner variable for which crystals grow isotropically to form simple and isotropic shape (spherical, cubic, etc.) [2]. Obviously, such approaches become questionable in the case the crystallization of organic products presenting anisotropic morphologies. In this work, we consider the transient behavior of crystals in a batch process with two-dimensional growth and nucleation; therefore, a two-dimensional population balance model needs to be employed which can be written as [3, 4]

where are two characteristic lengths for crystals, is the crystal size distribution, are the growth rate in direction, is a time-dependent solute concentration, is a time-dependent temperature, and is the creation/depletion rate and may include aggregation, nucleation, and breakage. It should be mentioned that here we assume the crystallizer is well mixed and the crystal size distribution is considered to be independent of the spatial coordinates. PBE (1) represents many of the crystals in the pharmaceuticals, photographic, and other industries for which their growth is associated with the change of two internal coordinates.

Since the PBEs are strongly nonlinear and the expression of , , and is sometimes very complex, it is hard to find the analytical solutions in the real crystallization operation. That is why numerical simulation becomes popular in solving PBEs.

So far, there has been plenty of numerical works on solving PBEs, and the available methods can be mainly divided into the following five categories [4]: (1) the method of moments [3, 5]; (2) the method of characteristics [6, 7]; (3) the method of weight residuals/orthogonal collocation [8]; (4) the Monte Carlo method [9, 10], and (5) the finite difference schemes/discrete population balances [11, 12]. The brief discussion and references about these methods can refer to the works [4, 13, 14].

Since the PBEs are hyperbolic equations (convective-reactive equations), numerical diffusion and dispersion are caused when ordinary schemes of discrete methods are used. That is why high resolution schemes are needed. The high resolution schemes, which were originally developed for gas dynamics, are now widely used in solving PBEs. Puel et al. [2] and Ma and Wang [15] used the method of classes to solve PBEs; they reported that the advantage of this method was the accurate simulation of the nucleation and the main drawback was the high computational time. Ma et al. [4, 16] and Gunawan et al. [17] adapted the high resolution of finite volume method (HR schemes) to solve PBEs, and different schemes were compared by Qamar et al. [18]. HR schemes are the most popular method for solving PBEs, in which different flux limiter functions were constructed to reduce the numerical diffusion and to capture the sharp peaks and discontinuities of the crystal size distribution. Based on the high resolution of finite volume method, moving mesh technique was introduced by Qamar et al. [14]. They showed the accuracy and efficiency of these adaptive schemes in resolving sharp peaks of crystal size distribution. Majumder et al. [13] developed Lattice Boltzmann method (LBM) for solving PBEs. They reported that LBM was at least as accurate as HR method but cheaper in computational cost. Hermanto et al. [19] used weighted essentially nonoscillatory (WENO) method to solve one-dimensional PBEs; they found that WENO method provided much higher order accuracy than other previously considered methods for PBEs.

Since the crystal size distribution in crystallization can be very sharp and sometimes discontinue, higher accurate schemes are needed especially when higher accurate distribution is desired. In this paper, WENO method is introduced to solve the two-dimensional PBEs. The higher accuracy of this method has been numerical proofed by Hermanto et al. [19]. However, they only deal with one-dimensional case. The two-dimensional PBEs are more complex but represent a large class of crystals that grow anisotropically.

This paper is organized as follows: in Section 2, the detailed procedure of WENO method for solving two-dimensional PBEs is presented; in Section 3, several numerical test problems are given to show the advantage of using WENO method; and finally the conclusions are given in Section 4.

2. WENO Method

WENO method for simulation fluid dynamics has received significant attention in the past two decades [20]. This method is developed based on the essentially nonoscillatory (ENO) method. It improved the accuracy of the original ENO method to the optimal order in smooth regions while maintaining the essentially nonoscillatory property near discontinuities [21]. One of the main advantages of WENO is its high accuracy. In this section, we present the detailed implement of WENO method in solving PBEs.

The two-dimensional PBE (1) can be notationally simplified by

where , , , and . The computational domain is assumed as , and the cells are expressed as , , , where is the total cell number in direction and is the total cell number in direction. Set the center of cell as , then , . In this work, uniform meshes are applied, which means, the internal mesh sizes are and .

To construct the WENO scheme, there are two different views, one is based on finite volume version and the other is based on finite difference version in space. Here, we use the finite difference WENO schemes constructed by Jiang and Shu [21].

In the finite difference WENO schemes, the computational value on the grid nodes (viz., ) is known. Equation (2) is discretized in the conservative form in space and leads to

Here, is approximate to and is approximate to , with the numerical flux on the cell boundaries.

To construct the and , it is important to use the upwind idea in the following way [21]:

To find , set

and to find , set

Defining the smoothness [21]

and the weights [21]

the 5th order WENO scheme in direction is finally written as [21]

Here, is a small positive number to prevent the denominators from becoming zero. In this work, is set to .

Similarly, to find , set

and to find , set

Defining the smoothness

and the weights

the 5th order WENO scheme in direction is written as

Equations (9) and (14) are 5th order accuracy in space. Take the , for example, it covers six cells: , , , , , and and three stencils: , , and . is a 5th order approximation to which provides the smallest truncation error on these six-cell stencils. Furthermore, it is a convex combination of , , and . If one of the stencils is chosen according to certain conditions, it is reduced to ENO scheme which only provides third order accuracy [21].

For time discretization, we use the third order TVD (total variation nonincreasing) Runge-Kutta schemes [21]. Equation (2) can be written as

where . The third order Runge-Kutta schemes can be expressed as

with and .

3. Numerical Examples

In order to show the high resolution of the WENO method, we consider several test problems for PBEs in crystallization. We compare the solutions of WENO method with HR method, which is a well-established high resolution method. The corresponding error norms are defined as

where is the exact solution.

The HR scheme for solving PBE (1) is written as [4, 17, 18]

when is the time step, are the space step in direction and direction, respectively. The expression of and is defined as follows:

where , , and .

3.1. Validity of the Simulation

In order to show that the 5th order WENO scheme has been implemented correctly, we now consider the following two-dimensional linear equation:

With the initial condition , , and the boundary conditions are assumed to be periodic. The analytical solution of this problem is .

Table 1 shows the accuracy of the WENO scheme at . It is notable that present WENO scheme works well in this smooth case and the -order and -order are around 5. Therefore, the validity of the numerical method is proved.

Table 1: Accuracy of WENO scheme in 2D linear equation.

-error
-order
-error
-order

2/5
0.479449
—
0.127868
—

1/5
0.0241227
4.31
0.00682424
4.23

1/10
0.000739797
5.03
0.000215145
4.99

1/20

5.08

5.11

1/40

5.10

4.94

Example 1 (size-independent growth (smooth distribution)). We first consider a homogenous PBE with size-independent growth which is expressed as

with the initial distribution set as . This is a Gaussian-type distribution which is very usual in the real crystallization processes. Here, the computational domain is taken as and the growth rate of crystals is set as . The analytical solution of this problem is

The computational domain is discretized into cells with . Figure 1 gives a comparison of the simulated CSDs with the analytical CSDs at . It is obvious that the result obtained by WENO agrees with the analytical distribution very well, while the HR method shows some diffusion as the peak is lower than the analytical distribution. Table 2 presents the and errors for HR and WENO method. It can be seen that WENO method provides less errors.

Table 2: Total errors in the schemes.

Method
HR
WENO

-error
0.00308458

-error
0.0107827
0.000131731

(a)

(b)

(c)

Figure 1: Final distribution comparison for Example 1: (a) exact, (b) HR, and (c) WENO.

Example 2 (size-independent growth (nonsmooth distribution)). Since most of initial particles distributions in the real processes can be very sharp, it is necessary to analyze the performance of WENO for nonsmooth distribution. The homogenous PBE (21) is considered here, and the initial data for CSDs is given by

The computational domain is taken as , and the growth terms are taken as . The analytical solution is the same as Example 1 which is given as (22). Numerical results obtained with grid size and time interval are shown in Figure 2 at the final time . It can be seen that in this nonsmooth distribution case, WENO method also shows higher accuracy and less diffusion. Table 3 compares the and errors for HR and WENO method. One can clearly see that WENO method provides less error. It is satisfied that WENO method does not show oscillation for capturing the sharp front. This is very important since the real crystallization processes always involve propagation of sharp fronts occurring due to nucleation.

Table 3: Total errors in the schemes.

Method
HR
WENO

-error
0.0115517
0.00635183

-error
0.0626252
0.0423699

(a)

(b)

(c)

Figure 2: Final distribution comparison for Example 2: (a) exact, (b) HR, and (c) WENO.

Example 3 (size-dependent growth). Next, we consider a homogenous PBE with size-dependent growth which is expressed as follows:

Here, we consider that the growth rate is dependent linearly on size

where , are the constants and the initial distribution is

where is the mean size of the distribution and is the total initial number of particles. The analytical solution is

Parameters in this problem are listed in Table 4. Figure 3 gives the comparison for final CSDs with different methods. Again results of WENO method have better resolution. This is also proofed with the corresponding errors which are listed in Table 5.

Table 4: Parameters for Example 3.

Description
Symbols
Value
Unit

Total -size distribution

0.1
μm

Total -size distribution

0.1
μm

Mesh size

10⁻³
μm

Mesh size

10⁻³
μm

Mean size of the distribution

10⁻²
μm

Total initial number of particles

1
Dimensionless

Growth rate

0.1
μm/s

Growth rate

0.2
μm/s

Simulation time

4.0
s

Table 5: Total errors in the schemes.

Method
HR
WENO

-error
0.000550266
0.000135244

-error
0.0195435
0.00982604

(a)

(b)

(c)

Figure 3: Final distribution comparison for Example 3: (a) exact, (b) HR, and (c) WENO.

Example 4 (2D batch crystallization: size-independent growth). We now consider a 2D batch crystallization example. Assume that the nucleation takes place at negligible crystal size; therefore, the corresponding PBE can be written as [4]

For potassium dihydrogen phosphate (KDP), the crystal shape is shown in Figure 4, where and are the width and length of the crystal, respectively. It is obvious that the volume of a single crystal is [4]

The growth rates are assumed to be independent of and and follow the simple power law function of the supersaturation [4]:

where is the relative supersaturation, is the solute concentration, and is the saturated solute concentration. The expression for is given as [4]

where is the temperature and follows the linear decreasing profile

The nucleation rate term is given as [4, 17]

Here, is generally assumed to be related to the stirring power, and is the total volume of the crystal [1]:

The mass balance of a batch system for the solute of which the crystalline phase is being assembled is [1]

where . For KDP, this equation is

The simulated initial distribution is given as [4]

Parameters in this problem are listed in Table 6. Since this problem has no analytical solution, the “exact” solution was obtained by WENO scheme on fine grid (). Figure 5 gives the final CSDs comparison. It is obvious that both HR method and WENO method exhibit diffusion on capturing the size distribution with crystals growing from nuclei. In addition, the performance of these two methods on predicting the distributions associated with crystals growing from seeds is similar. By comparison, WENO method shows better accuracy than HR method in this case.

Table 6: Parameters for Example 4.

Description
Symbols
Value
Unit

Total -size distribution

500
μm

Total -size distribution

900
μm

Mesh size

2
μm

Mesh size

2
μm

-Growth exponent

1.48
Dimensionless

-Growth rate constant

12.21
μm/s

y-Growth exponent

1.74
Dimensionless

y-Growth rate constant

100.75
μm/s

Nucleation rate constant

1/μm³s

Nucleation exponent

2.04
Dimensionless

Density of crystals

g/μm³

Figure 4: Characteristic lengths of KDP crystals.

(a)

(b)

(c)

Figure 5: Final distribution comparison for Example 4: (a) exact, (b) HR, and (c) WENO.

Example 5 (2D batch crystallization: size-dependent growth). We now consider a 2D batch crystallization example with growth rate dependent on characteristic lengths. KDP crystals are also used. The PBE is similar as (28) in Example 4 except for growth rates of and are size-dependent.

Here, the growth rates are given as [17]

and the temperature is assumed to follow the expression [17]

The saturated solute concentration , nucleation rate , and solute concentration are the same as Example 4. The kinetic parameters of , , , , , , and are also the same as Example 4 which are listed in Table 6.

The initial distribution is taken as [14]

and the computational zone is assumed as .

The final crystal size distribution is shown in Figure 6. Also, the “exact” solution is obtained by WENO method on fine grid (), while the other results are obtained with grid size . From Figure 6, it is clear that WENO method has better resolution.

(a)

(b)

(c)

Figure 6: Final distribution comparison for Example 5: (a) exact, (b) HR, and (c) WENO.

4. Conclusion

This article presents the application of WENO method for solving two-dimensional PBEs in crystallization processes. Numerical examples show that WENO method has higher resolution than HR method including the cases of PBEs with nucleation rate and size dependent/independent growth rate. Therefore, WENO method is recommended in crystallization simulation when the crystal size distributions are sharp and higher accuracy is desired.

Acknowledgments

The financial supports provided by the Postdoctoral Science Foundation of China (2012M521403), Henan Scientific and Technological Research Project (no. 122102210198), and Key Scientific and Technological Research Project of Department of Education of Henan Province (nos. 2011B110014 and 12B110006) are fully acknowledged.

References

C. Borchert and K. Sundmacher, “Morphology evolution of crystal populations: modeling and observation analysis,” Chemical Engineering Science, vol. 70, pp. 87–98, 2012.

F. Puel, G. Févotte, and J. P. Klein, “Simulation and analysis of industrial crystallization processes through multidimensional population balance equations. Part 1: a resolution algorithm based on the method of classes,” Chemical Engineering Science, vol. 58, no. 16, pp. 3715–3727, 2003.

H. Briesen, “Simulation of crystal size and shape by means of a reduced two-dimensional population balance model,” Chemical Engineering Science, vol. 61, no. 1, pp. 104–112, 2006.

D. L. Ma, D. K. Tafti, and R. D. Braatz, “High-resolution simulation of multidimensional crystal growth,” Industrial and Engineering Chemistry Research, vol. 41, no. 25, pp. 6217–6223, 2002.

C. B. B. Costa, M. R. W. Maciel, and R. M. Filho, “Considerations on the crystallization modeling: population balance solution,” Computers and Chemical Engineering, vol. 31, no. 3, pp. 206–218, 2007.

S. Qamar, I. Angelov, M. P. Elsner, A. Ashfaq, A. Seidel-Morgenstern, and G. Warnecke, “Numerical approximations of a population balance model for coupled batch preferential crystallizers,” Applied Numerical Mathematics, vol. 59, no. 3-4, pp. 739–753, 2009.

S. Qamar, A. Ashfaq, I. A. Angelov, M. P. Elsner, G. Warnecke, and A. Seidel-Morgenstern, “Numerical solutions of population balance models in preferential crystallization,” Chemical Engineering Science, vol. 63, no. 5, pp. 1342–1352, 2008.

M. Nicmanis and M. J. Hounslow, “Finite-element methods for steady-state population balance equations,” AIChE Journal, vol. 44, no. 10, pp. 2258–2272, 1998.

D. Meimaroglou and C. Kiparissides, “Monte Carlo simulation for the solution of the bi-variate dynamic population balance equation in batch particulate systems,” Chemical Engineering Science, vol. 62, no. 18–20, pp. 5295–5299, 2007.

H. Zhao, A. Maisels, T. Matsoukas, and C. Zheng, “Analysis of four Monte Carlo methods for the solution of population balances in dispersed systems,” Powder Technology, vol. 173, no. 1, pp. 38–50, 2007.

M. J. Hounslow, “Discretized population balance for continuous systems at steady state,” AIChE Journal, vol. 36, no. 1, pp. 106–116, 1990.

S. Kumar and D. Ramkrishna, “On the solution of population balance equations by discretization—I. A fixed pivot technique,” Chemical Engineering Science, vol. 51, no. 8, pp. 1311–1332, 1996.

A. Majumder, V. Kariwala, S. Ansumali, and A. Rajendran, “Lattice Boltzmann method for multi-dimensional population balance models in crystallization,” Chemical Engineering Science, vol. 70, pp. 121–134, 2012.

S. Qamar, A. Ashfaq, G. Warnecke, I. Angelov, M. P. Elsner, and A. Seidel-Morgenstern, “Adaptive high-resolution schemes for multidimensional population balances in crystallization processes,” Computers and Chemical Engineering, vol. 31, no. 10, pp. 1296–1311, 2007.

C. Y. Ma and X. Z. Wang, “Crystal growth rate dispersion modeling using morphological population balance,” AIChE Journal, vol. 54, no. 9, pp. 2321–2334, 2008.

D. L. Ma, D. K. Tafti, and R. D. Braatz, “Optimal control and simulation of multidimensional crystallization processes,” Computers and Chemical Engineering, vol. 26, no. 7-8, pp. 1103–1116, 2002.

R. Gunawan, I. Fusman, and R. D. Braatz, “High resolution algorithms for multidimensional population balance equations,” AIChE Journal, vol. 50, no. 11, pp. 2738–2749, 2004.

S. Qamar, M. P. Elsner, I. A. Angelov, G. Warnecke, and A. Seidel-Morgenstern, “A comparative study of high resolution schemes for solving population balances in crystallization,” Computers and Chemical Engineering, vol. 30, no. 6-7, pp. 1119–1131, 2006.

M. W. Hermanto, R. D. Braatz, and M.-S. Chiu, “High-order simulation of polymorphic crystallization using weighted essentially nonoscillatory methods,” AIChE Journal, vol. 55, no. 1, pp. 122–131, 2009.

Y. Xing and C.-W. Shu, “High order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms,” Journal of Computational Physics, vol. 214, no. 2, pp. 567–598, 2006.

G.-S. Jiang and C.-W. Shu, “Efficient implementation of weighted ENO schemes,” Journal of Computational Physics, vol. 126, no. 1, pp. 202–228, 1996.