Next Article in Journal
Effects of Ultrasonic Treatment on Grain Refinement and Gas Removal in Magnesium Alloys
Next Article in Special Issue
Characterization and Formation Mechanism of Ag2MoO4 Crystals via Precipitation Method: Influence of Experimental Parameters and Crystal Morphology
Previous Article in Journal
Growth, Structure, and Electrical Properties of AgNbO3 Antiferroelectric Single Crystal
Previous Article in Special Issue
Effects of Dipentaerythritol and Cellulose as Additives on the Morphology of Pentaerythritol Crystals
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Euler’s First-Order Explicit Method–Peridynamic Differential Operator for Solving Two-Dimensional Population Balance Equations in Crystallization

1
School of Mathematics & Statistics, Henan University of Science & Technology, Luoyang 471023, China
2
Longmen Laboratory, Luoyang 471023, China
*
Author to whom correspondence should be addressed.
Crystals 2024, 14(3), 234; https://doi.org/10.3390/cryst14030234
Submission received: 1 February 2024 / Revised: 22 February 2024 / Accepted: 25 February 2024 / Published: 28 February 2024
(This article belongs to the Special Issue Crystallization Process and Simulation Calculation, Second Edition)

Abstract

:
The population balance equations (PBEs) serve as the primary governing equations for simulating the crystallization process. Two-dimensional (2D) PBEs pertain to crystals that exhibit anisotropic growth, which is characterized by changes in two internal coordinates. Because PBEs are the hyperbolic equations, it becomes imperative to establish a high-resolution scheme to reduce numerical diffusion and numerical dispersion, thereby ensuring accurate crystal size distribution. This paper uses Euler’s first-order explicit (EE) method–Peridynamic Differential Operator (PDDO) to solve 2D PBE, namely, the EE method for discretizing the time derivative and the PDDO for discretizing the internal crystal-size derivative. Five examples, including size-independent growth with smooth and non-smooth distributions, size-dependent growth, nucleation, and size-independent/dependent growth for batch crystallization are considered. The results show that the EE–PDDO method is more accurate than the HR method and that it is as good as the fifth-order Weighted Essential Non-Oscillatory (WENO) method in solving 2D PBE. This study extends the EE–PDDO method to the simulation of 2D PBE, and the advantages of the EE-PDDO method in dealing with discontinuous and sharp front problems are demonstrated.

1. Introduction

Crystallization operation is widely used in chemical and pharmaceutical industries. The study of crystallization kinetics forms the foundation for developing crystallization theory and offers guidance in the design of crystallizers. The population balance equation (PBE) is a widely accepted crystallization kinetics model. Therefore, it is of great significance to solve the PBE accurately. Generally, the population balance model only considers one internal variable, which describes how the crystal grows isotropically and forms simple isotropic shapes (spherical, cubic, etc.) [1]. However, for the anisotropic, crystalline, organic products, the population balance model, which has two internal variables, is needed. Our objective is to present the high-resolution scheme for solving two-dimensional (2D) PBEs.
At present, a lot of numerical methods have been proposed for solving PBEs. They can be mainly classified into five categories [2,3,4]: (1) the method of moments [5,6]; (2) the method of characteristic [7,8]; (3) the method of weight residuals/orthogonal collocation [9]; (4) the Monte Carlo method [10,11]; and (5) the finite-difference scheme/discrete population balances [12,13]. Nguyen et al. [14] used the Extended Quadrature Method of Moments (EQMOM) method and Two-Size Moment (TSM) method to solve PBE. They reported that compared with the sectional method, the EQMOM and TSM showed their ability to describe accurately the fine-particle population with a much lower number of variables. Qamar et al. [7] proposed the method of characteristics to solve the population balance model in two coupled priority crystallizers. They found that the model becomes more complex and the solution exhibits strong discontinuity when PBE is coupled with computational fluid dynamics. Ruan [15] introduced the Chebyshev spectral collocation method to solve PBE. They suggested that the diffusive term should be added in cases of size-independent growth in which the PBE is a convection equation. Ma et al. [16] used a Monte Carlo simulation to solve PBEs and analyzed the factors that affect the size distribution of bubbles on distillation trays. For the discrete methods, Qamar and collaborators used high-resolution, finite-volume schemes [17] and moving grid technology [4] to solve PBEs. They reported that the schemes did a good job when capturing discontinuities in crystal size distribution. Ruan et al. [18,19] introduced the fifth-order Weighted Essential Non-Oscillatory (WENO) method to solve one-dimensional (1D) and 2D PBEs. They found that the fifth-order WENO method is valid for solving PBEs and that accurate crystal size distributions can be obtained. In our previous work [20], Euler’s first-order explicit (EE) method–Peridynamic Differential Operator (PDDO) was used to solve 1D PBEs. Compared with the second-order upwind [21] and HR-van [22] methods, the EE–PDDO method shows high accuracy in predicting crystal size distributions.
In this study, we attempt to extend the EE–PDDO method to the simulation of 2D PBEs. The PDDO method is an effective numerical method for solving differential equations that has been developed in recent years by Madenci et al. [23,24,25] based on the Peridynamics (PD) theory [26,27,28]. The PDDO is a nonlocal differential operator, and the main advantage of the PDDO method is that it can deal with the discontinuities or singularities without additional processing. Therefore, based on these advantages, there has been some work using the PDDO method in engineering. For example, Dorduncu [29] used the PDDO method to analyze the stress of sandwich plates with functionally graded cores; Gao et al. [30] used the PDDO method in the simulation of incompressible fluid with heat and mass transfer; and Li et al. [31] used the PDDO method to analyze of thermodynamic problems of functionally graded materials. Because the crystal size distributions in the crystallization process are very sharp and sometimes discontinuous, the PDDO method is suitable for solving PBEs.
The paper is organized as follows: in Section 2, the PBE and EE–PDDO method are presented, together with the stability of the EE–PDDO method; in Section 3, five examples, including size-independent growth with smooth and non-smooth distributions, size-dependent growth, nucleation, and-size-independent/dependent growth for batch crystallization, are considered, and the results of the EE–PDDO method are compared with those of the HR method and the fifth-order WENO method; finally, the conclusions are drawn in Section 4.

2. Background Theory

2.1. Population Balance Equation

