**1. Introduction**

Aerosols are dispersed systems of solid or liquid particles suspended in air. Aerosols are created during various natural processes and industrial operations [1], such as clouds, smoke particles from forest fires, sand and dust storms, volcanic dust, spores and seeds from plant life, fluidized catalysts, spray drying of fluids, atomization of liquid fuels, fly ash from stacks, soot in combustion flames, etc. Aerosol evolution is very complex, involving various physical

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

and chemical processes. The most popular description of aerosol evolution is using the population balance equation (PBE) [2], which is also called the general dynamic equation in aerosol field [3]. The PBE is a convection‐diffusion type equation with source terms taking into account the particle formation and depletion from various aerosol dynamics, mainly nuclea‐ tion, surface growth (condensation), and coagulation.

Analytical solutions to the PBE are available for only a few specific cases [4–7]. Mostly, numerical means are resorted to solve the PBE. There are generally three typical methods for solving the PBE, i.e., direct discretization, method of moments, and stochastic methods. Direct discretization may be in the particle size space, such as discretization by sections [8] or finite difference/finite element [9] and may also be in the functional space, such as those in the collocation method [10], Galerkin method [11], and weighted residual method [12]. Among the various discretizations, the sectional method pioneered by Gelbard et al. [8] is the most popular. There are quite a few sectional methods developed aiming to alleviate the nonposi‐ tiveness and nonconservation problems found in earlier works [13], and to obtain better efficiency and flexibility, such as the fixed/moving pivot method [14, 15] and cell average method [16, 17].

Other than solving the discretized PBE, a widely used alternative approach is the method of moments, which solves a group of moments equations derived from the PBE. These moments (generally several low order moments) provide important information on the particle size distribution (PSD) function, i.e., number density, volume fraction, polydispersity, etc. How‐ ever, owing to the nonlinearity of the PBE, the governing equations for lower order moments generally contain higher order moments, which are not closed after cutting off at a certain order of moment. Various closure methods have been developed, such as unimodal log‐normal approximation [18], quadrature/direct quadrature [19, 20], interpolative closure [21], Taylor expansion [22], etc. 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. Another issue closely tied to method of moments is the realizability problem, when the moments sequence fails to find its corresponding distribution function [23, 24]. Among all the drawbacks, the inability to provide the PSD is the most severe, when the PSD is highly needed.

Stochastic simulation (or Monte Carlo simulation) does not solve the PBE, but to mimic the evolution of aerosol particles through a stochastic particle system. Based on the Marcus‐ Lushnikov process modeling of coagulation [29], various stochastic methods have been developed to simulate the general aerosol evolution. It is very difficult to categorize precisely so many diversified methods. According to the evolution of time, there is time‐driven or event‐ driven Monte Carlo [30, 31]. According to the different ways to weight the numerical particles, there are the mass flow algorithms [32, 33], differentially weighted algorithm [34], and the more general weighted flow algorithms [35]. Also, there are various types of methods to improve the simulation efficiency, such as the majorant kernel method [36] and the leaping method [37].

Among the three typical methods, stochastic simulation is the most flexible for accommodating complex physical models for various aerosol dynamics, and is easy to implement without special numerical issues, which are often encountered in the other two methods, e.g., nonpo‐ sitiveness and nonconservation problems in direct discretization, realizability problem in method of moments. However, stochastic simulation generally has convergence rate propor‐ tional to , where is the number of total numerical particles used. The simulation cost (CPU time and memory) is usually proportional to 2 [25]. Hence, stochastic simulation is not numerically efficient when pursuing a high precision solution. Recently, Zhou et al. [26] developed an operator splitting Monte Carlo (OSMC) method, which combines stochastic simulation for coagulation process and deterministic integration for nucleation and surface growth (condensation). Stochastic simulation for coagulation is much more efficient than directly solving the Smoluchowski integro‐differential equation (the coagulation model equation) with traditional discretization [27]. On the other hand, deterministic integration for nucleation and surface growth processes is far more efficient than stochastic simulations. The stochastic simulation and the deterministic integration in OSMC are synthesized under the framework of operator splitting. The accuracy and efficiency of the OSMC have been demon‐ strated through various testing cases [26]. Especially, the OSMC has been implemented with Message Passing Interface (MPI) to achieve very promising parallel simulation efficiency [26]. Here, the impact of various operator splitting schemes on the accuracy and efficiency of the OSMC is addressed.

