Parallel Monte Carlo Simulation of Aerosol Dynamics

Abstract

A highly efficient Monte Carlo (MC) algorithm is developed for the numerical simulation of aerosol dynamics, that is, nucleation, surface growth, and coagulation. Nucleation and surface growth are handled with deterministic means, while coagulation is simulated with a stochastic method (Marcus-Lushnikov stochastic process). Operator splitting techniques are used to synthesize the deterministic and stochastic parts in the algorithm. The algorithm is parallelized using the Message Passing Interface (MPI). The parallel computing efficiency is investigated through numerical examples. Near 60% parallel efficiency is achieved for the maximum testing case with 3.7 million MC particles running on 93 parallel computing nodes. The algorithm is verified through simulating various testing cases and comparing the simulation results with available analytical and/or other numerical solutions. Generally, it is found that only small number (hundreds or thousands) of MC particles is necessary to accurately predict the aerosol particle number density, volume fraction, and so forth, that is, low order moments of the Particle Size Distribution (PSD) function. Accurately predicting the high order moments of the PSD needs to dramatically increase the number of MC particles.

1. Introduction

Population balance equations (PBEs), describing particulate entities conservation, have applications in many branches, such as aerosol dynamics [1, 2], polymerization [3], and so forth. Here, the discussion mainly focuses on aerosol dynamics. Aerosol refers to a colloid suspension of fine solid particles or liquid droplets, which can be found in clouds, air pollution such as smog and smoke, soot in combustion flames, cement dust, and so forth. Analytical solutions to the PBEs are available for only a few specific cases [4, 5]. The most employed general method and its variants for solving the PBEs is based on dividing the particle size domain into sections as developed by Gelbard et al. [6]. Usually, these direct discretization methods work well for either coagulation or condensation dominated aerosol dynamic processes, but they encounter various problems [7] when dealing with combined nucleation, coagulation and condensation, for example, numerical error in specifying a given distribution in a section, non-positiveness and non-conservation of mass and number concentration due to numerical diffusion.

Other than solving the discretized PBEs, a widely used alternative approach has been the method of moments, which solves a group of moments equations derived from the PBEs. These moments (generally several low order moments) provide important information on the PSD function, that is, number density, volume fraction, polydispersity, and so forth. However, owing to the nonlinearity of the PBEs, the governing equations for lower order moments contain higher order moments, which are not closed after cutting off. Various closure methods have been developed, such as unimodal log-normal approximation [8], quadrature/direct quadrature [9, 10], interpolative closure [11], Taylor expansion [12 –14], and so forth. Although the method of moments is computationally highly efficient, the complexity of the closure makes it difficult to accommodate complex physical models for aerosol dynamics with high flexibility.

Monte Carlo (MC) simulation does not “solve” the PBEs, but to mimic the evolution of aerosol particles through a stochastic particle system. Modeling coagulation by stochastic particle systems was introduced by several authors [15 –17], which is usually called Direct Simulation Algorithm (DSA). Eibeck and Wanger [18] proposed a stochastic algorithm that models the weighted Smoluchowki's coagulation equation, which exhibits lower variance than the DSA for predicting the tail of the PSD. There are modifications to the DSA aiming for higher computational efficiency. The majorant kernel method [19, 20] adopts an acceptance-rejection technique to speed up the simulation. The increased efficiency highly relies on finding a “good” majorant kernel. The so-called τ leaping method [21, 22] is an approximate accelerated stochastic simulation, which enables multiple coagulation events at a time. The extension of the DSA to model general aerosol dynamics (nucleation, surface growth, etc. also included) goes to diverse directions, time or event driven, constant volume or number MC simulations [23 –25]. Most of these algorithms share the same idea [26] to randomly choose a aerosol dynamic process (either nucleation, or coagulation, or surface growth) at a time to determine the aerosol particles evolution. Other than selecting all processes randomly, Patterson et al. [27] proposed the linear process deferment algorithm, which separates the surface reaction from other dynamic processes, and models the surface reaction independently for some steps, and then includes the other processes. Their idea to accelerate the simulation is similar to the τ leaping method in some sense.

In this paper we propose a parallel hybrid algorithm of stochastic simulation and deterministic integration for solving the PBE, which is called Operator Splitting Monte Carlo (OSMC). OSMC adopts operator splitting techniques (similar idea is also found in Patterson et al. [27]) to separate coagulation from other processes. Coagulation is simulated by the DSA, but all other processes are modeled by deterministic integration of their dynamic equations. Stochastic simulation is much more efficient than directly solving the Smoluchowski integro-differential equation with traditional discretization [23]. However, using stochastic method for nucleation and surface growth is far more expensive than their corresponding deterministic solution methods. The combination of stochastic and deterministic methods makes OSMC very efficient and flexible for aerosol dynamics simulation. OSMC is used to simulate a series of well-chosen testing cases. Detailed comparison with analytical and/or other numerical solutions is provided. The paper first introduces the OSMC algorithm, then discusses the simulation results of the testing cases.

2. Methodology

2.1. Governing Equations

The evolution of the PSD function, n (v,t), is described by the PBE (also called the general dynamic equation) [1]

\frac{\partial n}{\partial t} + \nabla \cdot n \vec{u} = \nabla \cdot D \nabla n + {[\frac{\partial n}{\partial t}]}_{nucl} + {[\frac{\partial n}{\partial t}]}_{growth} + {[\frac{\partial n}{\partial t}]}_{coag},

(1)