The 2D PBE [2,5,32,33] can be written as follows:
f ( r 1 , r 2 , t ) t + j = 1 2 { G j ( r 1 , r 2 , c ( t ) , T ( t ) ) f ( r 1 , r 2 , t ) } r j = h ( f ( r 1 , r 2 , t ) , c ( t ) , T ( t ) )
where r 1 , r 2 are two characteristic lengths for crystals, f ( r 1 , r 2 , t ) is the crystal size distribution function, G 1 , G 2 are the growth rate in the r 1 , r 2 direction, c is a time-dependent solute concentration, T is a time-dependent temperature, and h is the source term, including crystal aggregation, nucleation, and breakage. This study assumes that the mold is well-mixed and that the crystal size distribution is independent of spatial coordinates. Because the PBEs have strong nonlinearity and the expressions of G 1 , G 2 and h are sometimes very complex, to find the analytical solutions is difficult. Numerical simulation becomes a powerful tool in solving PBEs.

2.2. The PDDO Method

Silling and collaborators [34] proposed the PD theory. As shown in Figure 1, in an interaction domain H x , point x and point x interact with each other. Here, ξ = x x represents the initial relative position between the two points. The interaction domain H x is also called the family of point x . Each point has its own family members. The size and shape of each family can be different (as shown in Figure 1, the size and shape of H x and H x are different). All of the family members in the interaction domain H x interact with the point x , and the degree of nonlocal interaction among the family members is represented by the weight function ω ( ξ ) . The commonly used shape of interaction domain H x in 2D is the square at which the point x is at the center. The size of the interaction domain H x is determined by the horizon δ . The horizon has the following expression: δ = m Δ x , where m is the integer parameter and Δ x is the grid size. One can refer to Ref. [24] for more details.
Based on the PD theory, the PDDO method was introduced by Madenci et al. [24]. We now consider the Taylor’s series expansion of f ( x ) = f ( x + ξ ) in a 2D scalar function, which is
f ( x + ξ ) = n 1 = 0 N n 2 = 0 N n 1 1 n 1 ! n 2 ! ξ 1 n 1 ξ 2 n 2 n 1 + n 2 f ( x ) x 1 n 1 x 2 n 2 + R ( N , x )
where R ( N , x ) is the remainder. The orthogonal function g N p 1 p 2 ( ξ ) is introduced, which yields
p 1 + p 2 f ( x ) x 1 p 1 x 2 p 2 = H x f ( x + ξ ) g N p 1 p 2 ( ξ ) d V
in which p i ( i = 1 , 2 ) represents the partial derivative order with respect to variable x i . The orthogonality of PD function g N p 1 p 2 ( ξ ) can be expressed as:
1 n 1 ! n 2 ! H x ξ 1 n 1 ξ 2 n 2 g N p 1 p 2 ( ξ ) d V = δ n 1 p 1 δ n 2 p 2
where n i = 0 , , N , δ n p is the Kronecker symbol. The PD function can be constructed as
g N p 1 p 2 ( ξ ) = q 1 = 0 N q 2 = 0 N q 1 a q 1 q 2 p 1 p 2 ω q 1 q 2 ( ξ ) ξ 1 q 1 ξ 2 q 2
where ω q 1 q 2 ( ξ ) is the weight function related to the term ξ q .
The unknown coefficients a q 1 q 2 p 1 p 2 can be derived from the following equation:
q 1 = 0 N q 2 = 0 N q 1 A ( n 1 n 2 ) ( q 1 q 2 ) a q 1 q 2 p 1 p 2 = b n 1 n 2 p 1 p 2
where q i = 0 , , N . The coefficient matrix is defined as
A ( n 1 n 2 ) ( q 1 q 2 ) = H x ω q 1 q 2 ( ξ ) ξ 1 n 1 + q 1 ξ 2 n 2 + q 2 d V
The right-hand term is
b n 1 n 2 p 1 p 2 = n 1 ! n 2 ! δ n 1 p 1 δ n 2 p 2
For N = 1, the coefficient matrix A , the matrix of unknown coefficients a , and the vector b in Equation (6) can be written as [35]:
A = H x ω q 1 q 2 ( ξ ) 1 ξ 1 ξ 2 ξ 1 ξ 1 2 ξ 1 ξ 2 ξ 2 ξ 1 ξ 2 ξ 2 2 d V
a = a 00 00 a 00 10 a 00 01 a 10 00 a 10 10 a 10 01 a 01 00 a 01 10 a 01 01
and
b = 1 0 0 0 1 0 0 0 1
Determination of the unknown coefficients through a = A 1 b establishes the PD functions g N p 1 p 2 ( ξ ) .

2.3. Application of the EE–PDDO Method in 2D PBE

The computational domain is rectangular in two internal coordinates, which is [ a r 1 , b r 1 ] × [ a r 2 , b r 2 ] . We use a rectangular grid ( r 1 i , r 2 ( j ) ) , where r 1 i = a r 1 + ( i 1 2 ) Δ r 1 , ( i = 1 , 2 , , M r 1 ) and r 2 j = a r 2 + ( j 1 2 ) Δ r 2 , ( j = 1 , 2 , , M r 2 ) . The step sizes are Δ r 1 = ( b r 1 a r 1 ) / M r 1 and Δ r 2 = ( b r 2 a r 2 ) / M r 2 . The time points at which we will compute the solution are labeled as t 0 , t 1 , , t M t . We use a uniform time step Δ t ; therefore, the time points are t k = k Δ t , ( k = 0 , 1 , , M t ) . The time interval is 0 t t e n d . Therefore, Δ t and M t are connected through the equation Δ t = t e n d / M t .
Substituting the point ( r 1 i , r 2 j , t k ) into Equation (1), one obtains
f ( r 1 i , r 2 j , t k ) t + { G 1 ( r 1 i , r 2 j , c ( t k ) , T ( t k ) ) f ( r 1 i , r 2 j , t k ) } r 1 + { G 2 ( r 1 i , r 2 j , c ( t k ) , T ( t k ) ) f ( r 1 i , r 2 j , t k ) } r 2 = h ( f ( r 1 i , r 2 j , t k ) , c ( t k ) , T ( t k ) )
We use the EE method for discretizing the time derivative and obtain the following equation:
f i , j k + 1 = f i , j k + Δ t h i , j k G 1 f r 1 r 1 = r 1 i , t = t k f i , j k G 1 r 1 r 1 = r 1 i , t = t k G 2 f r 2 r 2 = r 2 j , t = t k f i , j k G 2 r 2 r 2 = r 2 j , t = t k
Here, f ( r 1 i , r 2 j , t k ) , G 1 ( r 1 i , r 2 j , c ( t k ) , T ( t k ) ) , G 2 ( r 1 i , r 2 j , c ( t k ) , T ( t k ) ) , and h ( f ( r 1 i , r 2 j , t k ) , c ( t k ) , T ( t k ) ) are abbreviated as f i , j k , G 1 , G 2 , h i , j k , respectively. The crystal growth rate G 1 ( r 1 , r 2 , c ( t ) , T ( t ) ) , G 2 ( r 1 , r 2 , c ( t ) , T ( t ) ) , and the source term h ( f ( r 1 , r 2 , t ) , c ( t ) , T ( t ) ) are usually given in specific cases. To discretize the internal crystal-size derivatives, f r 1 r 1 = r 1 i , t = t k and f r 2 r 2 = r 2 j , t = t k become important.
Using the PDDO method for discretizing the internal variables and taking the polynomial degree N = 1 , we obtain the following equation:
f r 1 r 1 = r 1 i , t = t k = y = 1 N ( i ) f k ( r 1 ( y ) ) g 1 10 ( ξ 1 ( y , i ) , ξ 2 ( y , i ) ) Δ A y f r 2 r 2 = r 2 j , t = t k = y = 1 N ( j ) f k ( r 2 ( y ) ) g 1 01 ( ξ 1 ( y , j ) , ξ 2 ( y , j ) ) Δ A y
where ξ 1 ( y , i ) = r 1 ( y ) r 1 ( i ) , ξ 2 ( y , j ) = r 2 ( y ) r 2 ( j ) , Δ A y = Δ r 1 Δ r 2 , and the PD function is constructed according to Equation (4).
When solving the hyperbolic equations, the direction of characteristics must be considered. Madenci et al. [21] reported an upwind weight function for solving hyperbolic equations. Here, the crystal growth rate is considered to be non-negative, which means G 1 k 0 and G 2 k 0 . We construct the following weight function (Figure 2):
ω ( ξ 1 , ξ 2 ) = e 4 ( ξ 1 σ 1 ) 2 e 4 ( ξ 2 σ 2 ) 2 ξ 1 , ξ 2 0 0 e l s e w h e r e .
or
ω ( ξ 1 , ξ 2 ) = 1 ξ 1 , ξ 2 0 0 e l s e w h e r e .

