A Fast Monte Carlo Method Based on an Acceptance-Rejection Scheme for Particle Coagulation

A fast acceptance-rejection scheme that can boost the performance of Monte Carlo methods for particle coagulation has been designed and validated in this article. The central idea behind this scheme is to perform a fast evaluation of the maximum coagulation rate by making full use of the information of those particle pairs that are already generated to determine the desired coagulation pair. To this end, numerical experiments are carried out to establish a connection between the average coagulation rate obtained from information of those already generated particle pairs and the maximum coagulation rate. The correctness of the proposed scheme is then validated by comparing it to some established methods, adopting a physically realistic Brownian coagulation kernel. The computational efficiency of the fast AR scheme is finally measured.


INTRODUCTION
The dynamics of dispersed particle systems is of interest in many areas of science and engineering including combustion, aerosol reactor engineering, atmospheric aerosols and chemical engineering (Friedlander, 1997;Ramkrishna, 2000).The evolution of a dispersed phase is governed by the population balance equation, which can account for all processes that generate, modify and remove particles from the population, such as coagulation, breakage, nucleation, condensation/evaporation.Among those mechanisms, coagulation is of primary importance with regard to the understanding of the behavior, handling and treatment of particle systems.The particle coagulation process is described by the general dynamic equation (GDE) (Friedlander, 1997) where n(v,t) is the particle size distribution at time t and β (u,v) is the coagulation rate for two particles with volume u and v.The first term on the right-hand of above equation is accounting for the formation of a particle with volume v due to the coagulation event between a particle of volume u and a particle of volume (v -u); the coefficient 1/2 is introduced since collisions are counted twice in the integral.The second term indicates the disappearance of a particle with volume v as a result of collisions with particles of other sizes.
Different methods are available at this point to solve above integro-differential equation, including sectional methods (Landgrebe and Pratsinis, 1990;Wu and Biswas, 1998;Jeong and Choi, 2001;Mitrakos et al., 2007;Yamamoto, 2012), methods of moment (Park and Lee, 2001;Terry et al., 2001;Yamamoto, 2004;Yu et al., 2008) and Monte Carlo methods (Garcia et al., 1987;Liffman, 1990;Simith and Matsoukas, 1998;Kruis et al., 2000;Lee and Matsoukas, 2000;Efendiev and Zachariah, 2002;Maisels et al., 2004;Zhao and Kruis, 2008).All these methods have different precision, algorithm complexity and computational efficiency.To be more specific, sectional methods take moderate computational time, yet their sectional representations often result in fairly complicated algorithms.Methods of moment are computationally efficient, however usually require knowledge about the shape of particle size distribution which may experience dramatic changes during the simulation, especially when another mechanism occurs simultaneously such as nucleation.Unlike those deterministic schemes, Monte Carlo methods solve Eq. ( 1) by simulating it as a stochastic process and estimating the interesting characteristics based on the 'data' from the computer experiment.The advantages of Monte Carlo methods include: • The discrete character embedded in Monte Carlo methods makes them the ideal choice for those inherently discrete processes, such as coagulation.• The methods do not require any a prior information of the particle size distribution.• Each individual particle is treated as a single object, which can have several properties such as diameter, charge, composition, morphology, etc. • The numerical algorithms are generic simple, robust and easy to code.A major limitation to Monte Carlo methods has been recognized as their heavy computing cost.This is because Monte Carlo methods are essentially a type of stochastic simulation method with its computational accuracy being proportional to the square root of the number of samples involved in simulation.That is, in order to increase the accuracy of Monte Carlo calculation by a factor of 10, the number of samples should be increased by a factor of 100.In addition, the limitations of current CPU speed and memory capacity have also posed a considerable obstacle to their applications.For example, Monte Carlo methods for particle dynamics are presently limited to some 10 5 -10 7 simulation particles in a single simulation.In this situation one has to take some effective measures to improve the computational performance of Monte Carlo algorithms, which may include using CPU-based high performance cluster (HPC) (Kruis et al., 2012), employing graphic processing units (GPU), designing better Monte Carlo algorithms with higher computing efficiency or jointly using the previous methods.
Two most widely used schemes by Monte Carlo methods to handle particle coagulations are the so-called Inverse (Nanbu, 1980) and acceptance-rejection (AR) scheme (Garcia et al., 1987;Efendiev and Zachariah, 2002).The fundamental difference between the two schemes lies in the different way they use to randomly choose the particle pair to be coagulated.Inverse method chooses a desirable particle pair by generating a random number between 0 and 1, and comparing it in turn with a normalized cumulative sum (prefix sum) of a sequence of coagulation rates.Cumulative sums are trivial to compute in a sequential way, but at the expense of considerable computer time as coagulation rate associated with a particle might change after each coagulation event.However, with Inverse scheme one can always find a particle pair after finite-times (up to the number of simulation particles) of attempts.Alternatively, AR scheme uses another simpler, more straightforward criterion to choose the desirable particles: it randomly picks a pair of particles, computes then the corresponding coagulation probability (based on the selected particle pair) and checks to see whether or not to accept it.AR scheme offers manyfold advantages in practical applications: first of all, particle selection depends only upon a ratio of the coagulation rate, calculated from the selected particle pair, to the maximum coagulation rate currently, which leads to a much easier implementation of Monte Carlo code.Secondly, the criterion for selecting particles is not only simply but sufficiently flexible; e.g., Monte Carlo method based on AR scheme still works even the foregoing maximum coagulation rate has been replaced with some larger values (of course, the probability of acceptance might drop).This feature essentially allows AR scheme to bypass the direct calculation of the 'true' maximum coagulation rate (in general computer time consuming).Thirdly, Monte Carlo method using AR scheme features an inherent parallelism.This is because choosing arbitrary particle pairs (to be coagulated) with AR scheme is uncorrelated, and hence can be done in parallel.This inherent feature about AR scheme allows Monte Carlo method that uses it to be easily implemented on some parallel architectures such as on Graphic Processing Units (GPU) (NVIDIA, 2008), in which a large number of lightweighted threads that can be executed simultaneously are available.
Perhaps a major difficulty that hampers Monte Carlo methods use AR scheme is the evaluation of the maximum coagulation rate.The reasons come from several sources: (1) The maximum coagulation rate can hardly be obtained explicitly due to the rather complexity of the coagulation kernel, i.e., an analytical evaluation is almost impossible; (2) Obtaining the maximum coagulation rate numerically is also computationally expensive because of the (simulation particle) number square dependence; (3) The situation might get worse as the evaluation of the maximum coagulation rate must be done repeatedly as particle size distribution keeps changing in the course of coagulation.In order to overcome this problem so as to improve the computing efficiency of Monte Carlo methods, a fast AR scheme has been designed and validated in this work.The basic idea used by the new AR strategy is to take advantage of the information of those particle pairs being generated for finding the desirable particle pair, and in such a way get around the direct calculation of the maximum coagulation rate.
The synopsis of this paper is outlined as follows.Section 2 contains a brief review of the Monte Carlo method adopted in the present work as well as an introduction to the general procedure of Monte Carlo method using AR scheme.Section 3 formulates in detail the idea behind the fast AR scheme.A validation study and an analysis of the efficiency improvement for the proposed AR scheme are presented in Section 4. Section 5 concludes the paper.

