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

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] 𝜕𝑓𝑟₁,𝑟₂,𝑡+𝜕𝑡²_𝑗=1𝜕𝐺_𝑗𝑟₁,𝑟₂,𝑐(𝑡),𝑇(𝑡)𝑓𝑟₁,𝑟₂,𝑡𝜕𝑟_𝑗=ℎ𝑓𝑟₁,𝑟₂,𝑡,𝑐(𝑡),𝑇(𝑡), (1) 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 𝜕𝑓+𝑓𝜕𝑡𝜕𝑝+𝑓𝜕𝑥𝜕𝑞𝜕𝑦=ℎ, (2) where 𝑝=𝐺₁𝑓, 𝑞=𝐺₂𝑓, 𝑥=𝑟₁, and 𝑦=𝑟₂. The computational domain is assumed as [𝑎,𝑏]×[𝑐,𝑑], and the cells are expressed as 𝐼_𝑖𝑗=[𝑥_𝑖−1/2,𝑥_𝑖+1/2]×[𝑦_𝑗−1/2,𝑦_𝑗+1/2], 1≤𝑖≤𝑁_𝑥, 1≤𝑗≤𝑁_𝑦, where 𝑁_𝑥 is the total cell number in 𝑥 direction and 𝑁_𝑦 is the total cell number in 𝑦 direction. Set the center of cell 𝐼_𝑖𝑗 as (𝑥_𝑖,𝑦_𝑗), then 𝑥_𝑖=(1/2)(𝑥_𝑖−1/2+𝑥_𝑖+1/2), 𝑦_𝑖=(1/2)(𝑦_𝑗−1/2+𝑦_𝑗+1/2). In this work, uniform meshes are applied, which means, the internal mesh sizes are Δ𝑥=𝑥_𝑖+1/2−𝑥_𝑖−1/2 and Δ𝑦=𝑦_𝑗+1/2−𝑦_𝑗−1/2.

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 𝑑𝑓𝑑𝑡^𝑖,𝑗+1𝑝Δ𝑥^{𝑖+1/2,𝑗}−𝑝_{𝑖−1/2,𝑗}+1𝑞Δ𝑦^{𝑖,𝑗+1/2}−𝑞_{𝑖,𝑗−1/2}=ℎ_𝑖,𝑗. (3) Here, (𝑝_{𝑖+1/2,𝑗}−𝑝_{𝑖−1/2,𝑗})/Δ𝑥 is approximate to 𝜕𝑝/𝜕𝑥|_(𝑖,𝑗) and (𝑞_{𝑖,𝑗+1/2}−𝑞_{𝑖,𝑗−1/2})/Δ𝑦 is approximate to 𝜕𝑞/𝜕𝑦|_(𝑖,𝑗), with 𝑝_{𝑖+1/2,𝑗},𝑝_{𝑖−1/2,𝑗},𝑞_{𝑖,𝑗+1/2},𝑞_{𝑖,𝑗−1/2} the numerical flux on the cell boundaries.