where $\vec{u}$ is the gas flow velocity, and D is the diffusion coefficient. The diffusion term is usually negligible owing to large Schmidt number for aerosol particles (Under standard atmospheric conditions, spherical particles with a diameter of 10 nm (resp. 100 nm) have a Schmidt number equal to 290 (resp. 2.2 × 10⁴) [1, page 34]). After neglecting the diffusion term, and introducing the Lagrangian time

t (\vec{x}) = \int_{0}^{\vec{x}} ‍ \frac{d {\vec{x}}^{'}}{u ({\vec{x}}^{'})},

(2)

that is, integrating the time that a fluid parcel travel to the spatial position $\vec{x}$ in a Lagrangian point of view, the new equation renders

\frac{\partial n}{\partial t} = {[\frac{\partial n}{\partial t}]}_{nucl} + {[\frac{\partial n}{\partial t}]}_{growth} + {[\frac{\partial n}{\partial t}]}_{coag} .

(3)

This work focuses on solving (3) with the OSMC algorithm. The convection of aerosol in the flow field is not considered here. There are many practically interested conditions, where the flow field can be considered as steady and one dimensional, for example, burner stabilized flames and counter flow diffusion flames. Then it is very easy to integrate (2) to build up a one-to-one mapping between the time and spacial coordinate. Hence the convection can be readily assimilated in the model equation (3). For general cases of aerosol evolution in three dimensional flows, it is possible to use the Lagrangian particle scheme to simulate the aerosol evolution along Lagrangian trajectories [28, 29]. Then (3) still applies along the Lagrangian trajectories. So the model equation (3) and the subsequent OSMC algorithm have very broad application.

The aerosol dynamic processes on the right hand side describe the interactions among molecules (or cluster of molecules) and particles. Nucleation is the process that dozens or hundreds of molecules form a stable critical nucleus (particle (Droplets and particles are treated as synonyms in this article)). It is usually modeled by the classic Becker-Döring theory [30] or its variant self-consistent correction theory [31]. Surface growth includes physical condensation and surface chemical reactions, which involves interactions between gas phase molecules and an aerosol particle. Condensation rate in the free molecule regime (i.e., when aerosol particles are smaller than the mean free path of the gas) is proportional to the surface area of particles, and it is proportional to the diameter in the continuum regime (i.e., when aerosol particles are bigger than the mean free path of the gas) [1]. An approximate interpolation formula for the entire range has been proposed by Fuchs and Sutugin [32]. Surface chemical reactions are too broad to address. A typical example in soot formation in combustion flames is the hydrogen-abstraction/carbon-addition soot growth mechanism [33]. Coagulation is the process that two particles coalesce to form a bigger particle. Here, coagulation may also refer to aggregation. The coagulation dynamics is described by the well-known Smoluchowski's equation

{[\frac{\partial n}{\partial t}]}_{coag} = \frac{1}{2} \int_{0}^{v} ‍ β (v, \tilde{v}) n (\tilde{v}) n (v - \tilde{v}) d \tilde{v} \begin{matrix} - \int_{0}^{\infty} ‍ β (v, \tilde{v}) n (v) n (\tilde{v}) d \tilde{v} . \end{matrix}

(4)

The collision kernel function β(v,u) describes the rate at which particles of size v coagulate with particles of size u. Analytical solutions for the Smoluchowski's equation are known [34] when β has one of the following three forms, (i) β(v,u) = 1, (ii) β(v,u) = v + u, and (iii) β(v,u) = vu. Practically, more complex physical models for β are developed. When particle diameter is smaller than the mean free path of air (the free molecule regime), β can be expressed as [1]

β (v, u) = {(\frac{6}{π})}^{2 / 3} {(\frac{π k_{B} T}{2 ρ})}^{1 / 2} {(\frac{1}{v} + \frac{1}{u})}^{1 / 2} {(v^{1 / 3} + u^{1 / 3})}^{2},

(5)

where v and u are the volume of the (spherical) colliders, T is the temperature, k_B is the Boltzmann's constant, and ρ is the density of particles.

2.2. Operator Splitting

The governing Equation (3) contains very different physical processes. Operator splitting is an efficient and powerful method to solve such evolution equations [35]. In operator splitting method, the time interval [0,t] is descritized as usual, and let δt_i denote one time step. The solution n(t) can be constructed by classic methods, such as the forward Euler method. Other than to integrate all dynamic processes together in one step (denoted as exp⁡(δtX)), operator splitting method separates the integration to multiple steps, such as

\exp (δ t X) = \exp (δ t X_{d}) \exp (δ t X_{s}) + 𝒪 (δ t^{2})

(6a)

= \exp (δ t X_{s}) \exp (δ t X_{d}) + 𝒪 (δ t^{2})

(6b)

= \exp (\frac{1}{2} δ t X_{d}) \exp (δ t X_{s}) \exp (\frac{1}{2} δ t X_{d}) + 𝒪 (δ t^{3})

(6c)

= \exp (\frac{1}{2} δ t X_{s}) \exp (δ t X_{d}) \exp (\frac{1}{2} δ t X_{s}) + 𝒪 (δ t^{3}),

(6d)

where X_d denotes nucleation and surface growth processes, which is solved by deterministic integration method, and X_s denotes coagulation, which is modeled by the Marcus-Lushnikov stochastic process [36]. Equations (6a) and (6b) are first order Lie schemes (the global error at t is 𝒪(δt), while is one order lower than the local error in a step), (6c) and (6d) second order Strang schemes. Figure 1 gives a schematic view of the operator splitting techniques.

Figure 1:

Operator splitting schemes. (a) first order Lie scheme; (b) second order Strang scheme. Inverting the order of X_d and X_s gives the other Lie and Strang schemes, respectively. X_s denotes MC simulation for coagulation, X_d denotes deterministic integration of the other aerosol dynamics.

