Euler’s First-Order Explicit Method–Peridynamic Differential Operator for Solving Population Balance Equations of the Crystallization Process

1  Introduction

Crystallization is a low-energy-consumption process that yields products of high purity. This process is employed in various fields, such as chemical industries, medicine, and biotechnology, and involves finer aspects, such as the nucleation and growth of crystals. The population balance equation (PBE) is a mathematical model that describes the crystal size distribution. An accurate solution of the PBE is important to understand and optimize the crystallization process for the improvement of product quality and minimization of production costs.

At present, a large number of numerical results exist in the study of the PBE. They can be roughly divided into three categories based on the numerical method involved: (i) the Monte Carlo method [1,2], (ii) the method of moments [3,4], and (iii) the discrete method [5–7]. The Monte Carlo method uses a large number of random samplings to solve the PBE and can be directly used to analyze multidimensional cases. Thus, the method has several applications in solving the PBE and two-phase flows [8,9]. However, this method is computationally expensive, is prone to generating noise, and hinders coupling with other equations. The method of moments is used to solve the equations for the various moments of the size-distribution function. It has been reviewed by Li et al. [4]. This method is computationally efficient, but reconstructing the size-distribution function from the moments is challenging [10]. Common examples of the method of moments are the quadrature method of moments [11,12] and direct quadrature method of moments [13]. The discrete method directly solves the PBE with the size-distribution function. This method discretizes the internal coordinates and obtains the size-distribution function directly by applying one of the following methods: the finite volume method [14,15], finite element method [16,17], collocation method [18,19], or lattice Boltzmann method [20]. This study focuses on the discrete method for solving the PBE. Since the PBE is a convection-reaction equation in mathematics, its strong hyperbolic property challenges the numerical method. Due to the strong hyperbolic nature of the equation, numerical schemes inevitably cause numerical dissipation and dispersion problems when solving such problems. Traditionally, this problem is solved by using a scheme called the upwind scheme. When calculating derivatives, it only calls more points from the upwind direction of the flow. Besides, the crystallization simulation also contain sharp crystal size distributions and discontinuities. Thus, the high resolution method are required.

The peridynamic differential operator (PDDO) method, a nonlocal differential operator [21–23], has been developed in recent years based on the concept of the peridynamics (PD) theory [24,25]. It connects local partial derivatives and nonlocal integrals through Taylor series expansions and the property of orthogonal functions. The PDDO method can be used to solve differential equations and derive smooth functions or scattered data from discontinuities or singular points. Thus, this method has been successfully applied to the study of stress analysis [26], heat and mass transfer [27], and thermodynamic problems [28]. Recently, Baryt et al. [29] modified the weight function of the PDDO based on the characteristics of hyperbolic equations and obtained a generalized PDDO scheme to accurately solve these equations. Herein, we attempt to apply this method to solve the PBE. The reason we use the PDDO method for dealing with the internal crystal-size derivatives in PBE is that the PDDO method has the ability to deal with the discontinuity in mathematical concepts. Besides, the PDDO can determine any arbitrary order of derivatives regardless of the presence of jump discontinuities or singularities. Also, the PDDO is free of the requirement of symmetric kernels; thus, it eliminates the necessity of ghost particles near the boundary.

In this study, the Euler’s first-order explicit (EE) method and PDDO method are combined (EE–PDDO) to solve the one-dimensional (1D)-PBE, and four cases are presented: (i) size-independent growth, (ii) size-dependent growth in a batch process, (iii) nucleation and size-independent growth, and (iv) nucleation and size-dependent growth in a continuous process. The effectiveness and high accuracy of the EE–PDDO method are verified by comparing it to the analytical solution and other numerical results. The stability of the EE–PDDO method, the weight function in the PDDO scheme, and the optimal time step size are discussed.

2  Population Balance Equation and Numerical Schemes

2.1 Population Balance Equation

The 1D-PBE is expressed as Eq. (1) [14,30].

∂f(L,t)∂t+∂{G(L,t)f(L,t)}∂L=h(L,t,f)(1)

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 and is positive, h(L,t,f) is a source term that includes the effects of aggregation, nucleation, and breakage. Since Eq. (1) is a convection-reaction equation in mathematics, it brings challenges to numerical simulation.

2.2 PDDO Method

