Numerical Error Analysis of Solvers Using Cumulative Number Distribution with Volume-Ratio Grid Spacing in Initially Particle-Free Nucleation-Condensation Systems

The numerical accuracy of solvers has not been fully investigated in initially particle-free nucleation-condensation systems, which are important in modeling particle generation in the free-molecular regime. The numerical error analysis of several schemes using the cumulative number function is conducted on the basis of Williams’ analytical solution and the Courant– Friedrich–Lewy (CFL) condition in volume-ratio grid spacing of the particle size. In this system, because the CFL number (obtained from the particle growth rate, grid size, and time step size) depends strongly on the particle size, we need to check the CFL condition (determining the numerical stability and accuracy) over a whole size range before the simulation. This has not been fully discussed in aerosol modeling. The present study introduces a new unique constant (DCFL) throughout the simulation in the volume-ratio grid spacing. In the sensitivity experiments, we can find the similarities of the simulated size distribution and its numerical error, which are determined by the DCFL. We can thus easily check the CFL condition, and choose an appropriate time step size, to satisfy the required accuracy in the nucleation-condensation module of the program.


INTRODUCTION
Analytical solutions (e.g., Smoluchowski, 1918;Ramabhadran et al., 1976;Seinfeld and Pandis, 1998) are used as exact values in benchmark tests of solvers in modeling condensation and coagulation systems (e.g., Jacobson and Turco, 1995;Mitrakos et al., 2007;Mohs and Bowman, 2011).Williams (1983) found analytical solutions for aerosol dynamics systems including both nucleation and condensation processes (hereafter termed "the nucleationcondensation system").In particular, initially particle-free nucleation and condensation are important in modeling particle generation.In the nucleation-condensation system, the nucleation rate is a source term using Dirac's delta function at the critical size, and the total particle number and volume increase over time.As a sharp size distribution of particles develops, numerical overshooting is likely to occur.However, the benchmark test has not yet been fully applied to the initially particle-free nucleation-condensation system.Thus, we need to estimate the numerical error of the solvers in detail, to ensure the conservation properties in the simulation of the aerosol dynamics.
The condensation is mathematically the same as the advection equation.In numerical fluid dynamics, the Courant-Friedrich-Lewy (CFL) condition (e.g., Ferziger and Meric, 1996) is introduced to discuss the numerical stability and accuracy.Chock and Winkler (2000) and Nguyen and Dabdub (2002) used constant CFL numbers in the benchmark tests of their condensation/evaporation solvers.However, in volume-ratio grid spacing of the particle size, the CFL condition has not been fully discussed in estimating the numerical error of the solver, because the CFL number rapidly changes with an increasing particle size.
In the present work, the exact solution of Williams (1983) and the CFL condition are applied to a benchmark test for evaluating the accuracy of several schemes using the cumulative number function.The use of the cumulative number has the merits that (1) it is difficult to overshoot even for an extremely sharp distribution of the particle size and (2) it is easy to apply the mathematical constraints of the conservation of the total particle number and the monotonicity of the cumulative number (Yamamoto, 2004).Because it is mathematically easy to extend the nucleation to the condensation equation system, the present work considers one program module, into which both nucleation and condensation processes are implemented.In this case, we need to modify the boundary condition at the critical size of nucleation, at which the cumulative number changes with time.At the present stage, the influence of the time dependent boundary on the cumulative number function has yet to be investigated in initially particle-free nucleationcondensation systems.Thus, we need to investigate the conservative properties of total particle number and volume simulated in this system.
In this article, the analytical solution is reviewed and the physical quantities for the numerical error estimation are defined in the next section.The numerical methods are described in the third section and appendix.The similarities of the size distributions and numerical errors are discussed using a unique constant D CFL in the fourth section.Here we find that D CFL is useful for determining an appropriate time step size in the model.