2.4. Stability Analysis of the EE–PDDO Method

We analyze the stability with the simplest form of PBE
f t + G 1 f r 1 + G 2 f r 2 = 0
Here, we ignore the impact of the source term and set the crystal growth rates G 1 , G 2 as constants. The corresponding EE–PDDO formula is as follows:
f i , j k + 1 = f i , j k Δ t [ G 1 y = 1 N ( i ) f k ( r 1 ( y ) ) g 1 10 ( ξ 1 ( y , i ) , ξ 2 ( y , i ) ) Δ A y + G 2 y = 1 N ( j ) f k ( r 2 ( y ) ) g 1 01 ( ξ 1 ( y , j ) , ξ 2 ( y , j ) ) Δ A y ]
Taking the horizon as δ 1 = Δ r 1 , δ 2 = Δ r 2 , Equation (18) can be written as Equation (19).
f i , j k + 1 = f i , j k Δ t G 1 Δ r 1 Δ r 2 q = 1 1 f i 1 , j + q k g 1 10 ( Δ r 1 , q Δ r 2 ) + Δ r 1 Δ r 2 q = 1 1 f i , j + q k g 1 10 ( 0 , q Δ r 2 ) + Δ r 1 Δ r 2 q = 1 1 f i + 1 , j + q k g 1 10 ( Δ r 1 , q Δ r 2 ) Δ t G 2 Δ r 1 Δ r 2 q = 1 1 f i + q , j 1 k g 1 01 ( q Δ r 1 , Δ r 2 ) + Δ r 1 Δ r 2 q = 1 1 f i + q , j k g 1 01 ( q Δ r 1 , 0 ) + Δ r 1 Δ r 2 q = 1 1 f i + q , j + 1 k g 1 01 ( q Δ r 1 , Δ r 2 )
By using Equations (9)–(11) and (15), we obtain the coefficients in PD function as follows:
a 00 10 = 0.9820 Δ r 1   2 Δ r 2 a 10 10 = 54.6013 Δ r 1   3 Δ r 2 a 01 10 = 0.0013 Δ r 1   2 Δ r 2   2 a 00 01 = 0.9820 Δ r 1 Δ r 2   2 a 10 01 = 0.0013 Δ r 1   2 Δ r 2   2 a 01 01 = 54.6013 Δ r 1 Δ r 2   3
Under the assumption of Δ r 1 = Δ r 2 = Δ r , G 1 = G 2 , we develop the numerical scheme of Equation (17) as
f i , j k + 1 = f i , j k λ [ 0.036 f i 1 , j 1 k 0.9634 f i 1 , j k 0.9634 f i , j 1 k + 1.968 f i , j k ]
where λ = G 1 Δ t Δ r is the Courant number.
The Fourier analysis [36] is used to explore the stability. In particular, the solution of Equation (21) has the form f i , j k + 1 = W k + 1 e α r 1   ( i ) I e β r 2   ( j ) I (where I = 1 and α , β are positive constants). Therefore, the amplification factor can be obtained as follows:
κ = W k + 1 W k = 1 λ [ 0.036 cos ( ( α + β ) Δ r ) + 0.036 I sin ( ( α + β ) Δ r ) 0.9634 ( cos ( α Δ r ) + cos ( β Δ r ) ) + 0.9634 I ( sin ( α Δ r ) + sin ( β Δ r ) ) + 1.968 ]
The stability requirement is κ 2 1 . Thus, the stability condition is shown in Equation (23).
0 < λ 0.5183
Similarly, taking the horizon as δ 1 = Δ r 1 , δ 2 = Δ r 2 , by using Equations (9)–(11) and (16), we obtain the coefficients in PD function as follows:
a 00 10 = 1 2 Δ r 1   2 Δ r 2 a 10 10 = 1 Δ r 1   3 Δ r 2 a 01 10 = 0 a 00 01 = 1 2 Δ r 1 Δ r 2   2 a 10 01 = 0 a 01 01 = 1 Δ r 1 Δ r 2   3
The corresponding numerical scheme is
f i , j k + 1 = f i , j k λ [ f i 1 , j 1 k + 0 f i 1 , j k + 0 f i , j 1 k + f i , j k ]
and the amplification factor is
κ = 1 λ + λ [ cos ( ( α + β ) Δ r ) I sin ( ( α + β ) Δ r ) ]
The stability requirement is κ 2 1 . Thus, the stability condition is shown in Equation (27).
0 < λ 1

2.5. Other Numerical Schemes

