Expressions for the collision frequency function in the free molecular regime (FMR) and continuum regime (CR). [64]
Coagulation describes the process by which particles that are suspended in a fluid come into contact and coalesce or adhere to one another. Particles interact because of Brownian (thermal) motion, hydrodynamics, and/or gravitational, Van der Waals, Coulomb, or other forces. During this work, only coagulation due to Brownian motion is considered as this is the dominant mechanism for particle interactions in aerosols. Coagulation due to Brownian motion has been extensively studied and expressions for the rate of coagulation are well known. Furthermore, these expressions have been validated experimentally. [64]
During the development of these expressions it is assumed that upon collision, particles stick together and form a third single spherical particle whose volume is equal to the sum of the original two, i.e., the assumptions are:
Spherical particles. This assumption implies that the rate of sintering or coalescence is instantaneous compared to the rate of collision. This is a good assumption for very small particles, as shown in Appendix B.
Accommodation or sticking coefficient of unity. Although little is known quantitatively about the accommodation coefficient of aerosol particles, their low kinetic energies make bounceoff unlikely.
The rate of collision (Z_{ij}) between particles of size υ_{i} and υ_{j} is given by Z_{ij} = β_{ij} n_{i} n_{j}, where n_{i} is the number of particles of size υ_{i1} < υ < υ_{i} per unit volume, and β_{ij} is the collision frequency function. Table 1 gives expression for β_{ij} in the free molecular regime (FMR) and in the continuum regime (CR), when coagulation is due to Brownian motion.
For a monodisperse aerosol with N particles per unit volume, the collision frequency simplifies to Z = ΒN^{2}/2. Then, the characteristic time for coagulation τ_{coll} = N/Z = 2/βN, where τ_{coll} is the average elapsed time between successive collisions for a particle of size υ. It can also be interpreted as the time required to halve the particle concentration of an initially monodisperse aerosol by coagulation, or the time to double the volume of the particles. Table 1 shows expressions for τ_{coll}.






FRM  

CR  

