A Majorant Kernel-Based Monte Carlo Method for Particle Population Balance Modeling

A computationally efficient Monte Carlo method using majorant kernel scheme for particle coagulation is presented in this note. The key is to use majorant kernel in the evaluation of the maximum coagulation rate, which is computationally time consuming, in hope of reducing the computational cost. The proposed scheme is verified by means of a deterministic sectional method for Brownian collision kernel. The computational efficiency of the scheme proposed has also been measured by comparing with sectional method and Monte Carlo methods using other schemes


INTRODUCTION
The growth and disintegration of dispersed particle systems are of interests in several areas of scientific and engineering including combustion, aerosols, chemical engineering (Friedlander, 1997;Ramkrishna, 2000), etc.The evolution of a dispersed phase can be governed by the population balance equation, which accounts for all the processes, such as coagulation, breakage, nucleation, condensation/ evaporation, etc., that generate and remove particles from the population.Among those mechanisms, coagulation is one of the most common particle growth mechanisms and is often decisive for the evolution of particle size distribution.Mathematically, the particle coagulation process can be governed by the following population balance equation (PBE) where n(v p , t) is the particle size distribution (PSD) at time t, β(u, v) is the coagulation rate for two particles with volume u and v.The first term on the right-hand side describes the rate of production of particles of 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 to avoid counting collisions twice in the integral.The second term indicates the rate of disappearance of particles of volume v as a result of collisions with any other particle.Generally speaking, there are two classes methods suitable for addressing above PBE: (1) The deterministic methods such as sectional methods (Landgrebe and Pratsinis, 1990;Wu and Biswas, 1998;Jeong and Choi, 2001;Mitrakos et al., 2007), methods of moment (Park et al., 2001;Terry et al., 2001;Yamamoto, 2004;Yu et al., 2008;Yamamoto, 2012); (2) The non-deterministic methods such as Monte Carlo (MC) methods (Garcia et al., 1987;Liffman, 1992;Smith and Matsoukas, 1998;Kruis et al., 2000;Efendiev and Zachariah, 2002;Maisels et al., 2004;Zhao et al., 2009).Among all these methods, MC methods are particularly suitable and effective due to several reasons: (1) MC methods are essentially a statistical method, making them suitable for describing the inherently discrete processes such as coagulation; (2) MC methods exhibit great flexibility and scalability when new characteristic of particles needs to be incorporated, because they treat each simulation particle as an individual object having some properties such as sizes, positions, charges, etc.; (3) In compare with those deterministic methods, MC algorithms are generic simple and easy to code.Of course, MC methods are also not perfect.A major disadvantage that hampers most MC methods is perhaps their heavy computational cost.It is well known that the computational accuracy of MC methods is adversely proportional to the square root of the number of samples involved; put it another way, in order to increase the accuracy of MC calculation by a factor of 10, the number of samples should be increased by a factor of 100.
Effectively choosing collision partners is at the heart of most MCs while dealing with particle coagulation.Two widely used strategies by various MCs to date are what is called the Inverse scheme (Nanbu, 1980;Kruis et al., 2012;Wei and Kruis, 2013a) and acceptance-rejection (AR) scheme (Garcia et al., 1987;Wei and Kruis, 2013b).Among them AR is more straightforward, flexible, and much easier to be implemented.A typical rule used by AR scheme reads where r is a random number, r~U[0, 1], P c is the coagulation probability, β(u, v) is the calculated coagulation rate based on the selected particle pair with volumes u and v, and β max is the present maximum coagulation rate.The basic idea behind Eq. ( 2) is that given a fixed β max , the particle pair that leads to a larger coagulation rate tends to have a higher probability to be accepted.Although Eq. ( 2) is very straightforward, its implementation in practice suffers from a problem, namely the evaluation of β max is computationally expensive, as the cost for a successful AR attempt (i.e., for successfully finding a desirable particle pair) is of O(n s 2 ) (n s is the number of simulation particles); to make things worse, there are in general hundreds of thousands of such attempts in a single simulation.
A feasible scheme for quickly estimating β max is the smart bookkeeping technique from (Kruis et al., 2012).Although this technique was initially developed for MCs using Inverse scheme, it is also suitable for AR based MCs.The central idea behind this technique is that after each coagulation event, only the sizes of two chosen particles have experienced changes, hence the coagulation rate reevaluation can be carried out among fairly limited particle pairs; accordingly, the bookkeeping technique scans merely those particle pairs with the inclusion of at least one particle from the chosen coagulation pair, rather than the entire particle pairs.However, smart bookkeeping technique always attempts to yield an accurate β max , which might influence its computational efficiency.
Another effective method for evaluating β max is the constant scheme, which uses a prescribed, over-estimated β c max to replace the real β max (Garcia et al., 1987;Lee and Matsoukas, 2000).Since the evaluation of β c max needs to be done only once prior to the beginning of the simulation, considerable computer time can thus be saved.However, employing a constant β max is tantamount to using a large denominator in Eq. ( 2), and incurs several problems: (1) As simulation proceeds, coagulation probability P c is decreasing rapidly, leading to an excessively degenerate acceptance rate and requiring more attempts (i.e., more computer time) to pick a desirable particle pair; (2) There is no general recipe for how β c max is chosen.If a too large β c max is used, it will lead to a degenerate particle acceptance rate; if instead a small β c max is chosen, it has the danger that the denominator in Eq. (2) becomes smaller than the numerator, thus breaking the criterion down; (3) In practice the evaluation of β c max might need to be done many times, e.g., when simulating particle transport and coagulation within a computational domain composed of a number of cells (Kruis et al., 2012), max should be calculated in each cell at each time step.This note proposes a novel scheme for quickly computing the maximum coagulation rate, namely through the use of majorant kernel.The concept of majorant kernel was initially proposed by (Eibeck and Wagner, 2000) and (Goodson and Kraft, 2002) to improve the efficiency of the numerical treatment of the coagulation problems.In particular, they designed a Monte Carlo algorithm using fictitious jumps and majorant kernel (in place of the original kernel) to simulate particle coagulation.A majorant kernel has at least two attractive features that are important to the present work: (1) a majorant kernel with respect to the source coagulation kernel is always larger than its source kernel; (2) a majorant kernel can be constructed in such a way that it is quite close to its source kernel.Intrinsically and mathematically, the second feature has imposed an effective (upper) bound to the source kernel; in other words, the majorant kernel would follow its source kernel in the course of coagulation.This means that we can replace the β max in Eq. ( 2) with a majorant kernel and in this way reduce the heavy load associated with the calculation of β max .This is the central idea behind our scheme to be described below.
The rest of this article is outlined as follows.A short introduction to the MC method used is given in the next section.Next, the detailed construction of a majorant kernel, which is based on a widely used, physically realistic Brownian coagulation kernel, is presented.In what follows, numerical experiments are carried out for investigating computational accuracy and efficiency of the proposed scheme.Finally the article is concluded.