In this section, we focus on the discretization in the internal crystal-size derivative. In order to be comparable with various schemes, the time derivative is discretized using the EE method.

2.5.1. HR Scheme

The HR scheme for solving PBE Equation (1) can be written as [2,17,19,22]
f i , j k + 1 = f i , j k + Δ t h i , j k G 1 f i + 1 2 , j k G 1 f i 1 2 , j k Δ r 1 G 2 f i , j + 1 2 k G 2 f i , j 1 2 k Δ r 2
where
G 1 f i + 1 2 , j k = G 1 i + 1 2 , j k f i , j k + 1 2 φ r i + 1 2 k f i + 1 , j k f i , j k G 2 f i , j + 1 2 k = G 2 i , j + 1 2 k f i , j k + 1 2 φ r j + 1 2 k f i , j + 1 k f i , j k
in which r i + 1 2 k = f i , j k f i 1 , j k f i + 1 , j k f i , j k , r j + 1 2 k = f i , j k f i , j 1 k f i , j + 1 k f i , j k , φ ( r ) = r + r 1 + r , Δ t is the time step, and Δ r 1 , Δ r 2 are the spatial step size in r 1 , r 2 direction, respectively.

2.5.2. Fifth-Order WENO Scheme

We use a conservative form for discretizing Equation (1) in space and obtain [19]
f i , j k + 1 = f i , j k + Δ t h i , j k G 1 f r 1 ( r 1 ( i ) , r 2 ( j ) ) G 2 f r 2 ( r 1 ( i ) , r 2 ( j ) )
where G 1 f r 1 ( r 1 ( i ) , r 2 ( j ) ) and G 2 f r 2 ( r 1 ( i ) , r 2 ( j ) ) are constructed by the upwind idea. Taking the first one as an example,
( G 1 f ) / r 1 ( i , j ) = ( G 1 f r 1 ) i , j   ( G 1 ) i 0 ( G 1 f r 1 ) i , j +   ( G 1 ) i < 0
In the fifth-order WENO scheme, to determine the ( G 1 f r 1 ) i , j , it covers six cells and three stencils. We do not present the implementation of the WENO scheme here; one can refer to our previous work [19] for more details.

3. Numerical Experiments

In order to compare the accuracy of the numerical methods, we define the errors as Equation (31).
L 1 e r r o r = i = 1 N j = 1 M f i , j e f i , j Δ r 1 Δ r 2 L 2 e r r o r = i = 1 N j = 1 M f i , j e f i , j 2 Δ r 1 Δ r 2
where f i , j e is the exact solution.

3.1. Size-Independent Growth (Smooth Distribution)

We first consider a homogenous PBE with size-independent growth. In this case, crystal aggregation, nucleation, and breakage are neglected, which means h ( f ( r 1 , r 2 , t ) , c ( t ) , T ( t ) ) = 0 . The corresponding PBE is expressed as
f ( r 1 , r 2 , t ) t + [ G 1 f ] r 1 + [ G 2 f ] r 2 = 0
In this example, we set the computational domain as [ 0 , 1   μ m ] × [ 0 , 1   μ m ] , the crystal growth rate is G 1 = G 2 = 0.1   μ m / s , and the initial distribution is set as
f ( r 1 , r 2 , 0 ) = exp 100 × ( ( r 1 0.25 ) 2 + ( r 2 0.25 ) 2 )
The exact solution to this problem has the following form [19]:
f ( r 1 , r 2 , t ) = f ( r 1 G 1 t , r 2 G 2 t , 0 )
Figure 3 shows the comparison of crystal size distribution between the exact solutions and the numerical solutions obtained from different schemes at t = 5   s . The internal crystal-size step is set as Δ r 1 = Δ r 2 = 0.01 . The time step is chosen as Δ t = 0.001 in both the HR and fifth-order WENO schemes. For the EE–PDDO scheme, the solution is obtained by considering a uniform time step size of Δ t = Δ r 1 G 1 = 0.1 with the horizon δ = Δ r 1 ( m = 1 ) and the upwind weight function in Equation (16). From Figure 3, it is noticeable that the results obtained using the fifth-order WENO and PDDO schemes agree with the exact solution very well, while the HR scheme presents a lower peak and shows some diffusion. Table 1 shows the L 1 and L 2 errors calculated using different schemes. It is evident that the fifth-order WENO and PDDO methods have smaller errors, while HR has the biggest error. It should be mentioned that the EE–PDDO method varies a lot when using the different weight functions and integer parameter m. In this example, by comparing the errors of EE–PDDO for different m, the upwind weight function in Equation (16) is superior, and m = 1 is suggested. In the case of Δ t = Δ r 1 G 1 , which means λ = 1 , Equation (25) deduces to f i , j k + 1 = f i 1 , j 1 k . This implies that the solution propagates along the characteristics.

3.2. Size-Independent Growth (Non-Smooth Distribution)

We again consider a homogenous PBE with size-independent growth, and Equation (32) is used. This time, we consider the non-smooth distribution in initial size distribution, which is
f ( r 1 , r 2 , 0 ) = 1   0 . 1 r 1 , r 2 0.3 0   elsewhere .
The computational domain is set as [ 0 , 1   μ m ] × [ 0 , 1   μ m ] , and the crystal growth rate is chosen as G 1 = G 2 = 0.1   μ m / s . The solution is obtained by considering the internal crystal-size step Δ r 1 = Δ r 2 = 0.01 and the time step Δ t = 0.1 . The exact solution of this example is the same as Equation (34).
Figure 4 shows the comparison of crystal size distribution between the exact solutions and the numerical solutions obtained from different schemes at t = 6   s . In the implementation of the EE–PDDO method, the upwind weight function shown in Equation (16) is used. In this example, results obtained using HR and fifth-order WENO schemes show some dissipation, while the result obtained using the EE–PDDO method does not cause any deviation from the exact solution. The result indicates that the EE–PDDO method has high accuracy and less dissipation, and there is no oscillation at the sharp front. Table 2 shows the L 1 and L 2 errors calculated using the different schemes. It can be seen that the error of the EE–PDDO methods is smallest, while the HR method and fifth-order WENO method are comparable.
Table 3 shows the L 1 , L 2 errors and CPU time of the EE–PDDO scheme in different Δ t or different Courant number λ . The program runs on a personal computer with 6 cores and a 3.30 GHz processor. The results show that when the integer parameter m = 1 in the EE–PDDO scheme, the L 1 and L 2 errors reach the minimum at Δ t = Δ r 1 G 1 = 0.1 ( λ = 1 ), and the CPU time is the shortest. As mentioned in Section 3.1, when λ = 1 , the EE–PDDO scheme implies the solutions propagation along with the characteristics. Therefore, the EE–PDDO method perfectly replicates the exact solution.