The PD theory is a nonlocal theory proposed by Silling et al. [31]. Point x interacts with Point x′ within an interaction domain Hx as shown in Fig. 1. The relative position vector between these points is defined as ξ=x′−x. Each point has its interaction domain (family). The PD range of each point is known as the horizon and can be denoted by the integer parameter (m) and grid size Δx, namely, δ=mΔx. The interaction domain of points may have different sizes and shapes. The degree of nonlocal interaction between the points is specified by the weight function ω(ξ).

Figure 1: PD interaction domains for the discretized Points x and x′

Madenci et al. [24,25] proposed the PDDO method based on the PD theory. The N-order Taylor series expansions of scalar fields f(x′)=f(x+ξ) are expressed as Eq. (2).

f(x+ξ)=∑n=0N1n!ξndnf(x)dxn+R(N,x)(2)

where R(N,x) is the remainder of the N-order approximation. Multiplying each term in Eq. (2) by PD functions gNp(ξ) and integrating over the domain of interaction Hx forms Eq. (3).

dpf(x)dxp=∫Hxf(x+ξ)gNp(ξ)dξ(3)

The PD function must be orthogonal to each term in the Taylor series expansions as shown in Eq. (4).

∫HxξngNp(ξ)dV=n!δnp(4)

where n=0,…N, δnp is the Kronecker symbol. PD functions gNp(ξ) can be expressed using Eq. (5) [24,25].

gNp(ξ)=∑q=0Naqpωq(ξ)ξq(5)

where ωq(ξ) is the weight functions associated with each term ξq.

The unknown coefficient aqp in Eq. (5) can be obtained from Eq. (6).

∑q=0NAnqaqp=bnp(6)

where q=0,…,N. The coefficient matrix is obtained from Eqs. (7a) and (7b).

Anq=∫Hxωq(ξ)ξn+qdξ(7a)

bnp=n!δnp(7b)

The PDDO method recovers local differentiation as the family size Hx decreases or the number of terms in the function gNp(ξ) increases. Thus, the degree of nonlocality reduces by decreasing the family size and increasing the number of terms in the Taylor series expansions.

2.3 EE–PDDO Method for Solving the Population Balance Equation

We initially introduced the internal coordinate grid points for computing the solution and labeled them sequentially as L1,L2,⋯,LML. We confined our attention to a uniform grid with a step size ΔL. The formula for the internal crystal-size points is Li=a+(i−12)⋅ΔL,(i=1,2,⋯,ML). ΔL and ML are connected via ΔL=(b−a)/ML. We introduced the time points for computing the solution, which is labeled sequentially as t0,t1,⋯,tMt. We confine our attention to a uniform grid with a time step size Δt. Hence, the formula for the time points is tk=k⋅Δt,(k=0,1,⋯,Mt). Δt and Mt are connected via Δt=tend/Mt.

We evaluated the differential Eq. (1) at the point (Li,tk) to obtain Eq. (8).

∂f(Li,tk)∂t+∂{G(Li,tk)f(Li,tk)}∂L=h(Li,tk,f(Li,tk))(8)

The EE method was used for the time derivative to obtain Eq. (9).

fik+1=fik+Δt⋅[hik−Gik∂f∂L|L=Li,t=tk−fik∂G∂L|L=Li,t=tk](9)

where G(Li,tk), f(Li,tk), h(Li,tk,f(Li,tk)) are abbreviated as fik, Gik, hik, respectively. The crystal growth rate G(L,t) and the source term h(L,t,f) are usually known in Eq. (9). Thus, the problem is to find a method to discretize the internal crystal-size derivative ∂f∂L|L=Li,t=tk.

The PDDO method is used to discretize the internal crystal-size derivative. Taking the polynomial degree N = 1, we have Eq. (10).

∂f∂L|L=Li,t=tk=∑j=1N(i)fjk⋅g11(Lj−Li)ΔL(10)

where the PD function is constructed using Eq. (5).

The direction of characteristics must be considered when solving hyperbolic equations. Baryt et al. [29] constructed an upwind weight function according to the direction of the incoming flow, when Gik≥0 and the weight function is as described in Eq. (11a) (Fig. 2).

Figure 2: Diagram of the upwind weight function