2.3. MC Simulation of Coagulation

The stochastic algorithm of Gillespie [36] is chosen to simulate coagulation. In essence, the simulation algorithm is as follows:

Step 1. Specify initial values (volume, surface, diameter etc.) for N particles, determine the simulator size V = N/n₀ (where n0 is the initial number density), and evaluate the pair coagulation rate C_ij = β(v_i,v_j)/V, (i = 1,…, N–1, j = i + 1,…, N).

Step 2. Generate a random time τ at which two particles coagulate, where τ satisfies the Poisson distribution

P (τ) = C_{0} \exp^{- C_{0} τ} .

(7)

Here C0 is the summation of all pair coagulation rate, $C_{0} = \sum_{i = 1}^{N - 1} ‍ \sum_{j = 1}^{i + 1} ‍ C_{i j}$ .

Step 3. Randomly choose two particles to coagulate according to the distribution P(i,j) = C_ij/C₀.

Step 4. Update N (one particle less) and C_ij, and repeat Steps 2 and 3, until the accumulated random time τ is larger than the simulation time t_stop.

Gillespie [36] provided three methods to realize Step 3. The partial conditioning method therein is adopted here, which is usually the most efficient and easy to implement for parallel computing. The calculation of C0 has complexity 𝒪(N²) and is the most computationally expensive step in the algorithm. However, there is simple way to update C0 by substracting from its current value the contributions of the two particles that collided and adding the contributions from the newly formed particle. Those operations have complexity 𝒪(N) and have negligible costs compared to the computation of C0. It is worth pointing out that the evolution of C0 can be computed by the majorant kernel method [19, 20], which is numerically more efficient than computing C0 directly. However, in our algorithm C0 is computed directly only at the every beginning, then C0 is updated by the simple and efficient scheme. Hence the majorant kernel method may not outperform our simple updating scheme.

2.4. Deterministic Integration of Nucleation and Surface Growth

Nucleation generates particles with a specific size, which does not depend on the existing aerosol particles. Although there are various nucleation theories for various kinds of aerosol, such as droplets [30, 31], soot particles [37], nucleation rate is generally a function of temperature and nucleated vapor concentration (or polycyclic aromatic hydrocarbon for soot). For a given nucleation model, the nucleation term in (3) can be easily integrated to give the number density of nucleated particles. Surface growth may refer to physical condensation process or surface chemical reactions. Other than on the temperature and surrounding vapor concentration, surface growth rate usually also depends on the size of aerosol particle. In the OSMC algorithm, the size of aerosol particle (modeled by MC particles) by is known. Surface growth term can also be readily integrated.

A self-adaptive fifth order Runge-Kutta method ([38], rkqs subroutine) is used to integrate the dynamic equations for nucleation and surface growth within a discrete time step δt. This choice adds flexibility to control the integration error to a given lever. Numerical experience shows that one-step explicit Euler integration is usually sufficient, that is, X_d(t + δt) = X_d(t) + δt(dX_d/dt).

2.5. Full Algorithm Including All Aerosol Dynamics

Figure 2 shows the flowchart of OSMC including all aerosol dynamics. Only the first order operator splitting (6a) is sketched in the figure.

Figure 2:

Flowchart of the OSMC. Filled blocks denote stochastic steps. Here up-sampling is limited to the doubling method, hence no stochasticity is introduced.

At time t = 0, the simulator is initialized with given parameters, that is, number of MC particles N, simulator size V, and initial PSD. The simulator size is determined as V = N/n₀, where n0 is the initial number density. If simulation starts from an empty case, that is, n₀ = 0, then the simulator is initialized with N particles of the same size for simplicity (the real value of the size is immaterial to the simulation results), and the simulator size is set to a huge number, say V = 10¹⁰ m³, which renders a tiny particle number density to approximate the condition n₀ = 0. If simulation starts from a case with a specific PSD, then the size of the N simulation particles is randomly assigned to satisfy the initial PSD. There are various convenient ways to generate random number to satisfy a given distribution [39], although they are usually not necessary in Monte Carlo simulation of aerosol dynamics, since mostly the simulation starts from an empty case.

In the particles nucleation step, the simulator size V is adjusted to reflect the change of particle number density due to nucleation. Nascent nucleated particles are added to the pool of MC particles. If the total number of MC particles exceeds the maximum allowable value, then down-sampling is performed, that is, exceeding MC particles are randomly removed from the pool to satisfy the number constraint. And then every MC particle undergoes the surface growth process according to its growth rate (usually depends on its size).

The coagulation simulation process (those steps grouped in the dashed bounding box X_s in the flowchart Figure 2) shows how to implement the stochastic algorithm of Gillespie [36]. Updating coagulation kernel is to calculate the pair coagulation rate C_ij = β(v_i,v_j)/V, (i = 1,…, N–1, j = i + 1,…, N). The random coagulation time τ is generated according to (7). The comparison statement τ_sum + τ<δt is to judge whether the time for two particles to coagulate is still admissible within the discrete time step δt. Here τ_sum is the accumulated coagulation time in the coagulation step, which is initialized to zero before a coagulation step begins. If the comparison statement is true, two particles are selected (Step 3 in the above algorithm) to perform coagulation, and the number of MC particles decreases by one. The subsequent up-sampling step is to keep the number of MC particles above a given threshold to avoid severe stochastic error.

The deterministic integration step X_d and the stochastic simulation step X_s can be organized according to the operator splitting techniques discussed in Section 2.2. The effects of different operator splitting will be reported elsewhere. One general conclusion is that the second order Strang scheme (6c) is optimal in considering the computational efficiency and precision.