Differentially-weighted Monte Carlo method
A differentially weighted Monte Carlo method (Zhao and Kruis, 2008) has been chosen in the present work.This Monte Carlo scheme introduces a weight property for each simulation particle.Mathematically, the weight value associated with a simulation particle indicates the number of real particles being represented by the latter.Using weight makes it possible to simulate a practical system with a large number of particles with fairly limited simulation particles, say 10 3 -10 5 particles.The particle weight as well as the particle diameter are adjusted after each coagulation event to keep track of the particle size distribution over the full particle size spectrum and preserve the simulation history of each particle.More precisely, differentially weighted Monte Carlo method describes particle coagulation using the following modified kernel 2 max( , ) '( , ) ( , ) where β(i,j) is the collision kernel between particle i and j (β(i,j) is equal to that β(u,v) in Eq. ( 1)); w i and w j stand for the weight in associated with particle i and j, respectively.Note that coagulation kernel β(i,j) is symmetric, i.e., β(i,j) = β(j,i); but β'(i,j) is not as w i is in general not equal to w j .

Acceptance-rejection Scheme for Selecting Particle pairs
The cornerstone of Monte Carlo methods when dealing with coagulation problems is how to randomly but reasonably choose desirable coagulation pairs.AR scheme uses a quite straightforward but effective criterion to do this work, namely , max '( , ) ' where r is a random number uniformly distributed between 0 and 1; p i,j is the coagulation probability calculated from as a ratio of the coagulation rate, β'(i,j) (based on the chosen particle pair i and j), to the maximum coagulation rate, β' max .
The general process of selecting coagulation pair by AR scheme is quite simple.A particle pair is randomly chosen and the corresponding coagulation probability based on the selected particles is calculated to compare with a random number, as shown in Eq. (3).The selected pair is accepted if the inequality in above equation is satisfied; otherwise it is rejected.For the latter case a new particle pair are chosen at random and the same judgment is made.This selection process will be repeated until a coagulation pair is found.After determining a coagulation pair, the average waiting time associated with that event, τ , can be computed as where n s is the number of simulation particles.The complexity of computing τ scales as n s 2 , which is relatively high for large n s .In this situation one may reduce above equation to ' , 2 ( -1) Put it in another way, the denominator of Eq. ( 4) has been replaced with only one particle pair -the selected coagulation pair (k,l) (Garcia et al., 1987).As long as τ is available, simulation time, t, can be updated via and in this manner particle coagulation is proceeding.
A cumbersome problem that is faced by Monte Carlo method using AR scheme is the calculation of β' max due to the reasons described in the previous section.A simple, intuitive solution to this problem is to use a constant value, max C β , to replace β' max so as to avoid the later calculation (Garcia et al., 1987).Note, the prescribed max C β is in general over-estimated in order to avoid the danger that the numerator in Eq. (3) (i.e., the coagulation rate β'(i,j)) becomes larger than the denominator (i.e., β' max ).The main flaw of this scheme is that the calculated probability is decreasing dramatically along with the coagulation as the denominator is held fixed, hence degenerating the overall computing efficiency.Thus, this scheme is although feasible theoretically, but it is inefficient in practical applications.
An alternative strategy is to use the so-called bookkeeping technique (Kruis et al., 2000), which is conceptually fairly straightforward.Since after the occurrence of each coagulation event only one pair of particles (the selected pair) has changed their sizes, the new coagulation rate between arbitrary two particles, but involve at least one particle from the selected pair, will be altered; otherwise it remains unchanged.The bookkeeping technique just uses this fact to skip the calculations with respect to those unchanged coagulation rates, and only update those changed parts.In this way a considerable computer time can be saved.The bookkeeping technique can be employed by Monte Carlo methods which use either Inverse or AR scheme.