ω(ξ)={0ξ>0e−(2ξ/δ)2ξ≤0(11a)

Again, when Gik<0, the weight function is expressed as Eq. (11b).

ω(ξ)={e−(2ξ/δ)2ξ≥00ξ<0(11b)

2.4 Stability Analysis of the EE–PDDO Method

The simplest form of the PBE is shown in Eq. (12).

∂f∂t+G∂f∂L=0(12)

G is the growth rate and is positive, and the source term is ignored. The corresponding EE–PDDO form is expressed as Eq. (13).

fik+1=fik−G⋅Δt⋅∑j=1N(i)fjk⋅g11(Lj−Li)ΔL(13)

Taking the horizon δ=ΔL, Eq. (13) can be written as Eq. (14).

fik+1=fik−G⋅Δt⋅[fi−1k⋅g(−ΔL)ΔL+fik⋅g(0)ΔL+fi+1k⋅g(ΔL)ΔL](14)

We obtain the coefficients for the PD function using Eqs. (6), (7) and (11) as shown in Eq. (15).

{a01=1ΔL2a11=1+e4ΔL3(15)

The numerical scheme of Eq. (12) is obtained as shown in Eq. (16).

fik+1=(1−λ)fik+λfi−1k(16)

with λ=GΔtΔL as the Courant number.

We now use the Fourier analysis to show the stability [32]. Let fik=WkeIrLi (where I=−1 and r is the positive constant). Substituting it in Eq. (16), we obtained the amplification factor as shown in Eq. (17).

κ=Wk+1Wk=1−λ[1−cos⁡(rΔL)+Isin⁡(rΔL)](17)

The stability requirement is |κ|2≤1. So stability is determined by using Eq. (18).

0<λ≤1(18)

Similarly, when one set of the horizon is δ=2ΔL, the numerical scheme is expressed as Eq. (19).

fik+1=(1−0.9032λ)fik+0.8065λfi−1k+0.0968λfi−2k(19)

The stability condition is shown in Eq. (20).

0<λ≤min{−0.3872cos2⁡(rΔL)+1.6130cos⁡(rΔL)−2−0.3497cos2⁡(rΔL)−1.3007cos⁡(rΔL)+1.6504}≈1.2(20)

2.5 Other Numerical Schemes

Herein we discuss internal crystal-size discretization, for which the treatment of the time derivative is the same as that in the EE method.

2.5.1 Second-Order Upwind Scheme

In the calculation of ∂∂L(Gk(L)fk(L)), we used the second-order upwind scheme shown in Eq. (21) [29].

∂(Gk(L)fk(L))∂L|L=Li={3(Gf)ik−4(Gf)i−1k+(Gf)i−2k2ΔLGik≥03(Gf)i+2k−4(Gf)i+1k+(Gf)ik2ΔLGik<0(21)

2.5.2 HR-Van Scheme

In the calculation of ∂∂L(Gk(L)fk(L)), the HR-van scheme shown in Eq. (22) was employed [14].

∂(Gk(L)fk(L))∂L|L=Li={1ΔL[Gik(fik+12ϕ(ri)(fi+1k−fik))−Gi−1k(fi−1k+12ϕ(ri−1)(fik−fi−1k))]Gik≥01ΔL[Gi+1k(fi+1k+12ϕ(ri+1)(fi+2k−fi+1k))−Gik(fik+12ϕ(ri)(fi+1k−fik))]Gik<0(22)

where ri=fik−fi−1kfi+1k−fik and ϕ(r)=|r|+r1+|r| with ϕ(r) being the flux limiter proposed by Van [14].

3  Numerical Experiments

We compared the accuracy of each numerical method by defining the following error shown in Eq. (23).

L2−error=∑i=0ML(fie−fin)2ΔL/∑i=0ML|fie|ΔL(23)

where fie is the exact solution at the node Li.

3.1 Size-Independent of Growth in a Batch Process

The first example is the crystallization in a batch process where the crystal growth G is independent of crystal size L. We neglected aggregation, nucleation, and breakage, and set h(L,t,f)=0. The corresponding PBE is written as Eq. (12).

3.1.1 Single-Mode with Non-Smooth Distribution

The initial crystal-size function satisfies the non-smooth distribution shown in Eq. (24) [33].

f(L,0)={1×1010,10≤L≤300,else(24)

In the calculation, the domain is set as [0,100µm] and crystal growth rate is G=1.0µm/s. The exact solution of this example is f(L,t)=f(L−Gt,0). The other parameters for the calculation are ML=200, tend=60 s, and ΔL=0.5 µm.

Fig. 3 shows the comparison of the exact solution with the results obtained by four numerical schemes at t = 60 s. The time step is chosen as Δt = 0.05 s in the first-order upwind, second-order upwind and HR-van schemes owing to the differences in the stability condition. The time step is chosen as Δt=ΔL=0.5 with δ=ΔL (m = 1) and Δt=1.1ΔL=0.55 with δ=2ΔL (m = 2) in the EE–PDDO scheme. The selection of Δt = 0.05 will be discussed later. Fig. 3a shows the predicted crystal size distributions at the whole crystal size. Compared to the initial crystal-size distribution, the solution shifts to right with the velocity G. Fig. 3a shows that the crystal size distribution predicted by the first-order upwind scheme has great dissipation, while the crystal size distribution predicted by the second-order upwind scheme exhibits significant dispersion. In the next example, we will show the L2 error of these upwind schemes in which the second-order upwind scheme is slightly better than that of the first-order upwind scheme. Therefore, in the following examples, we will only analyze and discuss the results of the second-order upwind scheme. Fig. 3b shows the predicted crystal size distributions at the local crystal size. Fig. 3b shows that the results obtained by the HR-van scheme have been greatly improved, and the obtained crystal size distribution is in good agreement with the exact solution. However, there are some errors where the crystal size distribution changes greatly. For the EE–PDDO method, the result of m = 1 is a reproduction of the exact solution, while the result of m = 2 also shows a good agreement with the exact solution. Although, the result of m = 2 of the EE–PDDO method has some dissipation where the crystal size distribution changes greatly, the obtained solution does not show oscillation. For the PBEs, it is shown that the first-order upwind will generate numerical dissipation, while the second-order upwind will generate numerical dispersion near discontinuous or steep gradients. The HR-van method is better than these two because it uses the flux limiters. From a perspective of finite difference method, in order to improve accuracy, a wider range of templates can be used, such as ENO and WENO formats. In this study, the EE–PDDO method effectively handles the discontinuity by introducing the upwind weight function.

Figure 3: Comparison of analytical and numerical crystal size distribution. (a) In the whole computational zone (b) zoomed figure

Fig. 4 shows the evolution of L2 error of different numerical schemes with the time steps at t = 60 s. The EE–PDDO method with m = 2 has the larger stability region, which means Δt can be the larger. The second-order upwind scheme has the smaller stability region. The EE–PDDO method with m = 1 has the smaller error at Δt=ΔL=0.5. The EE–PDDO method with different m values has some dissipation with the decrease of Δt. The HR-van method has good results under different values of Δt.

Figure 4: Evolution of the L2 error of different numerical schemes with the time steps

The method of directly solving the crystal size distribution function is often considered to be too expensive to couple with Computational fluid dynamics (CFD) [18,34]. Therefore, we now focus on the CPU time for these methods. Fig. 5 shows the CPU time of three methods with different numbers of internal crystal-size points. Since the stability conditions of these methods are different, the time step is selected as mentioned here: (i) Δt=0.1ΔL for the second-order upwind scheme; (ii) Δt=0.5ΔL for the HR-van scheme; and (iii) Δt=ΔL for the EE–PDDO method. The computations are carried out on a personal computer with 32 GB of random access memory and a 16-core 3.4 GHz processor. In this example, the computational time by the EE–PDDO method with different m is generally the same, which is less than that of the second-order upwind scheme and HR-van scheme.

Figure 5: Evolution of the required CPU time with the number of internal crystal-size points for different schemes

Therefore, it is best to select the EE–PDDO method with m = 1 and Δt = ΔL in the case of size-independent growth from the perspective of accuracy and efficiency. When λ=1, Eq. (16) deduces to fik+1=fi−1k. This means that the solution propagates along the characteristic.

3.1.2 Multi-Modal with Non-Smooth Distribution

The initial crystal size function of the multi-modal case is assumed as Eq. (25) [33]:

f(L,0)={0,L≤2.01×1010,2<L≤100,10<L≤181010cos2⁡(π(L−26)/64),18<L≤340,34<L≤4210101−(L−50)2/64,42<L≤580,58<L≤6610e−(L−70)2/2σ2,66<L≤740,L>74(25)

The computational zone is set as [0,100 µm] and the crystal growth rate G is 0.1 µm/s. The exact solution of this example is f(L,t)=f(L−Gt,0).

The other parameters are set as shown here: ML=100, tend=100 s, ΔL=1 µm, and σ=0.778. The time step is chosen as mentioned here: Δt = ΔL in the first-order upwind scheme, second-order upwind scheme and the HR-van scheme, and Δt=ΔL/G in the EE–PDDO method.

Fig. 6 shows the comparison of the analytical solution with numerical solutions obtained by different schemes at t = 100 s. The result of the EE–PDDO method with m = 1 is also a reproduction of the exact solution, while the result of the EE–PDDO method with m = 2 is similar to that of the HR-van method. Table 1 lists the L2 error and the CPU time for different schemes. From the perspective of accuracy and efficiency, the EE–PDDO method with m = 1 is also the best choice. Also, as shown in Table 1, the accuracy of the second-order upwind is sightly better than that of the first-order upwind. Therefore, we shall not show the performance of the first-order upwind scheme in the following examples.

Figure 6: Comparison of analytical and numerical crystal size distribution

3.2 Size-Dependent Growth in a Batch Process

We have considered the example of G dependent on L in a batch process. Aggregation, nucleation, and breakage are not considered. Thus, h(L,t,f)=0. The relationship between G and L satisfies the relation G(L,t)=G0L, where G0 is a constant. Therefore, the PBE can be expressed as Eq. (26).

∂f∂t+G0L∂f∂L+fG0=0(26)

The initial size distribution satisfies Eq. (27).

f(L,0)=N0L¯exp⁡(−LL¯)(27)

where N0,L¯ are constants. The analytical solution for this example is as shown in Eq. (28) [14].

f(L,t)=N0L¯exp⁡(−LL¯e−G0t−G0t)(28)

We set the parameters as L¯=0.01 µm3, N0=1, G0=0.1 µm3/s, ML=1000, ΔL=0.001 µm, and Δt=0.001 s.

Fig. 7 shows the comparison of the numerical solutions of different schemes with the analytical solution at t = 4 s. Here, the normal weight in the EE–PDDO method is the weight function considered as ω(ξ)=e−(2ξ/δ)2. There is no significant difference between the solutions obtained by these numerical schemes in Fig. 7a. Compared to the initial crystal-size distribution, the solution shifts to the right, but instead of shifting, it exhibits an exponential decrease. Sightly differences of the crystal-size distribution are found on the left side with smaller crystal-size. Fig. 7b shows the comparison of numerical results in the vicinity of L=0 µm. The results of the second-order upwind scheme and HR-van scheme show some deviations near the left boundary. However, the results obtained by the EE–PDDO method show good agreement with the analytical solution irrespective of m and the weight function.

Figure 7: Comparison of analytical and numerical crystal size distribution. (a) In the whole computational zone (b) zoomed figure

Fig. 8 shows the L2 error with the number of internal crystal-size points for different schemes at t=4 s. The errors obtained by the second-order upwind and HR-van schemes are relatively large, and their convergent rates are low. The errors of the EE–PDDO scheme for different m values are almost the same; m = 1 is slightly better than that of m = 2. Fig. 8 shows that the results of the normal weight function of the EE–PDDO method is better than that of the upwind weight function since they have smaller errors and higher convergence rates. The convergence rates of the EE–PDDO method with the upwind weight function and normal weight function are 1 and 2, respectively.

Figure 8: Evolution of the L2 error with the number of spatial points for different schemes

The PBE becomes a convection-reaction equation in the size-dependent growth case due to the dependency of G and L. Thus, its hyperbolic property decreases. The equation can be solved using regular numerical schemes. Therefore, the EE–PDDO method with the normal weight function is the most appropriate choice for the simulation of size-dependent growth cases.

3.3 Nucleation and Size-Independent Growth in a Continuous Process

The PBE is expressed as Eq. (29) considering the nucleation and size-independent growth in a continuous process.

∂f∂t+G0∂f∂L=h(L,t,f)(29)

where the source term is h(L,t,f)=B0δ(L−L0), B0 is the constant, and δ(L) is the Dirac function. The analytical solution for this case is as shown in Eq. (30) [18,35].

f(L,t)=B0G0u(t−LG0)(30)

where u(x) is a unit step function and satisfies u(x)={1,x≥00,x<0.

The initial distribution is f(L0,0)=B0G0. The computational domain is set as [0,2 µm]. The other parameters are set as: B0=20/µm3/s, G0=1.8 µm/s, L0=5×10−4 µm, ΔL=0.01 µm, and Δt=ΔL/G0. The time steps for the second-order upwind and HR-van schemes are chosen as Δt = ΔL/10. The PBE is similar to that used in Section 3.1. The normal weight function leads to instability in the EE–PDDO method. Therefore, only the upwind weight function is discussed here.

Fig. 9 shows the comparison of analytical and numerical crystal-size distributions. The EE–PDDO method with m = 1 has a strong ability to accurately capture the discontinuous states. Although for m = 2 in EE–PDDO method, the solution is poorer than that of the HR-van method, the solution is non-oscillating. Table 2 lists the L2 errors and CPU time for different schemes. The advantages of the EE–PDDO method are high accuracy and efficiency.

Figure 9: Comparison of analytical and numerical crystal size distribution: (a) t = 0.5 s and (b) t = 1 s

3.4 Nucleation and Size-Dependent Growth in a Continuous Process

The corresponding PBE for the nucleation and size-dependent growth in a continuous process is given by Eq. (31).

∂f∂t+G∂f∂L+f∂G∂L=h(L,t,f)(31)

The source term is assumed as h(L,t,f)=B0δ(L−L0)−f(L)τ, where B0,τ are constants, and δ(L) is the Dirac function. The relationship between G and L is G(L)=G0(1+γL)z, where G0, γ, and z are constants. The analytical solution is shown in Eq. (32) [14].

limt→∞f(L)=B0G0(1+γL)−zexp⁡(1−(1+γL)1−zG0τγ(1−z))(32)

The initial distribution was assumed to be given by Eq. (32) and the simulation ended at t = 10000 s. The computational domain is set as [0, 4 μm]. The other parameters are set as shown here: B0=2×10−10/µm3/s, τ=100 s, G0=0.00168 µm/s, z=0.3, γ=1, ML=400, and Δt=0.01 s.

Fig. 10 shows the comparison between the analytical and numerical solutions obtained by the three schemes. Fig. 10a indicates that all the numerical results are in good agreement with the exact solutions except for the second-order upwind scheme. This deviation is clear in Fig. 10b.

Figure 10: Comparison of analytical and numerical crystal size distribution. (a) In the whole computational zone (b) zoomed figure

Table 3 listed the L2 errors and CPU time for different schemes. The errors of the HR-van scheme and EE–PDDO method with the upwind weight function of m = 1 are close. The corresponding results are poorer than those obtained by the EE–PDDO method with the normal weight function of m = 1. In the EE–PDDO method, the accuracy can be improved by increasing the order of polynomials. When the polynomial order is increased (N = 3), highly accurate numerical results are obtained. The accuracy of EE–PDDO method is determined by parameters. In general, increasing the value of N will improve the accuracy. However, from the viewpoint of numerical experiment, this parameter should be carefully taken. In this problem N=1 is a good choice. According to the suggestion of Madenci et al. [24], the range of m is N≤m≤N+2. We found that it is better when m=N. In this example, we appropriately increased the value of N, and achieved good results when N=3,m=3. However, the choice of N and m should be depended on the equation and experimental experience.

In this size-dependent growth case, the EE–PDDO method with the normal weight function is recommended.

4  Conclusions

In this paper, the EE–PDDO method for solving 1D-PBE is presented. We studied four different cases: size-independent growth, size-dependent growth in a batch process, nucleation and size-independent growth, and nucleation and size-dependent growth in a continuous process. These cases were solved by the EE–PDDO, second-order upwind, and HR-van methods. The stability of the EE–PDDO method, choice of the weight function in the PDDO method, and optimal time step were discussed. The conclusions are as follows:

(1)   The EE–PDDO method has high accuracy and efficiency in solving the 1D-PBE. This method is characterized by non-oscillation and high accuracy, especially in the discontinuous and sharp crystal size distribution.

(2)   The stability condition of the EE–PDDO method of the simplest PBE with the size-independent growth case is 0<λ≤1 when the polynomial degree N = 1 and horizon δ = ΔL (m = 1) for the upwind weight function. The stability condition is 0<λ≤1.2 for upwind weight function when the polynomial degree is N = 1 and the horizon is δ = 2ΔL (m = 2).

(3)   The upwind weight function should be used in the PDDO method to reach stability in the size-independent growth case. The time step size is recommended as Δt=ΔL/G or the Courant number λ=GΔtΔL=1. In the size-dependent growth case, results obtained with the normal weight in the PDDO scheme are better than that of the upwind weight function. Therefore, in this case, the normal weight in the PDDO is recommended.

We present the EE–PDDO method for solving 1D PBE. The detailed implementation and the corresponding techniques were discussed here. This method is also expected to be used for solving more complex crystallization processes.

Acknowledgement: The authors acknowledge the School of Mathematics and Statistics, Henan University of Science and Technology, for allowing the use of their high-performance computing resources.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: The authors confirm contribution to the paper as follows: study conception and design: Chunlei Ruan, Cengceng Dong; data collection: Kunfeng Liang, Zhijun Liu, Xinru Bao; analysis and interpretation of results: Chunlei Ruan, Cengceng Dong; draft manuscript preparation: Chunlei Ruan, Cengceng Dong. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available on request from the corresponding author upon reasonable request.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Tandon, P., Rosner, D. E. (1999). Monte Carlo simulation of particle aggregation and simultaneous restructuring. Journal of Colloid & Interface Science, 213(2), 273–286. [Google Scholar]

2. Meimaroglou, D., Kiparissides, C. (2007). Monte carlo simulation for the solution of the bi-variate dynamic population balance equation in batch particulate systems. Chemical Engineering Science, 62(18–20), 5295–5299. [Google Scholar]

3. Caliane, B., Maria, R., Rubens, M. (2007). Considerations on the crystallization modeling: Population balance solution. Computers & Chemical Engineering, 31(3), 206–218. [Google Scholar]

4. Li, D. Y., Li, Z. P., Gao, Z. M. (2019). Quadrature-based moment methods for the population balance equation: An algorithm review. Chinese Journal of Chemical Engineering, 27(3), 483–500. [Google Scholar]

5. Kumar, S., Ramkrishna, D. (1997). On the solution of population balance equations by discretization—I. A fixed pivot technique. Chemical Engineering Science, 52, 1311–1332. [Google Scholar]

6. Li, S., Liu, W. K. (1999). Reproducing kernel hierarchical partition of unity, part I—Formulation and theory. International Journal for Numerical Methods in Engineering, 45(3), 251–288. [Google Scholar]

7. Li, S., Liu, W. K. (1999). Reproducing kernel hierarchical partition of unity, part II—Applications. International Journal for Numerical Methods in Engineering, 45(3), 289–317. [Google Scholar]

8. Xu, Z., Zhao, H., Zheng, C. (2014). Fast Monte Carlo simulation for partice coagulation in population balance. Journal of Aerosol Science, 74(4), 11–25. [Google Scholar]

9. Xu, Z., Zhao, H., Zheng, C. (2015). Accelerating population balance-Monte Carlo simulation for coagulation dynamics from the markov jump model, stochastic algorithm and GPU parallel computing. Journal of Computational Physics, 281, 844–863. [Google Scholar]

10. John, V., Angelov, I., Oncul, A., Thevenin, D. (2007). Techniques for the reconstrucion of a distribution from a finite number of its moments. Chemical Engineering Science, 62(11), 2890–2904. [Google Scholar]

11. Lage, P. L. C. (2011). On the representation of QMOM as a weighted residual method—The dual quadrature method of generalized moments. Computers & Chemical Engineering, 35(11), 2186–2203. [Google Scholar]

12. Grosch, R., Briesen, H., Marquardt, M., Wulkow, M. (2007). Generalization and numerical investigation of QMOM. AICHE Journal, 53(1), 207–227. [Google Scholar]

13. Marchisio, D. L., Fox, R. O. (2005). Solution of population balance equations using the direct quadrature method of moments. Journal of Aerosol Science, 36(1), 43–73. [Google Scholar]

14. Gunawan, R., Fusman, I., Braatz, R. D. (2004). High resolution algorithms for multidimensional population balance equations. AICHE Journal, 50(11), 2738–2749. [Google Scholar]

15. Noor, S., Qamar, S. (2015). The space time CE/SE method for solving one-dimensional batch crystallization model with fines dissolution. Chinese Journal of Chemical Engineering, 23(2), 337–341. [Google Scholar]

16. Nicmanis, M., Hounslow, M. J. (1998). Finite-element methods for steady-state population balance equations. AICHE Journal, 44(10), 2258–2272. [Google Scholar]

17. Sandu, A. (2006). Piecewise polynomial solutions of aerosol dynamic equation. Aerosol Science and Technology, 40(4), 261–273. [Google Scholar]

18. Alzyod, S., Charton, S. (2020). A meshless radial basis method (RBM) for solving the detailed population balance equation. Chemical Engineering Science, 228(31), 1–32. [Google Scholar]

19. Ruan, C. (2019). Chebyshev spectral collocation method for population balance equation in crystallization. Mathematics, 7(317), 1–12. [Google Scholar]

20. Majumder, A., Kariwala, V., Ansumali, S., Rajendran, A. (2012). Lattice Boltzmann method for multi-dimensional population balance models in crystallization. Chemical Engineering Science, 70(18), 121–134. [Google Scholar]

21. Gunzburger, M., Lehoucq, R. B. (2010). A nonlocal vector calculus with application to nonlocal boundary value problems. Multiscale Modeling & Simulation, 8(5), 1581–1598. [Google Scholar]

22. Yan, J., Li, S., Kan, X., Zhang, A. M., Lai, X. (2020). Higher-order nonlocal theory of updated lagrangian particle hydrodynamics (ULPH) and simulations of multiphase flows. Computer Methods in Applied Mechanics and Engineering, 368(2), 113176. [Google Scholar]

23. Yu, H., Li, S. (2021). On approximation theory of nonlocal differential operators. International Journal for Numerical Methods in Engineering, 122(23), 6984–7012. [Google Scholar]

24. Madenci, E., Baryt, A., Dorduncu, M. (2019). Peridynamic differential operator for numerical analysis. Switzerland: Springer Nature Switzerland. [Google Scholar]

25. Madenci, E., Baryt, A., Futch, M. (2016). Peridynamic differential operator and its applications. Computer Methods in Applied Mechanics and Engineering, 304, 408–451. [Google Scholar]

26. Dorduncu, M. (2019). Stress analysis of sandwich plates with functionally graded cores using peridynamic differential operator and refined zigzag theory. Thin-Walled Structures, 146, 106468. [Google Scholar]

27. Gao, Y., Oterkus, S. (2019). Nonlocal modeling for fluid flow coupled with heat transfer by using peridynamic differential operator. Engineering Analysis with Boundary Elements, 105(2), 104–121. [Google Scholar]

28. Li, Z. Y., Huang, D., Xu, Y. P. (2021). Nonlocal steady-state thermoelastic analysis of functionally graded materials by using peridynamic differential operator. Applied Mathematical Modelling, 93, 294–313. [Google Scholar]

29. Baryt, A., Madenci, E., Haghighat, E. (2022). On the solution of hyperbolic equations using the peridynamic differential operator. Computer Methods in Applied Mechanics and Engineering, 391, 114574. [Google Scholar]

30. Ma, D. L., Tafti, D. K., Braatz, R. D. (2002). High-resolution simulation of multidimensional crystal growth. Industrial & Engineering Chemistry Research, 41(25), 6217–6223. [Google Scholar]

31. Silling, S. A., Epton, M., Weckner, O. (2007). Peridynamic states and constitutive modeling. Journal of Elasticity, 88(2), 151–184. [Google Scholar]

32. Mark, H. (2011). Introduction to numerical methods in differential equaitons. Berlin: Springer. [Google Scholar]

33. Qamar, S., Elsner, M. P., Angelov, I. A., Warnecke, G., Seidel-Morgenstern, A. (2006). A comparative study of high resolution schemes for solving population balances in crystallization. Computers & Chemical Engineering, 30(6–7), 1119–1131. [Google Scholar]

34. Cheng, J. C., Li, Q., Yang, C., Zhang, Y. Q., Mao, Z. S. (2018). CFD-PBE simulation of a bubble column in OpenFOAM. Chinese Journal of Chemical Engineering, 26(9), 1773–1784. [Google Scholar]

35. Hounslow, M. J. (1988). A discretized population balance for nucleation growth, and aggregation. AIChE Journal, 34, 1576–1589. [Google Scholar]