A few important points should be addressed in the implementation of the algorithm. Within a time step δt, nucleation dumps a number of new particles to the particulate simulation system in order to reflect the change of number density in physical space. The change of number density can be any real value, but the number of newly dumped particles, N_dump, must be integer. There is inevitably a round-off error for the number density nearly 1/N_dump. When N_dump is big, the round-off error is small. However, N_dump is directly proportional to the time step δt and the nucleation rate, which means that N_dump can be very small or even zero under the conditions of every small δt or nucleation rate. To avoid the severe round-off error under these conditions, a correction to N_dump is developed. First, set $N_{dump} = ⌊ I_{nucl} V δ t ⌋$ , where $⌊ \cdot ⌋$ denotes the integer floor function and I_nucl is the nucleation rate, and let $N_{rand} = I_{nucl} V δ t - ⌊ I_{nucl} V δ t ⌋$ (i.e., the decimal part). Next, take a random number, r, from the uniform distribution on the unit interval. If r<N_rand, then N_dump = N_dump + 1; else N_dump is kept unchanged. This stochastic correction aims to produce the right result in the long run or over larger number of repetitions of MC simulation. This correction method has proved very satisfactory in all sorts of simulations investigated here.

Another critical issue is how to control the number of MC particles during the simulation, since nucleation and coagulation may cause the number up or down to unappropriate values. There are a few popular methods in literature [40, 41], all of which are actually different ways of re-sampling the PSD during the simulation. Re-sampling contains up-sampling and down-sampling, either adding or removing particles. After re-sampling, the volume of the system V should also be rescaled accordingly to keep the particle number density unchanged before and after the re-sampling.

Smith and Matsoukas [41] used a constant number scheme for the simulation of coagulation, which dumps a new particle to the MC particle pool after every coagulation (when the particle number is decreased by one). Another popular scheme is the “topping up” method [40], which increases the number by a factor when it reaches the lower limit. A special case of the “topping up” is the doubling method, that is, when the particle number decreases by half, then duplicates all the particles. In this doubling method, no stochastic error is introduced. Maisels et al. [26] used a doubling-halving scheme to simulate nucleation and coagulation, when the particle number decreases to half of N, then double all the particles; when particle number increases to 2N, then remove every other particle. Maisels et al. [26] also extended the constant number scheme [41] to model nucleation and coagulation, and compared the results with those from the doubling-halving scheme. They found that both schemes give comparable results (with regard to accuracy, random fluctuation etc.), except that the fluctuation of particle volume fraction given by the constant number scheme is several orders higher than that by the doubling-halving scheme.

In our simulation, the number of MC particles is set within the range from N/2 to N. When the number decreases to N/2, all the particles are doubled. When the number exceeds N, it is down-sampled to N. This treatment has a few advantages over the constant number and the doubling-halving schemes. Obviously, it is computationally more economic that the doubling-halving scheme since it avoids particles increasing to 2N when computational cost quadruples. On another aspect, it has much less variance than the constant number scheme. In simulations including coagulation and nucleation, it is found that the constant number scheme shows very large variance when coagulation and nucleation rates are comparable. Because under this condition, large portion of the PSD would be randomly removed/duplicated incessantly. However, in the present doubling scheme, no re-sampling is needed when the particle number stays between N/2 and N. Therefore, much less uncertainty is introduced because of less frequent re-sampling.

During re-sampling it requires to randomly select m particles from existing N particles. There are two ways to select the m particles, either repeatable or non-repeatable. Repeatable selection is rather easy to implement as it only needs to generate m random integers all lying between 1 and N. However, in a testing MC simulation it is found that the variance of repeatable re-sampling is four times bigger than the non-repeatable case, under the worst condition m = N/2 (If m>N/2, the inverse selection can be adopted.). Non-repeatable selection requires all m particles are distinct. It is not efficient to repeatedly draw a particle randomly and check whether it has been selected before until m particles are chosen out. Durstenfeld [42] proposed a very smart and efficient algorithm for non-repeatable selection. We re-introduce the algorithm here, since it is very easy to follow. Given an integer sequence from 1 to N, first, randomly pick one number to permutate with the last one. Next, repeat the permutation in the fist step for the beginning N–1 numbers. And next for N–2, and so on for m steps. The last m numbers in the sequence are now randomly selected from the original sequence.

2.6. Parallelization

Computational storage and speed limit the applicability of MC to simulations with relatively small number of particles [43]. Parallel computing can greatly alleviate the storage and speed limits. For the simulation of aerosol dynamics, only coagulation requires to evaluate the interaction between two particles. Hence, it is only necessary to develop parallel scheme for coagulation. Other aerosol dynamics are treated with the same manner as in serial simulation, except that the information about the total particle number density and the simulator volume should be gathered and broadcast to every processor.

Section 2.3 has introduced the partial conditioning method of Gillespie [36], who also introduced other two methods, full conditioning and first-coalescence. But the partial conditioning method is usually the most computationally efficient, and the most convenient for parallel computing.

The Message Passing Interface (MPI) is used in the parallel computing. Every computing processor has N_p particles assigned initially, and there are totally N_CPU processors. The total particle number is N = N_p · N_CPU. Every processor needs to communicate with other N_CPU-1 processors. The total communication channels among the processors are (N_CPU-1)N_CPU/2 (since communication is reciprocal). However, in the current parallel algorithm a processor deals with only N_b = (N_CPU-1)/2 channels, not to communicate with all other N_CPU-1 processors. The communication between any two processors are accomplished through either the direct channel (if has) or indirect multi-channel. A bookkeeping table is created to record the communication structure, that is, sending to and receiving from which processor. Here, N_CPU is required to be an odd number, otherwise the computational load cannot be evenly distributed among the processors. This peculiarity would usually suffer a little waste of computational resource, since computer hardware is always constructed with even number processors.