The theory of coagulation is essentially a scheme for keeping account of particle collision as a function of particle size. The rate of formation of particles of size k by collision of particles of size i and j is given by ½
which is the dynamic equation for the discrete size spectrum when only coagulation is important. The solution of Equation (1) depends on the form of β_{ij}, which is determined by the kinetics of particle collision. In this chapter, the Monte Carlo method developed to solve Equation (1) is described. Initially, a general description of the Monte Carlo technique is presented and then the method is applied to solve the coagulation equation. Later, the model is validated by comparing results with analytical and numerical solutions and experimental observations.
1. Monte Carlo simulation
The Monte Carlo technique is a method of statistically sampling events to determine the average behavior of a system. [65] The name Monte Carlo was first applied to a class of mathematical methods by scientists working on the development of nuclear weapons in Los Alamos in the 1940s. The essence of those methods was the invention of games of chance whose behavior and outcome were used to study statistical and deterministic phenomena. Since then, the term Monte Carlo has been employed to describe a wide range of procedures involving the notion of sampling.
Historically, the first example of a computation by a Monte Carlo method was reported by Buffon in 1777, who described an unexpected method to calculate the value of π. [66] With the increasing availability of very high speed generalpurpose computers, many problems have become treatable and the Monte Carlo method has increased in popularity. The technique has been applied in many disciplines to simulate very diverse systems. The greatest successes of the Monte Carlo method have arisen in areas where the problem itself consists of some random process. It has been applied to simulate problems like chemical kinetics, diffusion, radiation, turbulence, and air pollution.
The solutions obtained through a Monte Carlo simulation are statistical in nature and subject to the laws of chance. This aspect is a drawback, but not a serious one. Convenience, ease, directness, and expressiveness of the method are important assets. Other general features of the method are:
Cyclic nature of the programming
Necessity of performing a large series of computations of a uniform type
Use of a comparatively small amount of memory for storage of intermediate results
Smoothing of errors
A distinction is sometimes made between
The
With the
1.1. Monte Carlo simulation of aerosol processes
The Monte Carlo method has also been used in the study of aerosol dynamics. Kaplan et al. [68] studied the agglomeration of chainlike combustion aerosols due to Brownian motion. They used the socalled
The basic algorithm is as follows: At the beginning of each simulation, monodisperse or lognormally distributed clusters are randomly distributed on the comers of a 2D or 3D lattice. During each time step, each cluster moves at random toward any of the 4 (2D) or 6 (3D) possible directions on the lattice. The probability of movement of a given cluster is dependent upon cluster size as predicted by the diffusion theory; smaller clusters have a proportionately higher probability of diffusion compared to large ones. The probabilities of movement are scaled by assuming the probability of movement of a single mer at, for example, 2000 K to be l. Clusters located on adjacent nodes have a probability of collision that is dependent upon the form of the collision frequency function. At the end of each time step, after each particle has been displaced, and all possible collisions have been checked, the particle size distribution is updated, and the simulation continues until the required time duration is reached. This model accounts only for coagulation, assumes cluster–cluster collisions to be irreversible, and neglects sintering.
Akhtar et al. [69]advanced this work and performed a Monte Carlo simulation of the gasphase coagulation and sintering of nanosized particles. Particle coagulation and growth were simulated on a 500 x 500 square lattice with periodic boundary conditions to reduce edge effects. During the sintering routine, a particle on any of the perimeter sites of the aggregate is picked at random and allowed to move to any neighboring unoccupied perimeter site on the same aggregate. The probability of acceptance of such movement is determined by the total potential of the new configuration. Configurations with lower total potential are more likely to be accepted. The number of times that a particle is selected (i.e., number of sintering steps) per coagulating step are determined by the ratio of sintering to growth rates.
This method has been successful in providing insight in the formation of agglomerates. However, the possibility of extending it to include coagulation and condensation processes is limited by the weakness of the method in computing the physical time, [69] which is critical in the case of processes occurring at two different timescales like the case of the encapsulation process where coagulation and condensation occurs simultaneously but at different rates.
Instead, coagulation can be simulated through a more straightforward and physical Monte Carlo approach developed by Husar, [55] based on the work of Schaad [67] and Sutherland. [70] In these approaches, position and velocity space are omitted and only collisions are simulated, not in the physical sense but as a transition in a Markov chain. The following section illustrates the concept of Markov chains, and the use of it to simulate coagulation processes.
1.2. The Markov process
A Markov process is a chain of events occurring in sequence with the condition that the probability of each subsequent event in the chain is not influenced by prior events. [71] The usual example is a drunk gentleman who begins a walk through a strange city. At each street comer that he reaches, he continues his walk by choosing completely at random one of the streets leading from the intersection. The history of this random walk is then a Markov chain, because his decision at any point is not influenced by where he has been.
Defining S_{i} = {x_{1}, x_{2},...,x_{p}} as the set of variables that fully describe the state of a system at time t_{i}, a Markov chain is defined as the finite set of states M{S_{1}, S_{2},...,S_{n}} through which the system evolves. At each of the discrete sequences of time t_{i}, the state S_{i} determines a set of conditional probabilities q_{i1}, q_{i2},..., q_{im}, where q_{ij} is the probability that the system which is in state Si, will be in state
In the example of the drunk gentleman in the strange city, the state or position of the gentleman is described by his {x, y} coordinate in the city at every discrete time. The collection of comers where he has been defines the Markov chain. Because of the randomness of his choice at each intersection, the transition probability at each comer is q_{ij} = q = l/4. However, if his choice were determined by his “feeling” of safety, affected for example by the luminosity of the street, or the number of people on the street, the transition probabilities would be different at each state (corner), i.e., q_{ij} ≠ q_{ij+1} ≠ q_{mn}, but his choice would be still unaffected by the places where he has been.
Because of the nature of the problem, it might be possible to simulate a sample walk by playing a roulette wheel with 100 sections. Before every spin, the roulette wheel is divided in 4 sectors, each corresponding to a possible direction, such that each sector has a number of uniformly spaced sections proportional to the sense of safety that the drunk gentleman feels in that specific direction. Starting at a given position, the gentleman’s most likely path and most likely position after certain time could be found by simulating a large number of histories, i.e., by running through the process many times. Another example of a Markov process is the coagulation of particles in aerosols.
2. Coagulation as a Markov process
Consider the evolution of an aerosol consisting initially of 5 particles of the same size (singlets). lf the number of singlets in the system at time t_{i} is denoted by n_{i1}, the number of doublets (particles made of two singlets) by n_{i2}, and so on, the state of the aerosol at the i^{th} time step (collision) is defined by the set Si = {n_{i1}, n_{i2},...,n_{is}}. Imposing the condition of irreversible collisions, the transition from the. initial state S_{1} = {5,0,0,0,0} to the final state S_{7} = {0,0,0,0,l} may proceed along any of the paths shown below.
Every state S_{i} defines a set of possible types of collisions, each with a probability that is completely uninfluenced by the history of the aerosol. Therefore, the evolution of the aerosol is properly defined as a Markov process.
At time t_{1}, only collisions among singlets are possible. Then q_{11} = 0 for j ≠ 1. Therefore, the aerosol evolves to the state S_{2}= {3, 1, 0, 0, 0}. At time t_{2}, collisions type 1 (among singlets) and type 2 (singles with doublets) will happen with probability q_{21} and q_{22}, respectively. In general, at time t_{i} a collision type j (particles made of psinglets with particles made of qsinglets) will occur with probability q_{ij.} The transition probability q_{ij} is proportional to the collision frequency Z_{pq} = β_{pq} n_{ip} n_{iq}. With the restriction:
then, the transition probability is given by:
The objective of the MarkovMonte Carlo (MarkovMC) technique is to determine the most likely state at a given time t_{i} by simulating the transitions S1 S_{i}, using random or statistical sampling of transitions. At every transition, the sampling of transitions can be done by, for example, aligning the collision probabilities q_{ij} on a 0–1 interval in an arbitrary order and then generating a random number R_{n} such that 0 < R_{n} < l. Rn identifies uniquely the types of collisions that have occurred, and consequently, the types of transitions that have happened.
Physically, the MarkovMC model simulates particle coagulation by reproducing the following algorithm. Time is divided into time steps. During each time step, the occurrence of every collision is registered and classified according to its type. When the information is normalized by Z, the total number of collisions that have occurred, it provides the probability of occurrence of each type of collision at the given time. This sampling process is Monte Carlo simulated by aligning the probabilities on a O–1 line and generating Z random numbers, each one identifying a type of collision. In other words, by generating Z random numbers and correlating that number to the type of collision being statistically chosen, the collisions that have occurred during the time step are reproduced. Assuming that the time step is short enough, the change in the particle size distribution due to the occurrence of the Z collisions is negligible, and the change can be implemented at the end of the time step with minor error. By decreasing the time step, the error decreases. In the limit, the time step is chosen such that the particle size distribution is updated after every collision.
2.1. Computation of time
Computation of the physical time elapsed between events (i.e., the time step) has been a serious drawback of different Monte Carlo approaches. Most of the past studies have been performed where quantification of time can be avoided. However, when coagulation and condensation occur simultaneously, computation of time is a critical issue.
Schaad [67] proposed to use wellknown solutions to calibrate the model with respect to time and then use those settings for the particular case of interest. Using the same approach, Kaplan and Gentry, [68] working with a Monte Carlo simulation to study agglomeration of nonspherical chainlike agglomerates in combustiongenerated.aerosols, choose the quantity (N_{o}/N1) as a measure by which to represent the rate of cluster growth, since coagulation theory for constant collision frequency functions predicts that this quantity is directly proportional to time. Clearly, these approaches offer reasonable approximations but are not exact.
With the Markovchain version of the Monte Carlo technique (i.e., the MarkovMC method), the average time between events is computed through the collision frequency function Z_{ij}. Here, an event is defined as a collision among particles of any type. In a unit time and unit volume, there are Z_{11} collisions among partic1es of size 1 and 1 occurring simultaneously with Z_{ij} collisions among particles of size i and j. In total, there are Z =
This implies that in a unit volume, on average there is one collision every 1/Z unit times. Then the average elapsed time between events is given by: [55]
If it were desired to simulate the process more rigorously the time increment should be selected from a
The coagulation equation is a timedependent mass conservation expression for aerosols. Then a numerical solution of this nonlinear integro differential equation can be obtained by updating the particle size distribution function after every MarkovMC simulated collision and by computing the elapsed time through Equation (4). The implementation of this method on a computer code is straightforward. The following section discusses the specific computational aspects involved during its implementation.
3. Description of the code
3.1. Discretization of the particle size range
The MarkovMC algorithm starts by dividing into