DIFFERENTIALLY WEIGHTED MC METHOD
A simplified differentially weighted MC method (Zhao et al., 2009) was adopted in the present work.This method assigns each simulation particle a weight value, which reflects the number of real particles being represented.While using this method the corresponding coagulation kernel takes the following form where w i and w j stand for the weight value associated with particle i and j, respectively; β'(d i , d j ) is the modified coagulation kernel due to the presence of particle weight.Note that hereafter the coagulation kernel β(d i , d j ) will be expressed in terms of particle sizes (i.e., not in terms of volumes like the β(u, v) in Eq. ( 1)).For a modified coagulation kernel, the previous AR criterion now becomes where β' max is the maximum coagulation rate estimated based on the kernel in Eq. ( 3).
For choosing a desirable coagulation pair, AR criterion starts from randomly selecting a pair of particles and then calculating their coagulation rate.In what follows, it checks on whether or not accept the chosen particle pair with Eq. (4); if the particle pair is accepted, update their sizes and weights; otherwise, pick another particle pair for a new examination.This procedure is repeated until a particle pair has been found successfully.Once a coagulation pair is eventually determined, the average waiting time associated with that event, 〈τ〉, can be evaluated via In a scenario with the number of simulation particles being large, one can use a single pair of particles to estimate 〈τ〉, and Eq. ( 5a) can be reduced to where (k, l) exactly identifies the chosen particle pair.Alternatively, one may increase the accuracy of Eq. (5b) by using the average of those particle pairs that were chosen but rejected for coagulation (Smith and Matsoukas, 1998); that is, where n r is the number of particle pairs being attempted to find a coagulation pair.As long as 〈τ〉 is obtained, the elapsed simulation time, t, can be updated as t → t +〈τ〉, and the coagulation procedure proceeds.