Every processor has a buffer array of size N_b. The element buffer of the array holds all the information received from a specific processor, which is determined in the bookkeeping table. In evaluating the matrix C_ij, every processor is assigned an equal size block from the matrix. The information for i and j particles is either both available locally, or one is available locally and the other could be find in the buffer sooner or later. The case of two alien particles has been ruled out during the construction of the bookkeeping table. This process is guaranteed to work, since all processors keep sending and receiving information (non-blocking communication) until all processors finish the evaluation of C_ij. After that, the global information (summation of all C_ij elements and the maximum element of C_ij) is gathered and broadcast to every processor.

The key point in this parallel algorithm is the use of the bookkeeping table, which evenly distributes computational load and launches communications simultaneously among processors. Its efficiency is tested in the next section.

3. Results and Discussion

3.1. Parallel Computing Efficiency

Generally, coagulation is the most computationally intensive part in all aerosol dynamics. Besides, only coagulation involving two particles interaction and requires message communication between computing processors in parallel computing. Therefore a testing case of coagulation in the free molecule regime is simulated by the OSMC to investigate the parallel computing efficiency under an eight-core workstation and a computer cluster. The testing case was originally used by Frenklach and Harris [44] to validate their moments method. Initially all particles are mono-dispersed with diameter 1.5 nm, and density 1800 kg·m⁻³. The particle number density is 10¹⁸ m⁻³. The temperature is kept constant 1800 K. The simulation time is t_stop = 4.5 ms.

The coagulation scheme (Section 2.3) requires N²/2 operations to initialize the symmetric matrix C_ij. Every coagulation event only needs 𝒪(N) operations. On the other hand, the average time during two subsequent coagulation is proportional to 1/N, which means that increasing N will cause coagulation events more often within the time interval [0,t_stop]. Hence the total operations should be proportional to N².

Table 1 compares the computational efficiency with different number of MC particles. First, it is found that the computation time is proportional to N² (with relative error less than 10%) when using the same number of computing processors, which agrees with the above analysis. This fact is used to estimate the serial computation time t_serial with large number of particles, which cannot be handled by a single processor, by scaling the computation time for a small serial simulation. The t_serial obtained in this way is then used to calibrate the efficiency of parallel computing. The efficiency is found insensitive to N_p (number of particles per processor) with fixed N_CPU (number of processors) and only marginal decreasing is observed when increasing N_p. While increasing N_CPU, the parallel efficiency decreases more significantly. Overall, the parallel efficiency is satisfactory. Under the workstation, it reaches 87% and 86% for the two parallel simulations, respectively. Under the computer cluster, it also reaches 87% for a half million particles simulation. For a 3.7 million simulation, it is 58%. Such efficiencies are promising.

Table 1:

Parallel simulations of coagulation in the free molecule regime under an eight-core workstation (WS) and a Beowulf class heterogeneous computer cluster (CC).

N_p (× 10³)	N _CPU	N = N_p·N_CPU(× 10³)	Time (s)	Speedup (t_serial/t_parallel)	Efficiency (speedup/N_CPU)

WS
8	1	8	8.6	1.00	1.00
8	7	56	68.7	6.12	0.87
32	1	32	125.8	1.00	1.00
32	7	224	1021.0	6.04	0.86
CC
20	1	20	63.1	1.00	1.00
20	31	620	2241.3	27.06	0.87
20	93	1860	9817.1	55.60	0.60
40	1	40	233.7	1.00	1.00
40	31	1240	8835.0	24.42	0.82
40	93	3720	37441.0	54.00	0.58

3.2. Free Molecule Regime Coagulation

The setup of free molecule regime coagulation testing case is introduced in Section 3.1, where the efficiency of parallel computation is discussed. In this section, the accuracy will be addressed.

There is no general analytical solution available for the free molecule regime coagulation. However, after introducing the self similarity transformation, the Smoluchowski's equation can be reduced to an ordinary integrodifferential equation, and the normalized PSD will be self-similar after enough long time [1]. Specifically, the total particle number density N_∞ is found to satisfy [45]

\frac{d N_{\infty}}{d t} = - \frac{a}{2} {(\frac{3}{4 π})}^{1 / 6} {(\frac{6 k_{B} T}{ρ})}^{1 / 2} ϕ^{1 / 6} N_{\infty}^{11 / 6} = - a A N_{\infty}^{11 / 6},

(8)

where the constant a is an integral function (Friedlander [1, (page 261)] gave a = 6.67 by numerical analysis.) of the similarity transformation introduced, ϕ is the volume fraction, and A is a short notation of the corresponding part. Since both T and ϕ are constant for this case, it is easy to integrate the equation to render

N_{\infty} (t) = {(\frac{6 N_{\infty}^{5 / 6} (0)}{6 + 5 t a A N_{\infty}^{5 / 6} (0)})}^{6 / 5} .

(9)

Figure 3 compares the number density obtained by the MC simulation, a sectional method [46], and two solutions evaluated according to (9) with the constant a = 5.80 and a = 6.67, respectively. The MC simulation and the sectional method give nearly identical solutions. The analytical solution with fitted constant (a = 5.80) also agrees very well with the numerical solutions. But the solution (a = 6.67) given by the numerical analysis of Friedlander [1] differs slightly from the others. It is worth pointing out that a = 6.67 is obtained after the PSD has reached self similarity state, but within the simulation time t_stop = 4.5 ms the PSD is far from self similarity (see discussion below). The simulation shows that a smaller constant a = 5.80 can be used to predict the evolution of the total number density before the PSD reaching the self-similarity state. It is also worth pointing out that the MC simulation uses only 50 particles. This reveals that the evolution of number density is insensitive to the precise form of the PSD, since simulations with N = 50 may not capture the PSD very accurately. This conclusion coincides with that for constant kernel coagulation. Simulations of coagulation in the continuum regime also comes to the same conclusion (results not shown in this paper). Consequently it seems that the MC simulation can always predict the number density very accurately with only a few dozens of particles with enough repetitions to limit the uncertainty.

