- Email: [email protected]
Brownian dynamic simulation for the aggregation of charged particles
Aerosol Science 32 (2001) 1369–1388
www.elsevier.com/locate/jaerosci
Brownian dynamic simulation for the aggregation of charged particles Hyungho Parka , Sangsoo Kima ;∗ , Hyuksang Changb a
Aerosol & Particle Technology Laboratory, Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, Kusong-Dong, Yusong-Gu, Taejon 305-701, South Korea b Environmental Aerosol Engineering Laboratory, Department of Environmental Engineering, Yeungnam University, Kyungsan 712-749, South Korea Received 17 October 2000; accepted 19 March 2001
Abstract A numerical technique based on Brownian dynamic simulation for the aggregation of charged particles in the free molecular regime is presented. The Langevin equation is used for tracking each particle making up an aggregate. A periodic boundary condition is used for calculation of the aggregation process in each cell with 500 primary particles of 16 nm in diameter. Particle motion is based on the thermal force and the electrostatic force. The electrostatic force on a particle in the simulation cell is considered as a sum of electrostatic forces from other particles in the original cell and its replicate cells. We assume that the electric charges accumulated on an aggregate are located on its center of mass, and aggregates are only charged with pre-charged primary particles. The morphological shape of aggregates is described in terms of the fractal dimension. The fractal dimension for the uncharged aggregate was Df = 1:761, and changed slightly for the various amounts of bipolar charge. However, in case of unipolar charge, the fractal dimension decreased from 1.641 to 1.537 with the increase of the average number of charges on the particles from 0.2 to 0.3 in initial states. During the early and middle stages of aggregation process, the average aggregate size in the bipolar charge state was larger than in the uncharged state, but was almost equal in the ;nal stage. On the other hand, in the unipolar charge state, the average size of c 2001 Elsevier an aggregate and the dispersion of particle volume decreased with increasing charge. Science Ltd. All rights reserved.
1. Introduction Nanoparticles have many properties that di>er from the corresponding bulk material, and these properties o>er many attractive applications such as in electric, optic and magnetic ∗
Corresponding author. Tel.: +82-42-869-3002; fax: +82-42-869-3210. E-mail address: [email protected] (S. Kim). c 2001 Elsevier Science Ltd. All rights reserved. 0021-8502/01/$ - see front matter PII: S 0 0 2 1 - 8 5 0 2 ( 0 1 ) 0 0 0 6 3 - 5
1370
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
Nomenclature d0 Df F kb ke L mp m0 n N0 N qi ; qj r r0 Rg t T u vavg v0 X
diameter of primary particles, m mass fractal dimension, dimensionless electrostatic force, kg m s−2 Boltzmann constant, kg m2 s−2 K −1 Coulomb constant, kg m3 s−2 C−2 length of side of simulation cell, m mass of aggregate, kg mass of primary particle, kg number of primary particles in aggregate, dimensionless total number of primary particles in simulation cell, dimensionless total number of aggregates in simulation cell, dimensionless charge quantities on ith and jth aggregates, C displacement of aggregate, m center of mass for aggregate, m radius of gyration for aggregate, m time, s gas temperature, K aggregate velocity, m s−1 volume of aggregate, m3 volume of primary particle, m3 random force due to collision between gas and particle, kg m s−2
Greek letters friction constant, s−1 electrostatic potential energy, kg m2 s−2 volume fraction of particles, dimensionless particle density, kg m−3 p unit time, s 1 dimensionless time, dimensionless Subscripts k components of rectangular coordinates, k = 1; 2; and 3
products (Kruis, Fissan, & Peled, 1998). The synthesis of nanoparticles is mostly done in a gas phase process. Especially, for making TiO2 or SiO2 particles, the Game process is adopted as an eHcient tool in many applications. In the Game synthesis of nanoparticles, the particles are produced by nucleation or chemical conversion, and growing by condensation or coagulation. These processes are very common to many aerosol processes. In case of thermally stable particles, the coagulation should be considered as a dominant mechanism in the particle growth
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1371
(Pratsinis, 1988; Hurd & Flower, 1988). Therefore, the thermally stable particles in the Game grow from spherical primary particles to aggregates. The aggregation characteristics are mainly a>ected by the morphology of aggregates (Julien & Bodet, 1987; Matsoukas & Friedlander, 1991; Vemury & Pratsinis, 1995). Because a fractal dimension can quantify the morphology of aggregates, the fractal dimension becomes a key variable for analyzing the aerosol dynamics. For the analysis of the aggregation process, two methods have been applied. One is the fractal modi;ed general dynamic equation and the other is the statistical approach. The former is more eHcient than the latter in saving CPU time in the calculation. The application of the fractal modi;ed general dynamic equation is limited, as the fractal dimension should be known in the calculation process. Therefore, the statistical approach is the only choice for the case of unknown fractal dimension. Some researchers have developed numerical simulation techniques to characterize the growth of aggregates. Using theoretical simulation, Meakin and Wasserman (1984) have obtained almost the same fractal dimension compared with the results of previous experiments. However, their method is not proper as an investigation tool for aerosol dynamics because of its size invariant assumption for the di>usion of a particle. Mountain, Mulholland, and Baum (1986) simulated the aggregation process of the aerosol by using Brownian dynamics. They obtained exactly not only the fractal dimension but also the size of aggregates. With the increase of new applications for nanoparticles recently, it has been important to control the size and the morphology of nanoparticles in each nanoparticle synthesis process. Among the numbers of control variables, the electrostatic force has been adopted as a means for precise control of the aggregation. Xiong, Pratsinis, and Mastrangelo (1992) simulated the growth of Game aerosols with the sectional aerosol model. In their results, if the aerosols were in a unipolar charge state, the particle size distribution was more uniform and the average size was smaller than that of an uncharged state. However, their results were obtained using the assumption of spherical particles. Therefore, a new simulation model is needed for explanation in realistic cases. In this study, we wanted to quantify the e>ect of electrostatic force on the morphology and growth of charged aggregates. Therefore, we introduced a new simulation model that could describe the electrostatic force between charged particles of in;nite number. This electrostatic simulation model was coupled with the Brownian dynamic simulation to deal with the aggregation of charged particles. We studied about not only fractal aggregates, but also spherical aggregates to show the e>ect of morphology in the aggregation of charged particles. On the other hand, we treated only a simple case of point charges located on the center of an aggregate to avoid some complexities of the calculation, although this method could be extended to some cases, where the charges were distributed in the volume or on the surface of an aggregate. 2. Theory 2.1. Power law The morphological state of the aggregates is de;ned as the fractal number by Forrest and Witten (1979) who found that the fractal dimensions of the Game aerosol are 1.7–1.9 in the
1372
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
synthesis of iron, zinc, and silicon oxide particles. They con;rmed that there was a power law relation between the characteristic length and the number of primary particles in each aggregate. Usually the radius of gyration of the aggregate is used as the characteristic length for the aggregate. The power law relation is expressed as f n ∼ RD g ;
(1)
where n is the number of primary particles in an aggregate, Rg is the radius of gyration, R2g = (1=n)(rj − r0 )2 ; and r0 is the center of mass de;ned as r0 = (1=n)rj . 2.2. Langevin equation Brownian dynamics describes the dynamic behavior of a particle whose mass and size are larger than those of the gas molecules they are immersed in. These particles are subject to stochastic collisions with the gas molecules, which lead to the random motion of the particles. It is assumed that a particle is in thermal equilibrium with a surrounding gas at the temperature of T and is subject to an external force F. In this case, the translational motion of a particle is described by the Langevin equation as d (mp uk ) = Fk − mp uk + Xk : dt
(2)
In Eq. (2), the LHS represents the inertia force of particle of mass mp , and the ;rst term on the RHS is systematic external force acting on the particle. This term includes van der Waals force, electrostatic force, magnetic force, image force, gravity, etc. The other forces except the electrostatic force can be negligible in the aggregation of charged nanoparticles because the magnitude of the other forces are much smaller than that of the electrostatic force. The second term in the RHS is the frictional drag of the Guid around the particle. is the friction constant, and its inverse −1 represents the particle relaxation time. In the course of its motion, the particle displaces the same length of its mean free path with the mean thermal velocity of the particle between each collision. The last term, Xk , is a random acceleration force caused by the random collisions between the particle and the surrounding gas molecules. The random force Xk is represented with the Gaussian random distribution function of which the mean and square mean are represented as Xk = 0;
Xk2 = 2mp kb T:
(3)
2.3. Solution of the Langevin equation Several methods have been used to integrate the Langevin equation (Turq, Lantelme, & Friedman, 1977). In this study, we use the method derived by Ermak and Buckholz (1980). The analytic solution of the Langevin equation is based on the assumption that a single particle is in a constant force ;eld within the integrating time step. The velocity and the displacement
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1373
are expressed as
F0; k (1 − e−t ) + B1 ; mp 1 F 2F0; k 1 − e−t + uk + u0; k − + 0; k t + B2 : −t mp mp 1+e
ukx = u0; k e−t +
(4)
rk = r0; k
(5)
As shown in Eqs. (4) and (5), the displacement is proportional to the new velocity, old velocity, old displacement and old force F0; k . B1 and B2 are random variables, which have the characteristic components of the Gaussian random distribution as B1 = 0;
kb T (1 − e−2t ); mp B2 = 0; 2kb T 1 − e−t 2 B2 = t − 2 ; mp 2 1 + e−t B12 =
B1 B2 = 0:
(6)
3. Numerical simulation 3.1. Particle charging In this simulation, we assumed that primary particles were charged with some additive ions. According to the comparison between the characteristic time for ion-particle and particle–particle collisions (Xiong et al., 1992), the collision time of ion-particle is shorter than that of particle– particle about four order of magnitude for 16 nm particles and some alkali ions. Therefore, it can be assumed that an aggregate is charged in the conditions of the pre-charged primary particles and a neutral gas. In this case, the charge quantity of an aggregate is the sum of elementary charges on the primary particles constituting the aggregate. Therefore, the number of charge on an aggregate is proportional to the volume of the aggregate. On the other hand, we set the value of initial charge on the primary particles using the results of Adachi, Kousaka, and Okuyama (1985). The charge quantity on the primary particle mainly can be obtained from the relation of the particle size, ion concentration, and resident time, etc. According to their results, small size particles, 16 nm in diameter, are charged partially at the given conditions. Therefore, we set average charge quantities to be maximum at 0.3 in the unipolar charged case. This value means that initially 30% of the primary particles individually carry a single charge. On the other hand, in the bipolar charging case, if the charge quantity is 0.3, it means that 15% of the particles carry a positive single charge and 15% of the particles carry a negative single charge. The bipolar charge quantity becomes large up to 3.0. Although this quantity seems an extreme case, we used this value to show some charging e>ects clearly.
1374
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
3.2. Initial condition To simulate the aggregation process, initial conditions for the velocity and the position of the primary particles have to be given. The positions of these particles are obtained with a uniform random number generator in 3-dimensions. On the other hand, it can be assumed that the primary particles experience only thermal force in the initial state. Therefore, we set the initial velocities with a normal random number generator that is characterized by the Gaussian distribution of zero mean and kb T=m0 variance. We calculated the trajectories of 500 primary particles in a simulation cell. The cell size was determined so that the particle volume fraction is 0.0087. Ten simulations for each initial condition were performed to obtain statistically meaningful average values (Mountain et al., 1986). The unit length was the same as the size of the primary particle d0 . The unit time, = (m0 d20 =kb T )1=2 , was ;xed as the time when a single sphere freely moved a distance equal to its diameter. In this work, we have ;xed the values of d0 , T , and p to be 16 nm, 1500 K and 2000 g m−3 , respectively. These conditions were obtained from the experimental conditions in a Game aerosol process. We performed the simulation for a small time step as 1=50 of the unit time to minimize the overlap between colliding particles. 3.3. Boundary condition The periodic boundary condition was used not only for maintaining a constant volume fraction but also for calculating the electrostatic force. The electrostatic force of the ith particle of the simulation cell is an in;nite summation of each force from the other charged particles in the original cell as well as in the replicated cells. Therefore, this force includes the electrostatic dispersion e>ect mentioned in the paper of Kasper (1980). The electrostatic force is expressed as Fi =
∞ j=1; j=i
ke
qi qj rˆij : rij2
(7)
We can obtain an accurate value when the force is calculated for suHciently many particles because this force is inversely proportional to the square of the distance. But as it is a very time consuming process, Brush, Sahlin, and Teller (1966) and Emark (1975) introduced a new technique for treating an in;nite number of pair potential. This potential includes the interaction with the images of the original cell in the replicated cells. Though in this procedure they could reduce the in;nite summation to a single summation over only the remaining ions in the original cell, their method can be applied only in a system of zero net charge. This constraint is a necessary one to get the ;nite potential. Therefore, it is impossible to apply this method directly to the case of a unipolar charged state. Fortunately, in this study, we need only the electrostatic force, not the whole potential. The force can be calculated directly from the gradient of the pair potential without accessing the pair potential. Though an in;nite number of particles a>ect the electrostatic force of a given particle, the force has a ;nite value because it is symmetrical in a macroscopic sense. In our simulation, electric charges locate in the center of mass of the aggregate and the dielectric constant of the
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1375
gas is uniform. By applying this method to the present system, the modi;ed electrostatic force by the in;nite number of particles can be represented as Fi = −∇i N * = − ke qi qj [∇ 1 (rij ) + ∇ 2 ( rij )];
(8)
j=i
where i is the potential energy of ith particle and N is the total number of aggregates in the original cell; qi and qj are the quantities of charges on the ith and jth aggregates. 1 and 2 are expressed as √ 1 (r) = erfc( r=L)=r − 1=L − Em =L; ∗
1 2( r ) = L * *
√ * * erfc( | r =L + l |)
l
*
*
| r =L + l |
*
+
*
exp(− l2 )cos(2 l · r =L)
l2
;
(9)
where L is the length of the side of cell; Em is the Madelung constant and its value, Em = * 2:83729479 for a simple cubic lattice (Kittel, 1967); l is a vector with integer components * and l is its magnitude; and ∗ denotes exclusion of the l = 0. To treat the equation of 2 , we have used the optimized expansion in Cubic harmonics, which was ;rst introduced by Hansen (1973). Consequently, their gradients are represented as √ 1 exp(− A) erfc( A) ∇ 1 (r)|k = − 2 2 xk eˆ k ; (10) + L A A3=2 * 1 1 ∇ 2 ( r )|k = − 2 # + xk + $xk3 + %xk7 eˆ k ; (11) L xk where #, , $, and % are described as #=−
2EFM ; A3
2 EGB − EJB 4EGB + 2 EFM − EKM + A2 A3 6EFM + 2 EDC − ECL 8EDC + + ; A4 A5
= 2 EH − EI +
$= −
4EG ; A2
%= −
8ED : A4
(12)
1376
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
The uppercase letters in Eq. (12) are de;ned as *
A = | r =L|2 = x12 + x22 + x32 ; B = x14 + x24 + x34 ; C = x18 + x28 + x38 ; D = d8 A4 + d10 A5 + d12 A6 ; E = exp(− A); F = c6 A3 + c8 A4 + c10 A5 + c12 A6 ; G = b4 A2 + b6 A3 + b8 A4 + b10 A5 ; H = a4 A2 + a6 A3 + a8 A4 + a10 A5 ; I = 4a4 A + 6a6 A2 + 8a8 A3 + 10a10 A4 ; J = 4b4 A + 6b6 A2 + 8b8 A3 + 10b10 A4 ; K = 6c6 A2 + 8c8 A3 + 10c10 A4 + 12c12 A5 ; L = 8d8 A3 + 10d10 A4 + 12d12 A5 ; M = x12 x22 x32 :
(13)
The coeHcients a1 ; : : : ; d12 are listed in Hansen (1973). Fig. 1 shows a comparison of the simple Coulomb force and the modi;ed electrostatic force. The values were computed along a diagonal axis. The simple Coulomb force was calculated just within only one cell. There is no di>erence if two particles are close together, but there is a relatively large di>erence as the distance between them becomes large. Because of the periodic boundary condition, the electrostatic force must be zero if two particles are located at the center and the edge of the simulation cell. The modi;ed electrostatic force satis;es this condition well. Therefore, our model for the electrostatic force of an in;nite number of charged particles can be regarded as a good model coinciding with the real system. 3.4. Calculation procedure We de;ne that an aggregate collides with the other if the distance between a particle in one aggregate and the other particle in other aggregate is located within the unit length, d0 . The momentum and the charge quantity are conserved between before and after the collision, and immediately the new-formed aggregate is in thermal equilibrium state. We assumed that the
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1377
Fig. 1. Comparison between modi;ed electrostatic force versus simple Coulomb force along the diagonal axis.
sticking probability was 1 after each collision. The aggregates grow with the repetition of the above sequence. The simulation procedure is summarized in Fig. 2. 4. Results 4.1. Sphere model of Brownian dynamic simulation A comparison between well-known general dynamic equation models and the Brownian dynamic simulation was performed to verify our simulation technique. It was assumed that colliding spherical particles coalesce instantaneously to form larger spherical particles. We de;ned this model as the sphere model. As shown in Fig. 3, in the early stage of aggregation a reduced number concentration is similar to that of a monodisperse aerosol. However, it becomes similar to the self-preserving case in the ;nal stage, where the reduced number concentration is directly related to the average volume of aggregates as N0 =N (t) − 1 = vavg =v0 − 1:
(14)
v0 is the volume of the primary particle 1 , calculated as (t=) × particle volume fraction, is the dimensionless time. As seen in this ;gure, the volume fraction does not a>ect the aggregation process only if the time axis is normalized with unit time and multiplied by the volume fraction. These results coincide with the previous studies (Mountain et al., 1986; Kaplan & Gentry, 1988). On the other hand, the growth of aggregates was studied for various charged states. Fig. 4 shows the variation of the reduced number concentration with the average number of charges on each primary particle, which varies from 0 to 0.3 in the unipolar charged states. The
1378
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
Fig. 2. Flow chart for Brownian dynamic simulation including electrostatic force.
reduced number concentration of the charged particles is almost the same as that of the uncharged particles in the early stage, but it becomes much smaller than that of the uncharged state in the ;nal stage. In the early stage, the repulsion force of the charged particles is smaller than the thermal force because aggregates have small charge quantities. However, in the middle or ;nal stage, the repulsion force becomes so strong as to overcome the thermal force because the aggregates have large charge quantities. Therefore, the aggregation is almost stopped under some typical reduced number concentrations for the aggregation of unipolar charged particles.
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1379
Fig. 3. Comparison of Brownian dynamic simulation versus general dynamic equations for instantaneously coalescing particles.
Fig. 4. E>ect of particle charge quantity on reduced number concentration for instantaneously coalescing particles.
4.2. Morphology of aggregate in fractal model The fractal model is almost the same as the sphere model except that there is no coalescence condition. Therefore, in this case, the morphology of the aggregates becomes a very important
1380
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
Fig. 5. Visualization of aggregates at N0 =N (t) = 99 for two di>erent charged cases, (a) initially uncharged case, n = 205, (b) 0.3 unipolar charged case, n = 91.
factor to determine the dynamic behavior of aggregates. The morphology of an aggregate is visualized for the case of the initially uncharged particles in Fig. 5(a). This aggregate, consisting of 205 primary particles, is obtained from a simulation based on 500 primary particles in a cube with edge length 31d0 . Its structure is quite open to the outer space and it is very similar to the structures of aggregates observed by many researchers (Hurd & Flower, 1988; Meakin & Wasserman, 1984; Mountain et al., 1986; Forrest & Witten, 1979). On the other hand, the structure of the aggregate in Fig. 5(b) is more open than that in Fig. 5(a). This aggregate, having 91 primary particles, was obtained from the same conditions of Fig. 5(a) except for the charge quantity. In this condition, the primary particles have 0.3 unipolar elementary charges. The e>ect of the charge quantity on primary particles will be quantitatively investigated in the following section.
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1381
Fig. 6. Radius of gyration, in units of d0 , versus number of primary particles in aggregates at N0 =N (t) = 99 for two di>erent charged cases, (a) uncharged case, (b) 0.3 unipolar charged case.
4.3. Fractal dimension of charged aggregates The morphology is quanti;ed with a fractal dimension. It is obtained from the relation between the number of primary particles and the radius of gyration, as represented in Eq. (1). If at least two aggregates are plotted and linked with a line in a log–log chart, the fractal dimension is obtained from the slope of the line. In this work, a large number of aggregates, which have more than four primary particles, are plotted to improve reliability. These data were obtained for N0 =N (t) = 99. The relation of Rg and n is almost linear in the log–log chart as observed in Fig. 6(a). In this case, the fractal dimension, Df , is 1.761 and its standard deviation is 0.060. The value for Df is in good agreement with the results of Mountain et al. (1986) on
1382
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
Table 1 Fractal dimension for various charged states; these data are obtained for N0 =N (t) = 99 Charged states
BC 3
BC 1
UN 0
UC 0.2
UC 0.3
Df Std Dev. of Df Proportional constant
1.756 0.046 5.473
1.732 0.044 6.069
1.761 0.060 5.658
1.641 0.052 6.122
1.537 0.043 5.998
Df = 1:7–1.9, Meakin and Wasserman’s results (1984) for Df = 1:75 ± 0:05 and Forrest and Witten’s experiments (1979) for Df = 1:7–1.9. According to our results, the power law shown in Eq. (1) is modi;ed as 1:761 Rg n = 5:658 : (15) d0 On the other hand, as shown in Fig. 6(b), if the primary particles are charged with unipolar charges, the fractal dimension becomes smaller than in the uncharged case. The fractal dimension, Df , is 1.537, and its standard deviation is 0.043 for the 0.3 unipolar-charged case. The e>ects of electrostatic force on the morphology have been summarized in Table 1, where “BC”, “UC”, and “UN” indicate bipolar, unipolar, and uncharged states; and the numeric values represent the average number of charges on the primary particles at the initial state. As observed in this table, the fractal dimensions vary little for the charge quantities in the bipolar states. However, in the unipolar states, the fractal dimensions decrease as the charge fractions increase. This tendency may be understood, as the charged small aggregate does not penetrate deeply into the large one having a large number of charges, because of the strong repulsion. In this work, the charges on the aggregate have been assumed concentrated on its center of mass. Therefore, the small aggregates are mainly stuck on the outside of large aggregates. Consequently, the structure of the aggregates in the unipolar charged states is more open than that of the uncharged state. On the other hand, the structures of the aggregates in the bipolar states are thought to be very similar to that of the uncharged state except the early stage of aggregation because the aggregates are neutralized immediately in performing a collision with an aggregate having opposite polarity. 4.4. Size distribution of aggregates Fig. 7 shows the number distribution of aggregates based on the volume for various charged states. These graphs obtained from 100 aggregates for each charged state. To test how many aggregates were needed for detailed number distribution, the average and the standard deviation of the distribution were obtained for the spherical aggregate as the number of samples increases from 50 to 200. As shown in Table 2, if the number of aggregates is larger than 100, the standard deviation converges on a stable value, however, for 50 aggregates, the standard deviation is more a>ected according to the selection of sampling group. Therefore, 100 aggregates were used to represent the number distribution of aggregates for optimizing the accuracy and calculation time.
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1383
Fig. 7. Aggregate number density for N0 =N (t) = 99 for (a) spherical aggregates, and (b) fractal aggregates. This histogram plotted in the interval of 50 primary particles. These data are obtained for N0 =N (t) = 99. Table 2 Optimal number of aggregates for stable number distribution; these data are obtained for N0 =N (t) = 99 No. of samples Average (v0 ) Std Dev. (v0 )
50 100 76.30
50 100 81.35
50 100 83.18
50 100 77.21
100 100 79.84
100 100 78.47
200 100 78.96
These datum were obtained when the 500 primary particles decreased to 5 aggregates, so the average size of aggregates was 100v0 and it was same to each other. Table 3(a) and (b) shows the average and the standard deviation of the number distribution for Fig. 7(a) and (b). In the uncharged case, the standard deviation is 160v0 , and this value is
1384
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
Table 3 Average volume and standard deviation of number distribution for (a) spherical aggregates, and (b) fractal aggregates; these data are obtained for N0 =N (t) = 99 (a) Spherical aggregates Charging state Average (v0 ) Std Dev. (v0 )
UN 100 79.84
BC 0.3 100 77.68
UC 0.1 100 41.46
UN 100 158.64
BC 0.3 100 147.04
UC 0.3 100 124.24
(b) Fractal aggregates Charging state Average (v0 ) Std Dev. (v0 )
much larger than that of the spherical case. This result comes from the di>erence of the collision cross-sectional area, which is larger for the fractal aggregate than the spherical aggregate. It is well consistent with the results of Wu and Friedlander (1993). On the other hand, if the primary particle is charged with ions, the volume distribution of aggregates is di>erent for the polarity of charge. In the unipolar charged case, the standard deviations of volume are 41:46v0 for the spherical aggregates, and 124:24v0 for the fractal aggregates, which are much smaller than that of the uncharged case. These results due to the repulsion force of the unipolar charged aggregates. As the aggregation proceeds, large aggregates having many charges are formed and strong electrostatic repulsion between aggregates is introduced. While small aggregates continue to collide among themselves and with larger ones, large aggregates rarely collide. Therefore, the distribution is sharper than the uncharged case. In the bipolar charged case, the standard deviation is 77:68v0 for the spherical aggregates, and 147:04v0 for the fractal aggregates, which are slightly less than that of the uncharged case. The bipolar charged particles collide frequently with each other more than uncharged particles due to their attraction force, and this force mainly a>ects the collision frequency in the early stage of aggregation. In the middle and ;nal stage, the aggregates are neutralized by the bipolar collision. Consequently, the small deviation of bipolar charged aggregates could be obtained because the number of middle size aggregates increases, however the number of small aggregates decreases slightly. 4.5. Growth of aggregates The growth of aggregates for various charged states has been simulated in Figs. 8 and 9. In Fig. 8, the horizontal axis represents the nondimensional time, and the vertical is the reduced number concentration, which represents the average volume. As shown in this graph, under the charge quantity 1, bipolar and uncharged cases show a very similar tendency. However, upper the charge quantity 1, the average volume of bipolar charged aggregate is larger than that of
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1385
Fig. 8. Reduced number concentration, N0 =N (t) − 1, versus reduced time, 1 ((t=) × volume fraction), for various bipolar charged states.
Fig. 9. Reduced number concentration, N0 =N (t) − 1, versus reduced time 1 ((t=) × volume fraction), for various unipolar charged states.
uncharged aggregate, especially in the early and middle stage of aggregation. This tendency is extinct in the ;nal stage. This can be understood from strong attraction between bipolar charged particles and the subsequent neutralization. There are two electrostatic forces in the bipolar case. One is the
1386
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
attraction force between opposite polarity charges, and the other is the repulsion force between the same polarity charges. Both of them cause aggregates to grow. Therefore, the reduced number concentration of aggregates increased mainly in the early stage of the aggregation, and this e>ect disappears immediately only if the bipolar collision is performed. This is shown clearly for the large charge quantity, over the charge quantity of 1, as shown in Fig. 8. Fig. 9 shows the growth rate of the unipolar charged particles. As the charge quantity increased, the reduced number concentration decreased. These results come from the electrostatic repulsion between the charged aggregates with same polarity. The curves start deviating from the uncharged case at typical values for the reduced number concentration. These values are about 2 and 3 for the two charge quantities of 0.3 and 0.2, where the average number of charges on the aggregates is 0.6. At this condition, we can infer that the electrostatic force is almost the same order of magnitude as the thermal force. There is another interesting result in Fig. 9. Aggregates formed from the fractal model for charged particles become larger than that of the sphere model for uncharged particles. This result can be explained as the e>ect of the competition of two e>ects. One is the fractal dimension, and the other is the electrostatic repulsion force. As the fractal dimension decreases, the collision cross section of an aggregate becomes large, and then the collision rate increases. On the other hand, as the charge quantity increases in the unipolar state, the repulsion force between colliding particles increases, and then the collision rate decreases. In our results, the collision rate is mainly a>ected by the decrease of the fractal dimension, not by the increase of the electrostatic repulsion force. Therefore, the fractal model should be considered in studying the aggregation of charged particles. 5. Conclusions A numerical technique for simulating the aggregation of charged particles was presented with a Brownian dynamic simulation in the free molecular regime. To treat the in;nite sum of electrostatic force, a modi;ed electrostatic force has been introduced. This model coincided with the periodic boundary condition very well. The fractal dimension for the uncharged case is Df = 1:761, and this value agreed with results from previous studies. In the bipolar state, the fractal dimension varied rarely as compared to the uncharged state, but the average size of aggregates was larger than the uncharged case in the early and middle stages of the aggregation process. In the unipolar case, the fractal dimension decreased from 1.641 to 1.537 and the average size of aggregates decreased as the average number of charges increased from 0.2 to 0.3. This calculation was carried under the assumption of the point charge concentrated at the center of aggregate. Therefore, our conclusions cannot be applied directly to other cases, where the charges are not concentrated at the center of an aggregate. The charge position is a>ected by the particle conductivity. For example, if the particle were a non-conductor, the charge would be ;xed on the original primary particle. If the particle were a conductor, the charge would tend to spread on the surface and its distribution is not uniform for the non-spherical aggregate case. They would be located in the positions minimizing the electrostatic potential in the aggregate. To ;nd out the position is a much time consuming process, and the position can be determined
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1387
only statistically. Therefore, it needs some additional models, and we will perform additional simulations and experiments to con;rm the e>ect of the charge location as future works. Acknowledgements S. Kim and H. Park express their gratitude for the support by a grant from the National Research Laboratory program of the Ministry of Science and Technology and the Brain Korea 21 program of the Ministry of Education. H. Chang expresses his gratitude for the support by the grant no. 98-02000-03-01-3 from the Basic Research Program of the Korea Science and Engineering Foundation. References Adachi, M., Kousaka, Y., & Okuyama, K. (1985). Unipolar and bipolar di>usion charging of ultra;ne aerosol particles. Journal of Aerosol Science, 16, 109–123. Brush, S. G., Sahlin, H. L., & Teller, E. (1966). Monte Carlo study of one-component plasma. I. The Journal of Chemical Physics, 45, 2102–2121. Emark, D. L. (1975). A computer simulation of charged particles in solution. I. Technique and equilibrium properties. The Journal of Chemical Physics, 62, 4189–4196. Ermak, D. L., & Buckholz, H. (1980). Numerical integration of the Langevin equation: Monte Carlo simulation. Journal of Computational Physics, 35, 169–182. Forrest, S. R., & Witten Jr., T. A. (1979). Long-range correlations in smoke-particle aggregates. Journal of Physics A: General Physics, 12, L109–L117. Hansen, J. P. (1973). Statistical mechanics of dense ionized matter. I. Equilibrium properties of the classical one-component plasma. Physical Review A, 8, 3096–3109. Hurd, A. J., & Flower, W. L. (1988). In situ growth and structure of fractal silica aggregates in a Game. Journal of Colloid and Interface Science, 122, 178–192. Julien, R., & Bodet, R. (1987). Aggregation and aggregates. Singapore: World Scienti;c. Kaplan, C. R., & Gentry, J. W. (1988). Agglomeration of chain-like combustion aerosols due to Brownian motion. Aerosol Science Technology, 8, 11–28. Kasper, G. (1980). Electrostatic dispersion of homopolar charged aerosols. Journal of Colloid and Interface Science, 81, 32–40. Kittel, C. (1967). Introduction to solid state physics (91pp.). New York: Wiley. Kruis, F. E., Fissan, H., & Peled, A. (1998). Synthesis of nanoparticles in the gas phase for electronic, optical and magnetic applications—a review. Journal of Aerosol Science, 29, 511–535. Matsoukas, T., & Friedlander, S. K. (1991). Dynamics of aerosol aggregate formation. Journal of Colloid and Interface Science, 146, 495–506. Meakin, P., & Wasserman, Z. R. (1984). Some universality properties associated with cluster–cluster aggregation model. Physics Letters A, 103, 337–341. Mountain, R. D., Mulholland, G. W., & Baum, H. (1986). Simulation of aerosol agglomeration in the free molecular and continuum Gow regimes. Journal of Colloid and Interface Science, 114, 67–81. Pratsinis, S. E. (1988). Flame aerosol synthesis of ceramic powders. Progress in Energy Combustion Science, 24, 197–219. Turq, P., Lantelme, F., & Friedman, H. L. (1977). Brownian dynamics: Its application to ionic solutions. The Journal of Chemical Physics, 66, 3039–3044. Vemury, S., & Pratsinis, S. E. (1995). Corona-assisted Game synthesis of ultra;ne Titania particles. Applied Physics Letters, 66, 3275–3277.
1388
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
Wu, M. K., & Friedlander, S. K. (1993). Enhanced power law agglomerate growth in the free molecule regime. Journal of Aerosol Science, 24, 273–282. Xiong, Y., Pratsinis, S. E., & Mastrangelo, S. V. R. (1992). The e>ect of ionic additives on aerosol coagulation. Journal of Colloid and Interface Science, 153, 106–117.
www.elsevier.com/locate/jaerosci
Brownian dynamic simulation for the aggregation of charged particles Hyungho Parka , Sangsoo Kima ;∗ , Hyuksang Changb a
Aerosol & Particle Technology Laboratory, Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, Kusong-Dong, Yusong-Gu, Taejon 305-701, South Korea b Environmental Aerosol Engineering Laboratory, Department of Environmental Engineering, Yeungnam University, Kyungsan 712-749, South Korea Received 17 October 2000; accepted 19 March 2001
Abstract A numerical technique based on Brownian dynamic simulation for the aggregation of charged particles in the free molecular regime is presented. The Langevin equation is used for tracking each particle making up an aggregate. A periodic boundary condition is used for calculation of the aggregation process in each cell with 500 primary particles of 16 nm in diameter. Particle motion is based on the thermal force and the electrostatic force. The electrostatic force on a particle in the simulation cell is considered as a sum of electrostatic forces from other particles in the original cell and its replicate cells. We assume that the electric charges accumulated on an aggregate are located on its center of mass, and aggregates are only charged with pre-charged primary particles. The morphological shape of aggregates is described in terms of the fractal dimension. The fractal dimension for the uncharged aggregate was Df = 1:761, and changed slightly for the various amounts of bipolar charge. However, in case of unipolar charge, the fractal dimension decreased from 1.641 to 1.537 with the increase of the average number of charges on the particles from 0.2 to 0.3 in initial states. During the early and middle stages of aggregation process, the average aggregate size in the bipolar charge state was larger than in the uncharged state, but was almost equal in the ;nal stage. On the other hand, in the unipolar charge state, the average size of c 2001 Elsevier an aggregate and the dispersion of particle volume decreased with increasing charge. Science Ltd. All rights reserved.
1. Introduction Nanoparticles have many properties that di>er from the corresponding bulk material, and these properties o>er many attractive applications such as in electric, optic and magnetic ∗
Corresponding author. Tel.: +82-42-869-3002; fax: +82-42-869-3210. E-mail address: [email protected] (S. Kim). c 2001 Elsevier Science Ltd. All rights reserved. 0021-8502/01/$ - see front matter PII: S 0 0 2 1 - 8 5 0 2 ( 0 1 ) 0 0 0 6 3 - 5
1370
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
Nomenclature d0 Df F kb ke L mp m0 n N0 N qi ; qj r r0 Rg t T u vavg v0 X
diameter of primary particles, m mass fractal dimension, dimensionless electrostatic force, kg m s−2 Boltzmann constant, kg m2 s−2 K −1 Coulomb constant, kg m3 s−2 C−2 length of side of simulation cell, m mass of aggregate, kg mass of primary particle, kg number of primary particles in aggregate, dimensionless total number of primary particles in simulation cell, dimensionless total number of aggregates in simulation cell, dimensionless charge quantities on ith and jth aggregates, C displacement of aggregate, m center of mass for aggregate, m radius of gyration for aggregate, m time, s gas temperature, K aggregate velocity, m s−1 volume of aggregate, m3 volume of primary particle, m3 random force due to collision between gas and particle, kg m s−2
Greek letters friction constant, s−1 electrostatic potential energy, kg m2 s−2 volume fraction of particles, dimensionless particle density, kg m−3 p unit time, s 1 dimensionless time, dimensionless Subscripts k components of rectangular coordinates, k = 1; 2; and 3
products (Kruis, Fissan, & Peled, 1998). The synthesis of nanoparticles is mostly done in a gas phase process. Especially, for making TiO2 or SiO2 particles, the Game process is adopted as an eHcient tool in many applications. In the Game synthesis of nanoparticles, the particles are produced by nucleation or chemical conversion, and growing by condensation or coagulation. These processes are very common to many aerosol processes. In case of thermally stable particles, the coagulation should be considered as a dominant mechanism in the particle growth
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1371
(Pratsinis, 1988; Hurd & Flower, 1988). Therefore, the thermally stable particles in the Game grow from spherical primary particles to aggregates. The aggregation characteristics are mainly a>ected by the morphology of aggregates (Julien & Bodet, 1987; Matsoukas & Friedlander, 1991; Vemury & Pratsinis, 1995). Because a fractal dimension can quantify the morphology of aggregates, the fractal dimension becomes a key variable for analyzing the aerosol dynamics. For the analysis of the aggregation process, two methods have been applied. One is the fractal modi;ed general dynamic equation and the other is the statistical approach. The former is more eHcient than the latter in saving CPU time in the calculation. The application of the fractal modi;ed general dynamic equation is limited, as the fractal dimension should be known in the calculation process. Therefore, the statistical approach is the only choice for the case of unknown fractal dimension. Some researchers have developed numerical simulation techniques to characterize the growth of aggregates. Using theoretical simulation, Meakin and Wasserman (1984) have obtained almost the same fractal dimension compared with the results of previous experiments. However, their method is not proper as an investigation tool for aerosol dynamics because of its size invariant assumption for the di>usion of a particle. Mountain, Mulholland, and Baum (1986) simulated the aggregation process of the aerosol by using Brownian dynamics. They obtained exactly not only the fractal dimension but also the size of aggregates. With the increase of new applications for nanoparticles recently, it has been important to control the size and the morphology of nanoparticles in each nanoparticle synthesis process. Among the numbers of control variables, the electrostatic force has been adopted as a means for precise control of the aggregation. Xiong, Pratsinis, and Mastrangelo (1992) simulated the growth of Game aerosols with the sectional aerosol model. In their results, if the aerosols were in a unipolar charge state, the particle size distribution was more uniform and the average size was smaller than that of an uncharged state. However, their results were obtained using the assumption of spherical particles. Therefore, a new simulation model is needed for explanation in realistic cases. In this study, we wanted to quantify the e>ect of electrostatic force on the morphology and growth of charged aggregates. Therefore, we introduced a new simulation model that could describe the electrostatic force between charged particles of in;nite number. This electrostatic simulation model was coupled with the Brownian dynamic simulation to deal with the aggregation of charged particles. We studied about not only fractal aggregates, but also spherical aggregates to show the e>ect of morphology in the aggregation of charged particles. On the other hand, we treated only a simple case of point charges located on the center of an aggregate to avoid some complexities of the calculation, although this method could be extended to some cases, where the charges were distributed in the volume or on the surface of an aggregate. 2. Theory 2.1. Power law The morphological state of the aggregates is de;ned as the fractal number by Forrest and Witten (1979) who found that the fractal dimensions of the Game aerosol are 1.7–1.9 in the
1372
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
synthesis of iron, zinc, and silicon oxide particles. They con;rmed that there was a power law relation between the characteristic length and the number of primary particles in each aggregate. Usually the radius of gyration of the aggregate is used as the characteristic length for the aggregate. The power law relation is expressed as f n ∼ RD g ;
(1)
where n is the number of primary particles in an aggregate, Rg is the radius of gyration, R2g = (1=n)(rj − r0 )2 ; and r0 is the center of mass de;ned as r0 = (1=n)rj . 2.2. Langevin equation Brownian dynamics describes the dynamic behavior of a particle whose mass and size are larger than those of the gas molecules they are immersed in. These particles are subject to stochastic collisions with the gas molecules, which lead to the random motion of the particles. It is assumed that a particle is in thermal equilibrium with a surrounding gas at the temperature of T and is subject to an external force F. In this case, the translational motion of a particle is described by the Langevin equation as d (mp uk ) = Fk − mp uk + Xk : dt
(2)
In Eq. (2), the LHS represents the inertia force of particle of mass mp , and the ;rst term on the RHS is systematic external force acting on the particle. This term includes van der Waals force, electrostatic force, magnetic force, image force, gravity, etc. The other forces except the electrostatic force can be negligible in the aggregation of charged nanoparticles because the magnitude of the other forces are much smaller than that of the electrostatic force. The second term in the RHS is the frictional drag of the Guid around the particle. is the friction constant, and its inverse −1 represents the particle relaxation time. In the course of its motion, the particle displaces the same length of its mean free path with the mean thermal velocity of the particle between each collision. The last term, Xk , is a random acceleration force caused by the random collisions between the particle and the surrounding gas molecules. The random force Xk is represented with the Gaussian random distribution function of which the mean and square mean are represented as Xk = 0;
Xk2 = 2mp kb T:
(3)
2.3. Solution of the Langevin equation Several methods have been used to integrate the Langevin equation (Turq, Lantelme, & Friedman, 1977). In this study, we use the method derived by Ermak and Buckholz (1980). The analytic solution of the Langevin equation is based on the assumption that a single particle is in a constant force ;eld within the integrating time step. The velocity and the displacement
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1373
are expressed as
F0; k (1 − e−t ) + B1 ; mp 1 F 2F0; k 1 − e−t + uk + u0; k − + 0; k t + B2 : −t mp mp 1+e
ukx = u0; k e−t +
(4)
rk = r0; k
(5)
As shown in Eqs. (4) and (5), the displacement is proportional to the new velocity, old velocity, old displacement and old force F0; k . B1 and B2 are random variables, which have the characteristic components of the Gaussian random distribution as B1 = 0;
kb T (1 − e−2t ); mp B2 = 0; 2kb T 1 − e−t 2 B2 = t − 2 ; mp 2 1 + e−t B12 =
B1 B2 = 0:
(6)
3. Numerical simulation 3.1. Particle charging In this simulation, we assumed that primary particles were charged with some additive ions. According to the comparison between the characteristic time for ion-particle and particle–particle collisions (Xiong et al., 1992), the collision time of ion-particle is shorter than that of particle– particle about four order of magnitude for 16 nm particles and some alkali ions. Therefore, it can be assumed that an aggregate is charged in the conditions of the pre-charged primary particles and a neutral gas. In this case, the charge quantity of an aggregate is the sum of elementary charges on the primary particles constituting the aggregate. Therefore, the number of charge on an aggregate is proportional to the volume of the aggregate. On the other hand, we set the value of initial charge on the primary particles using the results of Adachi, Kousaka, and Okuyama (1985). The charge quantity on the primary particle mainly can be obtained from the relation of the particle size, ion concentration, and resident time, etc. According to their results, small size particles, 16 nm in diameter, are charged partially at the given conditions. Therefore, we set average charge quantities to be maximum at 0.3 in the unipolar charged case. This value means that initially 30% of the primary particles individually carry a single charge. On the other hand, in the bipolar charging case, if the charge quantity is 0.3, it means that 15% of the particles carry a positive single charge and 15% of the particles carry a negative single charge. The bipolar charge quantity becomes large up to 3.0. Although this quantity seems an extreme case, we used this value to show some charging e>ects clearly.
1374
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
3.2. Initial condition To simulate the aggregation process, initial conditions for the velocity and the position of the primary particles have to be given. The positions of these particles are obtained with a uniform random number generator in 3-dimensions. On the other hand, it can be assumed that the primary particles experience only thermal force in the initial state. Therefore, we set the initial velocities with a normal random number generator that is characterized by the Gaussian distribution of zero mean and kb T=m0 variance. We calculated the trajectories of 500 primary particles in a simulation cell. The cell size was determined so that the particle volume fraction is 0.0087. Ten simulations for each initial condition were performed to obtain statistically meaningful average values (Mountain et al., 1986). The unit length was the same as the size of the primary particle d0 . The unit time, = (m0 d20 =kb T )1=2 , was ;xed as the time when a single sphere freely moved a distance equal to its diameter. In this work, we have ;xed the values of d0 , T , and p to be 16 nm, 1500 K and 2000 g m−3 , respectively. These conditions were obtained from the experimental conditions in a Game aerosol process. We performed the simulation for a small time step as 1=50 of the unit time to minimize the overlap between colliding particles. 3.3. Boundary condition The periodic boundary condition was used not only for maintaining a constant volume fraction but also for calculating the electrostatic force. The electrostatic force of the ith particle of the simulation cell is an in;nite summation of each force from the other charged particles in the original cell as well as in the replicated cells. Therefore, this force includes the electrostatic dispersion e>ect mentioned in the paper of Kasper (1980). The electrostatic force is expressed as Fi =
∞ j=1; j=i
ke
qi qj rˆij : rij2
(7)
We can obtain an accurate value when the force is calculated for suHciently many particles because this force is inversely proportional to the square of the distance. But as it is a very time consuming process, Brush, Sahlin, and Teller (1966) and Emark (1975) introduced a new technique for treating an in;nite number of pair potential. This potential includes the interaction with the images of the original cell in the replicated cells. Though in this procedure they could reduce the in;nite summation to a single summation over only the remaining ions in the original cell, their method can be applied only in a system of zero net charge. This constraint is a necessary one to get the ;nite potential. Therefore, it is impossible to apply this method directly to the case of a unipolar charged state. Fortunately, in this study, we need only the electrostatic force, not the whole potential. The force can be calculated directly from the gradient of the pair potential without accessing the pair potential. Though an in;nite number of particles a>ect the electrostatic force of a given particle, the force has a ;nite value because it is symmetrical in a macroscopic sense. In our simulation, electric charges locate in the center of mass of the aggregate and the dielectric constant of the
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1375
gas is uniform. By applying this method to the present system, the modi;ed electrostatic force by the in;nite number of particles can be represented as Fi = −∇i N * = − ke qi qj [∇ 1 (rij ) + ∇ 2 ( rij )];
(8)
j=i
where i is the potential energy of ith particle and N is the total number of aggregates in the original cell; qi and qj are the quantities of charges on the ith and jth aggregates. 1 and 2 are expressed as √ 1 (r) = erfc( r=L)=r − 1=L − Em =L; ∗
1 2( r ) = L * *
√ * * erfc( | r =L + l |)
l
*
*
| r =L + l |
*
+
*
exp(− l2 )cos(2 l · r =L)
l2
;
(9)
where L is the length of the side of cell; Em is the Madelung constant and its value, Em = * 2:83729479 for a simple cubic lattice (Kittel, 1967); l is a vector with integer components * and l is its magnitude; and ∗ denotes exclusion of the l = 0. To treat the equation of 2 , we have used the optimized expansion in Cubic harmonics, which was ;rst introduced by Hansen (1973). Consequently, their gradients are represented as √ 1 exp(− A) erfc( A) ∇ 1 (r)|k = − 2 2 xk eˆ k ; (10) + L A A3=2 * 1 1 ∇ 2 ( r )|k = − 2 # + xk + $xk3 + %xk7 eˆ k ; (11) L xk where #, , $, and % are described as #=−
2EFM ; A3
2 EGB − EJB 4EGB + 2 EFM − EKM + A2 A3 6EFM + 2 EDC − ECL 8EDC + + ; A4 A5
= 2 EH − EI +
$= −
4EG ; A2
%= −
8ED : A4
(12)
1376
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
The uppercase letters in Eq. (12) are de;ned as *
A = | r =L|2 = x12 + x22 + x32 ; B = x14 + x24 + x34 ; C = x18 + x28 + x38 ; D = d8 A4 + d10 A5 + d12 A6 ; E = exp(− A); F = c6 A3 + c8 A4 + c10 A5 + c12 A6 ; G = b4 A2 + b6 A3 + b8 A4 + b10 A5 ; H = a4 A2 + a6 A3 + a8 A4 + a10 A5 ; I = 4a4 A + 6a6 A2 + 8a8 A3 + 10a10 A4 ; J = 4b4 A + 6b6 A2 + 8b8 A3 + 10b10 A4 ; K = 6c6 A2 + 8c8 A3 + 10c10 A4 + 12c12 A5 ; L = 8d8 A3 + 10d10 A4 + 12d12 A5 ; M = x12 x22 x32 :
(13)
The coeHcients a1 ; : : : ; d12 are listed in Hansen (1973). Fig. 1 shows a comparison of the simple Coulomb force and the modi;ed electrostatic force. The values were computed along a diagonal axis. The simple Coulomb force was calculated just within only one cell. There is no di>erence if two particles are close together, but there is a relatively large di>erence as the distance between them becomes large. Because of the periodic boundary condition, the electrostatic force must be zero if two particles are located at the center and the edge of the simulation cell. The modi;ed electrostatic force satis;es this condition well. Therefore, our model for the electrostatic force of an in;nite number of charged particles can be regarded as a good model coinciding with the real system. 3.4. Calculation procedure We de;ne that an aggregate collides with the other if the distance between a particle in one aggregate and the other particle in other aggregate is located within the unit length, d0 . The momentum and the charge quantity are conserved between before and after the collision, and immediately the new-formed aggregate is in thermal equilibrium state. We assumed that the
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1377
Fig. 1. Comparison between modi;ed electrostatic force versus simple Coulomb force along the diagonal axis.
sticking probability was 1 after each collision. The aggregates grow with the repetition of the above sequence. The simulation procedure is summarized in Fig. 2. 4. Results 4.1. Sphere model of Brownian dynamic simulation A comparison between well-known general dynamic equation models and the Brownian dynamic simulation was performed to verify our simulation technique. It was assumed that colliding spherical particles coalesce instantaneously to form larger spherical particles. We de;ned this model as the sphere model. As shown in Fig. 3, in the early stage of aggregation a reduced number concentration is similar to that of a monodisperse aerosol. However, it becomes similar to the self-preserving case in the ;nal stage, where the reduced number concentration is directly related to the average volume of aggregates as N0 =N (t) − 1 = vavg =v0 − 1:
(14)
v0 is the volume of the primary particle 1 , calculated as (t=) × particle volume fraction, is the dimensionless time. As seen in this ;gure, the volume fraction does not a>ect the aggregation process only if the time axis is normalized with unit time and multiplied by the volume fraction. These results coincide with the previous studies (Mountain et al., 1986; Kaplan & Gentry, 1988). On the other hand, the growth of aggregates was studied for various charged states. Fig. 4 shows the variation of the reduced number concentration with the average number of charges on each primary particle, which varies from 0 to 0.3 in the unipolar charged states. The
1378
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
Fig. 2. Flow chart for Brownian dynamic simulation including electrostatic force.
reduced number concentration of the charged particles is almost the same as that of the uncharged particles in the early stage, but it becomes much smaller than that of the uncharged state in the ;nal stage. In the early stage, the repulsion force of the charged particles is smaller than the thermal force because aggregates have small charge quantities. However, in the middle or ;nal stage, the repulsion force becomes so strong as to overcome the thermal force because the aggregates have large charge quantities. Therefore, the aggregation is almost stopped under some typical reduced number concentrations for the aggregation of unipolar charged particles.
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1379
Fig. 3. Comparison of Brownian dynamic simulation versus general dynamic equations for instantaneously coalescing particles.
Fig. 4. E>ect of particle charge quantity on reduced number concentration for instantaneously coalescing particles.
4.2. Morphology of aggregate in fractal model The fractal model is almost the same as the sphere model except that there is no coalescence condition. Therefore, in this case, the morphology of the aggregates becomes a very important
1380
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
Fig. 5. Visualization of aggregates at N0 =N (t) = 99 for two di>erent charged cases, (a) initially uncharged case, n = 205, (b) 0.3 unipolar charged case, n = 91.
factor to determine the dynamic behavior of aggregates. The morphology of an aggregate is visualized for the case of the initially uncharged particles in Fig. 5(a). This aggregate, consisting of 205 primary particles, is obtained from a simulation based on 500 primary particles in a cube with edge length 31d0 . Its structure is quite open to the outer space and it is very similar to the structures of aggregates observed by many researchers (Hurd & Flower, 1988; Meakin & Wasserman, 1984; Mountain et al., 1986; Forrest & Witten, 1979). On the other hand, the structure of the aggregate in Fig. 5(b) is more open than that in Fig. 5(a). This aggregate, having 91 primary particles, was obtained from the same conditions of Fig. 5(a) except for the charge quantity. In this condition, the primary particles have 0.3 unipolar elementary charges. The e>ect of the charge quantity on primary particles will be quantitatively investigated in the following section.
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1381
Fig. 6. Radius of gyration, in units of d0 , versus number of primary particles in aggregates at N0 =N (t) = 99 for two di>erent charged cases, (a) uncharged case, (b) 0.3 unipolar charged case.
4.3. Fractal dimension of charged aggregates The morphology is quanti;ed with a fractal dimension. It is obtained from the relation between the number of primary particles and the radius of gyration, as represented in Eq. (1). If at least two aggregates are plotted and linked with a line in a log–log chart, the fractal dimension is obtained from the slope of the line. In this work, a large number of aggregates, which have more than four primary particles, are plotted to improve reliability. These data were obtained for N0 =N (t) = 99. The relation of Rg and n is almost linear in the log–log chart as observed in Fig. 6(a). In this case, the fractal dimension, Df , is 1.761 and its standard deviation is 0.060. The value for Df is in good agreement with the results of Mountain et al. (1986) on
1382
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
Table 1 Fractal dimension for various charged states; these data are obtained for N0 =N (t) = 99 Charged states
BC 3
BC 1
UN 0
UC 0.2
UC 0.3
Df Std Dev. of Df Proportional constant
1.756 0.046 5.473
1.732 0.044 6.069
1.761 0.060 5.658
1.641 0.052 6.122
1.537 0.043 5.998
Df = 1:7–1.9, Meakin and Wasserman’s results (1984) for Df = 1:75 ± 0:05 and Forrest and Witten’s experiments (1979) for Df = 1:7–1.9. According to our results, the power law shown in Eq. (1) is modi;ed as 1:761 Rg n = 5:658 : (15) d0 On the other hand, as shown in Fig. 6(b), if the primary particles are charged with unipolar charges, the fractal dimension becomes smaller than in the uncharged case. The fractal dimension, Df , is 1.537, and its standard deviation is 0.043 for the 0.3 unipolar-charged case. The e>ects of electrostatic force on the morphology have been summarized in Table 1, where “BC”, “UC”, and “UN” indicate bipolar, unipolar, and uncharged states; and the numeric values represent the average number of charges on the primary particles at the initial state. As observed in this table, the fractal dimensions vary little for the charge quantities in the bipolar states. However, in the unipolar states, the fractal dimensions decrease as the charge fractions increase. This tendency may be understood, as the charged small aggregate does not penetrate deeply into the large one having a large number of charges, because of the strong repulsion. In this work, the charges on the aggregate have been assumed concentrated on its center of mass. Therefore, the small aggregates are mainly stuck on the outside of large aggregates. Consequently, the structure of the aggregates in the unipolar charged states is more open than that of the uncharged state. On the other hand, the structures of the aggregates in the bipolar states are thought to be very similar to that of the uncharged state except the early stage of aggregation because the aggregates are neutralized immediately in performing a collision with an aggregate having opposite polarity. 4.4. Size distribution of aggregates Fig. 7 shows the number distribution of aggregates based on the volume for various charged states. These graphs obtained from 100 aggregates for each charged state. To test how many aggregates were needed for detailed number distribution, the average and the standard deviation of the distribution were obtained for the spherical aggregate as the number of samples increases from 50 to 200. As shown in Table 2, if the number of aggregates is larger than 100, the standard deviation converges on a stable value, however, for 50 aggregates, the standard deviation is more a>ected according to the selection of sampling group. Therefore, 100 aggregates were used to represent the number distribution of aggregates for optimizing the accuracy and calculation time.
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1383
Fig. 7. Aggregate number density for N0 =N (t) = 99 for (a) spherical aggregates, and (b) fractal aggregates. This histogram plotted in the interval of 50 primary particles. These data are obtained for N0 =N (t) = 99. Table 2 Optimal number of aggregates for stable number distribution; these data are obtained for N0 =N (t) = 99 No. of samples Average (v0 ) Std Dev. (v0 )
50 100 76.30
50 100 81.35
50 100 83.18
50 100 77.21
100 100 79.84
100 100 78.47
200 100 78.96
These datum were obtained when the 500 primary particles decreased to 5 aggregates, so the average size of aggregates was 100v0 and it was same to each other. Table 3(a) and (b) shows the average and the standard deviation of the number distribution for Fig. 7(a) and (b). In the uncharged case, the standard deviation is 160v0 , and this value is
1384
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
Table 3 Average volume and standard deviation of number distribution for (a) spherical aggregates, and (b) fractal aggregates; these data are obtained for N0 =N (t) = 99 (a) Spherical aggregates Charging state Average (v0 ) Std Dev. (v0 )
UN 100 79.84
BC 0.3 100 77.68
UC 0.1 100 41.46
UN 100 158.64
BC 0.3 100 147.04
UC 0.3 100 124.24
(b) Fractal aggregates Charging state Average (v0 ) Std Dev. (v0 )
much larger than that of the spherical case. This result comes from the di>erence of the collision cross-sectional area, which is larger for the fractal aggregate than the spherical aggregate. It is well consistent with the results of Wu and Friedlander (1993). On the other hand, if the primary particle is charged with ions, the volume distribution of aggregates is di>erent for the polarity of charge. In the unipolar charged case, the standard deviations of volume are 41:46v0 for the spherical aggregates, and 124:24v0 for the fractal aggregates, which are much smaller than that of the uncharged case. These results due to the repulsion force of the unipolar charged aggregates. As the aggregation proceeds, large aggregates having many charges are formed and strong electrostatic repulsion between aggregates is introduced. While small aggregates continue to collide among themselves and with larger ones, large aggregates rarely collide. Therefore, the distribution is sharper than the uncharged case. In the bipolar charged case, the standard deviation is 77:68v0 for the spherical aggregates, and 147:04v0 for the fractal aggregates, which are slightly less than that of the uncharged case. The bipolar charged particles collide frequently with each other more than uncharged particles due to their attraction force, and this force mainly a>ects the collision frequency in the early stage of aggregation. In the middle and ;nal stage, the aggregates are neutralized by the bipolar collision. Consequently, the small deviation of bipolar charged aggregates could be obtained because the number of middle size aggregates increases, however the number of small aggregates decreases slightly. 4.5. Growth of aggregates The growth of aggregates for various charged states has been simulated in Figs. 8 and 9. In Fig. 8, the horizontal axis represents the nondimensional time, and the vertical is the reduced number concentration, which represents the average volume. As shown in this graph, under the charge quantity 1, bipolar and uncharged cases show a very similar tendency. However, upper the charge quantity 1, the average volume of bipolar charged aggregate is larger than that of
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1385
Fig. 8. Reduced number concentration, N0 =N (t) − 1, versus reduced time, 1 ((t=) × volume fraction), for various bipolar charged states.
Fig. 9. Reduced number concentration, N0 =N (t) − 1, versus reduced time 1 ((t=) × volume fraction), for various unipolar charged states.
uncharged aggregate, especially in the early and middle stage of aggregation. This tendency is extinct in the ;nal stage. This can be understood from strong attraction between bipolar charged particles and the subsequent neutralization. There are two electrostatic forces in the bipolar case. One is the
1386
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
attraction force between opposite polarity charges, and the other is the repulsion force between the same polarity charges. Both of them cause aggregates to grow. Therefore, the reduced number concentration of aggregates increased mainly in the early stage of the aggregation, and this e>ect disappears immediately only if the bipolar collision is performed. This is shown clearly for the large charge quantity, over the charge quantity of 1, as shown in Fig. 8. Fig. 9 shows the growth rate of the unipolar charged particles. As the charge quantity increased, the reduced number concentration decreased. These results come from the electrostatic repulsion between the charged aggregates with same polarity. The curves start deviating from the uncharged case at typical values for the reduced number concentration. These values are about 2 and 3 for the two charge quantities of 0.3 and 0.2, where the average number of charges on the aggregates is 0.6. At this condition, we can infer that the electrostatic force is almost the same order of magnitude as the thermal force. There is another interesting result in Fig. 9. Aggregates formed from the fractal model for charged particles become larger than that of the sphere model for uncharged particles. This result can be explained as the e>ect of the competition of two e>ects. One is the fractal dimension, and the other is the electrostatic repulsion force. As the fractal dimension decreases, the collision cross section of an aggregate becomes large, and then the collision rate increases. On the other hand, as the charge quantity increases in the unipolar state, the repulsion force between colliding particles increases, and then the collision rate decreases. In our results, the collision rate is mainly a>ected by the decrease of the fractal dimension, not by the increase of the electrostatic repulsion force. Therefore, the fractal model should be considered in studying the aggregation of charged particles. 5. Conclusions A numerical technique for simulating the aggregation of charged particles was presented with a Brownian dynamic simulation in the free molecular regime. To treat the in;nite sum of electrostatic force, a modi;ed electrostatic force has been introduced. This model coincided with the periodic boundary condition very well. The fractal dimension for the uncharged case is Df = 1:761, and this value agreed with results from previous studies. In the bipolar state, the fractal dimension varied rarely as compared to the uncharged state, but the average size of aggregates was larger than the uncharged case in the early and middle stages of the aggregation process. In the unipolar case, the fractal dimension decreased from 1.641 to 1.537 and the average size of aggregates decreased as the average number of charges increased from 0.2 to 0.3. This calculation was carried under the assumption of the point charge concentrated at the center of aggregate. Therefore, our conclusions cannot be applied directly to other cases, where the charges are not concentrated at the center of an aggregate. The charge position is a>ected by the particle conductivity. For example, if the particle were a non-conductor, the charge would be ;xed on the original primary particle. If the particle were a conductor, the charge would tend to spread on the surface and its distribution is not uniform for the non-spherical aggregate case. They would be located in the positions minimizing the electrostatic potential in the aggregate. To ;nd out the position is a much time consuming process, and the position can be determined
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
1387
only statistically. Therefore, it needs some additional models, and we will perform additional simulations and experiments to con;rm the e>ect of the charge location as future works. Acknowledgements S. Kim and H. Park express their gratitude for the support by a grant from the National Research Laboratory program of the Ministry of Science and Technology and the Brain Korea 21 program of the Ministry of Education. H. Chang expresses his gratitude for the support by the grant no. 98-02000-03-01-3 from the Basic Research Program of the Korea Science and Engineering Foundation. References Adachi, M., Kousaka, Y., & Okuyama, K. (1985). Unipolar and bipolar di>usion charging of ultra;ne aerosol particles. Journal of Aerosol Science, 16, 109–123. Brush, S. G., Sahlin, H. L., & Teller, E. (1966). Monte Carlo study of one-component plasma. I. The Journal of Chemical Physics, 45, 2102–2121. Emark, D. L. (1975). A computer simulation of charged particles in solution. I. Technique and equilibrium properties. The Journal of Chemical Physics, 62, 4189–4196. Ermak, D. L., & Buckholz, H. (1980). Numerical integration of the Langevin equation: Monte Carlo simulation. Journal of Computational Physics, 35, 169–182. Forrest, S. R., & Witten Jr., T. A. (1979). Long-range correlations in smoke-particle aggregates. Journal of Physics A: General Physics, 12, L109–L117. Hansen, J. P. (1973). Statistical mechanics of dense ionized matter. I. Equilibrium properties of the classical one-component plasma. Physical Review A, 8, 3096–3109. Hurd, A. J., & Flower, W. L. (1988). In situ growth and structure of fractal silica aggregates in a Game. Journal of Colloid and Interface Science, 122, 178–192. Julien, R., & Bodet, R. (1987). Aggregation and aggregates. Singapore: World Scienti;c. Kaplan, C. R., & Gentry, J. W. (1988). Agglomeration of chain-like combustion aerosols due to Brownian motion. Aerosol Science Technology, 8, 11–28. Kasper, G. (1980). Electrostatic dispersion of homopolar charged aerosols. Journal of Colloid and Interface Science, 81, 32–40. Kittel, C. (1967). Introduction to solid state physics (91pp.). New York: Wiley. Kruis, F. E., Fissan, H., & Peled, A. (1998). Synthesis of nanoparticles in the gas phase for electronic, optical and magnetic applications—a review. Journal of Aerosol Science, 29, 511–535. Matsoukas, T., & Friedlander, S. K. (1991). Dynamics of aerosol aggregate formation. Journal of Colloid and Interface Science, 146, 495–506. Meakin, P., & Wasserman, Z. R. (1984). Some universality properties associated with cluster–cluster aggregation model. Physics Letters A, 103, 337–341. Mountain, R. D., Mulholland, G. W., & Baum, H. (1986). Simulation of aerosol agglomeration in the free molecular and continuum Gow regimes. Journal of Colloid and Interface Science, 114, 67–81. Pratsinis, S. E. (1988). Flame aerosol synthesis of ceramic powders. Progress in Energy Combustion Science, 24, 197–219. Turq, P., Lantelme, F., & Friedman, H. L. (1977). Brownian dynamics: Its application to ionic solutions. The Journal of Chemical Physics, 66, 3039–3044. Vemury, S., & Pratsinis, S. E. (1995). Corona-assisted Game synthesis of ultra;ne Titania particles. Applied Physics Letters, 66, 3275–3277.
1388
H. Park et al. / Aerosol Science 32 (2001) 1369–1388
Wu, M. K., & Friedlander, S. K. (1993). Enhanced power law agglomerate growth in the free molecule regime. Journal of Aerosol Science, 24, 273–282. Xiong, Y., Pratsinis, S. E., & Mastrangelo, S. V. R. (1992). The e>ect of ionic additives on aerosol coagulation. Journal of Colloid and Interface Science, 153, 106–117.