To construct the 𝜕𝑝/𝜕𝑥|_(𝑖,𝑗) and 𝜕𝑞/𝜕𝑦|_(𝑖,𝑗), it is important to use the upwind idea in the following way [21]: 𝜕𝑝𝜕𝑥_(𝑖,𝑗)=𝑝_𝑥⁻_𝑖,𝑗𝐺_1𝑖𝑝≥0_{𝑥+𝑖,𝑗}𝐺_1𝑖<0,𝜕𝑞𝜕𝑦_(𝑖,𝑗)=𝑞_𝑦⁻_𝑖,𝑗𝐺_2𝑗𝑞≥0_{𝑦+𝑖,𝑗}𝐺_2𝑗<0. (4) To find (𝑝_𝑥)⁻_𝑖,𝑗, set 𝑢₁=𝑝_{𝑖−2,𝑗}−𝑝_{𝑖−3,𝑗}Δ𝑥,𝑢₂=𝑝_{𝑖−1,𝑗}−𝑝_{𝑖−2,𝑗},𝑢Δ𝑥₃=𝑝_𝑖,𝑗−𝑝_{𝑖−1,𝑗}Δ𝑥,𝑢₄=𝑝_𝑖+1,𝑗−𝑝_𝑖,𝑗,𝑢Δ𝑥₅=𝑝_𝑖+2,𝑗−𝑝_𝑖+1,𝑗,Δ𝑥 (5) and to find (𝑝_𝑥)⁺_𝑖,𝑗, set 𝑢₁=𝑝_𝑖+3,𝑗−𝑝_𝑖+2,𝑗Δ𝑥,𝑢₂=𝑝_𝑖+2,𝑗−𝑝_𝑖+1,𝑗,𝑢Δ𝑥₃=𝑝_𝑖+1,𝑗−𝑝_𝑖,𝑗Δ𝑥,𝑢₄=𝑝_𝑖,𝑗−𝑝_{𝑖−1,𝑗},𝑢Δ𝑥₅=𝑝_{𝑖−1,𝑗}−𝑝_{𝑖−2,𝑗}.Δ𝑥 (6) Defining the smoothness [21] 𝐼𝑆₁=13𝑢12¹−2𝑢₂+𝑢₃²+14𝑢₁−4𝑢₂+3𝑢₃²,𝐼𝑆₂=13𝑢12²−2𝑢₃+𝑢₄²+14𝑢₂−𝑢₄²,𝐼𝑆₃=13𝑢12³−2𝑢₄+𝑢₅²+143𝑢₃−4𝑢₄+𝑢₅² (7) and the weights [21] 𝐼𝜔₁=𝑎₁𝑎₁+𝑎₂+𝑎₃,𝐼𝜔₂=𝑎₂𝑎₁+𝑎₂+𝑎₃,𝐼𝜔₃=𝑎₃𝑎₁+𝑎₂+𝑎₃,𝑎₁=1110𝐼𝑆₁+𝜀²,𝑎₂=6110𝐼𝑆₂+𝜀²,𝑎₃=3110𝐼𝑆₃+𝜀², (8) the 5th order WENO scheme in 𝑥 direction is finally written as [21] 𝑝_𝑥^±_𝑖,𝑗=𝐼𝜔₁𝑢₁3−7𝑢₂6+11𝑢₃6+𝐼𝜔₂−𝑢₂6+5𝑢₃6+𝑢₄3+𝐼𝜔₃𝑢₃3+5𝑢₄6−𝑢₅6. (9) Here, 𝜀 is a small positive number to prevent the denominators from becoming zero. In this work, 𝜀 is set to 10⁻⁶.

Similarly, to find (𝑞_𝑦)⁻_𝑖,𝑗, set 𝑣₁=𝑞_{𝑖,𝑗−2}−𝑞_{𝑖,𝑗−3}Δ𝑦,𝑣₂=𝑞_{𝑖,𝑗−1}−𝑞_{𝑖,𝑗−2},𝑣Δ𝑦₃=𝑞_𝑖,𝑗−𝑞_{𝑖,𝑗−1}Δ𝑦,𝑣₄=𝑞_𝑖,𝑗+1−𝑞_𝑖,𝑗,𝑣Δ𝑦₅=𝑞_𝑖,𝑗+2−𝑞_𝑖,𝑗+1,Δ𝑦 (10) and to find (𝑞_𝑦)⁺_𝑖,𝑗, set 𝑣₁=𝑞_𝑖,𝑗+3−𝑞_𝑖,𝑗+2Δ𝑦,𝑣₂=𝑞_𝑖,𝑗+2−𝑞_𝑖,𝑗+1,𝑣Δ𝑦₃=𝑞_𝑖,𝑗+1−𝑞_𝑖,𝑗Δ𝑦,𝑣₄=𝑞_𝑖,𝑗−𝑞_{𝑖,𝑗−1},𝑣Δ𝑦₅=𝑞_{𝑖,𝑗−1}−𝑞_{𝑖,𝑗−2}.Δ𝑦 (11) Defining the smoothness 𝐽𝑆₁=13𝑣12¹−2𝑣₂+𝑣₃²+14𝑣₁−4𝑣₂+3𝑣₃²,𝐽𝑆₂=13𝑣12²−2𝑣₃+𝑣₄²+14𝑣₂−𝑣₄²,𝐽𝑆₃=13𝑣12³−2𝑣₄+𝑣₅²+143𝑣₃−4𝑣₄+𝑣₅² (12) and the weights 𝐽𝜔₁=𝑎₁𝑎₁+𝑎₂+𝑎₃,𝐽𝜔₂=𝑎₂𝑎₁+𝑎₂+𝑎₃,𝐽𝜔₃=𝑎₃𝑎₁+𝑎₂+𝑎₃,𝑎₁=1110𝐽𝑆₁+𝜀²,𝑎₂=6110𝐽𝑆₂+𝜀²,𝑎₃=3110𝐽𝑆₃+𝜀², (13) the 5th order WENO scheme in 𝑦 direction is written as 𝑞_𝑦^±_𝑖,𝑗=𝐽𝜔₁𝑣₁3−7𝑣₂6+11𝑣₃6+𝐽𝜔₂−𝑣₂6+5𝑣₃6+𝑣₄3+𝐽𝜔₃𝑣₃3+5𝑣₄6−𝑣₅6. (14)

