Chebyshev Spectral Collocation Method for Population Balance Equation in Crystallization

The population balance equation (PBE) is the main governing equation for modeling dynamic crystallization behavior. In the view of mathematics, PBE is a convection–reaction equation whose strong hyperbolic property may challenge numerical methods. In order to weaken the hyperbolic property of PBE, a diffusive term was added in this work. Here, the Chebyshev spectral collocation method was introduced to solve the PBE and to achieve accurate crystal size distribution (CSD). Three numerical examples are presented, namely size-independent growth, size-dependent growth in a batch process, and with nucleation, and size-dependent growth in a continuous process. Through comparing the results with the numerical results obtained via the second-order upwind method and the HR-van method, the high accuracy of Chebyshev spectral collocation method was proven. Moreover, the diffusive term is also discussed in three numerical examples. The results show that, in the case of size-independent growth (PBE is a convection equation), the diffusive term should be added, and the coefficient of the diffusive term is recommended as 2G × 10−3 to G × 10−2, where G is the crystal growth rate.


Introduction
Crystallization operations are widely used in chemical and pharmaceutical industries.Crystallization kinetics is the basis for the development of crystallization theory.Simultaneously, studies on crystalline kinetics also guide the design of crystallizers and the quality control of crystalline products [1].The population balance equation (PBE) is a widely accepted equation to model the dynamic crystallization behavior [2].Therefore, it is of great significance to accurately solve the PBE and to achieve accurate crystal size distribution (CSD).
Currently, there are a many numerical results on the study of PBE.In terms of numerical methods, they can be roughly divided into the following five categories [1,2]: (1) the method of moments [3]; (2) the method of characteristics [4,5]; (3) the method of weight residuals/orthogonal collocation [6,7]; (4) the Monte Carlo method [8,9]; (5) the finite-difference schemes/discrete population balances [10].The review by Ramkrishna [2] is helpful in this regard.The method of moments approximates the CSD through its moments.The predictions under certain conditions are closed.However, the moment closure conditions are violated for more complex systems.The method of characteristics tries reducing the partial differential equations (PDEs) to ordinary differential equations (ODEs) using lines of characteristics.The approach is efficient with simple physics, but it is not suitable for complex physics.The method of weight residuals/orthogonal collocation approximates the distribution by the basis function; the main weakness is the dependence on the choice of basis function.The Monte Carlo method tracks the histories of individual particles.It has an incomparable degree of omnipotence.However, its computational cost is relatively expensive.The finite-difference/volume methods are approximated by finite-difference/volume schemes.They can be applied to complex cases.However, they usually require a high number of grid points in order to achieve desirable accuracy.
Recently, a large number of numerical methods with high accuracy were constructed to solve PBE, e.g., the HR-van method with a flux limiter [1,11], the Lattice Boltzmann method [12], the weighted essentially non-oscillatory (WENO) method [13], the spectral method [14,15], etc.Among these, the spectral method has the advantage of being able to accurately solve the PBE using a lower number of grid points.It was also found that the spectral collocation method is more efficient than other spectral methods [14] and finite-volume schemes [15].
Because the PBE is a convection-reaction equation in mathematics, its strong hyperbolic property may challenge the numerical methods.The numerical methods may gain unstable results, which are not applicable for physics.In order to weaken the hyperbolic property of PBE, a diffusive term is added herein.Although the diffusive term makes the numerical solution deviate from the exact solution to some extent, it guarantees the stability of the numerical results.
In this paper, Chebyshev spectral collocation is presented to solve the one-dimensional PBE.The CSD of three types of crystallization in batch crystallizers and continuous crystallizers were calculated.By comparing the analytical results and other numerical results, the validity and high accuracy of the Chebyshev spectral collocation method were verified.Furthermore, the effects of the diffusive term are also discussed using three numerical examples.