Figure 3:

Evolution of the number density for coagulation in the free molecule regime. The program code for the sectional method is from Prakash et al. [46], enough number of bins has been used to render the solution presented here. a = 5.80 and a = 6.67 are the constants in the similarity transformation solution (9) corresponding to numerical fitting (current) and numerical analysis in Friedlander [1].

The self-similar solution of the PSD (after reaching the self-similar form) means that if the PSD at one time is known, then the PSD at any other time is also known through the self-similar transformation. Figure 4 shows the PSDs at two different times. Both times are larger than the time required to reach the self-similar form, t_SP, which is estimated as t_SP = 0.02's according to [1, Equation (7.78)]. With N = 50, the simulated PSDs exhibit good self-similar property, and also agree with that from a larger simulation of N = 3200. However, it is clear that the simulation with small N cannot correctly predict the number density of large particles (η>10 in Figure 4), which agrees with the findings in the simulations of constant kernel coagulation.

Figure 4:

Self-similar PSD (histogram) for coagulation in the free molecule regime. The results are obtained by averaging over 10000 repeated MC simulations. ϕ is the volume fraction. η is dimensionless particle volume normalized by the average volume (ϕ/N_∞). ψ is the normalized particle density function.

3.3. Free Molecule Regime Coagulation and Constant Rate Nucleation

This testing case has the same setting as introduced in Section 3.1 for coagulation, and it also includes nucleation with constant rate of 1.0 × 10²⁰ m⁻³·s⁻¹. Nucleated particles are assumed to have the same diameter 1.5 nm as that of the initial particles.

Figure 5 shows the simulation results with different N. Solutions from the sectional method [46] are also given as reference. All these numerical results, that is, number density, volume fraction and second moment and so forth, agree with each other very well. The second moment given by the sectional method is slightly bigger than those from the MC simulations. It is believed that the sectional method has lower accuracy than the MC simulation with N = 2000. Since increasing N further also gives the same solution as N = 2000. But it is hard to increase the accuracy of the sectional method used here to very high level [46].

Figure 5:

Results for free molecule regime coagulation and constant rate nucleation. (a) number density, (b) average diameter, (c) total particle volume, and (d) second moment.

3.4. Condensation and Coagulation

