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

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.


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]  ( 1 , where  1 ,  2 are two characteristic lengths for crystals,  is the crystal size distribution,  1 ,  2 are the growth rate in  1 ,  2 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  1 ,  2 , 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.
Since the PBEs are hyperbolic equations (convectivereactive 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 twodimensional 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 twodimensional 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.

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  =  1 ,  =  2 ,  =  1 , and  =  2 .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 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].
To construct the /| (,) and /| (,) , it is important to use the upwind idea in the following way [21]: 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 10 −6 .

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, Δ 1 , Δ 2 are the space step in  1 direction and  2 direction, respectively.The expression of ( 1 )   + , and ( 2 )  , + is defined as follows: ( 1 )
Table 1 shows the accuracy of the WENO scheme at  = 2.0.It is notable that present WENO scheme works well in this smooth case and the  1 -order and  2 -order are around 5. Therefore, the validity of the numerical method is proved.
Here, the computational domain is taken as [0, 1 m] × [0, 1 m] and the growth rate of crystals is set as  1 =  2 = 0.1 m/s.The analytical solution of this problem is The computational domain is discretized into 100 × 100 cells with Δ 1 = Δ 2 = 0.01.shows some diffusion as the peak is lower than the analytical distribution.Table 2 presents the  1 and  2 errors for HR and WENO method.It can be seen that WENO method provides less errors.
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 [0, 1 m] × [0, 1 m], and the growth terms are taken as  1 =  2 = 0.1 m/s.The analytical solution is the same as Example 1 which is given as (22).Numerical results obtained with grid size Δ 1 = Δ 2 = 0.01 and time interval Δ = 0.01 are shown in Figure 2 at the final time  = 6 s.It can be seen that in this nonsmooth distribution case, WENO method also shows higher accuracy  and less diffusion.Table 3 compares the  1 and  2 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.
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  0 ,   0 are the constants and the initial distribution is where  is the mean size of the distribution and  0 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.
Example 4 (2D batch crystallization: size-independent growth).We now consider a 2D batch crystallization For potassium dihydrogen phosphate (KDP), the crystal shape is shown in Figure 4, where  1 and  2 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  1 and  2 and follow the simple power law function of the supersaturation [4]: where  = ()/ sat () − 1 is the relative supersaturation,  is the solute concentration, and  sat is the saturated solute concentration.The expression for  sat is given as [4]  sat () = 9.3027 × 10 −5  2 − 9.7629 × 10 −5  + 0.2087, where  is the temperature and follows the linear decreasing profile The nucleation rate term is given as [4,17]  0 ( () ,  ()) =      () .
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 where G = ( 1 ,  2 ).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 (1500 × 2700).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.
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  1 and  2 are size-dependent.
Here, the growth rates are given as [17]  1 (

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

Table 2 :
Total errors in the schemes.

Table 3 :
Total errors in the schemes.

Table 5 :
Total errors in the schemes.