One-Dimensional PBE
PBE is a differential equation which describes the evolution of a population of particles.In the crystallization process, it can reflect the crystal nucleation, growth, agglomeration, and breakage [2].
One-dimensional PBE [1,2,11] can be written as follows: where f (L, t) is the crystal size distribution function, L is the crystal size variable, t is the time, G(L, t) is the growth rate of crystals, and h(L, t, f ) is the source term and may include the aggregation, nucleation, breakage, etc.Here, crystal size variable L is typically the characteristic length, volume, or mass, but it can also represent age, composition, and other characteristics of an entity in a distribution [1].The growth rate G(L, t) can be a function of size and other variables, such as the temperature and the concentration of chemical species in solution [1].The source term h(L, t, f ) is a creation/depletion rate, and it can be a function of other variables including the distribution, which occurs in nucleation processes resulting from particle-particle interactions and in agglomeration processes [1].Since many of these processes involve integrals, PBE is considered an integro-differential equation describing the complex crystallization processes [2].Since Equation ( 1) is a convection-reaction equation, its strong hyperbolic property may challenge the numerical methods.In order to weaken the hyperbolic property of Equation (1), a diffusive term Here, ε is a small positive number.Theoretically, the closer ε is to 0, the more accurate the result is.However, in numerical experiments, when ε approaches 0, the effect of the diffusive term becomes smaller, which leads to an unstable numerical solution.Therefore, in this paper, the value of ε was determined by numerical tests.

Chebyshev Spectral Collocation Method
The Chebyshev spectral collocation method is a kind of collocation method [16,17].On the one hand, as a spectral method, it has high accuracy.On the other hand, as a collocation method, it transforms the derivative term into a differential matrix, which makes the equations simple and easy to solve.To date, this method was successfully applied to solve differential equations in the fields of fluid mechanics [18,19], quantum mechanics [20,21], etc. Equation ( 2) can be written as follows: where The fourth-order Runge-Kutta method [22] is used to discretize the time, leading to Here, ∆t is the time step, t n = n∆t is the time on the nth step (n = 0, 1, • • • ), f n is the value of f on t n and the initial distribution is written as f 0 .It is worth mentioning that the variable f n in Equation ( 4) is only a function of the spatial variable L. In this way, the space can be discretized using the spectral collocation method.
In the spectral collocation method, the solution of f n (L) can be written as follows [16,17]: where P k (L) is the Lagrange interpolation polynomial at the node and u k is the function value at interpolation node L k .
In the Chebyshev spectral collocation method, Chebyshev points are chosen as the interpolation nodes.In the standard calculation domain [−1, 1], the Chebyshev points are Lj = cos(jπ/N)(j = 0, 1, • • • N).As for the general calculation domain [a, b], the Chebyshev points can be obtained by the coordinate transformation of standard Chebyshev points, namely Here, we take the calculation of R(t n , f n (L)) in Equation (4) as an example to show the implementation of the Chebyshev spectral collocation method.In fact, the right-hand term R(t n , f n (L)) Here, we need to show the calculation of d f n (L) dL and d 2 f n (L) dL 2 terms.In the Chebyshev spectral collocation method [16,17], the derivative term d f n (L) dL can be transformed into a differential matrix as follows: Here, D N is a (N + 1) × (N + 1) Chebyshev derivative matrix, and u = (u 0 , u 1 , • • • , u N ) T is a N + 1 function value vector at the interpolation nodes.In the standard calculation domain [−1, 1], the Chebyshev derivative matrix DN satisfies the following properties [16,17]: where ( DN ) ij represents the element of the i + 1th row and j + 1th column in the Chebyshev matrix, and c i is the constant which is determined by the following equation: As for the general calculation domain [a, b], the Chebyshev derivative matrix D N can be calculated using the standard Chebyshev derivative matrix DN , which is Similarly, in the Chebyshev spectral collocation method [16,17], the second-order derivative term dL 2 can be written as follows: Here, D 2 N is the second-order Chebyshev derivative matrix.The relationship between the second-order Chebyshev derivative matrix D 2 N and the first-order Chebyshev derivative matrix D N is