THE FAST ACCEPTANCE-REJECTION SCHEME
This article presents a fast AR scheme that can estimate β' max very efficiently, hence improving the computational performance of Monte Carlo method.To illustrate the fundamental idea of the scheme, one first notice a basic fact when using Monte Carlo method to handle particle coagulation: many times of attempts are often needed in order to find a desirable coagulation pair.Suppose that there are n r attempts being made before successfully selecting a coagulation pair.This means that n r particle pairs are already generated.Denote the average coagulation rate of those n r particle pairs by ˆsum β , then The idea behind the fast AR scheme is just to estimate an effective β' max from ˆsum β .For this purpose, a relation between β' max and ˆsum β should be established in the first place max ' where γ is a time-variant coefficient.Clearly, if γ is in prior known, then one may quickly evaluate β' max as now ˆsum β can be trivially obtained without much computational efforts.Unfortunately, how γ behaves with simulation time can hardly be obtained analytically, and one has to resort to some numerical means to determine it.For this purpose, one should first choose a coagulation kernel, e.g., a realistic Brownian kernel in the free-molecule regime (Kodas, 1999) where β fr (i,j) is the Brownian kernel; d i and d j are diameters with respect to particle i and j, respectively; K fr is the coagulation coefficient in the free-molecule regime, K fr = (3κ B T/ρ p ) 1/2 with T being the temperature, κ B the Boltzmann constant and ρ p the particle density.
For the free-molecule regime, the resulting size distributions are self-preserving and close to log-normal (Lai et al., 1972).Hence, unimodal log-normal distributions are generated where σ g is the geometric standard deviation and d g the geometric mean diameter.However, other mechanisms such as nucleation or transport can introduce new particles in the simulation volume.This can be simulated by considering a bimodal particle size distribution with the indices 1 and 2 referring to the two modes and N 2 /N 1 as the ratio of the number concentrations.By varying σ g , various unimodal particle size distributions can be simulated.Furthermore, by discretizing the specific particle size distribution, i.e., dividing the size distribution into limited size bins, β' max can be obtained by considering all possible combinations of different size bins; ˆsum β is calculated by picking randomly a large number of different particle pairs, each time from two different size bins, and testing with Eq. (3).At this point γ can be inferred using Eq.(7b) with ˆsum β and ˆsum β being known.Fig. 1 and Fig. 2 show how γ behaves as σ g changes for a unimodal and bimodal particle size distributions, respectively.Basically, γ increases along with σ g for the case of unimodal size distribution, and seems to be independent of d g .More precisely, γ grows from about 1 with σ g ≈ 1, to more than 40 with σ g ≈ 2.0.Larger values for σ g do not occur in freemolecular coagulation.The bimodal case is more complicated because γ is now dependent upon more parameters, like the ratio N 2 /N 1 , as shown in Fig. 2. As can be seen, N 2 /N 1 is not strongly influencing γ, but the broader size distribution (i.e., the bigger difference between 1 g d and 2 g d ) is leading to a pronounced larger γ.It can be seen from these figures that γ will not exceed a certain value for a given distribution.For instance, γ ≤ 40 is for the unimodal distribution with σ g = 2, and γ ≤ 33 for the bimodal distribution with σ g1 = σ g2 = 1.5.For simplicity, one can take a worst-case estimation of γ, and keep it constant throughout the simulation.For example, γ = 40 is adequate to deal with the coagulation in the interesting size regimes, but using some larger γ (γ > 40) is also feasible.As a matter of fact, γ in-between 40 and 200 has been chosen for the later simulation.
If rejected, return to step 2 and repeat step 2 to 4. 7. Update the diameter and weight of the chosen particles.8. Update the simulation time and check to see whether it becomes larger than some prescribed total computing time.If so, terminate the simulation; otherwise return to step 2 and repeat step 2 to 6.

