A New Analytical Solution for Solving the Smoluchowski Equation Due to Nanoparticle Brownian Coagulation for Non-Self-Preserving System

The Smoluchowski equation has become a fundamental equation in nanoparticle processes since it was proposed in 1917, whereas the achievement of its analytical solution remains a challenging issue. In this work, a new analytical solution, which is absolutely different from the conventional asymptotic solutions, is first proposed and verified for nonself-preserving nanoparticle systems in the free molecular regime. The Smoluchowski equation is first converted to the form of moment ordinary differential equations by the performance of Taylor expansion method of moments and subsequently resolved by the separate variable technique. In the derivative, a novel variable, g = m0m2/m1, where m0, m1 and m2 are the first three moments, is first revealed which can be treated as constant. Three specific models are proposed, two with a constant g (an Analytical Model with Constant g (AMC), and a Modified Analytical Model with Constant g (MAMC)), and another with varying g (a finite Analytical Model with Varying g (AMV)). The AMC model yields significant errors, while its modified version, i.e., the MAMC model, is able to produce highly reliable results. The AMV is verified to have the capability to solve the Smoluchowski equation with the same precision as the numerical method, but an iterative procedure has to be employed in the calculation.


INTRODUCTION
Brownian coagulation is regarded as the most important inter-particle mechanism modifying the size distribution of particles in processes involving nanoparticle synthesis, polymerization, aerosol, emulsification and flocculation (Fox, 2008;Anand et al., 2012;Buesser and Pratsinis, 2012;Menz et al., 2014).Because of Brownian coagulation, the particle systems are always unstable, even in the absence of turbulence (Friedlander, 2000;Crowe et al., 2011;Sitarski, 2012;Wei, 2013).The theory for Brownian coagulation was originally devised for liquids by Smoluchowski (1917) and was subsequently applied in other chemical process and environmental fields.The governing equation for Brownian coagulation has been established for nearly a century in its integral-differential form (Müller, 1928), which is usually called the Smoluchowski equation (SE), but the solution for it remains a challenging issue due to its nonlinear integral-differential property (Debry et al., 2003;Wang et al., 2007;Wei and Kruis, 2013).Currently, the process of Brownian coagulation receives much more attention in the newly emerging fields of aerosol science, cloud science, supercritical fluid processes and nanoparticle synthesis engineering (Garrick, 2011;Yu and Lin, 2010a;Buesser and Pratsinis, 2012), where the quick and reliable predication of the size distribution of particles has become a critical issue (Seipenbusch et al., 2008).In such conditions, a reliable analytical solution for solving SE due to Brownian coagulation becomes necessary.
The solution for SE due to Brownian coagulation has been widely studied by researchers for many years, including the analytical solution and numerical solution, which have been reviewed in (Kruis et al., 2000;Vogel et al., 2014).Although the numerical solution is better than the analytical solution in precision (Wang et al., 2007;Yuan and Fox, 2011;Zhao and Zheng, 2009;Wei and Kruis, 2013;Menz et al., 2014), it consumes much more computational expenses and thus may not be suited for the quick determination of particle properties (Claudotte et al., 2010;Rigopoulos, 2010).
The analytical solution was firstly used to solve the SE due to Brownian coagulation by Smoluchowski (1917), but this work is limited to monodisperse particles.In fact, in the atmospheric and process engineering conditions, most particles are polydisperse, rather than monodisperse (Friedlander, 2000;Kim and Kim, 2002;Khalili et al., 2010;Buesser and Pratsinis, 2012), and the Smoluchowski analytical solution is only an ideal model and cannot be applied in real conditions.In the recent studies on engineered nanoparticle exposure, the Smoluchowski analytical solution was further developed and extended from monodisperse systems to bidisperse systems, but it is still limited to its inability to characterize the evolution of size distributions (Seipenbusch et al., 2008).An analytical method for solving size-resolved coagulation has been developed based on the use of similarity transformation for the size distribution function (Swift and Friedlander, 1964;Friedlander and Wang, 1966).It is determined that the shapes of the distribution at different times are similar when reduced by a scale factor, and the particle must be at the self-preserving status (Swift and Friedlander, 1964).The self-preserving theory is considered to be among the most important theories in particle science in which the geometric standard deviation (GSD) of the size distribution is kept constant, whereas it is inconsistent with the property of some particles because the GSD may be an arbitrary value more than 1.0.Lee and Chen (1984), Lee et al. (1984) and Park et al. (1999) found that it is feasible to assume the size distribution to follow the lognormal distribution, and based on this assumption, his group successfully developed a series of analytical models for solving SE due to Brownian coagulation at different size regimes.However, as they acknowledged in their works, the formulas for obtaining total particle number concentration, geometric standard deviation or geometric mean diameter have to be coupled and dependent of each other, thus, iterative techniques have to solve the coupled equations.
Since proposed in 2008 (Yu et al., 2008), the Taylorexpansion method of moments (TEMOM) has been widely employed in studies due to its simple mathematical construction of ordinary differential equations and highly reliable precision (Yu andLin, 2009a, b, 2010b;Xie et al., 2012).It has also been developed in different versions with specific requirements (Lin and Chen, 2013;Chen et al., 2014).Similar to the solution of Lee et al. (1984), the asymptotic behavior of the TEMOM ordinary differential equations (TEMOM ODEs) has been investigated in (Xie and Wang, 2013;Xie, 2014).In recent studies on the TEMOM ODEs, a novel variable, g = m 0 m 2 /m 1 2 , where m 0 , m 1 and m 2 are the first three moments, it was revealed that they can be treated as constant in the derivation; thus, the analytical solution for SE can be directly achieved.In fact, in the studies of the TEMOM ODEs, it was found that the geometric standard deviation of particle number distribution can be characterized by the definition ln 2 (σ g ) = (1/9)ln(g), where σ g is the GSD of the size distribution, although no log-normal distribution was employed.For GSD, it has been confirmed that it is a constant of 1.345 from the TEMOM method and QMOM method (Yu et al., 2008) and 1.355 from log-normal method (Lee et al., 1984) in the free molecular regime.The variance of GSD with the Knudsen number is presented in Fig. 1.It needs to emphasize the TEMOM ODEs with three order Taylor expansion series is the best selection for deriving the analytical solution due to its both simple mathematical structure and high precision, which has been discussed in (Yu et al., 2008).The newly proposed solution has been verified with high precision for self-preserving particles, but its ability for non-self-preserving particles still needs to be verified.
Therefore, in this work, an alternative analytical solution for solving SE due to Brownian coagulation will be proposed and verified as the particle deviates from the self-preserving status.The investigated particles are limited to the free molecular regime; thus, the newly proposed solution is only valid with a Knudsen number of more than 10, as shown in Fig. 1.To validate the analytical solution with precision, the numerical method from the highly reliable Runge-Kutta algorithm (NM) for the same ODEs will be used as a reference.