Equations (9) and (14) are 5th order accuracy in space. Take the (𝑝_𝑥)⁻_𝑖,𝑗, for example, it covers six cells: 𝐼_{𝑖−3,𝑗}, 𝐼_{𝑖−2,𝑗}, 𝐼_{𝑖−1,𝑗}, 𝐼_𝑖,𝑗, 𝐼_𝑖+1,𝑗, and 𝐼_𝑖+2,𝑗 and three stencils: 𝑆⁰_𝑖,𝑗={𝐼_{𝑖−3,𝑗},𝐼_{𝑖−2,𝑗},𝐼_{𝑖−1,𝑗},𝐼_𝑖,𝑗}, 𝑆¹_𝑖,𝑗={𝐼_{𝑖−2,𝑗},𝐼_{𝑖−1,𝑗},𝐼_𝑖,𝑗,𝐼_𝑖+1,𝑗}, and 𝑆²_𝑖,𝑗={𝐼_{𝑖−1,𝑗},𝐼_𝑖,𝑗,𝐼_𝑖+1,𝑗,𝐼_𝑖+2,𝑗}. (𝑝_𝑥)⁻_𝑖,𝑗 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 𝑑𝑓𝑓,𝑑𝑡=𝐿 (15) where 𝐿(𝑓)=ℎ−(𝜕𝑝/𝜕𝑥)−(𝜕𝑞/𝜕𝑦). The third order Runge-Kutta schemes can be expressed as 𝑓⁽¹⁾=𝑓⁽⁰⁾+Δ𝑡𝐿𝑓⁽⁰⁾,𝑓⁽²⁾=𝑓⁽¹⁾+Δ𝑡4−3𝐿𝑓⁽⁰⁾+𝐿𝑓⁽¹⁾,𝑓⁽³⁾=𝑓⁽²⁾+Δ𝑡𝑓12−𝐿⁽⁰⁾−𝐿𝑓⁽¹⁾+8𝐿𝑓⁽²⁾, (16) with 𝑓⁽⁰⁾=𝑓^𝑛 and 𝑓⁽³⁾=𝑓^𝑛+1.

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 𝐿₁-error=^𝑁𝑀_𝑖=1𝑗=1𝑓^𝑒_𝑖,𝑗−𝑓_𝑖,𝑗Δ𝑟₁Δ𝑟₂,𝐿₂-error=^𝑁𝑀_𝑖=1𝑗=1𝑓^𝑒_𝑖,𝑗−𝑓_𝑖,𝑗²Δ𝑟₁Δ𝑟₂, (17) where 𝑓^𝑒_𝑖,𝑗 is the exact solution.