APPLICATIONS OF ANALYTICAL SOLUTION AND CFL CONDITION
The governing equation of the function for the particle here, the growth rate is simplified as where v S is the smallest volume of the system (i.e., the critical size of nucleation), v L is the largest volume of the system, and  is 2/3 for the free-molecule regime (e.g., Williams, 1983;Pratsinis, 1988).The  value is derived in the small-particle limit of K N >> 1, where K N is the Knudsen number for the precursor gas molecules in air.We consider the idealized conditions that (i) the initial particle distribution is zero for all particle sizes ( ( ,0) 0 ), (ii) the nucleation rate is constant (P = constant), and (iii) the free-molecule regime is considered ( = 2/3).The relation between time and the largest particle size (v L ) is obtained from Eq. ( 2), which is rewritten as From the integration of Eq. ( 3) with respect to time, This implies that the diameter of the largest particle ( v L 1/3 ) linearly increases with time at the rate of A/3.Because the total particle number linearly increases with time at the constant rate of P, the cumulative number distribution is shown in Fig. 1 (a).In this case, particles in the linear size distribution (line a-b) at t = ∆t undergo advection at a constant rate of A/3 (thick arrow) everywhere in the space of C -v 1/3 , and the linear distribution becomes that of line c-d at t = 2∆t.From the differentiation of the cumulative size distribution with respect to v, f(v, t) is written as (5) The final form of Eq. ( 5) is the same as the analytical solution of Williams (1983).The volume-weighted size distribution f(v, t)v and cumulative number function C(v, t) are rewritten as, in the size range from v S to v L .Because v L increases with time in Eq. ( 4), the domain of integration is time dependent in Eq. ( 7).Thus, although the integrand f(v, t) does not depend on time, the integral C(v, t) is time dependent.The similarities of the normalized size distributions are shown in Figs. 1(b) and (c).The normalized f(v, t)v increases to (v L /v S ) 1/3 along the curve, while the normalized C(v, t) becomes the time-independent similar form.For estimating the numerical error of solver, we introduce the following quantities, The exact solutions of the total particle number and volume are obtained from Eqs. (4), ( 8) and (9).Eq. ( 1) is rewritten as in the following coordinate system (Berry, 1967;Middleton and Brock ,1975) as J is defined as a logarithmic-scale measure of particle size.For b = 1, Eq. ( 11) has the volume-ratio form of a J-1 , where a is the volume ratio of v(J + 1)/v(J).For a = e (=  (v, t).The gray thick and dashed arrows indicate the growth rate of particles in this space (A/3) and nucleation rate (P), respectively.
2.71828…), Eq. ( 11) has the exponential form of exp[b(J -1)], where b is a logarithmic-scale increment of the volume grid.Both the volume-ratio and exponential forms are well used for grid spacing of the particle size in simulations of aerosol dynamics.In order to discuss the numerical accuracy of solver, the present work introduces the Courant-Friedrich-Lewy number The final form of Eq. ( 14) can be applied to the freemolecular regime.Under the CFL condition that C CFL is smaller than unity, the advection for one time step does not exceed one grid size.Generally, smaller C CFL leads to more stable and accurate solution for a given grid spacing.
In this case, as the particle size decreases, C CFL becomes larger.Thus C CFL has the largest value at the critical size.
Because C CFL strongly depends on the particle size, we newly introduce a unique constant D CFL throughout a given simulation: The characteristic length D CFL is determined by only a few experiment parameters (A, ∆t, a and b), and is not dependent on t and v.It is noted that C CFL is simply rewritten as C CFL = D CFL /D where D is the particle diameter.Because an inequality of D CFL < D satisfies the CFL condition (C CFL < 1), D CFL indicates the minimum diameter satisfying the CFL condition.