3.3. Size-Dependent Growth

We now consider the PBE in the case of size-dependent growth. Similarly, the aggregation and breakage are neglected. The corresponding PBE is as follows:
f ( r 1 , r 2 , t ) t + [ G 1 r 1 f ] r 1 + [ G 2 r 2 f ] r 2 = 0
The growth rate is assumed to be increased linearly with the crystal size, which means G 1 r 1 = G 0 r 1 and G 2 r 2 = G 0 r 2 , with G 0 and G 0 as the constants.
The initial size distribution satisfies the following equation:
f ( r 1 , r 2 , 0 ) = N 0 L ¯ exp r 1 + r 2 L ¯
Here, L ¯ is the average size of the distribution, and N 0 is the total number of initial particles.
The exact solution of this example is [19]
f ( r 1 , r 2 , t ) = N 0 L ¯ exp r 1 L ¯ e G 0 t G 0 t r 2 L ¯ e G 0 t G 0 t
The parameters used in this example are as follows: the computational domain is [ 0 , 0.1   μ m ] × [ 0 , 0.1   μ m ] , the step sizes are Δ r 1 = 10 3   μ m and Δ r 2 = 10 3   μ m , respectively, L ¯ = 10 2   μ m , N 0 = 1 , G 0 = 0.1   μ m / s , G 0 = 0.2   μ m / s , and the simulation ended at t = 4   s .
Here, due to the dependency of crystal growth rate and crystal size, the hyperbolic property of PBE decreases. Correspondingly, the PBE becomes a convection–reaction equation, which can be solved through a regular numerical scheme. In this example, in the implementation of the EE–PDDO method, the upwind weight function shown in Equation (15) and the corresponding normal weight function with ω ( ξ 1 , ξ 2 ) = e 4 ( ξ 1 σ 1 ) 2 e 4 ( ξ 2 σ 2 ) 2 are used. The comparison of the final crystal size distributions between the exact solution and the numerical solutions obtained using different schemes are shown in Figure 5. In Figure 5d, normal weight function is used in the EE–PDDO method. The reason will be given later. From Figure 5, it is evident that the results obtained through EE–PDDO agree with the exact solution very well. Both EE–PDDO and fifth-order WENO methods have good resolution. The corresponding errors listed in Table 4 also prove that.
Table 4 lists the L 1 and L 2 errors for predicting crystal size distribution with different schemes. The results shows that the errors of the EE–PDDO scheme with normal weight function are the smallest, while the errors of the EE–PDDO scheme with upwind weight function are bigger than those of the fifth-order WENO method but smaller than those of the HR method. This is expected, because the PBE is a convection–reaction equation in this example and the reaction term reduces the hyperbolic property of the PBE. In hyperbolic equations, the feature velocity determines the speed at which information is spread. In the case of small feature velocity or a wave spread speed slower than the mesh size, the upwind scheme may introduce too conservative numerical dissipation, resulting in a decrease in the accuracy of the numerical results. Therefore, in the case of size-dependent growth, normal weight function is superior to upwind weight function.

3.4. Nucleation and Size-Independent Growth for 2D Batch Crystallization

We now consider the example of nucleation and size-independent growth. Here, we examine the potassium dihydrogen phosphate (KDP) crystal, whose shape is shown in Figure 6. The corresponding PBE can be written as
f ( r 1 , r 2 , t ) t + [ G 1 c , T f ] r 1 + [ G 2 c , T f ] r 2 = B 0 ( c , T ) δ ( r 1 ) δ ( r 2 )
where r 1 and r 2 are the width and length of the crystal, respectively, and the volume of the single crystal is given by Equation (40) [2].
V c = 1 3 r 1 3 + ( r 2 r 1 ) r 1 2
The crystal growth rates are independent of r 1 and r 2 and have the following forms [2]:
G 1 ( c ( t ) , T ( t ) ) = k g 1 S g 1 G 2 ( c ( t ) , T ( t ) ) = k g 2 S g 2
where S = c ( t ) c s a t ( T ) 1 is the relative supersaturation, c is the solute concentration, and c s a t is the saturated solute concentration, which can be expressed as [2]
c s a t ( T ) = 9.3027 × 10 5 T 2 9.7629 × 10 5 T + 0.2087
Here, T is the temperature, and it decreases linearly with the form
T ( t ) ( ° C ) = 33 t 720
The nucleation rate term is [2,22]
B 0 c t , T t = k b S b V t
Here, b is the nucleation index, and the total volume of the crystal V t is given by Equation (45) [37].
V t = 0 f r 1 , r 2 , t V c r 1 , r 2 d r 1 d r 2
The solution concentration is
d c t d t = ρ c 0 f r 1 , r 2 , t × 2 G 1 r 1 r 2 r 1 2 + G 2 r 1 2 d r 1 d r 2
The initial crystal size distribution is set as
f ( r 1 , r 2 , 0 ) = 3.4786 × 10 4 × r 1 2 + r 2 2 + 0.1363609 r 1 + r 2 26.5486 180   r 1 , r 2 212 0 elsewhere
The parameters of this example are as follows: the computational domain is [ 0 , 700   μ m ] × [ 0 , 1000   μ m ] , the step sizes are Δ r 1 = 3   μ m and Δ r 2 = 3   μ m , respectively, g 1 = 1.48 , k g 1 = 12.21   μ m / s , g 2 = 1.74 , k g 2 = 100.75   μ m / s , k b = 7.49 × 10 8 / μ m 3 / s , b = 2.04 , and ρ c = 2.338 × 10 12 g / μ m 3 .
Figure 7 shows the comparison of the crystal size distribution obtained using different methods at t = 7200 s. There is no exact solution in this example. Therefore, we use the numerical solution of the fifth-order WENO method on very fine grids (2100 × 3000) as the “exact” solution. In the EE–PDDO method, the upwind weight function shown in Equation (15) is used. As can be seen in Figure 7a, there are mainly two distributions in the final crystal size distribution, one in which crystals grow from the seeds and another in which the crystals grow from nuclei. The crystals that grow from the seeds have an average length of about 800 µm and an average width of about 450 µm. They grow by gradually elongating. The size of new nuclei at the origin is small, and as time passes, the nucleated crystals grow larger. By comparing with the “exact” solution in Figure 7a, it is evident that the EE–PDDO method is the best for predicting the distribution related to crystals grown from nuclei, but is the worst for predicting the distribution related to crystals grown from seeds. In all, the EE–PDDO method exhibits higher accuracy than the HR method and the fifth-order WENO method.

