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

Authors: Chunlei Ruan; Kunfeng Liang; Xianjie Chang; Ling Zhang
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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)∂f(r1,r2,t)∂t+∑j=12∂{Gj(r1,r2,c(t),T(t))f(r1,r2,t)}∂rj=h(f(r1,r2,t),c(t),T(t)),
where r1,r2 are two characteristic lengths for crystals, f is the crystal size distribution, G1,G2 are the growth rate in r1,r2 direction, c is a time-dependent solute concentration, T is a time-dependent temperature, and h 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 ofG1, G2, and h 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 Section2, 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)∂f∂t+∂p(f)∂x+∂q(f)∂y=h,
where p=G1f, q=G2f, x=r1, and y=r2. The computational domain is assumed as [a,b]×[c,d], and the cells are expressed as Iij=[xi-1/2,xi+1/2]×[yj-1/2,yj+1/2], 1≤i≤Nx, 1≤j≤Ny, where Nx is the total cell number in x direction and Ny is the total cell number in y direction. Set the center of cell Iij as (xi,yj), then xi=(1/2)(xi-1/2+xi+1/2), yi=(1/2)(yj-1/2+yj+1/2). In this work, uniform meshes are applied, which means, the internal mesh sizes are Δx=xi+1/2-xi-1/2 and Δy=yj+1/2-yj-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 valuef on the grid nodes (xi,yj) (viz., fi,j) is known. Equation (2) is discretized in the conservative form in space and leads to
(3)ddtfi,j+1Δx(p~i+1/2,j-p~i-1/2,j)+1Δy(q~i,j+1/2-q~i,j-1/2)=hi,j.
Here, (p~i+1/2,j-p~i-1/2,j)/Δx is approximate to ∂p/∂x|(i,j) and (q~i,j+1/2-q~i,j-1/2)/Δy is approximate to ∂q/∂y|(i,j), with p~i+1/2,j,p~i-1/2,j,q~i,j+1/2,q~i,j-1/2 the numerical flux on the cell boundaries.To construct the∂p/∂x|(i,j) and ∂q/∂y|(i,j), it is important to use the upwind idea in the following way [21]:
(4)∂p∂x|(i,j)={(px)i,j-(G1)i≥0(px)i,j+(G1)i<0,∂q∂y|(i,j)={(qy)i,j-(G2)j≥0(qy)i,j+(G2)j<0.
To find (px)i,j-, set
(5)u1=pi-2,j-pi-3,jΔx,u2=pi-1,j-pi-2,jΔx,u3=pi,j-pi-1,jΔx,u4=pi+1,j-pi,jΔx,u5=pi+2,j-pi+1,jΔx,
and to find (px)i,j+, set
(6)u1=pi+3,j-pi+2,jΔx,u2=pi+2,j-pi+1,jΔx,u3=pi+1,j-pi,jΔx,u4=pi,j-pi-1,jΔx,u5=pi-1,j-pi-2,jΔx.
Defining the smoothness [21]
(7)IS1=1312(u1-2u2+u3)2+14(u1-4u2+3u3)2,IS2=1312(u2-2u3+u4)2+14(u2-u4)2,IS3=1312(u3-2u4+u5)2+14(3u3-4u4+u5)2
and the weights [21]
(8)Iω1=a1a1+a2+a3,Iω2=a2a1+a2+a3,Iω3=a3a1+a2+a3,a1=1101(IS1+ε)2,a2=6101(IS2+ε)2,a3=3101(IS3+ε)2,
the 5th order WENO scheme in x direction is finally written as [21]
(9)(px)i,j±=Iω1(u13-7u26+11u36)+Iω2(-u26+5u36+u43)+Iω3(u33+5u46-u56).
Here, ε is a small positive number to prevent the denominators from becoming zero. In this work, ε is set to 10-6.Similarly, to find(qy)i,j-, set
(10)v1=qi,j-2-qi,j-3Δy,v2=qi,j-1-qi,j-2Δy,v3=qi,j-qi,j-1Δy,v4=qi,j+1-qi,jΔy,v5=qi,j+2-qi,j+1Δy,
and to find (qy)i,j+, set
(11)v1=qi,j+3-qi,j+2Δy,v2=qi,j+2-qi,j+1Δy,v3=qi,j+1-qi,jΔy,v4=qi,j-qi,j-1Δy,v5=qi,j-1-qi,j-2Δy.
Defining the smoothness
(12)JS1=1312(v1-2v2+v3)2+14(v1-4v2+3v3)2,JS2=1312(v2-2v3+v4)2+14(v2-v4)2,JS3=1312(v3-2v4+v5)2+14(3v3-4v4+v5)2
and the weights
(13)Jω1=a1a1+a2+a3,Jω2=a2a1+a2+a3,Jω3=a3a1+a2+a3,a1=1101(JS1+ε)2,a2=6101(JS2+ε)2,a3=3101(JS3+ε)2,
the 5th order WENO scheme in y direction is written as
(14)(qy)i,j±=Jω1(v13-7v26+11v36)+Jω2(-v26+5v36+v43)+Jω3(v33+5v46-v56).Equations (9) and (14) are 5th order accuracy in space. Take the (px)i,j-, for example, it covers six cells: Ii-3,j, Ii-2,j, Ii-1,j, Ii,j, Ii+1,j, and Ii+2,j and three stencils: Si,j0={Ii-3,j,Ii-2,j,Ii-1,j,Ii,j}, Si,j1={Ii-2,j,Ii-1,j,Ii,j,Ii+1,j}, and Si,j2={Ii-1,j,Ii,j,Ii+1,j,Ii+2,j}. (px)i,j- is a 5th order approximation to ∂p/∂x|(i,j) which provides the smallest truncation error on these six-cell stencils. Furthermore, it is a convex combination of Si,j0, Si,j1, and Si,j2. 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)dfdt=L(f),
where L(f)=h-(∂p/∂x)-(∂q/∂y). The third order Runge-Kutta schemes can be expressed as
(16)f(1)=f(0)+ΔtL(f(0)),f(2)=f(1)+Δt4(-3L(f(0))+L(f(1))),f(3)=f(2)+Δt12(-L(f(0))-L(f(1))+8L(f(2))),
with f(0)=fn and f(3)=fn+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(17)L1-error=∑i=1N∑j=1M|fi,je-fi,j|Δr1Δr2,L2-error=∑i=1N∑j=1M(fi,je-fi,j)2Δr1Δr2,
where fi,je is the exact solution.The HR scheme for solving PBE (1) is written as [4, 17, 18]
(18)fi,jn+1=fi,jn+Δt×{hi,jn-[(G1f)i+,jn-(G1f)i-,jn]Δr1-[(G2f)i,j+n-(G2f)i,j-n]Δr2}
when Δt is the time step, Δr1,Δr2 are the space step in r1 direction and r2 direction, respectively. The expression of (G1f)i+,jn and (G2f)i,j+n is defined as follows:
(19)(G1f)i+,jn=(G1)i+,jn(fi,jn+12ϕ(ri+n)(fi+1,jn-fi,jn)),(G2f)i,j+n=(G2)i,j+n(fi,jn+12ϕ(rj+n)(fi,j+1n-fi,jn)),
where ri+n=(fi,jn-fi-1,jn)/(fi+1,jn-fi,jn), rj+n=(fi,jn-fi,j-1n)/(fi,j+1n-fi,jn), and ϕ(r)=(|r|+r)/(1+|r|).


### 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:(20)∂f(r1,r2,t)∂t+∂f∂r1+∂f∂r2=0.
With the initial condition f(r1,r2,0)=sin((π/2)(r1+r2)), -2≤r1≤2,-2≤r2≤2, and the boundary conditions are assumed to be periodic. The analytical solution of this problem is f(r1,r2,t)=sin((π/2)(r1+r2-2t)).Table1 shows the accuracy of the WENO scheme at t=2.0. It is notable that present WENO scheme works well in this smooth case and the L1-order and L2-order are around 5. Therefore, the validity of the numerical method is proved.Table 1
Accuracy of WENO scheme in 2D linear equation.
Δ
r
1
(
Δ
r
2
)
L
1-error
L
1-order
L
2-error
L
2-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
2.18444
e
-
005
5.08
6.24975
e
-
006
5.11
1/40
6.36434
e
-
007
5.10
2.03834
e
-
007
4.94Example 1 (size-independent growth (smooth distribution)).
We first consider a homogenous PBE with size-independent growth which is expressed as(21)∂f(r1,r2,t)∂t+∂[G1f]∂r1+∂[G2f]∂r2=0,
with the initial distribution set as f(r1,r2,0)=exp{-100×((r1-0.25)2+(r2-0.25)2)}. 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 G1=G2=0.1μm/s. The analytical solution of this problem is
(22)f(r1,r2,t)=f(r1-G1t,r2-G2t,0).
The computational domain is discretized into 100×100 cells with Δr1=Δr2=0.01. Figure 1 gives a comparison of the simulated CSDs with the analytical CSDs at t=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 L1 and L2 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
L
1-error
0.00308458
4.05432
e
-
005
L
2-error
0.0107827
0.000131731Final distribution comparison for Example1: (a) exact, (b) HR, and (c) WENO.
(a)
(b)
(c)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
(23)f(r1,r2,0)={10.1≤r1,r2≤0.30elsewhere.
The computational domain is taken as [0,1μm]×[0,1μm], and the growth terms are taken as G1=G2=0.1μm/s. The analytical solution is the same as Example 1 which is given as (22). Numerical results obtained with grid size Δr1=Δr2=0.01 and time interval Δt=0.01 are shown in Figure 2 at the final time t=6s. It can be seen that in this nonsmooth distribution case, WENO method also shows higher accuracy and less diffusion. Table 3 compares the L1 and L2 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
L
1-error
0.0115517
0.00635183
L
2-error
0.0626252
0.0423699Final distribution comparison for Example2: (a) exact, (b) HR, and (c) WENO.
(a)
(b)
(c)Example 3 (size-dependent growth).
Next, we consider a homogenous PBE with size-dependent growth which is expressed as follows:(24)∂f(r1,r2,t)∂t+∂[G1(r1)f]∂r1+∂[G2(r2)f]∂r2=0.
Here, we consider that the growth rate is dependent linearly on size
(25)G1(r1)=G0r1G2(r2)=G0′r2,
where G0, G0′ are the constants and the initial distribution is
(26)f(r1,r2,0)=N0L-exp(-r1+r2L-),
where L- is the mean size of the distribution and N0 is the total initial number of particles. The analytical solution is
(27)f(r1,r2,t)=N0L-exp(-r1L-e-G0t-G0t-r2L-e-G0′t-G0′t).
Parameters in this problem are listed in Table4. 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 Example3.
Description
Symbols
Value
Unit
Totalr1-size distribution
r
1
0.1
μm
Totalr2-size distribution
r
2
0.1
μm
Mesh size
Δ
r
1
10−3
μm
Mesh size
Δ
r
2
10−3
μm
Mean size of the distribution
L
-
10−2
μm
Total initial number of particles
N
0
1
Dimensionless
Growth rate
G
0
0.1
μm/s
Growth rate
G
0
′
0.2
μm/s
Simulation time
t
4.0
sTable 5
Total errors in the schemes.
Method
HR
WENO
L
1-error
0.000550266
0.000135244
L
2-error
0.0195435
0.00982604Final distribution comparison for Example3: (a) exact, (b) HR, and (c) WENO.
(a)
(b)
(c)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)∂f(r1,r2,t)∂t+∂[G1(c,T)f]∂r1+∂[G2(c,T)f]∂r2=B0(c,T)δ(r1)δ(r2).
For potassium dihydrogen phosphate (KDP), the crystal shape is shown in Figure 4, where r1 and r2 are the width and length of the crystal, respectively. It is obvious that the volume of a single crystal is [4]
(29)Vc=13r13+(r2-r1)r12.
The growth rates are assumed to be independent ofr1 and r2 and follow the simple power law function of the supersaturation [4]:
(30)G1(c(t),T(t))=kg1Sg1G2(c(t),T(t))=kg2Sg2,
where S=c(t)/csat(T)-1 is the relative supersaturation, c is the solute concentration, and csat is the saturated solute concentration. The expression for csat is given as [4]
(31)Csat(T)=9.3027×10-5T2-9.7629×10-5T+0.2087,
where T is the temperature and follows the linear decreasing profile
(32)T(t)(C°)=33-t720.
The nucleation rate term is given as [4, 17]
(33)B0(c(t),T(t))=kbSbV(t).
Here, b is generally assumed to be related to the stirring power, and V(t) is the total volume of the crystal [1]:
(34)V(t)=∬0∞f(r1,r2,t)Vc(r1,r2)dr1dr2.
The mass balance of a batch system for the solute of which the crystalline phase is being assembled is [1]
(35)dc(t)dt=-ρc∬0∞f(r1,r2,t)×(G·∇Vc)dr1dr2,
where G=(G1,G2). For KDP, this equation is
(36)dc(t)dt=-ρc∬0∞f(r1,r2,t)×(2G1(r1r2-r12)+G2r12)dr1dr2.
The simulated initial distribution is given as [4]
(37)f(r1,r2,0)={-3.4786×10-4×(r12+r22)+0.1363609(r1+r2)-26.5486180≤r1,r2≤2120elsewhere.
Parameters in this problem are listed in Table6. 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 Example4.
Description
Symbols
Value
Unit
Totalr1-size distribution
r
1
500
μm
Totalr2-size distribution
r
2
900
μm
Mesh size
Δ
r
1
2
μm
Mesh size
Δ
r
2
2
μm
r
1-Growth exponent
g
1
1.48
Dimensionless
r
1-Growth rate constant
k
g
1
12.21
μm/s
y-Growth exponent
g
2
1.74
Dimensionless
y-Growth rate constant
k
g
2
100.75
μm/s
Nucleation rate constant
k
b
7.49
×
10
-
8
1/μm3s
Nucleation exponent
b
2.04
Dimensionless
Density of crystals
ρ
c
2.338
×
10
-
12
g/μm3Figure 4
Characteristic lengths of KDP crystals.Final distribution comparison for Example4: (a) exact, (b) HR, and (c) WENO.
(a)
(b)
(c)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 G1 and G2 are size-dependent.
Here, the growth rates are given as [17]
(38)G1(r1,c(t),T(t))=kg1Sg10.1(1+0.6r1),G2(r2,c(t),T(t))=kg2Sg20.1(1+0.6r2),
and the temperature is assumed to follow the expression [17]
(39)T(t)(C°)=32-4(1-e-t/310).
The saturated solute concentration csat, nucleation rate B0, and solute concentration c are the same as Example 4. The kinetic parameters of g1, kg1, g2, kg2, kb, b, and ρc are also the same as Example 4 which are listed in Table 6.
The initial distribution is taken as [14]
(40)f(r1,r2,0)={-0.0034786(r12+r22)+1.363609(r1+r2)-26.548618≤r1,r2≤210elsewhere,
and the computational zone is assumed as [0,40μm]×[0,80μm].
The final crystal size distribution is shown in Figure6. Also, the “exact” solution is obtained by WENO method on fine grid (Δr1=Δr2=0.1), while the other results are obtained with grid size Δr1=Δr2=1. From Figure 6, it is clear that WENO method has better resolution.Final distribution comparison for Example5: (a) exact, (b) HR, and (c) WENO.
(a)
(b)
(c)


## 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:(20)∂f(r1,r2,t)∂t+∂f∂r1+∂f∂r2=0.
With the initial condition f(r1,r2,0)=sin((π/2)(r1+r2)), -2≤r1≤2,-2≤r2≤2, and the boundary conditions are assumed to be periodic. The analytical solution of this problem is f(r1,r2,t)=sin((π/2)(r1+r2-2t)).Table1 shows the accuracy of the WENO scheme at t=2.0. It is notable that present WENO scheme works well in this smooth case and the L1-order and L2-order are around 5. Therefore, the validity of the numerical method is proved.Table 1
Accuracy of WENO scheme in 2D linear equation.
Δ
r
1
(
Δ
r
2
)
L
1-error
L
1-order
L
2-error
L
2-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
2.18444
e
-
005
5.08
6.24975
e
-
006
5.11
1/40
6.36434
e
-
007
5.10
2.03834
e
-
007
4.94Example 1 (size-independent growth (smooth distribution)).
We first consider a homogenous PBE with size-independent growth which is expressed as(21)∂f(r1,r2,t)∂t+∂[G1f]∂r1+∂[G2f]∂r2=0,
with the initial distribution set as f(r1,r2,0)=exp{-100×((r1-0.25)2+(r2-0.25)2)}. 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 G1=G2=0.1μm/s. The analytical solution of this problem is
(22)f(r1,r2,t)=f(r1-G1t,r2-G2t,0).
The computational domain is discretized into 100×100 cells with Δr1=Δr2=0.01. Figure 1 gives a comparison of the simulated CSDs with the analytical CSDs at t=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 L1 and L2 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
L
1-error
0.00308458
4.05432
e
-
005
L
2-error
0.0107827
0.000131731Final distribution comparison for Example1: (a) exact, (b) HR, and (c) WENO.
(a)
(b)
(c)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
(23)f(r1,r2,0)={10.1≤r1,r2≤0.30elsewhere.
The computational domain is taken as [0,1μm]×[0,1μm], and the growth terms are taken as G1=G2=0.1μm/s. The analytical solution is the same as Example 1 which is given as (22). Numerical results obtained with grid size Δr1=Δr2=0.01 and time interval Δt=0.01 are shown in Figure 2 at the final time t=6s. It can be seen that in this nonsmooth distribution case, WENO method also shows higher accuracy and less diffusion. Table 3 compares the L1 and L2 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
L
1-error
0.0115517
0.00635183
L
2-error
0.0626252
0.0423699Final distribution comparison for Example2: (a) exact, (b) HR, and (c) WENO.
(a)
(b)
(c)Example 3 (size-dependent growth).
Next, we consider a homogenous PBE with size-dependent growth which is expressed as follows:(24)∂f(r1,r2,t)∂t+∂[G1(r1)f]∂r1+∂[G2(r2)f]∂r2=0.
Here, we consider that the growth rate is dependent linearly on size
(25)G1(r1)=G0r1G2(r2)=G0′r2,
where G0, G0′ are the constants and the initial distribution is
(26)f(r1,r2,0)=N0L-exp(-r1+r2L-),
where L- is the mean size of the distribution and N0 is the total initial number of particles. The analytical solution is
(27)f(r1,r2,t)=N0L-exp(-r1L-e-G0t-G0t-r2L-e-G0′t-G0′t).
Parameters in this problem are listed in Table4. 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 Example3.
Description
Symbols
Value
Unit
Totalr1-size distribution
r
1
0.1
μm
Totalr2-size distribution
r
2
0.1
μm
Mesh size
Δ
r
1
10−3
μm
Mesh size
Δ
r
2
10−3
μm
Mean size of the distribution
L
-
10−2
μm
Total initial number of particles
N
0
1
Dimensionless
Growth rate
G
0
0.1
μm/s
Growth rate
G
0
′
0.2
μm/s
Simulation time
t
4.0
sTable 5
Total errors in the schemes.
Method
HR
WENO
L
1-error
0.000550266
0.000135244
L
2-error
0.0195435
0.00982604Final distribution comparison for Example3: (a) exact, (b) HR, and (c) WENO.
(a)
(b)
(c)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)∂f(r1,r2,t)∂t+∂[G1(c,T)f]∂r1+∂[G2(c,T)f]∂r2=B0(c,T)δ(r1)δ(r2).
For potassium dihydrogen phosphate (KDP), the crystal shape is shown in Figure 4, where r1 and r2 are the width and length of the crystal, respectively. It is obvious that the volume of a single crystal is [4]
(29)Vc=13r13+(r2-r1)r12.
The growth rates are assumed to be independent ofr1 and r2 and follow the simple power law function of the supersaturation [4]:
(30)G1(c(t),T(t))=kg1Sg1G2(c(t),T(t))=kg2Sg2,
where S=c(t)/csat(T)-1 is the relative supersaturation, c is the solute concentration, and csat is the saturated solute concentration. The expression for csat is given as [4]
(31)Csat(T)=9.3027×10-5T2-9.7629×10-5T+0.2087,
where T is the temperature and follows the linear decreasing profile
(32)T(t)(C°)=33-t720.
The nucleation rate term is given as [4, 17]
(33)B0(c(t),T(t))=kbSbV(t).
Here, b is generally assumed to be related to the stirring power, and V(t) is the total volume of the crystal [1]:
(34)V(t)=∬0∞f(r1,r2,t)Vc(r1,r2)dr1dr2.
The mass balance of a batch system for the solute of which the crystalline phase is being assembled is [1]
(35)dc(t)dt=-ρc∬0∞f(r1,r2,t)×(G·∇Vc)dr1dr2,
where G=(G1,G2). For KDP, this equation is
(36)dc(t)dt=-ρc∬0∞f(r1,r2,t)×(2G1(r1r2-r12)+G2r12)dr1dr2.
The simulated initial distribution is given as [4]
(37)f(r1,r2,0)={-3.4786×10-4×(r12+r22)+0.1363609(r1+r2)-26.5486180≤r1,r2≤2120elsewhere.
Parameters in this problem are listed in Table6. 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 Example4.
Description
Symbols
Value
Unit
Totalr1-size distribution
r
1
500
μm
Totalr2-size distribution
r
2
900
μm
Mesh size
Δ
r
1
2
μm
Mesh size
Δ
r
2
2
μm
r
1-Growth exponent
g
1
1.48
Dimensionless
r
1-Growth rate constant
k
g
1
12.21
μm/s
y-Growth exponent
g
2
1.74
Dimensionless
y-Growth rate constant
k
g
2
100.75
μm/s
Nucleation rate constant
k
b
7.49
×
10
-
8
1/μm3s
Nucleation exponent
b
2.04
Dimensionless
Density of crystals
ρ
c
2.338
×
10
-
12
g/μm3Figure 4
Characteristic lengths of KDP crystals.Final distribution comparison for Example4: (a) exact, (b) HR, and (c) WENO.
(a)
(b)
(c)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 G1 and G2 are size-dependent.
Here, the growth rates are given as [17]
(38)G1(r1,c(t),T(t))=kg1Sg10.1(1+0.6r1),G2(r2,c(t),T(t))=kg2Sg20.1(1+0.6r2),
and the temperature is assumed to follow the expression [17]
(39)T(t)(C°)=32-4(1-e-t/310).
The saturated solute concentration csat, nucleation rate B0, and solute concentration c are the same as Example 4. The kinetic parameters of g1, kg1, g2, kg2, kb, b, and ρc are also the same as Example 4 which are listed in Table 6.
The initial distribution is taken as [14]
(40)f(r1,r2,0)={-0.0034786(r12+r22)+1.363609(r1+r2)-26.548618≤r1,r2≤210elsewhere,
and the computational zone is assumed as [0,40μm]×[0,80μm].
The final crystal size distribution is shown in Figure6. Also, the “exact” solution is obtained by WENO method on fine grid (Δr1=Δr2=0.1), while the other results are obtained with grid size Δr1=Δr2=1. From Figure 6, it is clear that WENO method has better resolution.Final distribution comparison for Example5: (a) exact, (b) HR, and (c) WENO.
(a)
(b)
(c)


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



---

*Source: 125128-2013-09-23.xml*