Analytical solutions for condensation and coagulation (Of course, the collision kernels are limited to the three types that the solutions to the Smoluchowski's equation are known.) are available for a few special cases [47]

α(v) = α₀v, β(v,u) = β₀,

α(v) = α₀v, β(v,u) = β₀·(v + u),

α(v) = α₀, β(v,u) = β₀,

where α(v) is the condensation rate, α₀ and β₀ are constants.

In case (iii) condensation and coagulation are independent with each other, which is a trivial case for the current MC method that intrinsically separates the two processes. Here, results for case (i) from MC simulation are reported. For simplicity, the parameters are set as α₀ = 1, β₀ = 1. The initial PSD is

n (v, 0) = \exp (- v),

(10)

which has total number density 1, and average volume 1.

Under these conditions, the solution for case (i) is [47]

n (v, t) = \frac{M_{0} {(t)}^{2}}{M_{1} (t)} \exp (- \frac{M_{0} (t)}{M_{1} (t)} v),

(11)

where the zeroth and first moments, M0 and M1, are

\begin{matrix} M_{0} (t) = \frac{2}{2 + t}, \\ M_{1} (t) = \exp (t) . \end{matrix}

(12)

Further, we derive the second and third moments here

M_{2} (t) = (2 + t) \exp (2 t),

(13)

M_{3} (t) = \frac{3}{2} {(2 + t)}^{2} \exp (3 t) .

(14)

Figure 6 compares the MC simulations (Strang1 scheme) of the condensation and coagulation case with the analytical solutions (12)–(14), which agree extremely well with each other. First order Lie1 and Lie2 schemes are found to give almost the same results as the Strang1 scheme (not shown). In this case, coagulation is independent of condensation, and condensation does not change the number density, hence the solution for the number density is the same as the case of pure constant rate coagulation (not shown here). Since the condensation rate is a linear function of particle volume, the evolution of volume fraction M1 is deterministic, which has increasing rate $α_{0} \sum_{i = 1}^{N} ‍ v_{i}$ . However, if condensation rate has nonlinear dependency on particle volume, this conclusion is generally not true, because the increasing rate may not be proportional to the total volume fraction any more. During the simulation time, coagulation alone increases the average particle volume by four folds (since the number density decreases to 0.25, see (13)). However, condensation makes particles grow at exponential rate. This case is condensation dominated, and different operator splitting schemes produce same results. It is also found that the simulation results are almost not affected by the time step δt chosen for this case.

Figure 6:

Evolution of moments for condensation and coagulation with Strang1 scheme, δt = 0.1. (a) M0, (b) M1, M2 and M3. Strang1 scheme, δt = 0.1.

4. Conclusions

A new hybrid algorithm combining stochastic simulation and deterministic method for solving the PBE for aerosol dynamics is proposed. Nucleation and surface growth are handled with deterministic means, while coagulation is simulated with a stochastic method (Marcus-Lushnikov stochastic process). Operator splitting techniques are used to synthesize the deterministic and stochastic parts in the algorithm. The algorithm is parallelized using the MPI.

Parallel simulation of aerosol coagulation in the free molecule regime shows that the parallel computing efficiency decreases with increasing the number of CPUs. While for fixed number of CPUs, increasing the number of MC particles does not affect the efficiency significantly. Near 60% parallel efficiency is achieved for the maximum testing case with 3.7 million MC particles running on 93 CPUs in parallel.

For verification purpose, the algorithm has been applied to a series of testing cases, that is, free molecule regime coagulation, free molecule regime coagulation + constant rate nucleation, and constant kernel coagulation and condensation. All the simulation results agree very well with known analytical solutions or numerical solutions from a sectional method [46]. Generally, it is found that only small number (hundreds or thousands) of MC particles is necessary to accurately predict the aerosol particle number density, volume fraction, and so forth, that is, low order moments of the PSD function. Accurately predicting the high order moments of the PSD needs to dramatically increase the number MC particles.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Footnotes

Acknowledgment

Kun Zhou would like to thank F. Bisetti for insightful discussion and valuable help. Kun Zhou is partially supported by the National Natural Science Foundation of China (Grant no. 11302110). Zhu He is supported by the National Key Technology R&D Program of Chian (Grant no. 2011BAK06B02).

References

Friedlander

, Smoke, Dust, and Haze: Fundamentals of Aerosol Dynamics, Oxford University Press, New York, NY, USA, 2nd edition, 2000.

Ramkrishna

, Population Balances: Theory and Applications to Particulate Systems in Engineering, Academic Press, San Diego, Calif, USA, 2000.

Kiparissides

, “Polymerization reactor modeling: a review of recent developments and future directions,” Chemical Engineering Science, vol. 51, no. 10, pp. 1637–1659, 1996.

Peterson

T. W.

Gelbard

, and Seinfeld

J. H.

, “Dynamics of source-reinforced, coagulating, and condensing aerosols,” Journal of Colloid and Interface Science, vol. 63, no. 3, pp. 426–445, 1978.

Gelbard

and Seinfeld

J. H.

, “The general dynamic equation for aerosols. Theory and application to aerosol formation and growth,” Journal of Colloid and Interface Science, vol. 68, no. 2, pp. 363–382, 1979.

Gelbard

Tambour

, and Seinfeld

J. H.

, “Sectional representations for simulating aerosol dynamics,” Journal of Colloid and Interface Science, vol. 76, no. 2, pp. 541–556, 1980.

Kim

Y. P.

and Seinfeld

J. H.

, “Simulation of multicomponent aerosol dynamics,” Journal of Colloid and Interface Science, vol. 149, no. 2, pp. 425–449, 1992.

Pratsinis

S. E.

, “Simultaneous nucleation, condensation, and coagulation in aerosol reactors,” Journal of Colloid and Interface Science, vol. 124, no. 2, pp. 416–427, 1988.

McGraw

, “Description of aerosol dynamics by the quadrature method of moments,” Aerosol Science and Technology, vol. 27, no. 2, pp. 255–265, 1997.

10.

Marchisio

D. L.

and Fox

R. O.

, “Solution of population balance equations using the direct quadrature method of moments,” Journal of Aerosol Science, vol. 36, no. 1, pp. 43–73, 2005.

11.

Frenklach

, “Method of moments with interpolative closure,” Chemical Engineering Science, vol. 57, no. 12, pp. 2229–2239, 2002.

12.

M. Z.

Lin

J. Z.

, and Chan

T. L.

, “A new moment method for solving the coagulation equation for particles in Brownian motion,” Aerosol Science and Technology, vol. 42, no. 9, pp. 705–713, 2008.

13.

M. Z.

and Lin

J. Z.

, “Taylor-expansion moment method for agglomerate coagulation due to Brownian motion in the entire size regime,” Journal of Aerosol Science, vol. 40, no. 6, pp. 549–562, 2009.

14.

M. Z.

Lin

J. Z.

Jin

H. H.

, and Jiang

, “The verification of the Taylor-expansion moment method for the nanoparticle coagulation in the entire size regime due to Brownian motion,” Journal of Nanoparticle Research, vol. 13, no. 5, pp. 2007–2020, 2011.

15.

Marcus

A. H.

, “Stochastic coalescence,” Technometrics, vol. 10, no. 1, pp. 133–148, 1968.

16.

Gillespie

, “The stochastic coalescence model for cloud droplet growth,” Journal of the Atmospheric Sciences, vol. 29, pp. 1496–1510, 1972.

17.

Lushnikov

A. A.

, “Certain new aspects of the coagulation theory,” Izvestiya, Atmospheric and Oceanic Physics, vol. 14, pp. 738–743, 1978.

18.

Eibeck

and Wagner

, “Stochastic particle approximations for Smoluchoskis coagulation equation,” Annals of Applied Probability, vol. 11, pp. 1137–1165, 2001.

19.

Eibeck

and Wagner

, “An efficient stochastic algorithm for studying coagulation dynamics and gelation phenomena,” SIAM Journal on Scientific Computing, vol. 22, no. 3, pp. 802–821, 2001.

20.

Goodson

and Kraft

, “An efficienct stochastic algorithm for simulating nano-particle dynamics,” Journal of Computational Physics, vol. 183, no. 1, pp. 210–232, 2002.

21.

Gillespie

D. T.

, “Approximate accelerated stochastic simulation of chemically reacting systems,” Journal of Chemical Physics, vol. 115, no. 4, pp. 1716–1733, 2001.

22.

Gillespie

D. T.

, “Stochastic simulation of chemical kinetics,” Annual Review of Physical Chemistry, vol. 58, pp. 35–55, 2007.

23.

Rajamani

Pate

W. T.

, and Kinneberg

D. J.

, “Time-driven and event-driven Monte Carlo simulations of liquid-liquid dispersions: a comparison,” Industrial and Engineering Chemistry Fundamentals, vol. 25, no. 4, pp. 746–752, 1986.

24.

Debry

Sportisse

, and Jourdain

, “A stochastic approach for the numerical simulation of the general dynamics equation for aerosols,” Journal of Computational Physics, vol. 184, no. 2, pp. 649–669, 2003.

25.

Zhao

and Zheng

, “A new event-driven constant-volume method for solution of the time evolution of particle size distribution,” Journal of Computational Physics, vol. 228, no. 5, pp. 1412–1428, 2009.

26.

Maisels

Einar Kruis

, and Fissan

, “Direct simulation Monte Carlo for simultaneous nucleation, coagulation, and surface growth in dispersed systems,” Chemical Engineering Science, vol. 59, no. 11, pp. 2231–2239, 2004.

27.

Patterson

R. I. A.

Singh

Balthasar

Kraft

, and Norris

J. R.

, “The linear process deferment algorithm: a new technique for solving population balance equations,” SIAM Journal on Scientific Computing, vol. 28, no. 1, pp. 303–320, 2006.

28.

Zhou

Attili

, and Bisetti

, “Direct numerical simulation of aerosol growth processes in a turbulent mixing layer,” in Proceedings of the 7th International Conference on Computational Fluid Dynamics (ICCFD '12), Big Island, Hawaii, 2012.

29.

Zhou

, “Monte Carlo simulation for soot dynamics,” Thermal Science, vol. 16, pp. 1491–1494, 2012.

30.

Becker

and Döring

, “Kinetische Behandlung der Keimbildung in übersättigten Dämpfern,” Annals of Physics, vol. 24, pp. 719–752, 1935.

31.

Girshick

S. L.

and Chiu

C.-P.

, “Kinetic nucleation theory: a new expression for the rate of homogeneous nucleation from an ideal supersaturated vapor,” Journal of Chemical Physics, vol. 93, no. 2, pp. 1273–1277, 1990.

32.

Fuchs

N. A.

and Sutugin

A. G.

, “High-dispersed aerosols,” in Topics in Current Aerosol Research, Hidy

G. M.

and Brock

J. R.

, Eds., Pergamon, New York, NY, USA, 1971.

33.

Frenklach

and Wang

, “Detailed modeling of soot particle nucleation and growth,” Proceedings of the Combustion Institute, vol. 23, no. 1, pp. 1559–1566, 1991.

34.

Aldous

D. J.

, “Deterministic and stochastic models for coalescence (aggregation and coagulation): a review of the mean-field theory for probabilists,” Bernoulli, vol. 5, no. 1, pp. 3–48, 1999.

35.

McLachlan

R. I.

and Quispel

G. R. W.

, “Splitting methods,” Acta Numerica, vol. 11, pp. 341–434, 2002.

36.

Gillespie

D. T.

, “An exact method for numerically simulating the stochastic coalescence process in a cloud,” Journal of the Atmospheric Sciences, vol. 32, no. 10, pp. 1977–1989, 1975.

37.

Appel

Bockhorn

, and Frenklach

, “Kinetic modeling of soot formation with detailed chemistry and physics: laminar premixed flames of C2 hydrocarbons,” Combustion and Flame, vol. 121, no. 1–2, pp. 122–136, 2000.

38.

Press

W. H.

Flannery

B. P.

Teukolsky

S. A.

, and Vetterling

W. T.

, Numerical Recipes in Fortran 77: The Art of Scientific Computing, Cambridge University Press, New York, NY, USA, 2nd edition, 1992.

39.

Fishman

G. S.

, Monte Carlo: Concepts, Algorithms, and Applications, Springer, New York, NY, USA, 1996.

40.

Liffman

, “A direct simulation Monte-Carlo method for cluster coagulation,” Journal of Computational Physics, vol. 100, no. 1, pp. 116–127, 1992.

41.

Smith

and Matsoukas

, “Constant-number Monte Carlo simulation of population balances,” Chemical Engineering Science, vol. 53, no. 9, pp. 1777–1786, 1998.

42.

Durstenfeld

, “Algorithm 235: random permutation,” Communications of the ACM, vol. 7, p. 420, 1964.

43.

Laurenzi

I. J.

Bartels

J. D.

, and Diamond

S. L.

, “A general algorithm for exact simulation of multicomponent aggregation processes,” Journal of Computational Physics, vol. 177, no. 2, pp. 418–449, 2002.

44.

Frenklach

and Harris

S. J.

, “Aerosol dynamics modeling using the method of moments,” Journal of Colloid and Interface Science, vol. 118, no. 1, pp. 252–261, 1987.

45.

Lai

F. S.

Friedlander

S. K.

Pich

, and Hidy

G. M.

, “The self-preserving particle size distribution for Brownian coagulation in the free-molecule regime,” Journal of Colloid and Interface Science, vol. 39, no. 2, pp. 395–405, 1972.

46.

Prakash

Bapat

A. P.

, and Zachariah

M. R.

, “A simple numerical algorithm and software for solution of nucleation, surface growth, and coagulation problems,” Aerosol Science and Technology, vol. 37, no. 11, pp. 892–898, 2003.

47.

Ramabhadran

T. E.

Peterson

T. W.

, and Seinfeld

J. H.

, “Dynamics of aerosol coagulation and condensation,” AIChE Journal, vol. 22, no. 5, pp. 840–851, 1976.