3.5. Nucleation and Size-Dependent Growth for 2D Batch Crystallization

We now consider the example of nucleation and size-dependent growth in 2D batch crystallization. The PBE is shown in Equation (39), and the growth rates G 1 and G 2 are size-dependent, which are given in Equation (48) [22].
G 1 ( r 1 , c ( t ) , T ( t ) ) = k g 1 S g 1 0.1 ( 1 + 0.6 r 1 ) G 2 ( r 2 , c ( t ) , T ( t ) ) = k g 2 S g 2 0.1 ( 1 + 0.6 r 2 )
Here, the temperature is expressed as follows [22]:
T ( t ) ( ° C ) = 32 4 ( 1 e t / 310 )
The initial condition is assumed as [4]
f ( r 1 , r 2 , 0 ) = 0.0034786 r 1 2 + r 2 2 + 1.363609 r 1 + r 2 26.5486 18   r 1 , r 2 21 0   elsewhere
The saturated solute concentration c s a t , nucleation rate B 0 , and solute concentration c are the same as in Section 3.4. The kinetic parameters of g 1 , k g 1 , g 2 , k g 2 , k b , b , and ρ c are also the same as in Section 3.4.
Figure 8 shows the comparison of crystal size distribution between the exact solutions and the numerical solutions obtained using different schemes at t = 160   s . The “exact” solution is obtained through the fifth-order WENO method on a fine spatial grid ( Δ r 1 = Δ r 2 = 0.1 ), and the other results are obtained with a spatial grid size of Δ r 1 = Δ r 2 = 1 . As illustrated in Figure 8d, the result of the EE–PDDO method is obtained by using the upwind weight function in Equation (15). Here, we choose the polynomial degree N = 2 and the integer parameter m = 2. From Figure 8, it is obvious that the results obtained using the fifth-order WENO and EE–PDDO methods agree with the exact solution very well, while the results obtained using the HR method show some diffusion. The accuracy of the EE–PDDO method is comparable to that of the fifth-order WENO method. However, the stencils used in EE–PDDO are less than those used in fifth-order WENO.

4. Conclusions

This study presents the EE–PDDO method to solve 2D PBEs in crystallization. Five examples are considered, namely, size-independent growth with smooth and non-smooth distributions, size-dependent growth, nucleation, and size-independent/dependent growth. The results of the EE–PDDO method are compared with the exact solution and other numerical methods. We draw the following conclusions:
(1)
The EE–PDDO method shows a strong ability to solve 2D PBEs. The accuracy of the EE–PDDO method is better than the HR method and is comparable to the fifth-order WENO method. In the case of size-independent growth with smooth and non-smooth distributions, results obtained using the EE–PDDO method perfectly replicates the exact solution. In the nucleation and size-dependent/independent growth cases, results obtained using the EE–PDDO method exhibit non-oscillation and high accuracy, especially in cases of sharp crystal size distribution.
(2)
The Fourier analysis is used to explore the stability of the EE–PDDO method with the simplest form of PBE. In the case of polynomial degree N = 1 and integer parameter m = 1 , the stability condition is 0 < λ 0.5183 for the upwind weight function in Equation (15), and the stability condition is 0 < λ 1 for the upwind weight function in Equation (16) with λ as the Courant number. The optimal time step size needs to be satisfied with λ = 1 for the size-independent growth with smooth and non-smooth distributions.
(3)
The weight functions of the EE–PDDO method are discussed for different examples. In the example of size-independent growth with smooth and non-smooth distributions, the upwind weight function in Equation (16) is recommended; in the example of size-dependent growth, the normal weight function of ω ( ξ 1 , ξ 2 ) = e 4 ( ξ 1 σ 1 ) 2 e 4 ( ξ 2 σ 2 ) 2 is recommended; in the example of nucleation and size-dependent/independent growth cases, the upwind weight function in Equation (15) is recommended.

Author Contributions

Methodology, C.R.; software, C.D.; validation, C.D. and C.R.; writing—original draft preparation, C.D.; writing—review and editing, C.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Key Science and Technology Project of Longmen Laboratory, grant number No. LMYLKT-001.

Data Availability Statement

