Nanofluidics and NanofluidsView this Special Issue
The Asymptotic Behavior of Particle Size Distribution Undergoing Brownian Coagulation Based on the Spline-Based Method and TEMOM Model
In this paper, the particle size distribution is reconstructed using finite moments based on a converted spline-based method, in which the number of linear system of equations to be solved reduced from 4m × 4m to (m + 3) × (m + 3) for (m + 1) nodes by using cubic spline compared to the original method. The results are verified by comparing with the reference firstly. Then coupling with the Taylor-series expansion moment method, the evolution of particle size distribution undergoing Brownian coagulation and its asymptotic behavior are investigated.
Particle size distribution (PSD) is one of the most important properties of aerosol particles, including transport, sedimentation, and so on . It is also of utmost interest in many industrial applications, such as powder preparation and particle synthesis [2, 3]. The evolution of the PSD undergoing different dynamic processes is usually described in the framework of population balance equations (PBEs) mathematically , which have a strong nonlinear structure in most cases and cannot be solved analytically. With a high-computational efficiency, the moment-based method has become a powerful tool for investigating aerosol microphysical processes, in which cases some statistical characteristics of the PSD, namely, the moments of the PSD, are obtained . However, the detailed information about the target PSD is out of reach. Theoretically, the PSD is equal to the moments of all orders. The proof about the uniqueness of reconstruction in the case all moments are known is given with an appropriate condition that the range of the PSD is in a finite interval . But in practice, only a finite number of moments are obtained.
Using a given number of moments to reconstruct the PSD is known as the finite-moment problem or inverse problem in mathematical analysis . Generally, this problem is distinguished between the three types for the monovariate case: the Hausdorff moment problem with the PSD supported on the closed interval [a, b], where [a, b] are the lower and upper limits of the domain of PSD; the Stieltjes moment problem with the PSD supported on [0, +∞); and the Hamburger moment problem with the PSD supported on (−∞, +∞) . Until now, there exist several frequently used reconstruction methods in the literature mainly for the Hausdorff moment problem, including but not limited to parameter-fitting method, Kernel density function-based method, maximum entropy method, and spline-based method. The parameter-fitting method is to assume the PSD as a simple function (i.e., log-normal distribution or gamma distribution), where the parameters in the function are determined by the given low-order moments . It is the fastest and easiest method but with drawbacks that need a priori knowledge about the solution and limited to simple shapes, even though a weighted sum of different simple functions can be used . The kernel density function-based method is a positivity-preserving representation and can be regarded as the development of parameter-fitting method, which approximates the PSD by a superposition of weighted kernel density functions . This method gives rise to an ill-posed problem for determining weights, and a large number of available moments are needed to ensure accuracy. Based on the maximization of the Shannon entropy or the minimization of the relative entropy from information theory, the maximum entropy method is a notable method which needs relatively less knowledge of the prior distribution or the number of moments compared to the previous two methods [8, 11]. With the advantage of no priori assumptions on the shape of the PSD as well as that the needed number of moments only depends on that of interpolation nodes, the spline-based method proposed by John et al.  has attracted some researchers’ attention, such as the investigation on particle aggregation and droplet coalescence [12, 13]. And an adaptive spline-based algorithm with a wider application for nonsmooth and multimodal distributions was developed later . More relevant research about the comparisons between these methods can be found in the literature [14, 15].
In this paper, we will use the spline-based method to reconstruct the PSD coupling with PBEs describing Brownian coagulation in the free molecule regime and continuum regime. Compared to the original method, the number of linear system of equations to be solved is significantly reduced through substituting the continuous conditions. The correctness of this new treatment is verified by comparing with the reference results in . Then with the moments obtained by the Taylor-series expansion moment method (TEMOM) , the evolution of PSD due to Brownian coagulation and its asymptotic behavior are investigated.
2. Theory and Modeling
2.1. Modeling of Spline-Based Method
In the original method, the support of the target PSD [a, b] is divided into m subintervals: a = x1 < x2 < ⋯ < xm+1 = b. In each subinterval, the PSD is approximated by a spline (piecewise polynomial) si(l)(x) of degree l; thus, there exist (l + 1)m unknowns. For cubic spline (l = 3), the splines s(l)(x) and their first and second derivatives are continuous at each node xi (i = 2, 3,…, m), which give 3(m − 1) conditions. With the smooth boundary conditions, which means s(l)(x) and their first and second derivatives are null at nodes x1 and xm+1, there still require m − 3 additional conditions, which have to be supplemented by the known moments. Then, a 4m × 4m ill-conditioned linear system is obtained. In order to improve the accuracy of calculation, the interval should be set as small as possible, which is controlled by tolred. For example, the last (or the first) subinterval is divided into n smaller subintervals: xm = xm1 < xm2 < ⋯ < xmn = xm+1; if the ratio of 2-norm of s(l)(xmn) to the maximum of s(l)(x) is less than tolred, the last node is reset as xm+1 = (xm + xm+1)/2. Furthermore, tolneg and tolsing are introduced to guarantee that the value of s(l)(x) is nonnegative. More detailed procedure is shown in .
In this paper, we use a converted ansatz for s(l)(x) to reduce the number of the linear system. For cubic spline, the second derivative in each node is set as Li, then can be written in the following form using linear interpolation:
Then, s(x) and their first derivatives can be gotten through integrating:where and are integral constants. With the continuity of the spline and their first derivatives at xi (i = 2, 3,…, m), we can getin which Δxi is the length of the ith subinterval and and are related to the left boundary conditions. Thus, the sum of the number of moments and boundary conditions needed to solve the equations is m + 3.
Usually, we consider that the value of PSD out of the support [a, b] is small enough and can be set as zero:and the first derivatives are denoted aswhere and are zero for smooth boundary conditions. Then, (3) can be simplified by substituting the left boundary conditions:
Together with the right boundary conditions, we can get the following formula:in which and are given as follows (i = 2, 3,…, m):
The kth order moment Mk of the PSD is defined as follows:
Thus, the kth order moment of s(x) isin which Ii are
In order to represent Li explicitly, (10) is arranged as follows:where
Now together with (7), a (m + 3) × (m + 3) linear system for Li, q1, and q2 is obtained. Next, we will discuss the number of moments that should be supplemented (note that, in this paper, all cases calculated with q1 and q2 are zero):(1)If q1 and q2 are unknowns, (m + 1) moments are needed to solve these equations.(2)If q1 and q2 are zero or any other constants given, (m − 1) moments are needed; if the value of at boundary is given (such as smooth conditions in , namely, L1 = Lm+1 = 0), (m − 3) moments are needed. And in this case, the order of the coefficient matrix is (m + 1) × (m + 1) or (m − 1) × (m − 1).(3)If q1 and q2 obey some relationships, for example, q1 = (s1(x2) − s1(x1))/Δx1, q2 = (sm(xm+1) − sm(xm))/Δxm, then L1 = −L2/2 and Lm+1 = −Lm/2 can be derived and (m − 1) moments are needed.
For quadratic spline, we can also get a (m + 1) × (m + 1) linear system in the same way by denoting thatwhere li is the first derivative in each node. The corresponding s(x) and linear system are as follows:where
2.2. Modeling of PBE and TEMOM
The population balance equation describing irreversible Brownian coagulation with continuous monovariable can be written as follows :where is the number density function of the particles with volume from υ to υ + dυ at time t and is the collision frequency function between particles with volume υ and υ1. In the free molecule and continuum regime, are represented separately asin which B1 = (3/4π)1/6(6kbT/ρp)1/2 and B2 = 2kbT/3μ, where kb is the Boltzmann constant, T is the temperature, ρp is the particle density, and μ is the gas viscosity.
With the definition of the kth order moment, Mk, (17) is transformed to a series of original differential equations by multiplying both sides with vk and then integrating over all particle sizes:
Using the Taylor-series expansion technology to approximate the collision frequency function and fractional moments, the moment equations are closed without any other artificial assumption [16, 18]. In the original TEMOM model, the first three moments can be obtained easily using the fourth-order Runge-Kutta method with M1 remaining constant due to the mass conservation requirement. The corresponding higher and fractional moments are as follows :where MC = M0M2/M12 is a dimensionless moment. Obviously, the reconstruction depends heavily on the reliability of known moments. Based on the log-normal size distribution assumption, the maximum relative error for Mk of this model is discussed by Xie , and the results demonstrate that the error of Mk for k ≤ 2 with a small standard deviation is acceptable. Furthermore, theoretical analysis of the PBE is feasible because of the relative simple form of this model [20, 21], and the explicit asymptotic solutions are as follows:and MC tends to a constant 2.200126847 or 2, respectively. Using the similarity transformation η = v/(M1/M0), the PSD can be arranged as follows:
According to the theory of self-preserving, ψ(η) does not change with time at a large t , and its moments only depend on k and MC:
3. Results and Discussions
One difficulty of the inverse problem is the ill-conditioned coefficient matrix of the linear system. Another is that the value of s(x) is nonnegative. By using the pseudoinverse routine, a least-squares solution of the linear system is obtained, in which the singular values smaller than tolsing are set as zero (see Remark 4.2 in John et al. ). Moreover, the parameter α is introduced to avoid large difference in the order of magnitude. In this paper, we will follow this treatment. Figure 1 shows the results of the reconstruction about Example 2.1 in  by using quadrature spline and cubic spline proposed in this paper. And the parameters tolred, tolneg, and tolsing are set as the same of those in the literature to maintain consistency. It should be noted that the tolerance values have an influence on the results [7, 14]. The great agreements with the references verify the validity of this new converted method. However, an underlying flaw is that only the continuity of s(x) is necessary in practice. Moreover, the sensitivity of tolsing to solution may increase when the number of the linear system sharply decreases. It can also be seen that some inflexion points appear with m increasing. This may be caused by the increasing condition number of the linear system.
The reconstruction of in the free molecule regime and continuum regime for different moments (m) by cubic spline is shown in Figure 2, where the references are from Lai et al.  and Friedlander and Wang . The initial interval is set as [1e − 5, 10], and the spacing of adjacent nodes are equidistant logarithmically. Both the left and right boundaries are adjusted according to the comparison results of ||s(ηmn)||2/max(s(η)) and tolred, with tolred = 1e − 2 for m = 6 and 1e − 4 for m = 7 or 8, respectively. Besides, tolred = −1e − 2 and the initial value of tolsing are set as 1e − 36. In addition to the first three integral moments m0, m1, and m2, the fractional moments m1/3, m2/3, m4/3, and m5/3 are chosen to reconstruct , for the reason that these fractional moments with volume-based variable are proportional to the integral moments with length-based variable. It can be seen that the results for m = 7 show relatively small differences compared to the references. Generally, the differences may be caused by two parts: one is the error of TEMOM model and the other is the error of spline-based method. Scaling the moments Mk and time t by Mk = MkM00(k−1)/M10k, τ = tB1M005/6M101/6, or τ = tB2M00, the evolution of dimensionless n(υ, t) for m = 7 at long time is presented in Figure 3 with the initial conditions given as = 1, = 1, and = 4/3. Obviously, the particle number decreases and the average volume increases with time advancing due to coagulation.
By establishing the ansatz s(x) on the basis of the continuity of second derivation, the number of linear ill-conditioned system can be reduced significantly from 4m × 4m to (m + 3) × (m + 3) for (m + 1) nodes by using cubic spline, although only the continuity of s(x) is necessary in practice. Then, coupling with the asymptotic solutions of TEMOM  and the theory of self-preserving, the evolution of the PSD due to Brownian coagulation in the free molecule regime and continuum regime and its asymptotic behavior are obtained easily.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
This work was supported by the National Natural Science Foundation of China with Grant nos. 50806023 and 11572138.
S. K. Friedlander, Smoke, Dust, and Haze: Fundamentals of Aerosol Dynamics, Oxford University Press, London, UK, 2nd edition, 2000.
E. Aamir, Z. K. Nagy, C. D. Rielly, T. Kleinert, and B. Judat, “Combined quadrature method of moments and method of characteristics approach for efficient solution of population balance models for dynamic modeling and crystal size distribution control of crystallization processes,” Industrial & Engineering Chemistry, Research, vol. 48, no. 18, pp. 8575–8584, 2009.View at: Publisher Site | Google Scholar