Brownian Collision Kernel
A widely-used, physically realistic Brownian collision kernel with interpolation formula of Fuchs is adopted in the present work, in which β(d i , d j ) takes the form (Jacobson, 2005) where c i and c j are the mean velocity of particle i and j, g i and g j are the Fuchs length of particle i and j, D i and D j are the Brownian diffusion coefficient associated with particle i and j.These parameters can be in turn expressed, taking particle i as an example, as where µ is the gas dynamical viscosity, m pi is the mass of particle i, A 1 , A 2 and A 3 are the coefficients of Cunningham correction factor, K n is the Knudsen number, K n = 2λ/d i where λ is the mean free path of surrounding air molecules, λ i is the mean free path of particle i, Note that both λ and µ are temperature-dependent.For example, µ can be obtained by using Sutherland-law (White, 1991) 3/2 0 0 0 where T 0 is the reference temperature, µ 0 is the reference viscosity at T 0 , S is Sutherland's constant.

Construction of Majorant Kernel
In this section we illustrate how to construct a majorant kernel,  (d i , d j ), given the coagulation kernel in Eq. (6a).First, it is known from the features of majorant kernel that where d min and d max represent the timely minimum and maximum particle diameter.As both d min and d max tend to change in coagulation, their values need to be updated from time to time.Obviously it takes a O(n s ) computer time to find d min and d max for n s simulation particles.
In order to construct  (d i , d j ), we have to first simplify the original kernel β(d i , d j ), as it is rather complicated.We rearrange β(d i , d j ) as Then we conclude from Eq. ( 8) that and Next we can construct a majorant kernel for each of the simplified kernel.We may follow the steps in Appendix B and obtain the following majorant kernel, 1  (d i , d j ), for with , where ρ p is the particle density. 1  (d i , d j ) in Eq. ( 10) becomes now a polynomial with each term being a power function in terms of particle diameter d i and d j .Similarly, one possible majorant kernel with respect to β 2 (d i , d j ), 2  (d i , d j ), can be obtained as (Appendix B) where K fm is the coagulation coefficient, Once majorant kernels 1  (d i , d j ) and 2  (d i , d j ) are available, an effective majorant kernel with respect to kernel β(d i , d j ) can be constructed as ˆˆ, min , , ,