METHODS
The present work estimates the numerical error of nucleation-condensation solvers using F(J, t) and the cumulative number function C(J, t).The base model is the same as Yamamoto (2004), in which the merits for the uses of cumulative number function and semi-Lagrangian scheme were described in detail.In the present case including the nucleation, the boundary condition in the J space is time dependent at the critical size.The cumulative number function C(J, t) and its differential with respect to J are ( ) respectively.The governing equation of C(J, t) is: in the condensation size region of J > 1.At the critical size of J = 1 (v(1) = v S ), we need to consider two cases.For evaporation (u(1) < 0 and For nucleation (u(1) > 0 and P > 0), Here u(J MAX )F(J MAX , t) represents the removal of particles from the system (1 ≤ J ≤ J MAX ) toward larger particles (J > J MAX ) via condensation, or the input from larger particles to the system via evaporation.C(J, t) has mathematical constraints (the total particle number and monotonicity): These two mathematical constraints are used to correct the overshoot/undershoot of F(J, t) and conserve the total particle number C(1, t).
The present study applies a cubic interpolated propagation advection (CIP) scheme (Yabe and Aoki, 1991;Yabe et al. 2001) to the semi-Lagrangian condensation method (CIP-SL, Yamamoto, 2004), because the CIP is one of the most accurate methods for simulating advection (which has essentially the same equation system as condensation).For comparison, several schemes used in numerical fluid dynamics are also applied to the advection procedure in our previous method (Appendix).

RESULTS AND DISCUSSIONS
The present work uses the parameters of a = 1.2, b = 1.0 and r S = 0.005 m for the grid spacing, A = 10 -8 cm/s and  = 2/3 for the growth rate, and P = 10 -2 #/cm 3 /s for the nucleation rate.The time step size ∆t of the time integration is set at 0.4 s.The size distribution starts from an initially aerosol-free state (F(J, 0) = C(J, 0) = 0) and is calculated to 1600 s (Case I B in Table 1).Fig. 2 shows the numerical and analytical solutions of f(v, t)v in this case.The cut-off front is formed at the largest size v L , and shifted toward larger size over time.The simple Euler's (ELR) and firstorder upward (UPW1) advection schemes lead to strongly noisy and diffusive size distributions, respectively.Although the third-order upward (UPW3) advection scheme (Kawamura and Kuwahara, 1984) somewhat improves the simulated size distribution, there remains the numerically diffusive problem around the cut-off front.In the cases of the CIP-SL and UPW3 schemes without the constraint of the cumulative number (gray circles and dashed curves in Fig. 2(b), respectively), undershoots to negative values occur near the front.Thus, it is difficult to simulate the size distribution around the cut-off of v L .On the other hand, the CIP-SL scheme with correction using the cumulative number distribution (Yamamoto, 2004) agrees well with the exact solution.This method completely prevents undershoot near the cut-off front, owing to the simple mathematical constraints of the cumulative number.
Fig. 3 shows the histories of the total particle number and volume and their relative errors in the CIP-SL and UPW3 schemes.The relative errors of the total particle number and volume are strictly defined as where the subscripts num and exa refer to numerical and exact solutions, respectively.Although the total particle number is well conserved because of the cumulative function constraint, the total volume is not exactly conserved in the both cases.dV/V has initially large values of ~10%, because the denominator V exc is very small in Eq. ( 24) at a few time steps from the start.The initial error rapidly decreases immediately after the start, and gradually approaches 0.5% in the CIP-SL scheme.Although it is difficult to simulate the particle size distribution in the presence of the cut-off front, the numerical solution better preserves the physical quantities.On the other hand, because the numerical diffusion occurs throughout the simulation using the UPW3 scheme, the asymptotic error of the total volume has large percentages of ~2%.The sensitivities of the numerical solution to grid spacing (a) and time step size (∆t) are now investigated.Because there are the numerical constraints (the artificial diffusion and CFL condition) for the UPW3 scheme, it is difficult to cover the wide parameter ranges for grid spacing and time step size.Thus, hereafter, we consider only the numerically stable CIP-SL method.Even when the intervals of the particle size and time step are large, the results are in good agreement with the exact solution (Fig. 4).In addition, several experiments are also conducted under numerical conditions that P varies from 10 -2 to 10 +1 #/cm 3 /s (Exp.III B ) and A varies from 10 -8 to 10 -10 cm/s (Exp.V B ).The numerical solutions (Fig. 5) have a similar form to the exact solutions in Figs. 2 and 4, though the distribution functions differ.As shown in Fig. 1(b), the size distribution normalized by PA -1 v S -1/3 is the same among Exps.I B , III B , and V B (Figs. 2 and 5).
Next, we discuss the time histories of the numerical errors for the sensitivity experiments listed in Table 1.The errors of the total volume rapidly decrease with time in the early stages of Exps.I, III and V (the left graph in Fig. 6).After 800 s (or 8000 s in Exp V), the error is time independent for large a (≥ 1.2).For a very fine grid (a = 1.05), the error gradually decreases over time after the initial rapid decrease.As the grid spacing decreases, the error also decreases.This is because the interpolation becomes more accurate with a decreasing grid interval.
In Exps.II, IV and VI (the right graph in Fig. 6), the errors of the total volume gradually decrease over time for a large time step size, while they rapidly decrease in the early stage for a small time step size.Finally, the errors approach 0.3% for small ∆t (0.1 s in Exps.II and IV and 10 s in Exp.VI).The accuracy increases with a decreasing time step size.In a simulation that requires high accuracy in the early stage, we should use a small time step size.However, an extremely small ∆t does not greatly improve accuracy, as the error approaches an asymptotic value (see the results of 0.1 and 0.4 s in Exp.II).Thus, we can use an appropriately large ∆t for long-term simulation, because the error gradually decreases toward an asymptotic value (~0.3%).
The histories of the errors are the same among Exps.I B , III B and V B (dashed curve in Fig. 6), although the experimental parameters are different.This implies that the numerical error has a similar pattern in the different experiments.The non-dimensional parameter C CFL is useful for discussing the error analysis in the same manner as numerical fluid dynamics.However, C CFL strongly depends on the particle size in this system, and thus the value at the cut-off front changes over time.On the other hand, D CFL is a constant and unique parameter throughout a given simulation, and indicates the minimum particle size satisfying the CFL condition (C CFL = D CFL /D < 1).Thus, the numerical errors became same among the experiments with a given D CFL in Table 1 (e.g., D CFL = 2.72 × 10 -4 μm for Exps.I B , III B and V B ).For the same time step size in Exps.I, III and V, the grid spacing parameters of a = 1.05, 1.2, 1.5, 1.8 and 2.1 correspond to the D CFL values of 1.02 × 10 -3 , 2.72 × 10 -4 , 1.22 × 10 -4 , 8.44 × 10 -5 and 6.69 × 10 -5 μm, respectively.As the grid spacing becomes coarser, the minimum particle size satisfying the CFL condition (D > D CFL ) decreases.
For the grid spacing of a = 1.2 (Exps.I, III and V), the sensitivity of the numerical error to time step size can be explained as the sensitivity to D CFL .Fig. 7(a) shows the CFL number plotted for each D CFL .C CFL strongly depends on the particle size, and the CFL condition is satisfied over the whole size range in all the simulations.When D CFL is introduced, the CFL condition is easily checked by only an inequality of D CFL < D. The numerical error is improved with a decreasing D CFL in Fig. 7(b), because smaller ∆t leads to smaller D CFL .The time histories of the numerical errors are represented as a unique function with respect to D CFL (Fig. 7(b)).The numerical errors are almost the same for D CFL smaller than 10 -4 .Thus, when the accuracy within 1% is required in the experiment of a = 1.2, we should choose a time step size satisfying the D CFL value of 10 -4 μm.D CFL is useful for determining the unique similarities of the size distributions and numerical errors, and thus for choosing an appropriate time step size.

CONCLUSIONS
The author conducted the numerical error analysis for solvers using cumulative number function in initially particlefree nucleation-condensation systems, based on Williams (1983)'s analytical solution and the CFL condition.Because C CFH strongly depends on particle size in the case of volumeratio grid spacing, we need to check the CFL condition before calculating the condensation (advection) term over the whole size range.This has not been fully discussed in aerosol modeling.In the present study, D CFL is newly introduced in the free-molecular regime.The CFL condition can be easily checked by an inequality of D CFL < D, because D CFL is a unique constant for a given simulation.
Around the cut-off front forming at the largest size, numerical diffusion and undershoot are likely to occur in the size distributions simulated by the high-accuracy methods (UPW3 and CIP-SL).Corrections using the cumulative number function prevent undershoot in the size distribution, and exactly conserve total particle number.The present study investigates the numerical error of the semi-Lagrangian solver using the CIP scheme (Yabe and Aoki, 1991) and simple mathematical corrections for preventing undershoot and exactly conserving total particle number (Yamamoto, 2004).The simulated total number and size distribution agree well with the exact ones.As the grid spacing and time step size are set at smaller values, the accuracy of the solver increases in the idealized nucleation-condensation system.Although the relative error of the total particle volume is not small at a few time steps from the start of simulation, it rapidly decreases immediately after the start and gradually approaches a small value.
Under the conditions of different P and A, we can see a similar form of the size distributions.In such cases, the numerical errors are also the same in the similar experiments Table 1.Experimental conditions.
The error analysis and its related similarity are expected to contribute to the inter-comparisons of different numerical schemes (Seigneur et al., 1986;Zhang et al., 1999).

ACKNOWLEDGMENTS
This work was partially funded by a Grant-in-Aid for Scientific Research (KAKENHI No. 23540514).
Step 1: Time Integration for J = 1 Advection from J < 1 is not considered, because particles smaller than the critical size (the smallest size) do not exist in the nucleation-condensation system.Thus, From this equation, the solution at time t + ∆t is Step 2: Numerical Procedure for J > 1 The CIP scheme (Yabe et al., 1991) is applied to the advection terms of Eqs. ( 10) and ( 18), because it is one of most accurate schemes.The governing equations are Eqs. ( 27) and ( 28) are separated into two steps using intermediate variables C * and F * : for the advection step, and: for the non-advection step.u(J)/J is approximated by the central difference.To obtain the intermediate variables using the semi-Lagrangian scheme with cubic interpolation, c J and f J are introduced: where x = J -u(J)∆t.For an integer j (j-th grid) satisfying the inequality j ≤ J -u(J)∆t < j + 1 (j ≤ x < j + 1): where The non-advection step is carried out after calculating C * (J) (= c J n+1 ) and F * (J) (= f J n+1 ).At the boundary J = 1, F(1, t + ∆t) is calculated as using the new values obtained in Step 1 and the values obtained for J > 1 in the non-advection step (Step 2).If C(J MAX , t + ∆t) > 0 owing to u(J MAX )F(J MAX , t) ≠ 0, we need to apply the transformation: C(J, t + ∆t) = C(J, t + ∆t) -C(J MAX , t + ∆t) (40) to preserve the total particle number in the size region between 1 and J MAX .The transformation does not affect F(J, t + ∆t) (= -dC/dJ).
For comparison, the cumulative number function in Eq. ( 27) is solved by various schemes.In the case of the thirdorder upwind (UPW3) advection scheme (Kawamura and Kuwahara, 1984), the advection term is represented as, For the first-order upwind (UPW1) advection scheme, In these cases, F(J, t) is obtained from -C(J + 1, t) + C(J, t).

Step 3: Corrections for Preventing Overshoot/Undershoot and Conserving Particle Number
After the numerical procedure in Step 2 is carried out in each time step, the mathematical constraints ( 21) and ( 22) are applied to the correction: C(J, t + ∆t) = C(1, t + ∆t) and F(J, t + ∆t) = 0 if J ≤ J TOT , (44) C(J, t + ∆t) = 0 and F(J, t + ∆t) = 0 if J ≥ J ZERO , (45) C(J, t + ∆t) = C(J + 1, t + ∆t) and F(J, t + ∆t) = 0 if C(J, t + ∆t) -C(J + 1, t + ∆t) < 0, (46) where J TOT is the largest integer J satisfying the mathematically incorrect inequality C(J, t + ∆t) ≥ C(1, t + ∆t), and J ZERO is the smallest integer J satisfying the mathematically incorrect inequality C(J, t + ∆t) ≤ 0. C(J ≤ J TOT , t + ∆t) is corrected for overshoot (because it is greater than the total particle number), while C(J ≥ J ZERO , t + ∆t) is corrected for undershoot (because it becomes negative).Furthermore, when F(J, t + ∆t) < 0, the size distribution is approximated by the finite difference of the corrected C(J, t + ∆t): F(J, t + ∆t) = -C(J + 1, t + ∆t) + C(J, t + ∆t) if F(J, t + ∆t) < 0 for J ≥ 1. (47) The above corrections ensure the conservation of the total particle number and the prevention of negative C(J, t + ∆t) and F(J, t + ∆t).

Fig. 1 .
Fig. 1.Analytical solutions of (a) cumulative number function C(v, t), (b) normalized volume-weighted size distribution f(v, t)v, and (c) normalized cumulative number function C(v, t).The gray thick and dashed arrows indicate the growth rate of particles in this space (A/3) and nucleation rate (P), respectively.

Fig. 2 .Fig. 3 .
Fig. 2. Numerical and analytical solutions of f(v, t)v at 0, 400, 800, 1200 and 1600 s for (a) ELR and UPW1 and (b) UPW3 and CIP-SL schemes.The analytical solutions are shown by gray solid curves.The solid circles, short-dashed, long-dashed, thick curves indicate the CIP-SL, UPW3, UPW1 and ELR schemes, respectively.The gray color shows the schemes without the constraint of the cumulative number.