The HR scheme for solving PBE (1) is written as [4, 17, 18] 𝑓^𝑛+1_𝑖,𝑗=𝑓^𝑛_𝑖,𝑗×ℎ+Δ𝑡_{𝑛𝑖,𝑗}−𝐺₁𝑓^𝑛_𝑖+,𝑗−𝐺₁𝑓^𝑛_𝑖−,𝑗Δ𝑟₁−𝐺₂𝑓^𝑛_𝑖,𝑗+−𝐺₂𝑓^𝑛_𝑖,𝑗−Δ𝑟₂ (18) when Δ𝑡 is the time step, Δ𝑟₁,Δ𝑟₂ are the space step in 𝑟₁ direction and 𝑟₂ direction, respectively. The expression of (𝐺₁𝑓)^𝑛_𝑖+,𝑗 and (𝐺₂𝑓)^𝑛_𝑖,𝑗+ is defined as follows: 𝐺₁𝑓^𝑛_𝑖+,𝑗=𝐺₁^𝑛_𝑖+,𝑗𝑓^𝑛_𝑖,𝑗+12𝜙𝑟^𝑛_𝑖⁺𝑓^𝑛_𝑖+1,𝑗−𝑓^𝑛_𝑖,𝑗,𝐺₂𝑓^𝑛_𝑖,𝑗+=𝐺₂^𝑛_𝑖,𝑗+𝑓^𝑛_𝑖,𝑗+12𝜙𝑟^𝑛_𝑗⁺𝑓^𝑛_𝑖,𝑗+1−𝑓^𝑛_𝑖,𝑗, (19) where 𝑟^𝑛_𝑖⁺=(𝑓^𝑛_𝑖,𝑗−𝑓^𝑛_{𝑖−1,𝑗})/(𝑓^𝑛_𝑖+1,𝑗−𝑓^𝑛_𝑖,𝑗), 𝑟^𝑛_𝑗⁺=(𝑓^𝑛_𝑖,𝑗−𝑓^𝑛_{𝑖,𝑗−1})/(𝑓^𝑛_𝑖,𝑗+1−𝑓^𝑛_𝑖,𝑗), and 𝜙(𝑟)=(|𝑟|+𝑟)/(1+|𝑟|).

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: 𝜕𝑓𝑟₁,𝑟₂,𝑡+𝜕𝑡𝜕𝑓𝜕𝑟₁+𝜕𝑓𝜕𝑟₂=0. (20) With the initial condition 𝑓(𝑟₁,𝑟₂,0)=sin((𝜋/2)(𝑟₁+𝑟₂)), −2≤𝑟₁≤2,−2≤𝑟₂≤2, and the boundary conditions are assumed to be periodic. The analytical solution of this problem is 𝑓(𝑟₁,𝑟₂,𝑡)=sin((𝜋/2)(𝑟₁+𝑟₂−2𝑡)).

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

Example 1 (size-independent growth (smooth distribution)). We first consider a homogenous PBE with size-independent growth which is expressed as 𝜕𝑓𝑟₁,𝑟₂,𝑡+𝜕𝐺𝜕𝑡¹𝑓𝜕𝑟₁+𝜕𝐺₂𝑓𝜕𝑟₂=0, (21) with the initial distribution set as 𝑓(𝑟₁,𝑟₂,0)=exp{−100×((𝑟₁−0.25)²+(𝑟₂−0.25)²)}. This is a Gaussian-type distribution which is very usual in the real crystallization processes. Here, the computational domain is taken as [0,1𝜇m]×[0,1𝜇m] and the growth rate of crystals is set as 𝐺₁=𝐺₂=0.1𝜇m/s. The analytical solution of this problem is 𝑓𝑟₁,𝑟₂,𝑡=𝑓𝑟₁−𝐺₁𝑡,𝑟₂−𝐺₂𝑡,0. (22) The computational domain is discretized into 100×100 cells with Δ𝑟₁=Δ𝑟₂=0.01. Figure 1 gives a comparison of the simulated CSDs with the analytical CSDs at 𝑡=5s. 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.