Discrete  
Logarithmic  
Geometric  
Linear  
3.2. The concept of parcels
In actual combustion aerosols the particle concentration is typically greater than l x l0 [18] particles/m3. To directly simulate that many collisions would be impracticable, even with advanced supercomputers. To address this issue the concept of parcel was developed. A parcel is defined as a bundle of
3.3. Particle Concentration
Employing the concept of parcels, the discrete version of the particle size distribution (n_{i}) is stored in a 1D vector (array). n_{i} is the number of particles per unit volume of size:
where the actual particle size r¡ is given by:
and M_{i} is the total mass of the particles in section i. Notice that by using Equation (6), it has been assumed that all the particles within the section have the same size. Also notice that by allowing the size of the particles within a section to dynamically assume their mass mean value, instead of assuming that it remains constant, the total mass of the aerosol is conserved during the simulation of the evolution of the aerosol.
3.4. Aligning of the transition probabilities
In the MarkovMC method, the sampling of transitions between states is done by aligning the collision probabilities q^{ij}’s on a 01 interval in an arbitrary order. The imp1ementation of this step is critical because in addition to align the transition probabi1ities, it should provide an easy and quick way to relate every transition probability q_{ij} with the respective type of collision (ij).
A 1D vector
where
3.5. General algorithm
The overall algorithm to perform the MarkovMC simulation is described in the flow chart shown in Figure 4. Initially, the type of discretization, the starting particle size distribution n_{i} and the stopping condition (e.g., time for simulation) are specified as inputs to the program. Then, the computation of the collision coefficient β_{ij} for all the possible types of collisions is performed. This information is stored in an
4. Validation
Being a probabilistic approach, the calculation of the absolute accuracy of the MarkovMC method is generally impossible. However, its validity can be judged by, for example, comparing results with those obtained by other methods or a model problem whose solution is known accurately. Consequently, in this section the evolution of an initially monodisperse aerosol, with collision frequency function independent of particle size (β_{ij} = β), is MarkovMC simulated and compared to the corresponding analytical solution to test the accuracy of the method. Later, the restriction β = constant will be removed, and the MarkovMC results will be compared with other numerical results and experimental observations.
4.1. Evolution of an initially monodisperse aerosol with β = constant
4.1.1. The Smoluchowski’s solution
One of the few analytical solutions to the GDE can be obtained for the case of aerosols subjected to the following restrictions:
Continuum regime
Collision frequency function β independent of particle size
Initially monodisperse size distribution
With these assumptions and setting r_{i} = r_{j}, the collision frequency function reduces to:
Substituting this expression into Equation (1), the coagulation equation reduces to:
Summing over all values of
Integrating once,
where N_{o} is the initial total number of particle per unit vo1ume, and τ_{coll} is the characteristic coagulation time given by:
Physically, the characteristic time for coagulation is the time needed for an initial monodisperse aerosol to reduce the total particle concentration to half of its initial value. To obtain an expression for n_{i} we proceed inductively, i.e., first an expression for n_{1} is obtained, then this result is used to obtain an expression for n_{2}, and so on. Doing so, the general solution for the concentration of particles of size r_{i} is:
Evolution of an initially monodisperse aerosol with constant collision frequency function were MarkovMC simulated. Results are shown in Figure 5 along with the analytical solution for N, n_{1}/N_{o}, n_{2}/N_{o}, and n_{3}/N_{o}.The agreement is excellent, demonstrating the validity of the MarkovMC model and the timescaling.
Figure 6 shows the error in the total number of particles per unit volume at different stages of the evolution. The error in all cases is below 1% for t/τ_{coll}<10, showing a good agreement between the analytical solution and the MarkovMC simulation. The error arises from the computation of the time steps through Equation (4). This error can be minimized up to any degree by increasing the statistics of the simulation, i.e., by increasing the number of parcels
The increase in error with time results because with every event (collision), the number of parcels available for simulation decreases by one. Consequently, the statistics of the simulation decreases with the number of events performed and, thus, the error in the computation of the time steps increases. The cumulative error in time can be maintained below a certain level without increasing the use of computational resources by dividing the entire simulation into several sequential runs. A run is stopped before the cumulative error in time is excessive. Before the next subsequent run, the number of parcels in each section is multiplied by 1O and the scale factor
To determine the appropriate initial number of parcel