RESULTS AND DISCUSSION
The proposed fast AR algorithm has been validated in this section by comparing its calculation results against benchmarks results from some established methods, including a discrete sectional method (deterministic) and a Monte Carlo method (non-deterministic).The sectional method (Lu, 1994) uses 200 sections and a sectional spacing of 1.12, whereas the Monte Carlo method utilizes Inverse scheme (Zhao and Kruis, 2008).A physically realistic Brownian kernel has been selected as the coagulation kernel.The validation simulations are carried out for particles in a wide size range, namely in the free-molecule and continuum regime.Brownian coagulation kernel for particles in the free-molecule regime is already given in Eq. ( 8); in the continuum regime it takes the following form (Jacobson, 2005) ( , ) ( ) where β c (d j , d j ) stands for the kernel, K c is the coagulation coefficient in the continuum regime, K c = 2κ B T/(3μ) with μ being the dynamic viscosity of the gas.Here μ is temperature-dependent and can be evaluated, according to Sutherland-law (White, 1991), as where T 0 is the reference temperature, μ 0 is the reference viscosity at T 0 , S is Sutherland's constant.C i (C j ) in Eq. ( 12) is the Cunningham slip correction factor (Kasten, 1968), which is a function of particle size d i (d j ) 1 2.498 0.84 exp 0.435 where λ is the mean free path of air molecules, amounting to λ 0 = 0.066 μm at T 0 = 293 K and atmospheric pressure p 0 .At temperature T and pressure p, λ should be corrected, using Eq. ( 13), by The simulation parameters are tabulated in Table 1.Particles are initially mono-disperse with a geometric standard deviation σ g = 1 and number concentration N 0 = 10 17 /m 3 .The simulation is carried out over a long period of time, ranging from 10 -3 τ to 10 3 τ with τ being the characteristic time of the coagulation.Note, τ takes on different expressions, namely τ fr and τ c (Kodas, 1999), and 0 1 2 depending upon particles in the free-molecule or continuum regime, respectively.V t in Eq. ( 16) is the total aerosol volume (m 3 /m 3 ) and can be calculated as V t = N 0 •π 3 0 d /6.As an example, τ fr is 0.0518 s for an aerosol containing 10 17 particles of 3 nm diameter, coagulating at 300 K and atmospheric pressure.All simulation calculations are performed on a desktop system equipped with an Intel Core ™ 2 Quad 2.66 GHz Processor, 8 GB RAM.
The calculated geometric mean diameter, number concentration and geometric standard deviation for particles in the free-molecule regime, averaged over a 100 simulation runs (n run = 100), are compared in Figs.[3][4][5]respectively; show the similar comparisons but for particles in the continuum regime.Clearly, the obtained results agree with the benchmark results very well, and clearly reflect the governing features of the particle coagulation behavior.As can be seen, particle size has witnessed a monotonically increase, while particle number concentration has undergone a monotonically decrease.For example, particles in the continuum regime have grown from 100 nm at t = 0 to about 1 μm at t = 1000 τ c , and their number concentration has decreased by at least three orders in magnitude within the same period of time.In addition, the   geometric standard deviation in the free-molecule regime, σ g , increases from a mono-disperse size distribution at the early phase of the coagulation, σ g = 1, to a self-preserving size distribution, σ g = 1.46, as shown in Fig. 5. Statistical errors are also measured to further reveal the quantitative deviation between the fast AR and the benchmark.At this point only the sectional method is chosen as the benchmark because it is a deterministic method and tends to yield a fairly accurate solution.The mean statistical error, ε mean , is defined in Eq. ( 18) using d g as an example, sec sec in which sec g d is the particle diameter from the sectional method, and d g calculated from fast AR is averaged over n run simulation runs.ε mean is eventually evaluated over n t different time points between 10 -3 τ and 10 3 τ.Figs. 9 and 10 show ε mean with respect to d g and σ g under different γ(40 -200) for particles in the free-molecule and continuum regimes, respectively.As can be seen from the figures, ε mean is rather low, being less than 0.15% in the free-molecule regime; whereas in the continuum regime ε mean is somewhat larger.This can be partially explained by the fact that the enlargement of particle size interval in the continuum regime (due to coagulation) is much larger than that in the freemolecule regime (see Fig. 3 and Fig. 6), and thus requires much more simulation particles to extract the useful information.
Fig. 11 compares the computer time used by Monte Carlo methods that use fast AR and Inverse scheme under different simulation particles, when simulating in the freemolecule regime within t = 1000 τ fr .Note that bookkeeping technique has been adopted by Monte Carlo method using Inverse scheme.It is quite clear from the figure that fast AR scheme performs much better than its counterpart; e.g., around 1s is needed by fast AR scheme, which is ten times   smaller than that by Inverse scheme, in order to complete the calculation within t = 1000 τ fr in the presence of 1000 simulation particles.Furthermore, some insight into the relation between the computer time (wall-clock time) and the number of simulation particles has been pictured in Fig. 12, which shows computer time, t c , scales with the particle number as for the fast AR scheme.Evidently, the computing complexity of fast AR follows basically a n s dependence upon the number of particles, much better than that of the 2 s n dependence of the Inverse scheme.This also explains why Monte Carlo method using fast AR has a much better performance.The remarkable efficiency of the fast AR scheme makes it particularly suitable for simulating particle systems with huge number (Tolke and Krafczyk, 2008).an acceptance-rejection scheme is a non-trivial task, due to the demand of a very intensive evaluation of the maximum coagulation rate.In this article a fast, effective acceptancerejection scheme that can effectively be used by Monte Carlo methods has been developed and validated.The central idea of this scheme is to perform a fast evaluation of the maximum coagulation rate by taking advantage of the information of those particle pairs that are already generated in the course of choosing the coagulation particle pair, and in this way avoid the direct, time-consuming calculation.For this purpose, numerical experiments are carried out to establish a connection between the information of those generated particle pairs and the maximum coagulation rate.Furthermore, the proposed fast AR has been justified by comparing it to some established methods using a physically realistic Brownian coagulation kernel.The high computational efficiency exhibited by the AR scheme proposed would make it the method of choice while employing Monte Carlo methods to real-world simulations, e.g., when coupling Monte Carlo methods with computational fluid dynamics (CFD).

Fig. 9 .
Fig. 9. ε mean with respect to d g and σ g for particles in the freemolecule regime.n run = 100, n t = 16.

Fig. 10 .
Fig. 10.ε mean with respect to d g and σ g for particles in the continuum regime.n run = 100, n t = 16.

Fig. 11 .
Fig. 11.Comparison of computer time between Inverse (solid line) and fast AR method (dashed line) for different number of simulation particles.

Fig. 12 .
Fig. 12. Computer time by AR and Inverse as a function of number of simulation particles within simulation time up to 1000 τ fr .n s = 1000.

Table 1 .
Parameters used in the simulation.