(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 𝑓𝑟₁,𝑟₂,0=10.1≤𝑟₁,𝑟₂≤0.30elsewhere. (23) The computational domain is taken as [0,1𝜇m]×[0,1𝜇m], and the growth terms are taken as 𝐺₁=𝐺₂=0.1𝜇m/s. The analytical solution is the same as Example 1 which is given as (22). Numerical results obtained with grid size Δ𝑟₁=Δ𝑟₂=0.01 and time interval Δ𝑡=0.01 are shown in Figure 2 at the final time 𝑡=6s. 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.

(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: 𝜕𝑓𝑟₁,𝑟₂,𝑡+𝜕𝐺𝜕𝑡¹𝑟₁𝑓𝜕𝑟₁+𝜕𝐺₂𝑟₂𝑓𝜕𝑟₂=0. (24) Here, we consider that the growth rate is dependent linearly on size 𝐺₁𝑟₁=𝐺₀𝑟₁𝐺₂𝑟₂=𝐺^′₀𝑟₂, (25) where 𝐺₀, 𝐺^′₀ are the constants and the initial distribution is 𝑓𝑟₁,𝑟₂=𝑁,0⁰𝐿−𝑟exp¹+𝑟₂𝐿, (26) where 𝐿 is the mean size of the distribution and 𝑁₀ is the total initial number of particles. The analytical solution is 𝑓𝑟₁,𝑟₂=𝑁,𝑡⁰𝐿−𝑟exp¹𝐿𝑒^−𝐺0𝑡−𝐺₀𝑟𝑡−²𝐿𝑒^{−𝐺′0𝑡}−𝐺^′₀𝑡. (27)
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.

Table 5 

Total errors in the schemes.

(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] 𝜕𝑓𝑟₁,𝑟₂,𝑡+𝜕𝐺𝜕𝑡¹(𝑐,𝑇)𝑓𝜕𝑟₁+𝜕𝐺₂(𝑐,𝑇)𝑓𝜕𝑟₂=𝐵₀(𝑐,𝑇)𝛿𝑟₁𝛿𝑟₂. (28) 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] 𝑉_𝑐=13𝑟³₁+𝑟₂−𝑟₁𝑟²₁. (29)
The growth rates are assumed to be independent of 𝑟₁ and 𝑟₂ and follow the simple power law function of the supersaturation [4]: 𝐺₁(𝑐(𝑡),𝑇(𝑡))=𝑘_𝑔1𝑆^𝑔1𝐺₂(𝑐(𝑡),𝑇(𝑡))=𝑘_𝑔2𝑆^𝑔2, (30) 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⁻⁵𝑇²−9.7629×10⁻⁵𝑇+0.2087, (31) where 𝑇 is the temperature and follows the linear decreasing profile 𝑇(𝑡)^∘C𝑡=33−.720 (32)
The nucleation rate term is given as [4, 17] 𝐵₀(𝑐(𝑡),𝑇(𝑡))=𝑘_𝑏𝑆^𝑏𝑉(𝑡). (33) Here, 𝑏 is generally assumed to be related to the stirring power, and 𝑉(𝑡) is the total volume of the crystal [1]: 𝑉(𝑡)=^∞₀𝑓𝑟₁,𝑟₂,𝑡𝑉_𝑐𝑟₁,𝑟₂𝑑𝑟₁𝑑𝑟₂. (34) The mass balance of a batch system for the solute of which the crystalline phase is being assembled is [1] 𝑑𝑐(𝑡)𝑑𝑡=−𝜌_𝑐^∞₀𝑓𝑟₁,𝑟₂,𝑡×𝐆⋅∇𝑉_𝑐𝑑𝑟₁𝑑𝑟₂, (35) where 𝐆=(𝐺₁,𝐺₂). For KDP, this equation is 𝑑𝑐(𝑡)𝑑𝑡=−𝜌_𝑐^∞₀𝑓𝑟₁,𝑟₂×,𝑡2𝐺₁𝑟₁𝑟₂−𝑟²₁+𝐺₂𝑟²₁𝑑𝑟₁𝑑𝑟₂. (36)
The simulated initial distribution is given as [4] 𝑓𝑟₁,𝑟₂,0=−3.4786×10⁻⁴×𝑟²₁+𝑟²₂+0.1363609𝑟₁+𝑟₂−26.5486180≤𝑟₁,𝑟₂≤2120elsewhere. (37)
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.

Table 6 

Parameters for Example 4.

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] 𝐺₁𝑟₁,𝑐(𝑡),𝑇(𝑡)=𝑘_𝑔1𝑆^𝑔10.11+0.6𝑟₁,𝐺₂𝑟₂,𝑐(𝑡),𝑇(𝑡)=𝑘_𝑔2𝑆^𝑔20.11+0.6𝑟₂, (38) and the temperature is assumed to follow the expression [17] 𝑇(𝑡)^∘C=32−41−𝑒^−𝑡/310. (39) The saturated solute concentration 𝑐_sat, nucleation rate 𝐵₀, and solute concentration 𝑐 are the same as Example 4. The kinetic parameters of 𝑔₁, 𝑘_𝑔1, 𝑔₂, 𝑘_𝑔2, 𝑘_𝑏, 𝑏, and 𝜌_𝑐 are also the same as Example 4 which are listed in Table 6.
The initial distribution is taken as [14] 𝑓𝑟₁,𝑟₂=𝑟,0−0.0034786²₁+𝑟²₂+1.363609𝑟₁+𝑟₂−26.548618≤𝑟₁,𝑟₂≤210elsewhere, (40) and the computational zone is assumed as [0,40𝜇m]×[0,80𝜇m].
The final crystal size distribution is shown in Figure 6. Also, the “exact” solution is obtained by WENO method on fine grid (Δ𝑟₁=Δ𝑟₂=0.1), while the other results are obtained with grid size Δ𝑟₁=Δ𝑟₂=1. 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.