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 bounce-off unlikely.
The rate of collision (Zij) between particles of size υi and υj is given by Zij = βij ni nj, where ni is the number of particles of size υi-1 < υ < υ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 = ΒN2/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 general-purpose 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 chain-like combustion aerosols due to Brownian motion. They used the so-called
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 gas-phase 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 Si = {x1, x2,...,xp} as the set of variables that fully describe the state of a system at time ti, a Markov chain is defined as the finite set of states M{S1, S2,...,Sn} through which the system evolves. At each of the discrete sequences of time ti, the state Si determines a set of conditional probabilities qi1, qi2,..., qim, where qij 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 qij = 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., qij ≠ qij+1 ≠ qmn, 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 ti is denoted by ni1, the number of doublets (particles made of two singlets) by ni2, and so on, the state of the aerosol at the ith time step (collision) is defined by the set Si = {ni1, ni2,...,nis}. Imposing the condition of irreversible collisions, the transition from the. initial state S1 = {5,0,0,0,0} to the final state S7 = {0,0,0,0,l} may proceed along any of the paths shown below.
Every state Si 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 t1, only collisions among singlets are possible. Then q11 = 0 for j ≠ 1. Therefore, the aerosol evolves to the state S2= {3, 1, 0, 0, 0}. At time t2, collisions type 1 (among singlets) and type 2 (singles with doublets) will happen with probability q21 and q22, respectively. In general, at time ti a collision type j (particles made of p-singlets with particles made of q-singlets) will occur with probability qij. The transition probability qij is proportional to the collision frequency Zpq = βpq nip niq. With the restriction:
then, the transition probability is given by:
The objective of the Markov-Monte Carlo (Markov-MC) technique is to determine the most likely state at a given time ti by simulating the transitions S1 Si, using random or statistical sampling of transitions. At every transition, the sampling of transitions can be done by, for example, aligning the collision probabilities qij on a 0–1 interval in an arbitrary order and then generating a random number Rn such that 0 < Rn < l. Rn identifies uniquely the types of collisions that have occurred, and consequently, the types of transitions that have happened.
Physically, the Markov-MC 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 well-known 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 combustion-generated.aerosols, choose the quantity (No/N-1) 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 Markov-chain version of the Monte Carlo technique (i.e., the Markov-MC method), the average time between events is computed through the collision frequency function Zij. Here, an event is defined as a collision among particles of any type. In a unit time and unit volume, there are Z11 collisions among partic1es of size 1 and 1 occurring simultaneously with Zij 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 time-dependent 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 Markov-MC 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 Markov-MC 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 (ni) is stored in a 1D vector (array). ni is the number of particles per unit volume of size:
where the actual particle size r¡ is given by:
and Mi 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 Markov-MC method, the sampling of transitions between states is done by aligning the collision probabilities qij’s on a 0-1 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 qij with the respective type of collision (i-j).
A 1D vector
where
3.5. General algorithm
The overall algorithm to perform the Markov-MC simulation is described in the flow chart shown in Figure 4. Initially, the type of discretization, the starting particle size distribution ni 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 Markov-MC 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 Markov-MC simulated and compared to the corresponding analytical solution to test the accuracy of the method. Later, the restriction β = constant will be removed, and the Markov-MC 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 ri = rj, 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 No 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 ni we proceed inductively, i.e., first an expression for n1 is obtained, then this result is used to obtain an expression for n2, and so on. Doing so, the general solution for the concentration of particles of size ri is:
Evolution of an initially monodisperse aerosol with constant collision frequency function were Markov-MC simulated. Results are shown in Figure 5 along with the analytical solution for N, n1/No, n2/No, and n3/No.The agreement is excellent, demonstrating the validity of the Markov-MC 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 Markov-MC 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] | 1x106 | 1.4954 | 447.1 | |
1x10 [13] | 1x105 | 0.3904 | 260.9 | |
1x10 [14] | 1x104 | 0.2713 | 6.8 | |
1x10 [15] | 1x103 | 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 106 parcels and increasing the resolution every 0.9 x l06 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 l05 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 Markov-MC 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 Markov-MC 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 Markov-MC results obtained for the evolution of an aerosol for three different random number sequences. Figure 7.b shows that the Markov-MC 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 Markov-MC 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 self-preserving form after long periods of time, and appears to be independent of the initial size distribution. The self-preserving distribution for the case of β = constant is exponential in form, whereas the self-preserving 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 self-preserving form develops after about 3 dimensionless coagulation times (3τcoll) and that the shape of the asymptotic distribution varies with the Knudsen number (Kn) but not with the initial conditions. Later, Lee et al. [74] assumed a lognormal type of self-preserving 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 Markov-MC 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 No = l x l0 [18] water particles/m3 of l nm in radii.
Figure 8 shows the evolution of the aerosol as predicted by the Markov-MC method. The results are presented in terms of nondimensional particle size distribution (ni/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 self-preserving distribution that seems to be lognormal.
To determine the type of self-preserving distribution predicted by the Markov-MC method, the obtained Markov-MC distributions were fitted with different exponential distribution functions. The goodness-of-fit 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 Markov-Monte Carlo (Markov-MC) 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 Markov-MC method appear when condensation processes are added into the model. This is the subject of the next chapter.