## **2. Methodology**

and chemical processes. The most popular description of aerosol evolution is using the population balance equation (PBE) [2], which is also called the general dynamic equation in aerosol field [3]. The PBE is a convection‐diffusion type equation with source terms taking into account the particle formation and depletion from various aerosol dynamics, mainly nuclea‐

Analytical solutions to the PBE are available for only a few specific cases [4–7]. Mostly, numerical means are resorted to solve the PBE. There are generally three typical methods for solving the PBE, i.e., direct discretization, method of moments, and stochastic methods. Direct discretization may be in the particle size space, such as discretization by sections [8] or finite difference/finite element [9] and may also be in the functional space, such as those in the collocation method [10], Galerkin method [11], and weighted residual method [12]. Among the various discretizations, the sectional method pioneered by Gelbard et al. [8] is the most popular. There are quite a few sectional methods developed aiming to alleviate the nonposi‐ tiveness and nonconservation problems found in earlier works [13], and to obtain better efficiency and flexibility, such as the fixed/moving pivot method [14, 15] and cell average

Other than solving the discretized PBE, a widely used alternative approach is the method of moments, which solves a group of moments equations derived from the PBE. These moments (generally several low order moments) provide important information on the particle size distribution (PSD) function, i.e., number density, volume fraction, polydispersity, etc. How‐ ever, owing to the nonlinearity of the PBE, the governing equations for lower order moments generally contain higher order moments, which are not closed after cutting off at a certain order of moment. Various closure methods have been developed, such as unimodal log‐normal approximation [18], quadrature/direct quadrature [19, 20], interpolative closure [21], Taylor expansion [22], etc. 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. Another issue closely tied to method of moments is the realizability problem, when the moments sequence fails to find its corresponding distribution function [23, 24]. Among all the drawbacks, the inability to provide the PSD is the most severe,

Stochastic simulation (or Monte Carlo simulation) does not solve the PBE, but to mimic the evolution of aerosol particles through a stochastic particle system. Based on the Marcus‐ Lushnikov process modeling of coagulation [29], various stochastic methods have been developed to simulate the general aerosol evolution. It is very difficult to categorize precisely so many diversified methods. According to the evolution of time, there is time‐driven or event‐ driven Monte Carlo [30, 31]. According to the different ways to weight the numerical particles, there are the mass flow algorithms [32, 33], differentially weighted algorithm [34], and the more general weighted flow algorithms [35]. Also, there are various types of methods to improve the simulation efficiency, such as the majorant kernel method [36] and the leaping

Among the three typical methods, stochastic simulation is the most flexible for accommodating complex physical models for various aerosol dynamics, and is easy to implement without

tion, surface growth (condensation), and coagulation.

method [16, 17].

2 Aerosols - Science and Case Studies

when the PSD is highly needed.

method [37].

#### **2.1. Governing equations**

Aerosol particles are generally described statistically in a mesoscopic scale through defining the PSD, ( , ,), and ( , ,) denotes the number density of aerosol particles at location at time with aerosol properties between and + . Here, is usually called an internal variable, which can be volume, surface, constituent, etc. For simplicity, it is assumed to be the particle size (volume or diameter specified by the context). It is worth pointing out that the OSMC method can also deal with the cases with a vector type , i.e., combined properties to describe particles. In fact, the current OSMC code treats the diameter, surface, and volume of an aerosol particle independently, which can easily model particles of various morphologies, such as particles with fractal dimensions. The evolution of the PSD function satisfies the PBE (also called the general dynamic equation) [3]

$$\frac{\partial \mathbf{n}}{\partial t} + \nabla \cdot \mathbf{n} \vec{u} = \nabla \cdot \mathbf{D} \nabla \mathbf{n} + I\_{\text{mc}} + G\_{\text{cond}} + C\_{\text{coag}} \tag{1}$$

where is the gas flow velocity, and is the diffusion coefficient. The source terms of *I*nuc, *G*cond , and *C*coag on the right hand side denote nucleation, surface growth (condensation), and coagulation, respectively. The determination of the convection velocity requires solving the