The original contributions presented in this study are included in the article; further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Puel, F.; Févotte, G.; Klein, J.P. Simulation and analysis of industrial crystallization processes through multidimensional population balance equations. Part 1: A resolution algorithm based on the method of classes. Chem. Eng. Sci. 2003, 58, 3715–3727. [Google Scholar] [CrossRef]
  2. Ma, D.L.; Tafti, D.K.; Braatz, R.D. High-resolution simulation of multidimensional crystal growth. Ind. Eng. Chem. Res. 2002, 41, 6217–6223. [Google Scholar] [CrossRef]
  3. Majumder, A.; Kariwala, V.; Ansumali, S.; Rajendran, A. Lattice Boltzmann method for multi-dimensional population balance models in crystallization. Chem. Eng. Sci. 2012, 70, 121–134. [Google Scholar] [CrossRef]
  4. Qamar, S.; Ashfaq, A.; Warnecke, G.; Angelov, I.; Elsner, M.; Seidel-Morgenstern, A. Adaptive high-resolution schemes for multidimensional population balances in crystallization processes. Comput. Chem. Eng. 2007, 31, 1296–1311. [Google Scholar] [CrossRef]
  5. Briesen, H. Simulation of crystal size and shape by means of a reduced two-dimensional population balance model. Chem. Eng. Sci. 2006, 61, 104–112. [Google Scholar] [CrossRef]
  6. Costa, C.B.B.; Maciel, M.R.W.; Filho, R.M. Considerations on the crystallization modeling: Population balance solution. Comput. Chem. Eng. 2007, 31, 206–218. [Google Scholar] [CrossRef]
  7. Qamar, S.; Angelov, I.; Elsner, M.; Ashfaq, A.; Seidel-Morgenstern, A.; Warnecke, G. Numerical approximations of a population balance model for coupled batch preferential crystallizers. Appl. Numer. Math. 2009, 59, 739–753. [Google Scholar] [CrossRef]
  8. Qamar, S.; Ashfaq, A.; Angelov, I.; Elsner, M.; Warnecke, G.; Seidel-Morgenstern, A. Numerical solutions of population balance models in preferential crystallization. Chem. Eng. Sci. 2008, 63, 1342–1352. [Google Scholar] [CrossRef]
  9. Nicmanis, M.; Hounslow, M.J. Finite-element methods for steady-state population balance equations. AIChE J. 1998, 44, 2258–2272. [Google Scholar] [CrossRef]
  10. Meimaroglou, D.; Kiparissides, C. Monte Carlo simulation for the solution of the bi-variate dynamic population balance equation in batch particulate systems. Chem. Eng. Sci. 2007, 62, 5295–5299. [Google Scholar] [CrossRef]
  11. Zhao, H.; Maisels, A.; Matsoukas, T.; Zheng, C. Analysis of four Monte Carlo methods for the solution of population balances in dispersed systems. Powder Technol. 2007, 173, 38–50. [Google Scholar] [CrossRef]
  12. Hounslow, M.J. A discretized population balance for continuous systems at steady state. AIChE J. 1990, 36, 106–116. [Google Scholar] [CrossRef]
  13. Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization—I. A fixed pivot technique. Chem. Eng. Sci. 1996, 51, 1311–1332. [Google Scholar] [CrossRef]
  14. Nguyen, T.T.; Laurent, F.; Fox, R.O.; Massot, M. Solution of population balance equations in applications with fine particles: Mathematical modeling and numerical schemes. J. Comput. Phys. 2016, 325, 129–156. [Google Scholar] [CrossRef]
  15. Ruan, C. Chebyshev Spectral Collocation Method for Population Balance Equation in Crystallization. Mathematics 2019, 7, 317. [Google Scholar] [CrossRef]
  16. Ma, Z.; Song, H.; Zhang, M. Using Monte Carlo methods to solve the particle number weighing equation and analyze the impact of various factors on bubble distribution. Chem. Eng. 2003, 31, 12–16+1. [Google Scholar]
  17. Qamar, S.; Elsner, M.; Angelov, I.; Warnecke, G.; Seidel-Morgenstern, A. A comparative study of high resolution schemes for solving population balances in crystallization. Comput. Chem. Eng. 2006, 30, 1119–1131. [Google Scholar] [CrossRef]
  18. Ruan, C.; Zhao, X.; Liang, K.; Chang, X. Weighted essentially non-oscillatory method for solving population balances in crystallization processes. In Proceedings of the 2013 International Conference on Advanced Mechatronic Systems, Luoyang, China, 25–27 September 2013. [Google Scholar]
  19. Ruan, C.; Liang, K.; Chang, X.; Zhang, L. Weighted Essentially Nonoscillatory method for two-dimensional population balance equations in crystallization. Math. Probl. Eng. 2013, 2013, 125128. [Google Scholar] [CrossRef]
  20. Ruan, C.; Dong, C.; Liang, K.; Liu, Z.; Bao, X. Euler’s First-Order Explicit Method–Peridynamic Differential Operator for Solving Population Balance Equations of the Crystallization Process. Comput. Model. Eng. Sci. 2024, 138, 3033–3049. [Google Scholar] [CrossRef]
  21. Bekar, A.C.; Madenci, E.; Haghighat, E. On the solution of hyperbolic equations using the peridynamic differential operator. Comput. Methods Appl. Mech. Eng. 2022, 391, 114574. [Google Scholar] [CrossRef]
  22. Gunawan, R.; Fusman, I.; Braatz, R.D. High resolution algorithms for multidimensional population balance equations. AIChE J. 2004, 50, 2738–2749. [Google Scholar] [CrossRef]
  23. Madenci, E.; Barut, A.; Futch, M. Peridynamic differential operator and its applications. Comput. Methods Appl. Mech. Eng. 2016, 304, 408–451. [Google Scholar] [CrossRef]
  24. Madenci, E.; Barut, A.; Dorduncu, M. Peridynamic Differential Operator for Numerical Analysis; Springer International Publishing: Berlin/Heidelberg, Germany, 2019; Volume 10, pp. 15–100. [Google Scholar]
  25. Madenci, E.; Roy, P.; Behera, D. Advances in Peridynamics; Springer Nature: Berlin/Heidelberg, Germany, 2022. [Google Scholar]
  26. Gunzburger, M.; Lehoucq, R.B. A nonlocal vector calculus with application to nonlocal boundary value problems. Multiscale Model. Simulation. 2010, 8, 1581–1598. [Google Scholar] [CrossRef]
  27. Yan, J.; Li, S.; Kan, X.; Zhang, A.M.; Lai, X. Higher-order nonlocal theory of Updated Lagrangian Particle Hydrodynamics (ULPH) and simulations of multiphase flows. Comput. Methods Appl. Mech. Eng. 2020, 368, 113176. [Google Scholar] [CrossRef]
  28. Yu, H.; Li, S. On approximation theory of nonlocal differential operators. Int. J. Numer. Methods Eng. 2021, 122, 6984–7012. [Google Scholar] [CrossRef]
  29. Dorduncu, M. Stress analysis of sandwich plates with functionally graded cores using peridynamic differential operator and refined zigzag theory. Thin-Walled Struct. 2020, 146, 106468. [Google Scholar] [CrossRef]
  30. Gao, Y.; Oterkus, S. Non-local modeling for fluid flow coupled with heat transfer by using peridynamic differential operator. Eng. Anal. Bound. Elem. 2019, 105, 104–121. [Google Scholar] [CrossRef]
  31. Li, Z.; Huang, D.; Xu, Y.; Yan, K. Nonlocal steady-state thermoelastic analysis of functionally graded materials by using peridynamic differential operator. Appl. Math. Model. 2021, 93, 294–313. [Google Scholar] [CrossRef]
  32. Inguva, P.K.; Braatz, R.D. Efficient numerical schemes for multidimensional population balance models. Comput. Chem. Eng. 2023, 170, 108095. [Google Scholar] [CrossRef]
  33. Lenka, M.; Bhoi, S.; Sarkar, D. Two-dimensional population balance modelling and validation of combined cooling and antisolvent crystallization of l-asparagine monohydrate. CrystEngComm. 2023, 25, 1424–1435. [Google Scholar] [CrossRef]
  34. Silling, S.A.; Epton, M.; Weckner, O. Peridynamic states and constitutive modeling. J. Elast. 2007, 88, 151–184. [Google Scholar] [CrossRef]
  35. Chang, H.; Chen, A.; Kareem, A.; Hu, L.; Ma, R. Peridynamic differential operator-based Eulerian particle method for 2D internal flows. Comput. Methods Appl. Mech. Eng. 2022, 392, 114568. [Google Scholar] [CrossRef]
  36. Mark, H. Introduction to Numerical Methods in Differential Equations; Springer New York: New York, NY, USA, 2007. [Google Scholar]
  37. Borchert, C.; Sundmacher, K. Morphology evolution of crystal populations: Modeling and observation analysis. Chem. Eng. Sci. 2012, 70, 87–98. [Google Scholar] [CrossRef]