Other Numerical Schemes
We mainly focus our attention on the numerical schemes of spatial discretization.The numerical method for time discretization is the classic Runge-Kutta method which is shown in Equation (4).
The diffusive term "ε ∂ 2 f (L, t) ∂L 2 " is no longer taken into account in the following two schemes.Thus, the right-hand term in Equation ( 4

Second Upwind Scheme
In the calculation of ∂ ∂L (G n (L) f n (L)), we use the following second-order upwind scheme [22]: Here, f n i represents the value of f n at node L i , G n i represents the value of G n at node L i , and ∆L represents the distance between two adjacent nodes.During our simulation, the nodes are taken as uniformly distributed nodes, rather than Chebyshev points in the Chebyshev spectral collocation method.

HR-Van Scheme
In the calculation of ∂ ∂L (G n (L) f n (L)), the following HR-van scheme is used [1]: where , with φ(r) being the flux limiter proposed by Van.Similarly, this scheme is implemented under uniformly distributed nodes rather than Chebyshev points.

Numerical Experiments
In order to compare the accuracy of each numerical method, we define the following errors: Here, f e k is the exact solution at node L k .

Size-Independent Growth in a Batch Process
The first example is the crystallization in a batch process with the crystal growth rate G independent of the crystal size L. Here, we do not consider the aggregation, nucleation, and breakage, and we set h(L, t, f ) = 0.The crystal growth rate is set to G = 1.0 µm/s and the initial CSD is given by with the following Gaussian distribution: with µ = 20 µm, σ = 3.The other parameters in the calculation are set to N = 200, dt = 0.001, ε = 0.005, t end = 60 s.The exact solution of this example is f (L, t) = f (L − Gt, 0).
Figure 1 shows the comparison of exact solution with the results of three numerical scheme sat t = 60 s.In order to clearly show the performances with different schemes, Figure 1b shows the predicted crystal size distributions at the local crystal size.It can be seen that second-order upwind scheme is the worst, which causes both numerical dispersion and diffusion; both the HR-van scheme and the Chebyshev spectral collocation method cause numerical diffusion, but the Chebyshev spectral collocation method causes less numerical diffusion.Table 1 shows the L 1 and L 2 errors with different numerical schemes.The error table clearly shows that the Chebyshev spectral collocation method has the smallest error.Therefore, the Chebyshev spectral collocation method has the highest accuracy among these three methods.We now discuss the effects of the diffusive term.Figure 2 shows the numerical results obtained using the Chebyshev spectral collocation method with different coefficients of diffusive term.As can been seen in Figure 2a, when no diffusive term is added ( 0 = ε ), an oscillation happens near the right boundary at s t 2 = .Therefore, in this case, we cannot obtain stable numerical solutions.
However, as shown in Figure 2b, when a diffusive term is added ( errors increase.This is consistent with our explanation for the diffusive term; the closer ε isto0, the more accurate the result is.However, when ε approaches 0, the effect of the diffusive term becomes smaller which leads to an unstable numerical solution.The coefficient of the diffusive term In the size-independent growth case, PBE is a convection equation ( displaysa strong hyperbolic property.The diffusive term should be added in order to obtain a stable and convergent numerical solution.Figure 3 shows the relationship between the optimal ε (with minimum 1 L and 2 L errors) and the growth rate G.It is clear that the optimal ε has a linear relationship with the growth rate G. Therefore, the coefficient of diffusive term ( ε ) is recommended as 2G× 10 −3 to G×10 −2 , where G is the crystal growth rate.We now discuss the effects of the diffusive term.Figure 2 shows the numerical results obtained using the Chebyshev spectral collocation method with different coefficients of diffusive term.As can been seen in Figure 2a, when no diffusive term is added (ε = 0), an oscillation happens near the right boundary at t = 2 s.Therefore, in this case, we cannot obtain stable numerical solutions.However, as shown in Figure 2b, when a diffusive term is added (ε = 0.002), a convergent solution is obtained.The CSD function shifts to the right with the velocity of G. Table 2 lists the L 1 and L 2 errors with different coefficients of the diffusive term.When ε ≤ 0.001, a convergent numerical solution cannot be obtained.With the increase of ε (ε ≥ 0.002), both L 1 and L 2 errors increase.This is consistent with our explanation for the diffusive term; the closer ε is to 0, the more accurate the result is.However, when ε approaches 0, the effect of the diffusive term becomes smaller which leads to an unstable numerical solution.The coefficient of the diffusive term (ε) with G = 1 is recommended as 2 × 10 −3 to 10 −2 .
In the size-independent growth case, PBE is a convection equation ( and displaysa strong hyperbolic property.The diffusive term should be added in order to obtain a stable and convergent numerical solution.Figure 3 shows the relationship between the optimal ε(with minimum L 1 and L 2 errors) and the growth rate G.It is clear that the optimal ε has a linear relationship with the growth rate G. Therefore, the coefficient of diffusive term (ε) is recommended as 2G × 10 −3 to G × 10 −2 , where G is the crystal growth rate.

Size-Dependent Growth in a Batch Process
The second example involves crystallization in a batch process with the crystal growth rate G

Size-Dependent Growth in a Batch Process
The second example involves crystallization in a batch process with the crystal growth rate G

Size-Dependent Growth in a Batch Process
The second example involves crystallization in a batch process with the crystal growth rate G dependent on the crystal size L. Similar to the first example, aggregation, nucleation, and breakage are not taken into account, which means h(L, t, f ) = 0. Crystal growth rate is supposed as the linear function of crystal size, namely G(L, t) = G 0 L, where G 0 is a constant.The initial CSD satisfies the following equation: where N 0 , L are constants.The analytical solution for this case is [1] as follows: Parameters in the simulation are taken as L = 0.01 µm 3 , N 0 = 1, G 0 = 0.1 µm 3 /s, ε = 0, N = 200, ∆t = 0.00001 s.The number of grids in space used for the second upwind method and HR-van method is 1000.
Figure 4 shows the comparison of analytical solutions with numerical solutions obtained by three numerical schemes at t = 4 s.In Figure 4a, there is no significant difference between the analytical solution and the numerical solutions.Figure 4b gives the comparison of the numerical results near the boundary L = 0 µm.From Figure 4b, it is obvious that Chebyshev spectral collocation method performs the best and its result is almost coincident with the analytical solution; however, the second-order upwind scheme and HR-van scheme have some deviations from the analytical solution near the boundary.Table 3 gives the comparison of L 1 and L 2 errors for different numerical schemes.From Table 3, it is clear that the Chebyshev spectral collocation method has the smallest error, whereas the HR-van scheme has a smaller error and the second-order upwind scheme has the biggest error.Thus, the Chebyshev spectral collocation method also has the highest accuracy among these three schemes.
Mathematics 2019, 7, x FOR PEER REVIEW 9 of 13 where L N , 0 are constants.The analytical solution for this case is [1] as follows: Parameters in the simulation are taken as . The number of grids in space used for the second upwind method and HR-van method is 1000.
Figure 4 shows the comparison of analytical solutions with numerical solutions obtained by three numerical schemes at s t 4 = .In Figure4a, there is no significant difference between the analytical solution and the numerical solutions.Figure 4b gives the comparison of the numerical results near the boundary 4b, it is obvious that Chebyshev spectral collocation method performs the best and its result is almost coincident with the analytical solution; however, the second-order upwind scheme and HR-van scheme have some deviations from the analytical solution near the boundary.Table 3 gives the comparison of 1 L and 2 L errors for different numerical schemes.From Table 3, it is clear that the Chebyshev spectral collocation method has the smallest error, whereas the HR-van scheme has a smaller error and the second-order upwind scheme has the biggest error.Thus, the Chebyshev spectral collocation method also has the highest accuracy among these three schemes.
(a) (b)  Figure 5 presents the evolution of CSD obtained using the Chebyshev spectral collocation method.As we can see in Figure 5, there is no oscillation during the computation.Therefore, the diffusive term is not necessary.This is not surprising.In this size-dependent growth case, the PBE is a convection-reaction equation ( ).The reactive term weakens the  Figure 5 presents the evolution of CSD obtained using the Chebyshev spectral collocation method.As we can see in Figure 5, there is no oscillation during the computation.Therefore, the diffusive term is not necessary.This is not surprising.In this size-dependent growth case, the PBE is a convection-reaction equation ( The reactive term weakens the hyperbolic property of PBE and makes it easy to solve.This is reflected in the decay term "exp(−G 0 t)" of the analytical solution in Equation (18).
hyperbolic property of PBE and makes it easy to solve.This is reflected in the decay term " ) exp( 0 t G − " of the analytical solution in Equation ( 18).

Nucleation and Size-Dependent Growth in a Continuous Process
The third example is the crystallization in a continuous process with the nucleation and crystal growth rate G dependent on the crystal size L .Here, the source item is assumed to τ δ ,where z G , , 0 γ are constants.The analytical solution for this case is [1] as follows: In order to compare the prediction effects of each numerical method, Equation ( 19) is given as the initial CSD and the simulation is ended at .The number of grids in space used for the second upwind method and HR-van method is 400.
Figure 6 shows the comparison of analytical solutions with numerical solutions of three numerical schemes.From Figure 6a and Figure6b, we can conclude that the result of the second-order upwind scheme is the worst since the predicted result has the biggest deviation from the analytical solution.Results of the HR-van method and the Chebyshev spectral collocation method are satisfactory.Figure6b shows that the results obtained using the Chebyshev spectral collocation method are slightly better than those of the HR-van method.Table 4 shows the comparison of 1 L and 2 L errors for three different numerical schemes.From the error comparison, we know that the Chebyshev spectral collocation method is the best one, whereas the HR-van scheme is the better one and the second-order upwind scheme is the worst one.Thus, in this case, the Chebyshev spectral collocation method also shows the highest accuracy.

Nucleation and Size-Dependent Growth in a Continuous Process
The third example is the crystallization in a continuous process with the nucleation and crystal growth rate G dependent on the crystal size L. Here, the source item is assumed to h(L, t, f ) = , where B 0 , τ are constants, and δ(L) is a Dirac δ function.The relationship between crystal growth rate and crystal size is G(L) = G 0 (1 + γL) z , where G 0 , γ, z are constants.The analytical solution for this case is [1] as follows: In order to compare the prediction effects of each numerical method, Equation ( 19) is given as the initial CSD and the simulation is ended at t = 10000 s.In our simulation, parameters are set to B 0 = 2 × 10 −10 /µm 3 /s, τ = 100 s, G 0 = 0.00168 µm/s, z = 0.3, γ = 1, ε = 0, N = 40, dt = 0.001.The number of grids in space used for the second upwind method and HR-van method is 400.
Figure 6 shows the comparison of analytical solutions with numerical solutions of three numerical schemes.From Figure 6a,b, we can conclude that the result of the second-order upwind scheme is the worst since the predicted result has the biggest deviation from the analytical solution.Results of the HR-van method and the Chebyshev spectral collocation method are satisfactory.Figure 6b shows that the results obtained using the Chebyshev spectral collocation method are slightly better than those of the HR-van method.Table 4 shows the comparison of L 1 and L 2 errors for three different numerical schemes.From the error comparison, we know that the Chebyshev spectral collocation method is the best one, whereas the HR-van scheme is the better one and the second-order upwind scheme is the worst one.Thus, in this case, the Chebyshev spectral collocation method also shows the highest accuracy.

Conclusions
In this paper, the Chebyshev spectral collocation method was introduced to solve the one-dimensional PBE.The detailed implementation and the related techniques of the Chebyshev

Conclusions
In this paper, the Chebyshev spectral collocation method was introduced to solve the one-dimensional PBE.The detailed implementation and the related techniques of the Chebyshev

Conclusions
In this paper, the Chebyshev spectral collocation method was introduced to solve the one-dimensional PBE.The detailed implementation and the related techniques of the Chebyshev spectral collocation method are given.Three examples are presented, namely size-independent growth, size-dependent growth in a batch process, and with nucleation, and size-dependent growth in a continuous process.The numerical results showed that the Chebyshev spectral collocation method is effective and highly accurate in solving the PBE.However, in the case of size-independent growth, when PBE is a convection equation, it displays a highly hyperbolic property.A diffusive term should be added in order to obtain a stable solution.
For the more complex crystallization processes, e.g., aggregation and breakage, PBE becomes an integro-differential equation.Other methods, e.g., the cell average method and the fixed pivot method, can be combined with the spectral collocation method to obtain accurate CSD [15].Moreover, realistic PBE problems involve multidimensional states.The spectral collocation method can be easily implemented for two-dimensional problems, enabling its use in multivariate systems.However, the spectral collocation method as a spectral method has its limitations in solving PBE, i.e., the computational zone should be simple and the distribution of CSD should be sufficiently smooth.Therefore, under such conditions, the Chebyshev spectral collocation method is recommended to solve PBE for which a high accuracy of CSD is needed.

Figure 1 .
Figure 1.Comparison for final crystal size distribution (CSD) with analytical and numerical results: (a) in the whole crystal;(b) zoomed figure.

1 L and 2 L
solution is obtained.The CSD function shifts to the right with the velocity of G .Table2lists the errors with different coefficients of the diffusive term.When numerical solution cannot be obtained.With the increase of ε (

Figure 1 .
Figure 1.Comparison for final crystal size distribution (CSD) with analytical and numerical results: (a) in the whole crystal; (b) zoomed figure.

Figure 2 .
Figure 2.Comparison for numerical results with different coefficients of the diffusive term: (a) 0 = ε

.G
dependent on the crystal size L .Similar to the first example, aggregation, nucleation, and breakage are not taken into account, which means Crystal growth rate is supposed as the linear function of crystal size, namely is a constant.The initial CSD satisfies the following equation:

Figure 2 .
Figure 2.Comparison for numerical results with different coefficients of the diffusive term: (a) 0 = ε

G
dependent on the crystal size L .Similar to the first example, aggregation, nucleation, and breakage are not taken into account, which means is a constant.The initial CSD satisfies the following equation:

Figure 4 .
Figure 4.Comparison for final CSD with analytical and numerical results: (a) in the whole crystal;(b) zoomed figure.

Figure 4 .
Figure 4. Comparison for final CSD with analytical and numerical results: (a) in the whole crystal; (b) zoomed figure.

Figure 5 .
Figure 5.Evolution of CSD obtained using the Chebyshev spectral collocation method.
δ function.The relationship between crystal growth rate and crystal size is

Figure 5 .
Figure 5. Evolution of CSD obtained using the Chebyshev spectral collocation method.

Figure 6 .
Figure 6.Comparison for final CSD with analytical and numerical results: (a) in the whole crystal; (b) zoomed figure.

Figure 7
Figure7displays the evolution of CSD obtained using the Chebyshev spectral collocation method.It can be seen that it does not result in an unstable solution in the case where no diffusive term is added.In this case, the PBE becomes

Figure 7 .
Figure 7.Evolution of CSD obtained using the Chebyshev spectral collocation method.

Figure 6 .
Figure 6.Comparison for final CSD with analytical and numerical results: (a) in the whole crystal; (b) zoomed figure.

Figure 7
Figure7displays the evolution of CSD obtained using the Chebyshev spectral collocation method.It can be seen that it does not result in an unstable solution in the case where no diffusive term is added.In this case, the PBE becomes∂ f ∂t + G ∂ f ∂L + f ∂G ∂L = B 0 δ(L) − f (L) τ .It is a convection-reaction equation.As shown in example 2, the reactive term weakens the hyperbolic property of PBE.Therefore, the diffusive term is not necessary.

Figure 6 .
Figure 6.Comparison for final CSD with analytical and numerical results: (a) in the whole crystal; (b) zoomed figure.

Figure 7
Figure7displays the evolution of CSD obtained using the Chebyshev spectral collocation method.It can be seen that it does not result in an unstable solution in the case where no diffusive term is added.In this case, the PBE becomes

Figure 7 .
Figure 7.Evolution of CSD obtained using the Chebyshev spectral collocation method.

Figure 7 .
Figure 7. Evolution of CSD obtained using the Chebyshev spectral collocation method.

Table 1 .
1L and 2 L errors for different numerical schemes.

Table 1 .
L 1 and L 2 errors for different numerical schemes.

Table 2 .
1 L and 2L errors for different coefficients of the diffusive term.

Table 2 .
L 1 and L 2 errors for different coefficients of the diffusive term.

Table 2 .
1 L and 2L errors for different coefficients of the diffusive term.

Table 3 .
1L and 2 L errors for different numerical schemes.

Table 3 .
L 1 and L 2 errors for different numerical schemes.

Table 4 .
1L and 2 L errors for different numerical schemes.

Table 4 .
L 1 and L 2 errors for different numerical schemes.

Table 4 .
1L and 2 L errors for different numerical schemes.