Maximum of the Majorant Kernel
With majorant kernel with respect to β(d i , d j ) being available, its maximum value over the size interval [d min , d max ], max  (d i , d j ), can be evaluated from Eq. ( 12) as where 1,max  (d i , d j ) and 2,max  (d i , d j ) stand for the maximum coagulation rate of kernel β 1 (d i , d j ) and β 2 (d i , d j ), respectively.
The most straightforward scheme for quickly evaluating {k},max  (k = 1, 2) is to put d min and d max directly into Eqs.( 10) and ( 11), respectively, and take Due to the symmetric of the majorant kernels the order of d min and d max in Eq. ( 14) is not important.A pictorial proof for justifying this scheme can be seen in Fig. 1, in which the spatial distribution of β(d i , d j ) over the size interval [1, 100] at T = 300 K is shown.As can be seen, the maximum coagulation rate peaks at (d i , d j ) = (1, 100) nm (or (100, 1) nm), where the largest size difference between the two particles is reached.
As long as a possible version of the majorant kernel of β(d i , d j ) is available, we can immediately construct a majorant Fig. 1.The implementation of MC method using mojorant kernel scheme for particle coagulation.Note that a conditional update strategy is adopted for updating the particle size interval.
kernel for the modified kernel in Eq. ( 3) as where w max is the maximum particle weight.The rightmost term in the last line of Eq. ( 15) gives a feasible majorant kernel with respect to the modified kernel in Eq. ( 3 is now obtained in an analytical way, it is much faster than directly evaluating β' min , which is achieved by comparing a fair amount of particle pairs in a discrete manner.

Algorithm Implementation
The flow chart for the implementation of the proposed method in this article is schematically sketched in Fig. 2, in which t sim is the total simulation time.The overall procedure is straightforward.A particular point which deserves comment is the update of the particle size interval [d min , d max ].At first glance it would seem that the simplest way to do this would be by comparing d min and d max with all simulation particles after each collision event.However, one may do this more smartly through the use of a conditional update strategy, as one or more successive coagulation events may not give rise to a change in d max, , as a result, the previous size interval remains unchanged.Clearly, the conditional update strategy is a good choice for reducing computational cost.The recipe for implementing conditional update is also straightforward: The present d max is compared with the sizes of newly formed particles after each collision event; if their sizes are found to be smaller than d max , the preceding size interval will be kept; otherwise all particles will be checked to find the new d min and d max .

Algorithm Verification
The computational accuracy of the majorant kernel scheme described above is verified by comparing with a (deterministic) sectional method (Lu, 1994) using 200 sections and a sectional spacing of 1.12.The following identical simulation parameters are adopted for both the benchmark and Monte Carlo scheme: initial particle diameter d 0 = 1 nm, number concentration N 0 = 10 17 /m 3 , particles density ρ p = 1000 kg/m 3 , initial geometric standard deviation of particle size distribution σ g = 1.In addition, total simulation time, t sim , is set to in the interval [10 -3 , 10 3 ]τ char with τ char being the characteristic coagulation time.For particles in the free-molecule regime τ char can be given by (Kodas, 1999): where V t is the aerosol volume fraction, V t = 3 0 0 / 6 d N  .Specifically, with the foregoing parameters τ char can be computed as 0.0897 s.The simulations were performed on a desktop system with an Intel Core™ 2 Quad 2.66 GHz Processor and 8 GB RAM.
The Monte Carlo codes were run for 10 times and the average results with respect to particle size, d p , and standard geometric deviation, σ g , are shown in Figs. 3 and  4, respectively.It is seen that d p has witnessed a monotonic growth from a low of 1nm at t ≈ 0 up to a high of nearly 20 nm at t = 1000 τ char , whereas σ g increases from around 1 to a peak value 1.46 after about t = 10 τ char , and then decreases to 1.43.The simulation results are in very good agreement with the benchmark results.The computational accuracy of the proposed scheme can further be assessed by ε mean , a calculation error defined as (using d p as an example)  with n r being the simulation runs.Given n r = 10 and n t = 9, the errors for d p and σ g under different particle number n s are shown in Fig. 5.The overall errors for both d p and σ g are less than 0.4%, indicating an extremely good agreement between the two results.An interesting point is that a larger number of simulation particles do not systematically lead to a smaller calculation error, which may be due to the statistical errors in the computing.

Computational Efficiency
Ideally, the best computational efficiency will be obtained when β max used in Eq. ( 2) takes its true value; however, the value calculated from the majorant kernel is always somewhat larger, i.e., the value used in practice is overestimated.To quantitatively measure the closeness between the overestimated and the true value of the maximum where the denominator β' max (d i , d j ) stands for the true value.
It is clear that ψ ≥ 1 and its value dynamically reflects the discrepancy between the estimated and the true maximum coagulation rate in the course of coagulation.The smaller ψ will lead to a better computational efficiency, as less attempts for finding coagulation particle pairs are needed.The evolution of ψ in a typical coagulation procedure is shown in Fig. 6.It can be seen that it is less than 9 over the entire time interval, indicating that the majorant kernel (Eq.( 12)) follows its original kernel (Eq.( 8)) very well.
The evolution of d max in the same coagulation process is shown in Fig. 7, where d max increases in a step-wise manner with simulation time is clearly seen.This fact implies that d max tends to remain unchanged within some time intervals, and no update to particle size interval is unnecessary within those time intervals.Indeed, this fact has justified the conditional update strategy described in the previous section.Fig. 8 compares computer time (wall-clock time) needed by the forgoing different schemes under particle number n s = 5000.It can be clearly seen that among all the nondeterministic MC schemes, smart bookkeeping technique has the lowest computational efficiency, while majorant kernelbased scheme exhibits the best computational performance, and their performance gap is fairly pronounced.For example, the majorant scheme is 10 2 times faster than bookkeeping technique even for a long simulation time up to t/τ char = 100.The reason is that for both smart bookkeeping and constant scheme, an initial calculation for determining the maximum coagulation rate is required.Although this operation needs to be done only once, however, it is computationally expensive due to the O(n s 2 ) running time.As a result, computer time used for this calculation is much longer than that for the real coagulation computing.Indeed, this is a bottleneck for MCs using smart bookkeeping or constant scheme in some practical calculations, e.g., when coupling MC with CFD to handle particle coagulation and transport simultaneously.Suppose that a reactor domain consists of a number of cells (Kruis et al., 2012).Then if MC using either bookkeeping or constant scheme is employed, the maximum coagulation rate has to be evaluated in each cell and in each time step, and the resultant computational cost would be formidable.In this situation, the majorant kernel scheme that can avoid the direct calculation of the maximum coagulation rate has an obvious edge.Also, Fig. 8 presents the computer time of the benchmark, the sectional method, as indicated by the thick line.Somewhat surprisingly, majorant kernel scheme with a moderate number of particle is still superior to the deterministic method.
Fig. 9 compares the computer time required by different MC schemes for accomplishing the simulation up to t sim = 1000 τ char , but with different number of simulation particles, i.e., n s = 10 3 -10 4 .As can be seen, the majorant scheme spends less than 6s in completing the calculation under the largest particle number (n s = 10000) and longest simulation  time (t sim = 1000 τ char ); in contrast, it takes more than 1200 s for the bookkeeping technique to accomplish the same task.

CONCLUSIONS
A fast, practical MC scheme using majorant kernel concept for particle coagulation has been presented in this note.With this scheme, the evaluation of the maximum coagulation rate, which is computationally time-consuming and a bottleneck for MC, can be performed very efficiently.Satisfactory computational accuracy of the scheme has been shown by comparing it with a sectional benchmark.The numerical experiments also show that in terms of computational efficiency, the scheme proposed in this note far outperforms the sectional method and MCs using other schemes.As a result, this scheme would be a worthwhile choice for various MCs applied to some computationally intensive calculations; e.g., when MCs are chosen to model particle dynamics in a reactor together with a CFD solver, as described above.

ACKNOWLEDGMENT
I would like to acknowledge helpful discussion with Prof. H. Zhao.

APPENDIX A: JUSTIFICATION OF TWO LEMMAS
Proof of Lemma 1.Without loss of generality, one may assume a ≥ b, and denote f = a/b, clearly f ≥ 1.Further, denote m as m = [m] + {m}, where [m] is the integer part of m and {m} is the fractional part.To prove lemma 1 one requires: The last line of above expression is obvious as for f ≥ 1, f {m} ≥ 1; lemma 1 hence holds.From lemma 1, one has since 1/ε > 1 for ε (0, 1); hence, lemma 2 holds true.

APPENDIX B: CONSTRUCTING A MAJORANT KERNEL OF β 1 (d i , d j ) AND β 2 (d i , d j )
We first prove g i ≤ λ i and g j ≤ λ j .Simplifying Eq. (6c) yields Multiplying the numerator and denominator of the first term of the right-hand side of above equation by a factor d i From lemma 2, we have Similarly, we can prove that g j ≤ λ j .Then, using these relationships and lemma 2, we can follow from Eq. (9a) that ) using a simple inequality relation: a•b/(a + b) < a or b,  a, b > 0. Denote the kernels on the right-hand side of Eqs.(9a) and (9b) as β 1 (d i , d j ) and β 2 (d i , d j ), respectively.

Fig. 2 .
Fig. 2. The shape of Brownian collision kernel for particles in size interval [1, 100] nm at a temperature T = 300 K.

Fig. 3 .
Fig. 3. Comparison of σ g between MC and sectional method as a function of dimensionless simulation time.

Fig. 4 .
Fig. 4. Comparison of d p between MC and sectional method as a function of dimensionless simulation time.

Fig. 5 .
Fig. 5. Calculation errors ε mean for d p and σ g as a function of number of simulation particles.

Fig. 8 .
Fig. 8.Comparison of computer times for MCs using various schemes as well as sectional method.

Fig. 9 .
Fig. 9. Particle number dependence of the computer time for MCs using various schemes.
c i and c j in Eq. (6b) into Eq.