THEORIES
The integral-differential SE was firstly proposed by Müller in 1928 based on the Smoluchowski groundbreaking work for coagulation process (Müller, 1928;Smoluchowski, 1917), and it takes the following form, , , where n(v, t)dv is the particle number whose volume is between v and v + dv at time t, and β(v, v') is the collision kernel for two particles of volumes v and v'.In the free molecular regime, the collision kernel was derived from gas kinetic theory, which is expressed by (Friedlander, 2000), where T is the gas temperature and ρ is the mass density of the particles.The Taylor-expansion method of moments (TEMOM) is introduced in solving Eq. ( 1) with the closure model for the k-th moment (Yu et al., 2008), where u is the Taylor expansion point.In this work, u is defined to be m 1 /m 0 .The TEMOM ordinary differential Fig. 1.Changes in the geometric standard deviation (GSD, σ g ) as a function of the Knudsen number (Kn) over the entire size range.λ is the mean free path of gas molecules, 68.410 nm, and r is the particle radius.
equations (TEMOM ODEs) take the following expression (Yu et al., 2011), In this equation, g = m 0 m 2 /m 1 2 , and the k-th moment is defined as The zeroth moment, m 0 , represents the particle number concentration, which decreases with the growth of the particle size.The first moment, m 1 , remains constant in Brownian coagulation process and is proportional to the total particle mass.The second moment, m 2 , is usually used as an index to characterize the total light scattered, which increases with the growth of the particle size and polydisperity (Settumba and Garrick, 2003).If Eq. ( 4) are further disposed by a dimensionless solution with m k = M k m k0 and m k0 = Nv g0 k , then they have the following expression, where τ 1 = tNB 1 (v g0 ) 1/6 , and N and v g0 are the initial total particle number and geometric mean size, respectively.The same dimensionless time was used in the analytical solution for solving SE with the assumption of a log-normal distribution (Lee et al., 1984).If g can be regarded as a constant, the analytical solutions for Eq. ( 5) should be, where . M 00 , M 10 and M 20 are the zeroth, first and second moments at initial time, separately.In the above derivation, the separation of variables technique was used.Eq. ( 6) exhibits the novel property in that the first three moments are decoupled, and especially both M 0 and M 2 are only functions of the dimensionless time τ 1 .Here, it needs to point out in the derivation only three order terms of Taylor expansion series are reserved.
Although high precision is expected to achieve with higher order terms of Taylor expansion series, it is not feasible here because higher order terms of Taylor expansion series leads to too much complicated TEMOM ODEs whose analytical solution cannot be obtained.In fact, as shown in the following work, three order terms of Taylor expansion series meet the requirement of solving Smoluchowski Equation in precision.

COMPUTATIONS
The calculations were all performed on an Intel (R) Core i7-3573U CPU 2.5 GHz computer with 4GB memory.The fourth order Runge-Kutta method with a fixed time step of 0.001 was used to solve the TEMOM ODEs in the numerical solution.For all numerical and analytical solutions, the total time was up to 100.All the programs codes were written with the C++ Programming language under the platform of Microsoft Visual Studio 2008 compiler.The calculation of the relative error for any variables follows the definition (Yu et al., 2008), where ϕ is the arbitrary variable and ϕ NM is the referenced numerical variable.In the present study, a criterion is proposed for assessing the ability of the analytical method in solving SE due to Brownian coagulation, i.e., the investigated analytical method and the reference numerical method are considered to have the same precision in solving SE if the absolute relative errors, |RE%|, are less than 0.050.
In the calculation, all initial dimensionless moments take the following expressions: , w g = 3lnσ g0 and σ g0 is the initial GSD.In this work, the determinations of first three moments are reached by three specific models, two with constant GSD (AMC model and MAMC model) and another with varying GSD (AMV model).In the AMC and MAMC, the solutions for the first three moments are decoupled, whereas in the AMV the solutions for the formula, g = M 0 (τ 1 )M 2 (τ 1 )/M 1 (τ 1 ) 2 , as well as for the first three moments, have to be coupled where an iterative calculation is required.In the AMV, g needs to be obtained first with the first three moments at the last time step (τ 1 -∆t).It is then used in the acquirement of first three moments at the current time step (τ 1 ).In this way, the time step size (∆t) should be small enough so that no significant errors arise when introducing g to the solution for M 0 , M 1 and M 2 shown in Eq. ( 5).In the present study, the time step for the AMV is fixed at 0.001, same as the time step for the NM.Thus, the AMV is actually a finite analytical model.

RESULT AND DISCUSSION
In the free molecular regime, the successful achievement of an analytical solution for the SE is attributed to the novel dynamical property of particles that the GSD of the size distribution can be treated as a constant (Lee et al., 1984;Vemury et al., 1994;Friedlander and Wang, 1966).This treatment is valid for self-preserving particulates since the GSD in the free molecular has been confirmed to be a constant (Lee et al., 1984).In this work, it needs to be verified whether the newly proposed analytical method can also be valid at the non-self-preserving status, or whether the analytical method is able to produce the result within a reasonable error range for non-self-preserving particles.

The Zero Moment (M 0 )
The governing equations for particles dominated solely by Brownian coagulation mechanism must meet the following requirement since the total particle number concentration decreases and the second moment increases with time (Yu et al., 2008), In the free molecular regime, the effective range of the TEMOM ODEs in terms of GSD is confirmed to be at σ g ∈ (1.000, 1.601) through solving Eq. ( 9).In the acquirement of the effective range of GSD, the relationship between σ g and g, i.e., ln 2 (σ g ) = (1/9)ln(g), is employed.Within this range, the TEMOM model is expected to solve SE due to Brownian coagulation with any initial size distribution.Because the TEMOM ODEs are solved by an analytical solution, the effectiveness needs to be verified by comparison with the existing reliable methods.In the present work, the numerical solution with the highly reliable 4th-order Runge-Kutta algorithm is selected as the reference for comparison (Yu et al., 2008).Figs.2(a) and 2(b) separately show the RE% of the AMV and the AMC to the NM for M 0 as the particle initially follows the different distributions, i.e., different initial GSD.The RE% are obtained through solving Eq. ( 7) shown in Section 3. Six different initial GSD and σ g0 are distributed over the entire effective range of the TEMOM ODEs in terms of GSD.In the performance of the AMC, the GSD is always the same as the initial value, whereas in the AMV, the GSD is updated at the end of each interaction as discussed in Section 3. To make it more convenient to be analyzed, the horizontal axis is scaled by a log operation for Fig. 2(a).The same operation will be executed for Fig. 4(a) below.
Fig. 2(a) shows that the distribution initially deviates much more from the self-preserving size distribution, i.e., σ g0 deviates much more from 1.345, which was considered to be the value of the self-preserving status in the free molecular regime (Yu et al., 2008), the absolute value of RE% achieves larger values.The absolute RE%s are always limited to a very small range of variance, and they all have the same trend of finally converging to zero for all investigated distributions.The maximum RE%, 0.031, is reached at a distribution with σ g0 = 1.600, once the particle evolves for a short dimensionless time, ~1.000.The maximum absolute value of RE% for each distribution occurs at nearly the same dimensionless time, ~1.000.More importantly, the RE% for all distributions over the entire evolution time are less than the value of the criterion, 0.050, which is proposed in Section 3. Based on the definition for the criterion shown in Section 3, the AMV can be considered to produce M 0 with the same precision as the NM.
Compared to the AMV, the AMC produces relatively larger errors for M 0 as the same six distributions are investigated, which is clearly shown in Fig. 2(b).In the figure, except for the self-preserving size distribution with σ g0 = 1.345, which is a special case in the free molecular regime (Yu et al., 2008), the absolute values of RE% for the other distributions are all more than 0.05 once the particles evolve for a short dimensionless time, i.e., ~1.000.It is clear that as the initial GSD deviates much more from 1.345, the larger absolute values of RE% are generated, which is the same as the AMV shown in Fig. 2(a).In Fig. 2(b), the RE% approaches 0.250 in the distribution with σ g0 = 1.600.In such a case, the results generated by the AMC cannot be used as reliable results for solving SE due to Brownian coagulation.However, a novel property of the AMC is revealed in Fig. 2(b), in that for any distribution, the value of RE% quickly approaches its maximum or minimum value at the initial stage (the dimensionless time is below ~10.000) and then maintains this value at the later stage, which means the relative errors of the AMC to the NM for producing M 0 are unchanged once the particle evolves for a short time.More importantly, this difference depends only on one quantity, i.e., the initial GSD.Therefore, the modification of the AMC can be achieved by including an additional term accounting for the relative errors of the AMC to the NM.To achieve this, the relationship between the RE% of the AMC to the NM and the initial GSD for different distributions needs to be determined.
To determine the relationship between the initial GSD and the RE%, many more initial distributions up to 24 are investigated by the performance of the NM calculation and the AMC calculation.The RE% of the AMC to the NM relevant to the initial GSD are shown in Fig. 3.The data for the initial GSD, as well as the corresponding RE%, have been attached online as Data-Fig.3. It shows that the variance of RE% for M 0 is limited in the range from -0.240 to 0.138.The relationship between the initial GSD and RE% can be represented by a suitable mathematical equation by performing the 6 th degree polynomial fitting technique.
In mathematics, the fitting function fits much more to the initial data if the polynomial degree is selected as a larger value.Here, the mathematical equation relating the initial GSD to RE% is, f(σ g0 ) free-m0 = 0.00018895ω 6 + 0.001492ω 5 + 0.0043408ω 4 + 0.0046002ω 3 -0.031742ω 2 -0.13938ω -0.012157 (10) with ω = (σ g0 -1.359)/ 0.16273.f(σ g0 ) free-m0 is the function of σ g0 representing the relationship between the initial GSD and the RE% over the entire effective regime of TEMOM ODEs.Correspondingly, the residuals in the performance of curve fitting are also shown in Fig. 3, whose absolute values are always less than 0.0001, indicating that the fitting Eq. ( 10) is high reliable.Therefore, the analytical solution for M 0 in Eq. ( 3) can be modified as following, The Second Moment (M 2 ) The second moment, M 2 , is an important quantity

Fig. 3.
The RE% of the AMC to the NM for M 0 relevant to the initial GSD, as well as the corresponding fitting curve, f(σ g0 ) free-m0 (a); the residuals produced in the performance of the 6-th degree polynomial fitting operation (b).characterizing particle polydispersity and particle mean size (Friedlander, 2000).Because only the Brownian coagulation mechanism is involved, the second moment always increases with the growth of the particle size and polydisperity (Settumba and Garrick, 2003).In this section, the same six distributions as in Section 4.1 are selected and investigated.The RE% of the AMV and the AMC to the NM for M 2 are shown in Figs.4(a) and 4(b), respectively.In Fig. 4(a), the absolute RE% for all distributions initially increases to each maximum value and then quickly decreases and converges to zero.The lag times to reach the maximum value for all distributions are nearly the same, ~1.000.Over the entire time range, the distribution with the initial GSD deviating more from the self-preserving value of 1.345 has larger RE%.The absolute RE%, 0.069, is achieved in its maximum value at the distribution with σ g0 = 1.600 and at dimensionless time, ~1.000, but it quickly decreases to be less than the value of criterion at the dimensionless time of ~1.855.This means that under the criterion, the AMV is able to produce M 2 with the same precision as the NM once the particle evolves for a short time.
The RE% of the AMC to the NM for M 2 are shown in Fig. 4(b).In this figure, the RE% quickly increases to a constant at dimensionless time of ~10.000 and later remains at this value at the later stage for all six distributions.This is similar to the zeroth moment shown in Fig. 2(b).An important conclusion can be drawn that the RE% of the AMC to the NM for M 2 are also unchanged once the particle evolves for a short time.Therefore, it is feasible to obtain the relationship between the initial GSD and the RE% for M 2 , which is shown in Fig. 5.The curve obtained by performing the 7th degree polynomial fitting technique, as well as the corresponding residuals, are also shown in this figure.The data for the initial GSD and the corresponding RE% have been attached online as Data-Fig.5.The residuals are always below 0.010, ensuring that the fitting equation is reliable, although it is not as high as in Fig. 3

in precision.
Here, the fitting equation obtained by the performance of the fitting technique is, f(σ g0 ) free-m2 = 0.0037563ϵ 7 + 0.010874ϵ 6 -0.016636ϵ 5 -0.088877ϵ 4 -0.12637ϵ 3 -0.10002ϵ 2 + 0.0075218ϵ -0.00032143 (12) with ϵ = (σ g0 -1.359)/0.16273.When Eq. ( 12) is introduced in the solution for M 2 shown in Eq. ( 6), the modified solution should be, Finally, the modified AMC, i.e., the MAMC, takes the following expression, The MAMC shown in Eq. ( 14) reserves the advantage of its original version in the decoupled solution for the first three moments, but it generates more precise results because of the additional modified terms accounting for the difference between the original AMC and the NM.Once the first moments are achieved by solving Eq. ( 14), the key quantities of any particle including the total particle number concentration, the total particle volume concentration, the geometric mean size, and the geometric standard deviation of size distribution, are also determined (Lee and Chen, 1984;Settumba and Garrick, 2003;Yu et al., 2008).

Evolution of the Size Distribution Relevant to the AMV
The GSD is the most important index characterizing the polydisperse property of particles (Lee et al., 1984), and thus, the ability to exactly trace the evolution of GSD with time is the criteria to distinguish good and bad models (Pratsinis, 1988).Figs.6(a) and 6(b) separately show the varied GSD with time for both AMV and NM, as well as the corresponding RE% of the AMV to the NM.It is found from Fig. 6(a) that the GSD from the AMV are nearly the same as the NM, although the NM is calculated by fourth order Runge-Kutta interactions, consuming a relatively higher computational cost.For the four investigated distributions shown in Fig. 6(a), the GSDs are found to quickly converge to a constant, which is consistent with the self-preserving theory first found by Friedlander and Wang (1966) and later confirmed in Lee et al. (1984) and Yu et al. (2008).In Fig. 6(b), similar to the solution for both M 0 and M 2 shown in Figs.2(b) and 3(b), the RE% of the AMV to the NM for GSD are calculated by solving Eq. ( 7).The variance of RE% for GSD with time has the same trend as M 0 and M 2 .The absolute RE% initially increases to approach the maximum and then decreases until it converges to zero.All distributions achieve their own maximum value for GSD at nearly the same time, ~1.000, which is consistent with that for M 0 and M 2 shown in Figs.2(a) and 4(a).It is also shown the RE%s for GSD are limited to a very small range of variance with the maximum of 0.015, which is less than the value of criterion proposed in this work shown in Section 3. The AMV is thus verified to be able to characterize the GSD with the same precision as the NM.
Based on the above analysis of the three key statistical quantities, i.e., M 0 , M 2 and σ g , the AMV is verified to be able to solve SE due to Brownian coagulation for non-selfpreserving particles.Thus it has strong potential to replace the NM due to the inherent advantage of the analytical solution.

Verification of the MAMC Model
In Sections 4.1 and 4.2, the MAMC has been proposed as the modified version of the AMC, but its ability for solving SE due to Brownian coagulation needs to be verified quantitatively.To validate the MAMC, the NM is selected to be the reference for comparison.In Fig. 7, the comparison of the results produced by the NM, the MAMC and the AMC for M 0 , M 2 and σ g is exhibited.As examples, two different initial distributions, i.e., σ g0 = 1.200 and 1.500, are selected.For both M 0 and M 2 , it is shown that the results  from the MAMC quickly converge to the NM as the particle evolves for dimensionless time, ~10.000, whereas the AMC deviates clearly from the NM.For GSD, the MAMC also converges to the NM after the same time.For the two investigated initial distributions, the GSDs are found to converge to the same value, 1.345, which is also the value of a particle system dominated by Brownian coagulation with the self-preserving size distribution status (Yu et al., 2008).
To verify whether the MAMC has the ability to solve SE due to Brownian coagulation over the entire effective GSD range, many more distributions are required to be investigated and tested.To achieve this, six representative distributions, which span the entire effective GSD range, are selected.
The RE%s of the MAMC to the NM for M 0 , M 2 and σ g with time are obtained by solving Eq. ( 7) and are presented in Fig. 8.For all six representative distributions, the variances of RE% for M 0 , M 2 and σ g with time exhibit the same trend, and ultimately, they all converge to zero.This result indicates that the differences between the NM and MAMC for the above three key statistical quantities, i.e., M 0 , M 2 and σ g , decrease with time until the difference vanishes.In the figure, the absolute values of RE% for M 0 , M 2 and σ g for each distribution are less than the criterion value, 0.050, once the dimensionless time of evolution is more than ~10.000.In addition, it is found that the closer to 1.345 the GSD is the less lag time it takes to approach the criterion.
In fact, the dimensionless time, τ 1 = 10.000,corresponds to the physical time, t = 1.115 × 10 -3 s, under the conditions in Section 3. In such a case, an important conclusion can be drawn that the MAMC has the same capability to solve SE due to Brownian coagulation as the NM once the particle evolves for a very short time.For longer times of evolution, there is higher precision of the MAMC.The MAMC has a significant advantage because it does not require interactive calculation and, more importantly, the acquirement of moments is decoupled.In conclusion, the MAMC has potential to replace the NM in processes involving Brownian coagulation due to its superior properties, and therefore, it is regarded as the ideal analytical solution for solving SE due to Brownian coagulation.

Fig. 2 .
Fig. 2. The comparison of relative errors of the AMV (a) and the AMC (b) to the NM with six different initial size distributions for M 0 .

Fig. 4 .
Fig. 4. The comparison of the relative errors of the AMV (a) and the AMC (b) to the NM with six different initial size distributions for M 2 .

Fig. 5 .
Fig. 5.The RE% of the AMC to the NM for M 2 relevant to the initial GSD, as well as the corresponding fitting curve, f(σ g0 ) free-m0 (a); the residuals produced in the performance of the 7 th degree polynomial fitting operation (b).

Fig. 6 .
Fig. 6.The comparison of GSD produced by the NM and the AMV with time (a) and the RE% of the AMV to the NM (b).