Figure 1. Interaction of peridynamic points x and x with arbitrary family size and shape.
Figure 1. Interaction of peridynamic points x and x with arbitrary family size and shape.
Crystals 14 00234 g001
Figure 2. Diagram of the upwind weight function.
Figure 2. Diagram of the upwind weight function.
Crystals 14 00234 g002
Figure 3. The comparison of crystal size distribution between the exact solution and the numerical solutions in size-independent growth with smooth distribution. (a) Exact; (b) HR; (c) Fifth-order WENO; (d) EE–PDDO.
Figure 3. The comparison of crystal size distribution between the exact solution and the numerical solutions in size-independent growth with smooth distribution. (a) Exact; (b) HR; (c) Fifth-order WENO; (d) EE–PDDO.
Crystals 14 00234 g003
Figure 4. The comparison of crystal size distribution between the exact solution and the numerical solutions in size-independent growth with non-smooth distribution. (a) Exact; (b) HR; (c) Fifth-order WENO; (d) EE–PDDO.
Figure 4. The comparison of crystal size distribution between the exact solution and the numerical solutions in size-independent growth with non-smooth distribution. (a) Exact; (b) HR; (c) Fifth-order WENO; (d) EE–PDDO.
Crystals 14 00234 g004
Figure 5. The comparison of crystal size distribution between the exact solution and the numerical solutions in size-dependent growth. (a) Exact; (b) HR; (c) Fifth-order WENO; (d) EE–PDDO.
Figure 5. The comparison of crystal size distribution between the exact solution and the numerical solutions in size-dependent growth. (a) Exact; (b) HR; (c) Fifth-order WENO; (d) EE–PDDO.
Crystals 14 00234 g005
Figure 6. Diagram of the KDP crystals.
Figure 6. Diagram of the KDP crystals.
Crystals 14 00234 g006
Figure 7. The comparison of crystal size distribution between the exact solution and the numerical solutions in nucleation and size-independent growth. (a) Exact; (b) HR; (c) Fifth-order WENO; (d) EE–PDDO.
Figure 7. The comparison of crystal size distribution between the exact solution and the numerical solutions in nucleation and size-independent growth. (a) Exact; (b) HR; (c) Fifth-order WENO; (d) EE–PDDO.
Crystals 14 00234 g007
Figure 8. The comparison of crystal size distribution between the exact solution and the numerical solutions in nucleation and size-dependent growth. (a) Exact; (b) HR; (c) Fifth-order WENO; (d) EE–PDDO.
Figure 8. The comparison of crystal size distribution between the exact solution and the numerical solutions in nucleation and size-dependent growth. (a) Exact; (b) HR; (c) Fifth-order WENO; (d) EE–PDDO.
Crystals 14 00234 g008
Table 1. L 1 errors and L 2 errors for predicting crystal size distribution with different schemes in size-independent growth with smooth distribution.
Table 1. L 1 errors and L 2 errors for predicting crystal size distribution with different schemes in size-independent growth with smooth distribution.
MethodHRFifth-Order WENOEE–PDDO, N = 1, m = 1, Upwind Weight (Equation (16))
L 1 -error3.08458 × 10−34.05432 × 10−53.40151 × 10−4
L 2 -error1.07827 × 10−21.31731 × 10−46.74667 × 10−4
Table 2. L 1 errors and L 2 errors for predicting crystal size distribution with different schemes in size-independent growth with non-smooth distribution.
Table 2. L 1 errors and L 2 errors for predicting crystal size distribution with different schemes in size-independent growth with non-smooth distribution.
MethodHRFifth-Order WENOEE–PDDO, N = 1, m = 1, Upwind Weight (Equation (16))
L 1 -error1.15517 × 10−26.35183 × 10−34.10000 × 10−3
L 2 -error6.26252 × 10−24.23699 × 10−26.40312 × 10−2
Table 3. Errors and CPU time of EE–PDDO scheme in different Δ t .
Table 3. Errors and CPU time of EE–PDDO scheme in different Δ t .
Δ t ( λ ) 0.001   ( λ = 0.01 ) 0.005   ( λ = 0.05 ) 0.01   ( λ = 0.1 ) 0.05   ( λ = 0.5 ) 0.1   ( λ = 1 )
L 1 -error3.69075 × 10−23.64740 × 10−23.59031 × 10−22.96320 × 10−24.10000 × 10−3
L 2 -error1.17296 × 10−11.16245 × 10−11.14873 × 10−11.00555 × 10−16.40312 × 10−2
CPU time(s)125.30350.53747.13135.81835.508
Table 4. L 1 errors and L 2 errors for predicting crystal size distribution with different schemes in size-dependent growth.
Table 4. L 1 errors and L 2 errors for predicting crystal size distribution with different schemes in size-dependent growth.
MethodHRFifth-Order WENOEE–PDDO, N = 1, m = 1
Upwind Weight (Equation (15))
EE–PDDO, N = 1, m = 1
Normal Weight
L 1 -error5.50266 × 10−41.35244 × 10−44.29756 × 10−41.00368 × 10−5
L 2 -error1.95435 × 10−29.82604 × 10−37.59137 × 10−31.73679 × 10−4
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Dong, C.; Ruan, C. Euler’s First-Order Explicit Method–Peridynamic Differential Operator for Solving Two-Dimensional Population Balance Equations in Crystallization. Crystals 2024, 14, 234. https://doi.org/10.3390/cryst14030234

AMA Style

Dong C, Ruan C. Euler’s First-Order Explicit Method–Peridynamic Differential Operator for Solving Two-Dimensional Population Balance Equations in Crystallization. Crystals. 2024; 14(3):234. https://doi.org/10.3390/cryst14030234

Chicago/Turabian Style

Dong, Cengceng, and Chunlei Ruan. 2024. "Euler’s First-Order Explicit Method–Peridynamic Differential Operator for Solving Two-Dimensional Population Balance Equations in Crystallization" Crystals 14, no. 3: 234. https://doi.org/10.3390/cryst14030234

APA Style

Dong, C., & Ruan, C. (2024). Euler’s First-Order Explicit Method–Peridynamic Differential Operator for Solving Two-Dimensional Population Balance Equations in Crystallization. Crystals, 14(3), 234. https://doi.org/10.3390/cryst14030234

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop