1. Introduction
The population balance equation (PBE), also known as the general dynamic equation, plays an important role in describing the evolution of particle size distributions for various atmospheric and engineering applications, including the airborne transmission of COVID-19 [
1], atmospheric aerosol formation and transport [
2,
3,
4], combustion-generated nanoparticles [
5,
6], polymerization and crystallization processes [
7,
8], and nuclear radioactive aerosols [
9]. The basic form of the population balance equation originates from the Smoluchowski equation, which describes the temporal evolution of a continuous particle size distribution (PSD) caused by coagulation [
10]. Because of complex integral–differential characteristics, various numerical methods have been proposed to solve the PBE. The numerical methods can be categorized into the sectional method [
11,
12], method of moments [
13,
14,
15], stochastic method [
16,
17,
18], and spectral method [
19,
20]. These methods are based on different principles and exhibit different accuracies and efficiencies. Detailed comparisons are provided in previous studies [
21,
22].
Among the existing numerical methods, many assume that the PSD can be described by a specific distribution function, and the parameters of this predefined distribution can be solved using moment equations. The typical distribution patterns include the well-known log-normal distribution [
23], the improved log-skew normal distribution [
24], the modal aerosol approach [
25], and weighted sum of number density functions with a log-normal kernel, gamma kernel, or beta kernel [
26,
27]. Although these methods achieve moment closure and well predict moments, they generally lack accuracy and stability in predicting the size distribution evolution. This is mainly because the PSD information is missing when transforming the PBE into a finite set of moment equations [
28]. Moreover, the predefined distributions are not sufficiently flexible to represent the unsolved particle size distribution.
In recent years, radial basis functions (RBFs) have been extensively applied in solving complex numerical problems involving ordinary differential equations (ODEs), partial differential equations (PDEs), and integral equations [
29,
30,
31]. Moreover, RBFs are popular kernel functions used in neural networks and have many applications in predicting atmospheric pollutant concentrations [
32,
33,
34,
35]. The basic principle of RBFs is their superior ability to approximate continuous functions. Therefore, in a previous study, RBFs were adopted to solve the PBE for particle coagulation [
36]. The number density function was approximated using two-dimensional Gaussian radial basis functions with volume and time as the variables, which provided an accurate description of PSDs over a period of time but were limited to predicting long time-period, self-preserving distributions. Subsequently, Alzyod and Charton [
37] proposed a meshless radial basis method to solve the PBE for particle growth, nucleation, coagulation, and breakage. They used multiquadric basis functions and transformed the PBE into ODEs, and, in another study, proposed an adaptive radial basis method to model the hydrodynamic behavior of liquid–liquid dispersions [
38]. However, these method formulations are difficult to use and validations of the method are only for simple coagulation cases, which cannot guarantee the accuracy of the method for more general coagulation problems. Therefore, the RBF methods proposed in previous studies are incomplete and need further improvements and more sufficient validations.
This study further develops a reliable and robust RBF method to solve the PBE for coagulation and provides an adequate validation of the method under various possible coagulation situations.
Section 2 presents the theory and numerical details of the method. The governing equations were derived based on the radial basis function approximation. The integral terms were obtained using the Gaussian quadrature, along with mathematical treatments that decreased the quadrature error. The final derived equations are a set of ODEs that can be easily solved using commercial ODE solvers.
Section 3 validates the RBF method for various coagulation cases, including constant kernel and sum kernel coagulations, realistic Brownian coagulation in the continuum and free-molecular regimes, self-preserving distributions, and multimodal coagulation. The predictions obtained using the RBF method were compared with the analytical solutions and those obtained using a sectional method to validate its accuracy and efficiency. Finally,
Section 4 summarizes the content of this study.
2. Materials and Methods
A radial basis function (RBF) is a function whose value changes with the distance from a center point. Typical radial basis functions are the Gaussian, multiquadric, and inverse multiquadric functions. In this section, the theory and numerical procedures of the radial basis function method are presented in detail.
2.1. Governing Equations
For particle coagulation, the governing PBE to describe the time evolution of particle size distributions is expressed in the following integral–differential form:
where
is the number density function for the particle volume
at time
,
is the birth term of particles with volume
caused by coagulation,
is the death term of particles with volume
caused by coagulation, and
is the coagulation kernel for two particles with volumes
and
, respectively.
The principle of the RBF method is to approximate the number density function
using the weighted sum of the radial basis functions. In our previous study [
36], both volume
and time
were regarded as the input variables into the two-dimensional Gaussian functions, thereby complicating the solving process. A convenient approach is to use only a one-dimensional Gaussian function for volume
and to make the coefficients time-variant, which can be expressed as follows:
where
is the
i-th coefficient that changes with time
,
is the
i-th Gaussian kernel located at the center point
, and
is the standard deviation of the
i-th Gaussian kernel. It should be noted that the logarithmic transformation from
to
is used in Equations (2) and (3) because aerosol particle sizes are usually distributed on logarithmic scales rather than on linear scales.
In addition to the normal RBF approximations, Equations (2) and (3) also have physical implications. The Gaussian function coupled with the logarithmic transformation is equivalent to the widely used log-normal distribution. Thus, the approximation of
in Equations (2) and (3) can be interpreted as the superposition of many log-normal sub-distributions, which is also consistent with the idea of a previous study [
39]. Each sub-distribution represents a flexible component; therefore, all sub-distributions have the ability to approximate distributions of any shape.
Substituting Equation (2) into Equation (1) and calculating the derivative yields:
According to the collocation method principle, for any given point
, Equation (4) is automatically satisfied. For simplicity and realizability, the collocation points are set to be the same as the center points of the Gaussian kernels. In this case, Equation (4) is transformed into the following
p equations:
where
.
Using the matrix representation for Equation (5) leads to the following governing equations for the coefficients
:
In Equation (6), the left array represents the time derivatives of
, the first term on the right-hand side is a constant matrix that depends only on the kernel function and selected center points, and the second term is an array of coagulation integrals. The calculations for these integrals are presented in
Section 2.2.
The initial conditions for Equation (6) is simply substituting
into Equation (2) and using the initial size distribution, as follows:
By applying the same procedure as that for Equation (6), the initial values of
can be calculated as follows:
Thus, Equation (6), together with the initial conditions in Equation (8), constitutes the governing equations for , which is relatively concise compared with the original PBE.
2.2. Quadrature Approximation
For the governing Equation (6), the remaining issue is to calculate the complex integrals for the birth term
and death term
. In our previous study [
36],
and
were calculated using the Gauss–Legendre quadrature and Gauss–Laguerre quadrature, respectively; however, the two quadratures may not accurately predict long-term coagulation, because the integral error increases with time in the long-term. To address this problem, this study adopts special mathematical treatments for the integrals, and the quadrature accuracy can be guaranteed even for long time-period, self-preserving predictions, as illustrated in
Section 3.3.
For the birth term
, the integral domain of
should be transformed into [–1, 1] to satisfy the Gauss–Legendre quadrature. To achieve this, the transformation
is used, and the original integral becomes:
Applying the Gauss–Legendre quadrature to Equation (9) yields the following final quadrature formula for
:
where
and
are the nodes and weights of the Gauss-Legendre quadrature, respectively, and
is the number of Gauss-Legendre quadrature points. Calculations of the nodes and weights for the quadrature are provided in the literature [
40].
For the death term
, the previous Gauss–Laguerre quadrature is not suitable for long-term coagulation calculations. Thus, the infinite integral interval [0, ∞] is first truncated into a finite interval
, which has a negligible impact on the integral results. Here,
is the smallest volume, and
is the largest volume of the computation domain. Because the entire particle volume range is relatively large in actual situations, the integral should be calculated under logarithmic scales to decrease the quadrature error; therefore:
To satisfy the Gauss–Legendre quadrature rule, the above integral interval
is transformed into [–1, 1] using the relation
, where
and
. Equation (11) then becomes:
In Equation (12), the values of can be calculated by substituting the obtained coefficients of into Equation (2).
According to Equations (10) and (12), the coagulation integrals and can be calculated conveniently using the Gauss–Legendre quadrature. The quadrature calculation only involves multiplications and additions that can be further integrated into the matrix computation, which can considerably decrease the numerical costs.
2.3. Numerical Procedures
Based on the RBF approximation in
Section 2.1 and the quadrature calculation in
Section 2.2, the original PBE is simplified into the ODEs of the coefficients
, which can be easily calculated using commercial ODE solvers. Numerically, there are two steps involved in solving the temporal evolution of an initial PSD for a given coagulation kernel.
The first step is to define numerical parameters for the calculations. The range of the particle volume
is denoted as
, and the range of the physical time is denoted as
. The selection of the center points in Equation (2) is the same as that used in our previous study [
36]. Thus, the center points are set to be equally distributed in logarithmic scales and can be expressed as follows:
where
is the total number of center points.
The collocation points are the same as the center points defined in Equation (13). The standard deviation
of each Gaussian function is taken as twice the distance between the neighboring center points. Detailed reasons for this setting are provided in a previous study [
41]. Since the center points are equally distributed, all Gaussian functions share the same standard deviation, which is calculated as follows:
The second step is to use an ODE solver to calculate the temporal evolution of the coefficients based on Equations (6) and (8). Before running the numerical program, the initial condition should be determined. The initial coefficient can be computed from the initial PSD based on Equation (8). This process is concise and robust for any continuous PSD.
The quadrature approximations for any given coagulation kernel were calculated using Equations (10) and (12), respectively. The number of quadrature points was determined by the desired accuracy for the numerical problem. All the program codes were implemented in MATLAB R2018b on a personal computer with an Intel Core i7-6700 CPU and 16 G RAM. The calculations were mainly performed using a matrix computation to achieve a fast computing efficiency. The ODE45 solver was used, which implements a Runge–Kutta method with a variable time step for efficient computation. Thus, the number of the time steps changes according to the actual condition of the numerical problem to be solved. All the codes and models will be open-sourced at a later stage.
3. Results
To validate the RBF method, a thorough analysis including seven different numerical cases is presented in this section. The detailed expressions of the coagulation kernel functions and the initial size distributions of these cases are presented in
Table 1.
Cases 1 and 2 were constant kernel and sum kernel coagulations, respectively. These two cases were chosen because they can be analytically solved according to a study by Ramabhadran et al. [
42]. The initial distribution was a normalized exponential distribution with
and
. For Cases 1 and 2, the RBF method was compared with the analytical solutions to validate its accuracy.
Cases 3 and 4 were Brownian coagulations under continuum and free-molecular regimes, respectively. The two coagulation kernels were selected because they are widely used to describe submicron particle coagulation dynamics. The initial size distribution was a normalized log-normal distribution with
,
, and
. As there are no accurate analytical solutions for Brownian coagulation, a highly reliable sectional method was used as a benchmark model [
12]. The sectional code, using a sectional spacing factor of 1.05, was already verified in our previous study [
24]. Thus, the RBF method was compared with the sectional method for Cases 3 and 4.
Cases 5 and 6 were used to verify the self-preserving predictions of the RBF method for Brownian coagulation under the continuum and free-molecular regimes. Because self-preserving distributions indicate a long-term asymptotic behavior, the two cases were selected to verify the accuracy of the RBF method for long-term coagulation simulation. The initial distribution was the same as those in Cases 3 and 4. Reliable predictions of the self-preserving distribution were obtained using the sectional method with a sectional spacing factor of 1.05.
Case 7 was a bimodal particle coagulation, which was more complex than the above unimodal cases. Bimodal size distributions are common in actual situations and were selected for complete validation. The initial size distribution was a combination of two normalized log-normal sub-distributions, with the parameters , , and , while , , and . In this case, the RBF method was again compared with the sectional method to validate its accuracy.
The numerical settings for these cases are listed in
Table 2. The particle volume domain depends on the minimum volume
and maximum volume
. The minimum volume
is determined from the truncated smallest volume of the initial size distribution. The maximum volume
is determined to ensure that the particle volume does not exceed the upper limit at a given coagulation time. Setting a very large value for the
is feasible, but this increases the computational costs and is not recommended. The physical time
is determined by the demand for the numerical problem. Five center points are generally sufficient to cover a unit logarithmic interval; therefore, the number of center points was set to
. The larger the number of center points, the higher the accuracy and computational costs. The number of quadrature points varies under different situations because the quadrature accuracy depends on the coagulation kernel and numerical case. In the following sections, the validation results of the RBF method for seven numerical cases are presented.
3.1. Constant Kernel and Sum Kernel Coagulation
Figure 1 shows a comparison of the RBF method and the analytical solution for predicting the temporal evolution of the PSD for the constant kernel coagulation (Case 1) and sum kernel coagulation (Case 2). As coagulation proceeds, the PSD becomes shorter and shifts to the right side. The reason for this is that the particle number concentration decreases and the average particle size increases as particles collide and coagulate to form bigger particles. At different coagulation times, the predictions of the RBF method were the same as the results of the analytical solutions, meaning that the RBF method could almost accurately predict the temporal evolution of the PSDs for the two cases. In addition, the present RBF method could provide an explicit formula of the number density function from Equation (2). Therefore, the RBF method has similar accuracy and features as analytical methods.
To further quantify the error between the predictions of the RBF method and the results of the analytical solutions, the root mean square error (RMSE) was used in this study to describe the deviation, which is expressed as follows:
According to Equation (15), the RMSE value represents the average deviation of the RBF method predictions from the analytical solutions within a given range.
Based on Equation (15),
Figure 2 shows the temporal evolution of the RMSE values for constant kernel coagulation (Case 1) and sum kernel coagulation (Case 2). For Case 1, the maximum RMSE value over the entire coagulation time was 5.61 × 10
−5. This is a considerably small value, which again proves the accuracy of the RBF method compared with the analytical solutions. Moreover, the RMSE values in
Figure 2a decreased with time, indicating that the RBF method maintained a high accuracy during the entire coagulation process. For Case 2, the maximum RMSE value was 3.42 × 10
−5 (
Figure 2b), which is also considerably small compared to the absolute values in
Figure 1b. The RMSE curve first increased with the coagulation time and then gradually plateaued. This indicates that the RBF method also maintained a high accuracy during the sum kernel coagulation. From the above comparisons, the RBF method can provide almost the same predictions as can the analytical solutions for constant kernel and sum kernel coagulations.
3.2. Brownian Coagulation under the Continuum and Free-Molecular Regimes
Figure 3 shows a comparison of the RBF and the sectional methods for predicting the temporal evolution of the PSD for Brownian coagulation under the continuum (Case 3) and free-molecular (Case 4) regimes. The evolution patterns of the PSDs in
Figure 3 are similar to those in
Figure 1 as the result of particle coagulation. At different coagulation times, high consistency was observed between the predictions of the RBF method and those of the sectional method under both cases. Thus, the RBF method can accurately predict the temporal evolution of the PSD for actual Brownian coagulation problems.
In addition to predicting the size distribution evolution, the moments of the particles are of great interest in many situations. For the present RBF method, the
k-th moment
can be directly obtained by integrating Equations (2) and (3) from 0 to ∞, as follows:
Thus, the RBF method can also be used to predict moments of arbitrary order based on the obtained coefficients .
Based on Equation (16),
Figure 4 further compares the temporal evolution of the first two moments,
and
, predicted using the RBF and the sectional methods for Cases 3 and 4, respectively.
is the total particle number concentration and
is the total particle volume. Because particle coagulation means that two particles collide to form one larger particle, the particle concentration should decrease during coagulation, while the particle volume remains constant. This is highly consistent with the results in
Figure 4, which show that the
curves decrease with time and
curves remain horizontal. The predictions of
and
using the RBF method are highly consistent with the results of the sectional method. This indicates that the RBF method can provide accurate predictions of the temporal evolution of the first two moments. Thus, the RBF method can also be used as a benchmark for predicting the moment evolution for Brownian coagulation in the continuum and free-molecular regimes.
In addition to the continuum and free-molecular regimes, atmospheric PSDs also cover the intermediate transition regime. As for the transition regime, the mathematical expression of the coagulation kernel is more complicated [
23]. Since the RBF method directly uses numerical quadrature for the coagulation integrals, the solution process does not depend on the specific form of the coagulation kernel, as shown in Equations (10) and (12). Thus, the present RBF method can be easily extended to solve the transition regime and even more complicated coagulation regimes.
3.3. Self-Preserving Distributions
Particles undergoing Brownian coagulation finally reach a self-preserving state, indicating a long-term asymptotic behavior. Self-preserving distributions are important characteristics that need to be correctly predicted. In our previous study [
36], the two-dimensional RBF method lacked the ability to predict the self-preserving distribution, owing to the RBF approximation for both the volume
and time
. For the present RBF method, only the volume
is pre-discretized, and the time
is driven by the time stepping of ODEs, allowing the RBF method to predict self-preserving distributions.
A self-preserving distribution is normally expressed using the dimensionless particle volume
and dimensionless number density function
, as follows:
where
is the mean particle volume calculated as
(
is the total particle volume, and
is the total particle number concentration). The values of the first two moments can be calculated using Equation (16).
Figure 5 compares the self-preserving distributions predicted using the RBF and the sectional methods for Brownian coagulation under the continuum (Case 5) and free-molecular (Case 6) regimes. The self-preserving distribution curves were obtained by applying the transformation in Equations (17) and (18) for the number density function over the long-term coagulation. The entire physical times were set to 1000 and 100 for Cases 5 and 6, respectively, which is regarded as approaching the long-term limit (
Table 2). For both the continuum and free-molecular regimes, the RBF method provided almost the same predictions for the self-preserving distribution compared with the results obtained using the sectional method. Thus, the present RBF method is verified to be reliable in predicting long-term coagulation behavior and is not limited by the physical coagulation time.
3.4. Multimodal Coagulation
Previous studies have focused only on unimodal particle coagulation problems. In many actual situations, particle sizes are described by bimodal or multimodal distributions [
25]. Thus, the RBF method must be validated for multimodal particle coagulation. Here, a bimodal coagulation case (Case 7) was selected as a representative, and a detailed description of this case is provided at the beginning of
Section 3.
Figure 6 shows the temporal evolution of an initial bimodal distribution predicted using the RBF and sectional methods for Brownian coagulation under the continuum regime. The peak of the smaller particles rapidly decreased, and the peak of the larger particles gradually shifted to the right. The evolution results indicate that smaller particles were scavenged by larger particles due to intra-particle coagulation. At different coagulation times (
= 0, 0.3, and 1.0), the RBF method provided almost the same predictions of PSDs as did the sectional method, and differences between the two methods were difficult to identify. Thus, the present RBF method can accurately predict bimodal coagulation evolution.
For cases of more than two modes, the RBF method can be smoothly extended be-cause RBFs have been proved to have the best approximation ability to any continuous functions [
43]; therefore, the RBF method should not be limited by the number of modes of an initial size distribution.
3.5. Computational Time
In addition to the accuracy of the RBF method, computational speeds are very important when evaluating numerical methods.
Table 3 lists the computational times of the RBF method and sectional method for the different numerical cases. For Cases 1 and 2, the computational times of the sectional method are not listed because the two cases were analytically solved as discussed in
Section 3.1. From this table, the computational times of the RBF method were less than one second or just a few seconds, which were at least 100 times smaller than those of the sectional method. Thus, the RBF method was highly efficient for the tested cases compared to the sectional method. Based on an in-depth analysis, the computational costs of the RBF method resulted mainly from two factors. One was the numerical solution of the ODEs, and the other was the quadrature calculation. This indicates that the computational time increases with an increase in the number of center points and quadrature points. In addition, the range of the physical time significantly and greatly affects the computational time.
4. Conclusions
This study presents a reliable and robust radial basis function (RBF) method for predicting the evolution of particle size distributions during coagulation. The basic idea of the RBF method is to approximate the number density function using the weighted sum of the radial basis functions with the coefficients . The Gaussian function serves as the basis function and is distributed evenly in a logarithmic volume space. The governing equations were derived based on the RBF approximation, after which the original PBE was transformed into the simple ODEs of . The integral terms were obtained using the Gauss–Legendre quadrature and could be efficiently calculated using a matrix. Thus, the RBF method solves the number density function by using ODEs and numerical quadratures.
Subsequently, the accuracy and efficiency of the RBF method was validated under various numerical cases. These cases covered many key issues in this field and could either be analytically solved (Cases 1 and 2) or numerically solved using a highly reliable sectional method (Cases 3–7). The comparison results show that the RBF method can provide nearly the same predictions as can analytical solutions for constant kernel and sum kernel coagulations. For realistic Brownian coagulation, the RBF method can accurately predict the temporal evolution of PSDs as well as the moments of particles. Moreover, the RBF method provides almost accurate predictions of the self-preserving distributions for Brownian coagulation; therefore, it is verified to be reliable for long-term coagulation simulations. In addition to unimodal cases, the RBF method can also solve bimodal coagulation problems, and the predictions are highly consistent with the results of the sectional method. The computational time of the RBF method is typically of the order of seconds for all cases, meaning that the RBF method is numerically efficient.
Based on these validations, the RBF method has a high accuracy and has a high computational efficiency. It directly provides a continuous number density function through an explicit formula, which has features similar to those of analytical solutions. Thus, the proposed method can be used as a reference method for aerosol coagulation simulations. Future work should focus on extending and validating the RBF method for more complex particle dynamics, including nucleation, coagulation, and surface growth.