fluid dynamics control equations, and the coupling of the fluid dynamics with the aerosol dynamics is another important research field. Here, the gas velocity is assumed known, using some techniques (often computational fluid dynamics). It is worth pointing out that the coupling between the fluid dynamics and the aerosol dynamics is generally handled by the operator splitting method [11, 28–30], i.e., to solve the fluid dynamics and the aerosol dynamics one by one with the feedback from the other dynamics temporarily freezing, and to iterate between these two dynamics to obtain a converged solution.

Another simplification in this work is to neglect the diffusion term in the PBE. For aerosol particles, the mass diffusion due to the number concentration gradient (∇⋅∇) is usually very small compared to the air viscous diffusion, which can be described by the Schmidt number (i.e., the ratio of the viscous diffusion to the mass diffusion). For example, under standard atmospheric conditions, a spherical particle with a diameter of 10 nm has a Schmidt number equal to 290, and a 100 nm particle has Schmidt number 2.2 × 10<sup>4</sup> [3]. So it is legitimate to neglect the aerosol mass diffusion under very general conditions.

With the assumptions that the gas velocity is known and the mass diffusion term is negli‐ gible, the PBE can be simplified as

$$\frac{\partial \mathfrak{m}}{\partial t^\*} = I\_{\text{nuc}} + G\_{\text{cond}} + \mathbb{C}\_{\text{coag}}.\tag{2}$$

The convection term on the left hand side of the original PBE (1) has been implicitly solved with the introduction of the Lagrangian time \*

$$\text{tr}^\*\left(\vec{\mathbf{x}}\right) = \int\_0^{\vec{\mathbf{x}}} \frac{d\vec{\mathbf{x}}\prime}{\vec{\mathbf{u}}\prime\prime\prime}\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag\tag$$

The transformation converts the Eulerian point of view in Eq. (1) to the Lagrangian point of view in Eq. (2), with the help of the assumed known velocity field . There are many practically interested conditions, where the flow field can be considered as steady and one‐dimensional, e.g., burner stabilized flames and counter flow diffusion flames. Then it is very easy to integrate Eq. (3) to build up a one‐to‐one mapping between the time and the spatial coordinate. Hence, the convection can be readily assimilated in the model Eq. (2). 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 [31, 32]. So the model Eq. (2) and the subsequent OSMC method may have very broad applications. Later on, \* is simply denoted as , since the Lagrangian point of view is not explicitly referred any longer.

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. Surface growth includes physical condensation and surface chemical reactions, which involves interactions between gas phase molecules and an aerosol particle. Surface chemical reactions will not be discussed here due to their vastness. Coagulation is the process that two particles coalesce to form a new one. More details about the aerosol dynamic processes are presented in Section 2.2.

**Figure 1.** Flowchart of the OSMC. Filled blocks denote steps where stochasticity arises. Adapted from Reference [26].

#### **2.2. Operator splitting Monte Carlo method**

fluid dynamics control equations, and the coupling of the fluid dynamics with the aerosol dynamics is another important research field. Here, the gas velocity is assumed known, using some techniques (often computational fluid dynamics). It is worth pointing out that the coupling between the fluid dynamics and the aerosol dynamics is generally handled by the operator splitting method [11, 28–30], i.e., to solve the fluid dynamics and the aerosol dynamics one by one with the feedback from the other dynamics temporarily freezing, and to iterate

Another simplification in this work is to neglect the diffusion term in the PBE. For aerosol particles, the mass diffusion due to the number concentration gradient (∇⋅∇) is usually very small compared to the air viscous diffusion, which can be described by the Schmidt number (i.e., the ratio of the viscous diffusion to the mass diffusion). For example, under standard atmospheric conditions, a spherical particle with a diameter of 10 nm has a Schmidt number equal to 290, and a 100 nm particle has Schmidt number 2.2 × 10<sup>4</sup> [3]. So it is legitimate to

With the assumptions that the gas velocity is known and the mass diffusion term is negli‐

The convection term on the left hand side of the original PBE (1) has been implicitly solved

*u x*

The transformation converts the Eulerian point of view in Eq. (1) to the Lagrangian point of view in Eq. (2), with the help of the assumed known velocity field . There are many practically interested conditions, where the flow field can be considered as steady and one‐dimensional, e.g., burner stabilized flames and counter flow diffusion flames. Then it is very easy to integrate Eq. (3) to build up a one‐to‐one mapping between the time and the spatial coordinate. Hence, the convection can be readily assimilated in the model Eq. (2). 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 [31, 32]. So the model Eq. (2) and the subsequent OSMC method may have very broad applications. Later on, \* is simply