1x10 [12]  1x10^{6}  1.4954  447.1  
1x10 [13]  1x10^{5}  0.3904  260.9  
1x10 [14]  1x10^{4}  0.2713  6.8  
1x10 [15]  1x10^{3}  0.2224  3.2 
In Table 3, the Central Processing Unit (CPU) time is the amount of system time required to run the job, i.e., the number of seconds devoted to system operations required for the program. The user CPU time is the amount of time that the program requires of the CPU.
Figure 6 and Table 3 show that working with Charney, the balance between accuracy and use of computational resources is obtained by starting with 1 x 10^{6} parcels and increasing the resolution every 0.9 x l0^{6} events. However, the required number of starting parcels depends on the machine used and the number of sections actively being used. For example, early stages of the evolution of the monodisperse aerosol require less computational operations than a fully developed one, since at the early stages most of the sections are empty, and the number of operations per event is proportional to the number of sections being actively used. On the other hand, working with slower machines, resolution needs to be sacrificed to obtain results within a reasonable amount of time. For example, to obtain results within an hour, the number of starting parcels need to be reduced to l x l0^{5} for
4.1.2. Effect of the types of sectionalization and random number sequence
Several runs were made to observe the effect of different types of sectionalization on the Markov MC results. For the same conditions, discrete, geometric, linear, and logarithmic sectionalization was used. In all the cases, the results were the same. The possibility of using different types of sectionalization according to the problem being analyzed constitutes another advantage of the MarkovMC method over similar codes like the sectional method, which was designed to work on a logarithm type of sectionalization. A similar type of test was performed to see the effect of the number of sections being used on the evolution of a coagulating aerosol.
The evolution of an initially normal distributed aerosol was observed using a logarithmic type of sectionalization. Several runs were made, varying the number of section per decade being used. Figure 7.a shows the results obtained. This figure shows that the number of sections per decade can be reduced up to 18 and the evolution of the aerosol as predicted by the MarkovMC model is basically the same.
Similarly, the effect of the random number sequence being used was studied. C++, as any other programming language, provides a standard library function to generate pseudo random numbers with a uniform distribution. The sequence of random numbers produced by this function is determined by a seed number. Varying this seed number, different sequences can be obtained. Figure 7.b shows the MarkovMC results obtained for the evolution of an aerosol for three different random number sequences. Figure 7.b shows that the MarkovMC results are independent of the sequence of random numbers used.
4.2. Comparison with numerical and experimental results
In general, the collision frequency function β is not a constant but depends on the sizes of the colliding particles and the nature of the aerosol. Under these conditions, there is not an exact analytical solution to the coagulation equation. Many numerical and experimental works have been performed to describe the behavior of coagulating aerosols. The main results are summarized here.
Later, the MarkovMC results obtained for an atmospheric aerosol with β(r) are qualitatively compared with these numerical and experimental results.
Friedlander et al. [52] stated that the particle size distribution in a coagulating aerosol approaches an asymptotic selfpreserving form after long periods of time, and appears to be independent of the initial size distribution. The selfpreserving distribution for the case of β = constant is exponential in form, whereas the selfpreserving distribution for Brownian coagulation appears to be lognormal. Thus, for Brownian coagulating aerosols, the state of the aerosol at any time could be characterized by only two parameters, the mean particle size and the geometric standard deviation σ_{g.}
The theory was later supported by Hidy [73] who through numerical experiments concluded that the selfpreserving form develops after about 3 dimensionless coagulation times (3τ_{coll}) and that the shape of the asymptotic distribution varies with the Knudsen number (K_{n}) but not with the initial conditions. Later, Lee et al. [74] assumed a lognormal type of selfpreserving distribution and found an approximate analytical solution for the coagulation equation by using the method of moments. The solution predicts an asymptotic value of σ_{g} = 1.38 in the free molecular regime. The solution also estimates the time to reach the asymptotic distribution in the range from 0 to 473τ_{coll} depending on the initial conditions.
Matsoukas and Friedlander, [75] studying the dynamics of the formation of metal oxides (MgO and ZnO) in a flat flame where salts of magnesium and zinc were introduced in the form of a dry aerosol, measured the particle size distribution of samples collected at different heights above the flame. They observed that the distribution collapses into a single curve when plotted in normalized form and concluded that the normalized distribution is, to a good approximation, lognormal with a geometric standard deviation σ_{g}
To validate the MarkovMC method, results for a coagulating atmospheric aerosol (i.e., water in air at room temperature) with β(r) were obtained and compared with these numerical and experimental results. Initially, the aerosol consists of N_{o} = l x l0 [18] water particles/m^{3} of l nm in radii.
Figure 8 shows the evolution of the aerosol as predicted by the MarkovMC method. The results are presented in terms of nondimensional particle size distribution (n_{i}/N) as a function of the nondimensional time (t/τ_{coll}). The characteristic time for collision at the initial conditions is used to nondimensionalize time. Figure 8 shows that the particle size distribution gradually broadens out with time and evolves toward a selfpreserving distribution that seems to be lognormal.
To determine the type of selfpreserving distribution predicted by the MarkovMC method, the obtained MarkovMC distributions were fitted with different exponential distribution functions. The goodnessoffit was evaluated by the
5. Overview
A Monte Carlo method has been developed to simulate particle collisions in aerosols. Collisions are simulated as transitions in a Markov chain. The MarkovMonte Carlo (MarkovMC) method has been validated by comparing results with analytical and numerical results and with experimental observations. The method can handle any number of sections and any type of sectionalization. Conditions to obtain good accuracy with minimum use of computational resources were determined. The convenience, ease, directness, and expressiveness are some of the advantages of the method. However, the unique advantages of the MarkovMC method appear when condensation processes are added into the model. This is the subject of the next chapter.