'

¶ (2)

<sup>r</sup> <sup>r</sup> <sup>r</sup> r r (3)

*<sup>n</sup> IG C*

*<sup>x</sup> dx t x*

denoted as , since the Lagrangian point of view is not explicitly referred any longer.

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. Surface growth includes physical

0

( )= , ( ') ò

\*

\* nuc cond coag = . ¶ + +

between these two dynamics to obtain a converged solution.

neglect the aerosol mass diffusion under very general conditions.

*t*

gible, the PBE can be simplified as

4 Aerosols - Science and Case Studies

with the introduction of the Lagrangian time \*

Operator splitting is an efficient and powerful method to solve such evolution equations [33] as Eq. (2). In an operator splitting method, the time evolution of () is still handled by the classical methods, such as the forward Euler method or Runge‐Kutta method. Other than to integrate all aerosol dynamic processes together in one time step , operator splitting separates the integration into multiple fractional steps. There are various ways to construct the multiple fractional steps, such as nuc cond coag or 1 2 nuc cond 1 2coag, where 1/2 denotes a half integration step. The flowchart (**Figure 1**) shows an example how the operator splitting scheme nuc cond coag is carried out. Details of the various models for the aerosol dynamics are described in the following subsections, and the error analyses for various operator splitting schemes are given in Section 3.

#### *2.2.1. Nucleation*

Nucleation is the gas‐to‐particle transition to generate nascent nanoparticle. There are various nucleation models, largely based on the classical nucleation theory developed in the 1920– 1940s [34]. In the classical nucleation theory, there are two key parameters for the nucleation modeling, i.e., the nucleation rate (how many particles are nucleated per unit volume, unit time), and the critical size, which is the minimum size of nucleated particles with stable nucleus. According to the classical Becker‐Döring theory [3, 35], the nucleation rate and critical diameter ∗ are

$$I = \frac{P\_v \chi\_v}{k\_B T} \sqrt{\frac{2\sigma}{\pi m}} \exp\left(-\frac{16\pi \sigma^3 m^2}{\Im(k\_B T)^3 \rho\_p^2 (\ln S)^2}\right) \tag{4}$$

$$d\_p^\* = \frac{4\sigma v\_w}{k\_B T \ln S}.\tag{5}$$

Here, is the vapor pressure, is the mole fraction of the vapor, is the molecular mass of the nucleation vapor, is the particle density, is the temperature, <sup>=</sup> /sat is the vapor saturation ratio, is the Boltzmann constant, is the surface tension, and <sup>=</sup> / is the molecular effective volume of the nucleation vapor. Hence, the nucleation term in Eq. (2) can be expressed with the help of the Dirac delta function as

$$I\_{\rm nuc} = \int I \mathcal{S}(d\_p - d\_p^\*) \mathbf{d}d\_{p'} \tag{6}$$

which means that only particles with the critical size are nucleated.

Within the OSMC, nucleation model is not limited to the classical nucleation theory. Any other nucleation model can be readily included, should only the nucleation rate and the critical size be provided.

#### *2.2.2. Surface growth (condensation)*

The functional form of the particle volume growth rate due to condensation depends on the Knudsen number Kn = 2/, i.e., the ratio between the air mean free path and the particle

radius /2. In the free molecular regime (Kn <sup>≥</sup> <sup>10</sup>) and the continuum regime (Kn <sup>≤</sup> 0.1), the condensational growth rates can be derived from the gas kinetic theory and the gradient diffusion model, respectively. In the transition regime (0.1 < Kn < 10), interpolation formulas are usually used, such as the harmonic mean of the free molecular and continuum regimes. The expressions for volume growth rate and the condensation term cond are [3]

$$\mathbf{G}\_v = \begin{cases} \frac{\pi d\_p^2 (P\_v - P\_p) \upsilon\_m}{(\mathbf{Z} \pi m k\_p T)^{1/2}}, & \mathbf{Kn} \ge 10 \\\\ \frac{2 \pi D d\_p (P\_v - P\_p) \upsilon\_m}{k\_p T}, & \mathbf{Kn} \le 0.1 \\\\ \text{harmonic mean}, & 0.1 \le \mathbf{Kn} \le 10 \end{cases} \tag{7}$$

$$
\mathcal{G}\_{\rm cond} = -\frac{\partial \langle \mathfrak{n} G\_v \rangle}{\partial v}. \tag{8}
$$

Here, is the equilibrium vapor pressure for particles of size , which is usually assumed to be the saturation pressure of the condensation vapor. For particles close to the critical size, the Kelvin correction is also required [3].

Within the OSMC, any other surface growth model can also be accommodated, should only the particle volume growth rate be provided.

#### *2.2.3. Stochastic simulation of coagulation*

half integration step. The flowchart (**Figure 1**) shows an example how the operator splitting

dynamics are described in the following subsections, and the error analyses for various

Nucleation is the gas‐to‐particle transition to generate nascent nanoparticle. There are various nucleation models, largely based on the classical nucleation theory developed in the 1920– 1940s [34]. In the classical nucleation theory, there are two key parameters for the nucleation modeling, i.e., the nucleation rate (how many particles are nucleated per unit volume, unit time), and the critical size, which is the minimum size of nucleated particles with stable nucleus. According to the classical Becker‐Döring theory [3, 35], the nucleation rate and

nuc cond coag is carried out. Details of the various models for the aerosol

3 2 32 2

(4)

(6)

(5)

æ ö ç ÷ è ø

ps

*m*

Here, is the vapor pressure, is the mole fraction of the vapor, is the molecular mass of the nucleation vapor, is the particle density, is the temperature, <sup>=</sup> /sat is the vapor saturation ratio, is the Boltzmann constant, is the surface tension, and <sup>=</sup> / is the molecular effective volume of the nucleation vapor. Hence, the nucleation term in Eq. (2) can

r

scheme

*2.2.1. Nucleation*

6 Aerosols - Science and Case Studies

critical diameter ∗

be provided.

*2.2.2. Surface growth (condensation)*

operator splitting schemes are given in Section 3.

are

*v v*

be expressed with the help of the Dirac delta function as

which means that only particles with the critical size are nucleated.

*B B p P x <sup>m</sup> <sup>I</sup> xp kT m kT S*

> *B <sup>v</sup> <sup>d</sup> kT S* <sup>4</sup> = . ln

*ppp I Id d d* nuc = ( )d , \* - ò d

Within the OSMC, nucleation model is not limited to the classical nucleation theory. Any other nucleation model can be readily included, should only the nucleation rate and the critical size

The functional form of the particle volume growth rate due to condensation depends on the Knudsen number Kn = 2/, i.e., the ratio between the air mean free path and the particle

s

*p*

\*

2 16 = e , 3( ) (ln )

s

p

> Coagulation is the process that two particles coalesce/aggregate to form a new particle. The coagulation dynamics is modeled by the well‐known Smoluchowski's equation [3]

$$\mathcal{L}\_{\text{coag}} = \frac{1}{2} \int\_0^v \beta(\upsilon, \upsilon - \tilde{\upsilon}) n(\tilde{\upsilon}) n(\upsilon - \tilde{\upsilon}) d\tilde{\upsilon} - \int\_0^v \beta(\upsilon, \tilde{\upsilon}) n(\upsilon) n(\tilde{\upsilon}) d\tilde{\upsilon}. \tag{9}$$

The coagulation kernel function (,) describes the rate at which particles of size coagulate with particles of size . The first term on the right‐hand side denotes the formation rate of particles with volume due to coagulation between particles and . (Volume is conserved before and after coagulation.) The half factor (1/2) is to correct the double‐counting effect. The second term on the right‐hand side describes the depletion rate of particles with volume due to their coagulation with any other particles. Practically, the coagulation kernel function (,) is very hard to find, except for simplified conditions. In the free molecule regime (Kn ≥ 10), has the form [3]

#### 8 Aerosols - Science and Case Studies

$$
\beta \beta (v, \mu) = \left(\frac{6}{\pi}\right)^{2/3} \left(\frac{\pi k\_B T}{2\rho\_p}\right)^{1/2} \left(\frac{1}{v} + \frac{1}{\mu}\right)^{1/2} \left(v^{1/3} + u^{1/3}\right)^2. \tag{10}
$$

In the continuum regime (Kn ≤ 0.1), takes the form [3]

$$\mathcal{J}(\boldsymbol{v},\boldsymbol{u}) = \frac{2k\_B T}{3\mu} \left( \frac{1}{\boldsymbol{v}^{1/3}} + \frac{1}{\boldsymbol{u}^{1/3}} \right) (\boldsymbol{v}^{1/3} + \boldsymbol{u}^{1/3}),\tag{11}$$

where is the dynamic viscosity of the gas‐phase. In the transition regime (0.1 < Kn < 10), interpolation formula taking into account the former two regimes are generally used.

Although the Smoluchowski's equation was developed one century ago, and it has received broad attention in diverse fields, no analytical solutions are known for the kernels discussed above. Analytical solutions known so far are only limited to the following three simple kernels (or their linear combination) [5], i.e., (i) (,)=1, (ii) (,) = + , and (iii) (,) = .

Other than to solve the Smoluchowski's equation directly through various numerical means, the coagulation process is simulated by the Marcus‐Lushnikov stochastic process in the OSMC. The Marcus‐Lushnikov process is the natural stochastic analogue of the finite discrete Smoluchowski coagulation equation. The well‐known review paper [5] provides a wide‐ ranging survey on the Marcus‐Lushnikov process and the Smoluchowski's equation. Here, we follow the work of Gillespie [25], and give an intuitive introduction to the stochastic simulation of coagulation.

In a virtual box of volume , there are "well‐mixed" particles. Any two particles may coalesce/aggregate to form a new one at a random time. The propensity of coalescence between two particles and is determined by the coagulation kernel <sup>=</sup> (,)/. Let (,,) denote the probability that the next coalescence will occur in the time interval (, + ) between and ( < ), then it is found [25]

$$P(\tau, i, j) = C\_{\downarrow} \exp\left[ -\sum\_{k=1}^{N-1} \sum\_{l=k+1}^{N} C\_{kl} \tau \right]. \tag{12}$$

In essence, the simulation is to randomly pick up two particles and to coalesce (forming a new particle) at the random time according to the joint probability density function (,,). Repeat the former process as time advances. During the simulation, the number of particles drops, and the particle number density drops also.

The key step in the simulation is to generate a random point (,,) according to the density function (,,). The random coagulation time satisfies a Poisson process (the exponential distribution), and it can be generated as

Operator Splitting Monte Carlo Method for Aerosol Dynamics http://dx.doi.org/10.5772/65140 9

$$
\sigma = \frac{1}{C\_o} \ln \left( \frac{1}{r\_1} \right) \tag{13}
$$

where 0 <sup>=</sup> ∑ = 1 1∑ <sup>=</sup> + 1 , and 1 is a random number picked from the uniform distribution in the unit interval. Eq. (13) gives the classical way to generate a random variate from the exponential distribution with parameter 0.

After selecting the time at which two particles coagulate, the next step is to randomly select two particles according to the probability function 2(,) = /0. Loosely speaking, 0 is the total coagulation probability between any two particles, and 2(,) is the fraction of the coagulation probability over 0. Gillespie [25] introduced two methods to choose and , the so‐called full conditioning and partial conditioning methods. In the full conditioning method, the selection of and is accomplished in two steps. During the first step, calculate the partial sum <sup>=</sup> ∑ <sup>=</sup> + 1 , ( = 1,…,), then combine the partial sums in sequential to form a range with the index from 1 to denoting the part for each partial sum, and then it is straightforward to build up a linear mapping between the range and the unit interval (0,1). Draw a random number from the uniform distribution, and check the random number landing on which corresponding part in the range, and then select particle by the part index. During the second step, combine , ( <sup>=</sup> + 1,…,) in sequential to form another range, and build up the linear mapping between the new range with the unit interval, then draw another random number to determine the particle in a similar way as in the first step. In the partial conditioning method, it is actually an acceptance‐rejection method to generate a bivariate random number (,) from the function 2(,). In a general acceptance‐rejection algorithm to

generate a vector random variable with density ( ), it needs to choose a density function ( ) and a constant 1 such that ( ) ( ). The algorithm proceeds as follows [36]:

**1.** Generate from density .

*B p k T v u v u*

 r

p

*<sup>B</sup> k T*

m

p

In the continuum regime (Kn ≤ 0.1), takes the form [3]

b

(, + ) between and ( < ), then it is found [25]

drops, and the particle number density drops also.

distribution), and it can be generated as

t

b

8 Aerosols - Science and Case Studies

of coagulation.

1/2 2/3 1/2 <sup>6</sup> 1 1 1/3 1/3 2 ( , )= ( ). <sup>2</sup> æ ö æö æ ö

> *v u v u v u*

interpolation formula taking into account the former two regimes are generally used.

1/3 1/3 <sup>2</sup> 1 1 ( , )= ( ), <sup>3</sup> æ ö ç ÷ + + è ø

where is the dynamic viscosity of the gas‐phase. In the transition regime (0.1 < Kn < 10),

Although the Smoluchowski's equation was developed one century ago, and it has received broad attention in diverse fields, no analytical solutions are known for the kernels discussed above. Analytical solutions known so far are only limited to the following three simple kernels (or their linear combination) [5], i.e., (i) (,)=1, (ii) (,) = + , and (iii) (,) = .

Other than to solve the Smoluchowski's equation directly through various numerical means, the coagulation process is simulated by the Marcus‐Lushnikov stochastic process in the OSMC. The Marcus‐Lushnikov process is the natural stochastic analogue of the finite discrete Smoluchowski coagulation equation. The well‐known review paper [5] provides a wide‐ ranging survey on the Marcus‐Lushnikov process and the Smoluchowski's equation. Here, we follow the work of Gillespie [25], and give an intuitive introduction to the stochastic simulation

In a virtual box of volume , there are "well‐mixed" particles. Any two particles may coalesce/aggregate to form a new one at a random time. The propensity of coalescence between two particles and is determined by the coagulation kernel <sup>=</sup> (,)/. Let (,,) denote the probability that the next coalescence will occur in the time interval

> *N N ij kl k lk*

+ é ù -ê ú ë û

 t

å å (12)

1 =1 = 1

In essence, the simulation is to randomly pick up two particles and to coalesce (forming a new particle) at the random time according to the joint probability density function (,,). Repeat the former process as time advances. During the simulation, the number of particles

The key step in the simulation is to generate a random point (,,) according to the density function (,,). The random coagulation time satisfies a Poisson process (the exponential

*P i j C xp C*

( , , )= e . -

ç÷ ç ÷ ç ÷ + + ç ÷ èø è ø è ø

*v u*

1/3 1/3

(10)

(11)


The density function ( ) can be chosen freely. However, the optimal choice is to make the interval in 2 as small as possible, so that there is a better chance that ≤ ( ) so as to accept as the required without too often rejections and repetitions of the whole process. In the acceptance‐rejection method [25] to generate and according to 2(,), ( ) is chosen as a uniform function, and is an upper bound of all , (preferably the smallest upper bound). It is argued [25] that when the values of are relatively comparable with each other (ideally constant) the particle conditioning method should be more efficient than the full condition‐ ing method.

Within the OSMC, both the full and particle conditioning methods have been implemented to choose from easily.

## **2.3. Full algorithm for the OSMC**

**Figure 1** shows the flowchart of OSMC including all aerosol dynamics. Only the first‐order operator splitting is sketched in the figure.

At time = 0, the simulator is initialized with given parameters, i.e., number of Monte Carlo (MC) particles , simulator size , and the initial PSD. The simulator size is determined as <sup>=</sup> /0, where 0 is the initial number density. If simulation starts from an empty case, i.e., 0 = 0, then the simulator is initialized with particles of uniform size (the real value of the size is immaterial to the simulation results), and the simulator size is set to a huge value, say = 1010, which renders a tiny particle number density to approximate the condition 0 = 0. If simulation starts from a case with a specific PSD, then the size of the simulation particles is assigned by a stochastic process to satisfy the initial PSD. There are various convenient ways to generate random number to satisfy a given distribution [37], 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 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 threshold value (which is set beforehand to balance the stochastic error and numerical cost), then down‐ sampling is performed, i.e., 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 in **Figure 1**) shows how to implement the stochastic algorithm of Gillespie [25]. Updating coagulation kernel is to calculate the pair coagulation rate <sup>=</sup> ( , )/, ( = 1,…, 1, = + 1,…,). The random coagulation time is generated according to Eq. (13). The com‐ parison statement sum <sup>+</sup> <sup>&</sup>lt; is to judge whether the time for two particles to coagulate is still admissible within the discrete time step . Here sum is the accumulated coagula‐ tion 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 to perform coagula‐ tion, 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.
