NBS

PUBLICATIONS A 11 IDE ( 1 b 3 3 A

NATL. INST OF STANDARDS & TECH R.I.C.

All 10261 6328

Mulholland, George W/Simulation of aeros QC100 .U56 NO. 86-3342 1986 V19 C.1 NBS-P Simulation of Aerosol

Agglomeration in the Free

Molecular and Continuum

Regimes

George W. Mulholland

Raymond D. Mountain

Howard Baum

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Fire Research

Gaithersburg, MD 20899

March 1986

Prepared for: Defense Nuclear Agency

'"'““irtment of Defense

lington, DC 20305

100

.1156

# 86-3542

1986

c.2

Flow

RESEARCH

INFORMATION CENTER

NBSIR 86-3342

SIMULATION OF AEROSOL

AGGLOMERATION IN THE FREE

MOLECULAR AND CONTINUUM FLOW

REGIMES

George W. Mulholland

Raymond D. Mountain

Howard Baum

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Fire Research

Gaithersburg, MD 20899

March 1986

Prepared for: Defense Nuclear Agency

Department of Defense

Washington, DC 20305

U.S. DEPARTMENT OF COMMERCE. Malcolm Baldrige, Secretary

NATIONAL BUREAU OF STANDARDS. Ernest Ambler, Director

TABLE OF CONTENTS

Page

List of Figures iv

Abstract 1

1. INTRODUCTION 2

2. COMPUTATIONAL TECHNIQUE 4

3. FRICTION FORCE 8

3.1 Free Molecular Regime 8

3.2 Continuum Regime 10

4. RESULTS 20

4.1 Structure 20

4.2 Kinetics 23

4.3 Cluster Size Distribution 30

5. COAGULATION EQUATION 35

6. COMPARISON WITH EXPERIMENTS 41

7. ACKNOWLEDGMENT 42

8. REFERENCES 43

t

-iii-
LIST OF FIGURES

Page

Figure 1.

Figure 2.

Figure 3.

Figure 4.

Planar projection of an agglomerate with 90 spheres ...

vs k. for 3t = 0.05 and p = 0.0167

g

2 R vs k for continuum flow for p = 0.0167

g

Reduced number concentration versus time steps for free

molecular flow

21

22

24

26

Figure 5. Reduced number concentration versus reduced time for

Bt = 0.05 28

Figure 6. Reduced number concentration versus t N for ..continuum

flow 29

Figure 7. The cluster size distribution is plotted for free

molecular simulation for N

Q /N(t) = 4.5 31

Figure 8. The cluster size distribution is plotted for N

Q /N(t) = 8.7

for the same simulation as Figure 7 32

Figure 9. The cluster size distribution is plotted for the continuum

simulation for N

Q /N(t) = 3 33

Figure 10. The cluster size distribution is plotted for N

Q /N(t) = 9 for the same simulation as Figure 9 34

-lv-
SIMULATION OF AEROSOL AGGLOMERATION IN THE FREE

MOLECULAR AND CONTINUUM FLOW REGIMES

George W. Mulholland, Raymond D. Mountain and Howard Baum

Abstract

The formation of high temperature aerosol agglomerates is simulated by

following the Langevin trajectory of each particle with the boundary condition

that the particles stick upon collision. Both the free molecular and

continuum flow are treated. A new derivation of the friction force of an

agglomerate in the continuum limit is developed based on the evaluation of the

surface momentum flux at the Oseen flow limit. The agglomerates can be

described as a fractal, at least in regard to power law relationship between

mass and size, with a dimensionality of 1.7-1. 9 independent of the flow

regime. The particle growth is shown to be much more rapid in the free

molecular regime than in the continuum. The global kinetics are shown to be

consistent with a similarity analysis of the coagulation equation with a

modified coagulation coefficient. Comparison between the simulation and

coagulation theory at small time suggests a slight fluctuation enhancement in

the free molecule case and a small-time enhancement of the coagulation rate at

high concentration for the continuum case.

Keywords: aerosol agglomeration, agglomerates. Brownian motion, Darcy's

law, fractal dimension, free molecular regime, Smoluchowski equation, soot.

- 1 -
1 . INTRODUCTION

Aerosol agglomerates are formed in high temperature processes including

soot formed in flames and in engines, silica produced through the combustion

of silane, and metal particles produced by vaporization of a metal as a result

of joule heating. This last method is known as the exploding wire technique

[1] for aerosol production. The particle growth processes include nucleation,

surface growth (or condensation), and coagulation. Initially the particles

may coalesce upon colliding and sticking, but eventually the particles retain

their structure upon collision and a cluster made up of individual spherules

evolves. It is the latter cluster growth process that is the subject of this

paper.

One of the first models of cluster growth was developed by Sutherland [2]

for studying floe structures in colloids. Cluster growth in this model is

treated geometrically as a series of random collisions between particles and

particle clusters. Meakin [3] and Kolb et al. [4] have developed clustercluster aggregation models based on random walk motion of the clusters. The

diffusion coefficient of the cluster is assumed to have a power law dependence

on the number of spheres in the cluster.

In three dimensions the structures found in these models are quite open

and the degree of openness has been characterized in terms of a fractal dimensionality, D, which can be defined by the following equation:

D m ~ r .

- 2 -
The quantity m represents the mass of the object and r its radius. For a

solid three-dimensional object, D=3. Meakin [3] has obtained D * 1.8 for

agglomerates. The fact that the dimensionality is less than three implies

that the average density decreases with increasing radius.

We use the term agglomeration to refer to the growth process in which

clusters collide and stick together. Our usage of the term agglomeration is

essentially identical with Meakin's cluster-cluster aggregation terminology.

This differs from the model originally developed by Void [5] in which

individual particles are added one at a time to a growing cluster. The Voldlike growth process we term aggregation. Witten and Sander [6] have developed

a diffusion limited aggregation model in which the individual particles are

added one at a time via Brownian motion to a stationary aggregate. Richter _e_t

al . [7 ] have used this model for describing soot growth and obtain a more

dense structure than is obtained with the agglomeration models. The Richter

et al . structure has a fractal dimension of about 2.4.

Our simulation differs from the others in that we follow the trajectory

of every cluster. This enables straightforward comparison with the kinetics

of a real system. The so-called Brownian dynamics method that we use places a

premium on the use of an efficient computer code because of the complex structures that develop. The general procedure for carrying out the simulation is

described in the next section, while the architecture of the computer code is

described in a second paper [8].

A necessary input to the simulation is the friction force acting on the

aggl°mera te. In section 3 we derive expressions for the friction force in the

- 3 -
free molecular limit, which is approached at flame conditions, and in the

continuum limit, which becomes more valid as large agglomerates are formed in

a rising smoke plume. Our continuum result is compared with the results

obtained by Meakin £t_ _al . [9] based on the Kirkwood-Riseman theory.

One of the principal results of our study is the characterization of the

agglomerate structure with a fractal dimension. The results of this simulation are compared with Meakin' s results [3] based on random walk of the

clusters on a lattice. The second area of major interest is the global

kinetics of the agglomeration process. The results from the simulation are

compared with a prediction based on a similarity solution of the coagulation

equations in section 5. The results are also compared with available experimental results on the structure of the agglomerate and kinetics in section 6.

2. COMPUTATIONAL TECHNIQUE

This simulation is based on the assumption that the agglomerating

particles move subject to the rules of Brownian motion plus the condition that

if two agglomerates "touch", as defined below, they stick and the resulting

agglomerate also rigidly diffuses according to the Brownian motion rule. The

"rule" for Brownian motion is that the dynamics of a particle are governed by

the Langevin equation

d(mv )/dt = - mBv + f ( 1 ) X XX

where vx

is the xth cartesian component of the velocity of the center of mass

of the agglomerate of mass m and f is a stochastic force statisfying

- 4 -
2 ~1

<f > = 20mk T. As can be seen from eq. (1), 6 represents the relaxation

x o time of the agglomerate. The quantity 0 is the controlling parameter for the

dynamics. In addition, the center of mass of the aggregate moves according to

dx/dt = v

XJ (2)

subject to periodic boundary conditions. In the work reported here, we have

used the Ermack and Buckholtz [10] solution of the Langevin equations given by

r=r + i (v + v 1 tanh (0h/2) + §.

o 0 o 2 (3)

and

+ -0h , £ v = v e + B,

o 1 (4)

for the positions, {r}

,

and velocities, {v}

,

at time t+h given positions {r } o and velocities {v^} at time t. The stochastic integrals ^

and satisfy

< B

x > = < B

2 > = < B

l

• B

2 > = 0, (5)

2 k R T

< B^ > = 3 — [1 - exp ( —2 0 h ) ]

,

( 6 )

and

< b

2

2 >

6 k

B T

m0

2

[ 0h - 2 tanh ( 0h/2 )

]

(7)

- 5 -
This solution assumes that the diffusing particles are in thermal

equilibrium with the background gas through which they move. It should be

noted that this solution is contained in the 1943 Reviews of Modern Physics

article by Chandresakhar [11]. The development in this article is also an

important advance in the application of stochastic methods to random processes

in condensed matter physics.

Initially, we used a program which was a straightforward adaptation of

the molecular dynamics codes which had been used by Mountain [12] in other

studies. As a result of lessons learned by applying some of the techniques

developed in computer science to the molecular dynamics problem [13], we

attempted to apply these techniques to the agglomeration problem. The result

was a reduction by a factor of 4 in the execution time of the simulation

program. A detailed description of the efficient, scalar program is described

in a second paper [8]. The basic idea implemented in this program is that

particles which are within unit distance of other particles "touch" and are

therefore within the same cluster and are identified as such. Then list

manipulation techniques are used to implement the sticking conditions and to

identify the newly formed clusters in updated lists.

It is our practice to generate ten sets of initial conditions so that

adequate averaging is possible. The initial conditions consist of a set of

coordinates for 500 spheres of mass m

Q and unit diameter a. The coordinates

are generated by a uniform random number generator which places the initial

positions of the particles in a cube of side L. The initial velocities are

obtained using a random number generator which produces normally distributed

numbers with unit variance so that the particles are in thermal equilibrium

- 6-
with the background gas through which they diffuse. In this way, the temperature, T, and hence the energy unit, kgT, for the simulation, is specified.

The unit of time T = (m a^/kT)^^ is the time required for a single sphere to

freely stream a distance equal to its diameter o. The particles move

according to equations (3) and (4) for a time interval h. Periodic boundary

conditions are used to eliminate the influence of the surface of the cube of

side L on the motion of the agglomerates. After each interval h, the system

is examined to see if any agglomeration events have occurred. If a new

agglomeration event has occurred, the lists containing cluster information are

updated and the size and the radius of gyration of the existing clusters are

determined. This process is repeated until a specified number of iterations

of the equations of motion have been made. As noted above, averages over ten

sets of initial conditions are made. This results in good statistics for the

smaller clusters but relatively poor statistics for the few larger clusters

generated.

It should be noted that the clusters are only allowed to translate. In

principle the clusters should also be allowed to rotate, subject to a rotational Brownian motion rule. We intend to examine this topic in a future

study. Meakin [14] has shown for a two-dimensional model that rotation in

"equilibrium" with translation does not significantly change the fractal

dimensionality of the resulting clusters.

- 7 -
3 . FRICTION FORCE

3.1 Free Molecular Regime

Soot particles in flames approach free molecular conditions because of

both the small particle size (~ 0.02 urn) and the increased mean free path of

the gas (0.2 - 0.3 Mm) at the flame temperature. In this limit, Epstein [15]

showed that the friction coefficient for a sphere, K, is proportional to the

2 particle cross section, ttR , and to the average thermal speed, v, of the host

gas molecules.

K = m 3 = % 6ttR2 o 3 P V

m v

g

( 8 )

The quantity mg refers to the mass of a gas molecule, P to the gas pressure,

and 5 to the surface accommodation factor which has a value of 1 for specular

reflection and a value of 1.444 for diffuse reflection.

We approximate the friction coefficient of a cluster made up of k spheres

as k times the friction coefficient of a single sphere. In making this

approximation we neglect the shielding effect of the other spheres, but for a

tenuous, low density agglomerate this is a reasonable first approximation. In

any event, it is an approximation that will ultimately require testing.

As pointed out earlier in the paper, the quantity 8~^

, which appears in

eq . (1), corresponds to the relaxation time of the agglomerate. The apparent

agglomerate mean free path, X is obtained from v 8 where v refers to thr

P P P average thermal speed of the particle. If X is large compared to the

P

-8 -
particle diameter, the particle behavior is free molecular. Equivalently, the

free molecular condition corresponds to the relaxation time, 8 , being long

compared to the time, t, to free stream a particle diameter; that is,

8t « 1 . (9)

From the definition of t and from eq. (8), we obtain the following result in

the free molecular limit

Bt

4 /2 6P

P k T

s B

/mm , o g

( 10 )

where m

Q is the mass of an individual sphere, p

g refers to its density, and mg

the mass of the gas molecules. Our simulations are carried out for Bt = 0.2,

which corresponds to an 8 nm particle radius for a 1500 K flame temperature,

and for Bt = 0.05, which corresponds to a 3 nm radius. The density of the

individual particles is taken to be 2.0 g/cnP and ambient pressure is assumed.

In the simulation it is assumed that when two particles collide and stick

together, the new aggregate moves with the thermally accommodated velocity.

The average interparticle spacing at the start of the simulation is in the

range 2.7 particle diameters to 5.8 particle diameters. In an actual aerosol,

the interparticle spacing is much greater, on the order of 100 or greater.

This greater spacing provides time for thermal accommodation with the host gas

before a newly formed agglomerate collides with another agglomerate. An

estimate of the momentum relaxation time can be obtained by integrating the

Langevin equation, eq . (1), to obtain the root mean particle velocity

2 1/2

< v (t) > , assuming the particle is initially at rest.

- 9 -
—2 3

>. 1 /2

) (ID

2 1/2 —1 The time, tm> at which < vx

(t) > increased to (1-e ) of its final value is

given by

tm

0.255

B

( 12 )

The corresponding free streaming distance traveled would be 1.3 and 5.1

particle diameters for Bx equal to 0.2 and 0.05, respectively. So it is seen

that thermal accommodation may not occur at the high concentration, low Bx

case, before a second collision occurs in the simulation. To avoid this nonphysical effect, the agglomerates are thermally accommodated after each collision in the simulation.

As the soot particles grow to a size approaching a micrometer in diameter

in the near ambient conditions of the rising smoke plume, the continuum limit

may become a more valid description for the friction force because of the

increased particle size and decreased mean free path. Continuum flow is valid

in the limit of the characteristic length scale of the structure being much

greater than the mean free path of the gas. We first derive a general expression for the friction force in the continuum limit in terms of a surface

integral involving the stress tensor and momentum flux. The flow field is

derived based on Stokes flow in the vicinity of the agglomerate and Oseen flow

far from the agglomerate. The solution technique is similar to the mat-hid

asymptotic expansion applied by Proudman and Pearson [16] to a solid sph«r<*.

3.2 Continuum Regime

- 10 -
The technique is also described by Van Dyke [17]. In our case we have a

porous agglomerate rather than a solid sphere. We treat the agglomerate as a

porous medium obeying Darcy's law in regard to flow within the agglomerate.

Felderhof [18] has treated an analogous problem of the friction coefficient of

a polymer by a technique that apparently does not involve Oseen flow.

Because previous analyses have not included the convective term, we

present a consistent hydrodynamic analysis of the friction force including the

convective term even though the result is similar to that of the previous

analyses [18]. For the reader interested in the final expression for the

friction force used in the simulations, he may go directly to eq. (45) without

loss of continuity.

We begin our analysis with the equation for the momentum flux for an

incompressible fluid with density p.

We consider steady state flow for which the time derivative in eq. (13)

vanishes. To obtain the friction force, we first integrate the right hand

side over the entire fluid volume contained between the complex agglomerate

(13)

surface, s a> and the outer surface, s Q , of radius r

Q

. Then using the divergence theorem, we obtain

9 Pv.v.n. dA - <j> Pv.v.n.dA = J> a n dA - so l k k sa l k k

r so ik k

( 14 )

where n^ refers to the outward normal to the surface and where

- 11 -
' 1

3v

i

3vk \ V * - P

6 ik + °ik

= - P

6 ik + 11 + W[) 05

At the surface of the agglomerate, the fluid velocity vanishes so the second

term on the LHS of eq. (14) vanishes. The second term on the RHS of eq . (14)

corresponds to the force on the particle; that is, to the friction force, F^»

that we wish to calculate.

F.

i

) o n, dA

sa lk k

(16)

Thus we obtain the following .general result for the friction force of an

agglomerate of arbitrary shape in terms of a surface integral.

F.

i

Pv .v. In. dA

l k

J k

(17)

To solve for F^ we must determine the pressure and velocity for large r.

For incompressible flow, the steady state Navier-Stokes equation can be

written as'

P (v • V) v+Vp-nAv=0 (18)

with

V • v = 0. (19)

First we consider the behavior close to the agglomerate for which the Reynolds

number, £, is small and the convective term may be neglected. Taking the

divergence of eq . (18) and making use of eq. (19), we obtain

- 12 -
V

2 p = 0 ( 20 )

A trial solution for p in terms of spherical coordinates is given by

p = p oo- v b cosQ x(r).

Substituting from eq. (21) into eq. (20), we find

r

2 X" + 2 rx' - 2 x = 0.

The solution to eq. (22) remaining finite as r + 00 is given by

X = r > R,

r

where R is the radius of the agglomerate.

We chose the fluid velocity far from the agglomerate to be v

Q

in the

direction. It is convenient to use Felderhof's [18] suggested form for v

spherical coordinates.

v(r) = v {4> ( r ) cos6e - [

<J> ( r ) + l- r <j> ' (r)l sin6 e.}

o r 2 0

This form insures that V • v(r) = 0. Substituting from eqs. (24) and (21)

for v and p into eq . (18) with the convective term dropped, we find

r

2

<j>" (r) + 4r<t>'(r) + r

2 x'(r) = 0

- 13 -

( 21 )

( 22 )

(23)

z

in

(24)

(25)
Substituting from eq. (23) for x and solving the differential equation, we

obtain

<j> = l- A+ B

r > R

r 3 r

Far from the body such that £r/R is of order

(18), the convective term, can not be neglected,

following reduced variables

> +

v* = v/v

o

(26)

unity, the first term in eq

It is convenient to use the

x! = ex./R

l l

Since there is no longer spherical symmetry with the convective term included,

the equations will be expressed in terms of cartesian coordinates. From eqs.

(23) and (26) , we obtain v| to order £ for the Stokes solution.

r! = 6 . _ + e l i3

A<5 i3

2Rr

Ax! x^

Rr'

3 (28)

The lowest order term in the pressure equation is given by

P

? -

£

2 Ax^

(29)

The form of the solution of the Stokes problem suggests the following

expansion for v' and p' for the solution of the Navie r-Stokes equation to

leading order

- 14 -
( 30 ) v! = 6 + £v.

1 i3 i

£2 (31)

P

Substituting eqs. (30) and (31) into eqs. (18) and (19), we obtain to lowest

order

This version of the Oseen flow problem was analyzed by Lamb [19]. The

solution is

By expanding the exponential term in eq . (34) for the limit of small r'

corresponding to small Reynolds number, £, one finds that eq. (34) is

consistent with the Stokes solution given in eq . (28). We note that taking

the divergence of eq . (18) still results in a harmonic equation for the

pressure. So eq . (29) is still valid for the pressure.

The convective component of the friction force is given by

(32)

(33)

(34)

(35)

Substituting from eqs. (30) and (34) into eq . (35) and expanding, we obtain
( 36 )

F.

l

o k.3 ^ 2 ^ ^

£ <b v,n_dA ~ pv & <p v.v. n. dA

i 3 o

y

l k k

The first term vanishes by symmetry, the second term vanishes since fluid

is incompressible, and the last term vanishes as l/r0

* Substituting for v^

from eq. (34) and carrying out the integration in cylindrical coordinates

suggested by the cylindrically symmetric wake structure, we obtain

F = F = 0 x y

We find that the stress tensor component of F

3 vanishes in the limit of large

r , so eq. (37) gives the total contribution to the friction force.

We point out that Felderhof's result [18] is essentially identical to eq

( 37 ), however, in his case the entire contribution is from the stress tensor

term. He does not include the convective term found in eq . (18) even for the

integration over all space. It is a surprising coincidence that the two

expressions are identical.

To determine the coefficient A, we model the flow inside the agglomerate

as flow through a porous object satisfying Darcy's law.

F 0 = F = 4ttAv n

3 z o (37)

o

VpK + — v = 0

<

( 38 )

- 16 -
The quantity < is the hydrodynamic permeability. Again we assume the fluid to

be incompressible and find the following results for X and within the porous

sphere.

X = Cr r < R (39)

<|> = C< r < R (40)

The coefficient A appearing in eq. (37) is obtained by requiring continuity of

pressure and of radial and angular velocity at the surface of the porous

sphere. We obtain

A 7 • (41)

2/3 + k/R

So from eq. (37) we obtain the following expression for the friction

force

F =

z

6 u R v n o

, 3 *

1 + 2 2 R

(42)

In the impermeable sphere limit k/R2 + 0 the Stokes formula for the friction

2 force is recovered. For the case where k/R >> 1, called the free draining

limit for polymers, we find

F = % "nR3 V /k.

z j o (43)

There are not adequate experimental data to determine which of these limits is more nearly appropriate for aerosol agglomerates. We suspect that the limit

2 </R + 0 to be the correct limit for the reasons discussed below.

- 17 -
Assuming < to be inversely proportional to the agglomerate density and

taking R to be the radius of gyration of the agglomerate, we find that

< R

^1 g

So provided that the fractal dimension D > 1, which is the case as we show

below for agglomerate growth, then the term

y

k/R^ can be neglected for large

r . Felderhof and Deutch [20] derived an inverse density dependence for < in

g

the related problem of the friction coefficient for polymers. In our simulations, lacking specific information on k, we drop the second term in the

denominator and use the expression,

K = m3 = 6tttir . (44)

g

This Stokes-Einstein type expression suggests that the fluid is trapped in the

agglomerate and that the entire assembly moves as if it were a compact sphere

of radius Rg.

For the simulations carried out in this study, we assume the primary

sphere has a radius of 0.5 Um, P =2, and T = 298 K at ambient pressure.

s From eq . (44) and the definition of x, we obtain the following:

3x = 6 PR

g

6it \ 1 / 2 _= 84

6nR,

st =

6tt \ 1/2 167 ( R

g

/a ) P ok T

s B

single sphere

agglomerate

(45)

(46)

3x = 55 doublet (47)

- 18 -
The quantities a refers to the primary particle diameter and Rg the radius of

gyration. The value of 8t for the doublet was taken to be the orientation

averaged value of a spheroid with major axis 2 a and minor axis a [21]. This

value agrees within 10% with the measured friction coefficient based on a

sedimentation experiment [22]. For agglomerates with three or more spheres,

equation (46) is used.

Meakin et al . [9] have calculated the friction force in the continuum

limit based on the Kirkwood-Riseman theory [23]. The hydrodynamic interactions in this theory are the far field limits where a particle is treated as

a force field from a point in regard to its effect on other particles. For a

specified configuration of spheres in the cluster, the friction coefficient

was calculated and expressed in the form

K = 6rnR , h

where R^ is the hydrodynamic radius. Calculations were performed for agglomerates with 50 to 350 particles. The quantity Rh was fitted to a power law

relationship of the form R^ ~ k^ , where k refers to the number of spheres in

the agglomerate and the exponent Y has a value of 0.544 _+ 0.014. The radius

of gyration, Rg , was also found to have a power law relationship of. . the form

R ~ k

a with exponent a, 0.554 +_ 0.038, within the error limits of the

oexponent y. So we see that Rg and R^ have a similar dependence on k and

consequently the expressions for k obtained by our macroscopic hydrodynamic

description of the flow, eq . (44), and obtained by the detailed flow field

analysis by Meakin _et_ a^. are similar. There are no data available to test

the theories.

- 19 -
The values of St are obviously much greater in the continuum than in the

free molecular limit. In the continuum, the agglomerate diffuses only a small

fraction of a particle diameter before its direction of motion is significantly changed whereas the free molecular particles move several particle

diameters before changing their direction. Also in the free molecular case St

is independent of the number of spheres in the agglomerate k, while St

decreases with increasing k in the continuum limit.

4 . RESULTS

4.1 Structure

The structure of the agglomerate is quite open as indicated in Fig. 1 for

a planar projection of the structure. This particular simulation is for the

free molecular limit with St = 0.2, but similar structures are also obtained

in the continuum limit. A quantitative measure of the structure is obtained

from the dependence of the radius of gyration of the agglomerate, Rg » on the

number of spheres in the cluster, k. The radius of gyration is related to the

second moment of the mass distribution and is defined by

o k 9 mR = 2 m. r . , (48)

g i= l 1 1

where r^ is the distance from the center of mass of the agglomerate to the ith

sphere and m is the total mass of the agglomerate.

In Fig. 2, Rg i s plotted versus k for the case 8t = 0.05 and density

P = 0.0167. The quantity P is defined as the number of particles per volume

- 20-
21 -

Figure

1.

Planar

projection

of

an

agglomerate

with

90

spheres

is

displayed

for

a

simulation

based

on

500

spheres

randomly

distributed

in

a

cube

with

edge

length

llo

(8

t

=

0.2).
T3

a; .c

0) CO bO CD

CD

T3

Id

CO * a> 0) >

ca £ m A o <u cu u A •

U-i (0 4j <0 0 Id E V CD tao CO B 4J «r^ •U o CO

•rH

H 5 • e ac » a; O 3 5iC CU c oc

eo

u o • c a. —* •H <U cu c

id u-

«_i

O

ao 4) CO

a; c > •3 c

•vH Id 0 « 3 O c CO O c. CO o 0) 'O c

•H i- 1) 0 -D 4J V ;

u c X ~ — u U — 0

>N COu

GO • U-i

u

X

U-i o vO a> i-> 0 —1 CJ u o 3 X

CO 0) • u 3 o o 3 U- a C- 3 3 n X

CO 3 V u a si 0 u 2

•J

4-» 3 -3 73 —J CO 3 “3

CO CO

C

o- “ 3 X 3 u 0£ 'O

a; a;

O 1)

u > uu x II 0 u 3 73 0

3*h 7) .

CO caU3 1)

U 3 u 3 P C. >

rsivu3

SO

•H

u»

22
where volume is in units of a

3

. Over the cluster size range from about 5 to

200, R has a power law dependence on k given by

g

R = Ak

g

a (49)

where a = 0.55 and A = .42 for R in units of a . This power law dependence

O g

is consistent with a fractal dimensionality D

a given by [24]

(50)

which in this case is 1.8

of the quantity itself.

than about 50 arises from the relative infrequency of any one large agglomerate out of ten repeat simulations. We have also performed simulations

for Bt = 0.2 and for P = 0.05, 0.0167, and 0.005. In all cases the dimensionality inferred from the plots is in the range 1.7-1. 9. Our results for are

We find that over the agglomerate size range accessible to our simulations, 5 to 200, the fractal dimension in the continuum fluid mechanical

regime is essentially identical with the results found in the free molecular

simulations. Meakin [25] arrived at a similar conclusion that the size dependence of the agglomerate diffusion coefficient did not affect the fractal

dimension based on Monte Carlo simulations in two dimensions. The data for

the longest run in the continuum flow regime are given in Fig. 3, which corresponds to P = 0.0167. This fit of the data yields

in good agreement with Meakin* s results [3] for (1.75) based on the dependence of R? on k for Brownian trajectories of clusters on a lattice.

o

- 23 -
o

-24 -

Figure

3.

kg

vs k

for

continuum

flow

for

p

=

0.0167

with

the

dashed

line

corresponding

to

a

fractal

dimension

of

1.70.
R = 0.38 k

g

0.59 (51)

with a fractal exponent of 1.7.

4.2 Kinetics

A major interest in our study is the kinetics of cluster agglomeration.

A global measure of the kinetics is the time dependence of the total number of

clusters, N, defined by

N = E n. (52)

N

The quantity

N(t)

- 1 is convenient for displaying the results of the simulations. First, coagulation theory for a constant coagulation coefficient

predicts that this reduced number concentration is proportional to the product

of the initial number concentration and time. Secondly, the reduced number

concentration is simply related to the change in the average particle volume

N

o N(t) (53)

The results of the free molecular simulation for both Bx = 0.05 and 0.2

are contained in Fig. 4. As the concentration decreases from 0.05 to 0.0167

to 0.005, more time is required for agglomeration to take place. In all cases

we see that the slope increases with time from a value of about 1 at early

time to a value greater than 2 at the time when the total number of cluster

has decreased by a factor of about 20. Provided Bx is small enough so that

-25 -
omOx

00

03o

CMoX

00

CMoo

00o - 26 -

TIME

STEPS

Figure

4.

Reduced

number

concentration

versus

time

steps

for

free

molecular

flow

(0t=

0.05

o, 0t=

0.2 [))

at

three

concentrat

ions, P=

0.05,

0.0167,

and

0.005.
the free molecular limit applies, the kinetics should be insensitive to the

value of St. The fact that the results for Qr = 0.2 lie slightly below the

results for Bt = 0.05 suggests that the agglomerate trajectory is not quite

free molecular in this case. The condition for free molecular behavior

depends not only on the size of the individual sphere but also on the size of

the agglomerate.

In an attempt to correlate the kinetic data, we plot the data versus a

reduced time, t^, in Fig. 5

t = (time steps) (p) (54)

As pointed out above, coagulation theory for constant coagulation coefficient

N

predicts that the reduced number concentration, — 1, is proportional to

the product of the initial number concentration times time. We see that all

of the data are reduced to one curve. For long time the sl«ope increases to

about 2.6, though because of the limited data it is not certain that an asymptotic power law behavior has been reached.

The kinetics of the agglomeration process is considerably different in

the continuum regime. The slope increases only slightly with time up to a

value of about 1.15 compared to 2.6 in the free molecular regime. The use of

the reduced time variable t^ does not collapse all of the data to one curve

(see Fig. 6) as was the case for the free molecular simulations. Also, as

shown in Fig. 6, the small time behavior for the simulations seems to diverge

from the coagulation theory prediction instead of converging to it. This is

surprising since one would expect coagulation theory to be valid at small time

- 27 -
(DN

100

0.1 I I i i I 0.05 0.15 0.5 1.5 5

Figure 5. Reduced number concentration versus reduced time (time steps x for 8t = 0.05 and P = 0.05 (o), p = 0.0167 (), and

P = 0.005 (A). The dashed line corresponds to the caLculiti-Mi m

Lai et al . and the solid line to coagulation theory for constant:

coagulation.

-28 -
(ON

0.5 1.5 5 15 50

Figure 6. Reduced number concentration versus t^ for continuum flow with

P = 0.05 () and P = 0.0167 (0). The solid line corresponds to

the prediction based on the coagulation theory and the dashed

lines to the slopes for long time (1.15) and short time (0.85).

-29 -
when there are relatively few agglomerates. In the next section an explanation for this result will be given.

4.3 Cluster Size Distribution

The information on the cluster size distribution is limited by the

relatively small number of particles in the simulation. The cluster size

distributions are plotted in Figs. 7 and 8 for times corresponding to approximately a 5 fold and 9 fold drop in the initial number concentration. The free

molecular solution of the coagulation equation for coalescing droplets

obtained by Lai et al. [26] is consistent with the data. Even the

Smoluchowski solution [27] for the case of size independent coagulation

coefficient is within the scatter of the "data". The major qualitative

conclusion is that the cluster size distribution is more nearly exponential in

character than power law for large cluster sizes.

The cluster size distributions for the continuum regime are plotted in

Figs. 9 and 10 for times corresponding to approximately a 3 fold and 9 foLd

drop in the initial number concentration. The continuum solution of the

coagulation equation obtained by Hidy [28], as well as the Smoluchowski solution [27], are in reasonable agreement with the simulations. So, while there

is no a priori reason to expect agreement between the cluster agglomeration

simulations and the solution of the coagulation equation, we do in fact see

agreement within the uncertainty of the simulations for clusters of size 30

and less. We must await much larger simulations to better define the cluster

size distribution especially in the large agglomerate region.

- 30 -
I- I I I I I I I 1—1 I % jpt* I

1 3 5 10 20 30 50

k

Figure 7. The cluster size distribution is plotted for free molecular

simulation (St = 0.2, p = 0.05 [0] and St = 0.05, p = 0.05 [•])

for NQ /N(t) = 4.5. The solid curve corresponds to the solution of

the coagulation equation by Lai et al.

-31 -
k

Figure 8. The cluster size distribution is plotted for N

G /N(t) = 8.7 for the same simulations as Fig. 7. The solid line corresponds to Che Lai

et a l » solution and the dashed curve to the Sraoluchowski solution

to the coagulation equation.

-32 -
1 3 5 10 20 30 50

k

Figure 9. The cluster size distribution is plotted for the continuum

simulation (p = 0.0167) for N

Q /N(t) = 3. The solid line

corresponds to the Hidy solution for the coagulation equation.

- 33 -
k

Figure 10. The cluster size distribution is plotted for N

Q /N(t) = 1 for the same simulation as Fig. 9. The solid and dashed lines correspond

to the Hidy and Sraoluchowski solutions of the coagulation

equation, respectively.

- 34 -
5. COAGULATION EQUATION

In this section, we present the results of coagulation theory for

comparison with the simulations. We also use a similarity transform of the

coagulation equation to determine the long time kinetics for an "agglomeration"

type coagulation coefficient.

The evolution of the size distribution function n(v,t) for coalescing

droplets of volume v is described by the coagulation equation:

3n(v,t)/3t = j 1 K (V , o v - v ) n (v* , t) n (v - v’ , t) dv'

- n (v, t) / K (v, v’) n (v', t) dv'

o (55)

The coagulation coefficient K (v, v') for the free molecular (fm) and

continuum (c) regime are given by

* (v, v) . mn

( v

i/3 v

1 / 3

) 2 (i + A,)" 2 V V

fm (56)

2k T

B_ f 1/3 ,1/3^ /

1 K (v, v' ) = [v + v

2kgT

3n

(r + r') (- + ij-),

r x

J c (57)

were r and r' are the radii of spheres with volume v and v' , respectively,

For small time, most of the particles are individual spheres so that

treating K as a constant in eq . (55) is a good approximation. Integrating

-3 5-
both sides of eq. (55) over all particle sizes and solving the resulting

ordinary differential equation for N(t), we obtain

N(t)/N = (l + hr K N t)" 1

. o l o

( 58 )

The expressions for K in the two limits are obtained by setting v = v' .

K = 4 /I

6rkgT\l/2

fm (59)

K

8 k

B T

3n c (60)

It is convenient for comparison with the simulations to express eq. (58) in

the following form:

N

o /N(x

1 ) - 1 - \

K 3 ij,

a (61)

where = (t/t) p. In the two limits, we obtain

N /N (t J ol - 1 = 2 /tt t,

1 fm (62)

4ttt,

VN (t

i

J - 1 - -ft (63)

In the free molecular limit, the simulation results at small time for the

reduced number concentration are about 20 % greater than predicted by eq . (62)

(see Fig. 5). This increase may result from fluctuations arising in the

simulation which are not taken into account in the coagulation equation. '•>-

have also performed simulations without agglomeration in which a pair of

- 36 -
colliding particles is replaced with a single particle of mass m

Q to eliminate

any size dependence of the coagulation process. We still obtain slightly more

rapid growth compared to the prediction of the coagulation equation with

constant K.

As pointed out in the previous section, the simulation and the prediction

based on the coagulation equation for the continuum regime seem to be diverging rather than converging at small time (Fig. 6). Also the results of the

simulation depend on P while the solution of the coagulation equation is

independent of P. The cause of the discrepancy is a small-time enhancement in

the coagulation rate for very high particle concentration [29,30]. The

derivation of eq. (60) involves the calculation of the rate at which particles

arrive at the surface of a sphere, F(t). For equal size spheres of diameter

p, the solution of the diffusion equation with absorbing boundary conditions

[31] is given by

where D is the particle diffusion coefficient. The transient term is dropped

in deriving eq . (60); however, this approximation is not valid at small times.

For example, for tN

= 1.66, which corresponds to 0.05 s for p = 0.0167 and

0.0167 s for P = 0.05, we obtain sizeable transient terms for a 1 urn diameter

sphere in air

(64)

a 0.52

1/2 P = 0.0167

( TTDt )

a 0.90

1/2 P = 0.05

(tfDt)

- 37 -
This transient effect is the major cause for the difference in the small time

behavior for the simulation and the coagulation equation. At lower concentrations the effect is less pronounced since the time required for coagulation is

longer. For aerosols, P would seldom be greater than 10“^, and under these

conditions the transient terra is negligible. We are not able to simulate such

dilute conditions because of the long computation time required. As an aside,

values of p as large as 10“*- exist in colloidal systems, so in this case

deviations of the sort seen in Fig. 6 are to be expected.

We next consider the solution of the coagulation equation for an

"agglomerate" coagulation coefficient to compare with the results of the

simulation. Assuming K(v, v') to be symmetric, it can be shown that

GO 00 = / / K (v, v') n (v) n (v') dvdv ' . (65)

o o

In the free molecular limit, the coagulation coefficient consists of a

collision cross section term multiplying a velocity terra. The cross section

of a low density agglomerate will be much greater than a compact sphere. The

expression for the coagulation coefficient, eq . (56), is modified to apply to

an agglomerate by replacing the v^ 2 term, which is related to the collision

cross section, with the radius of gyration, RCT

» We then make use of the

fractal character of the agglomerate, eq . (49), to express R T in terms of

O

v^/D. So the v dependence of K for the agglomerate can be expressed as

K (v, v' ) ~ (v 1/D + v’ 1/D

) 2 (i + -L-) l/2

. (bh)

V V

38-
To solve for N(t) it is convenient to use the similarity variables ip, the

dimensionless size distribution function, and v, the dimensionless particle

volume

if, =

n- V

N(t)

2

(67)

vN(t)

V

( 68 )

The volume concentration V is conserved during particle coagulation.

Substituting eq. (66) into eq. (65) and expressing n and v in terms of p and

v, we obtain

dN

dt

5_ 2 2 1_ N

2 '

° V

D

' 2

A. (69)

The quantity A includes dimensionless integrals involving various moments of

t( v )« Integrating eq . (69), we obtain

N(t)

where p = (j

-

-1

(70)

(71)

For D = 1.80, we obtain p = 2.57 which is in good agreement with the

simulation result of 2.6. The value of p is very sensitive to the fractal

exponent. For D set equal to 3 , we recover the free molecular solution of Lai

et al . [26], p = 6/5.

- 39 -
Next we consider the continuum limit. In eq . (57) the first term

involving radius relates to the collision radius, whereas the second relates

to the diffusion coefficient. For the case of an agglomerate, we replace both

radii with the radius of gyration but for different reasons. In the first

case it is because the characteristic collision radius of an agglomerate is

R . In the second case the diffusion coefficient of an agglomerate is

O

inversely proportional to the friction force, which is proportional to R for

o an agglomerate. Using the relationship between Rg and v, we obtain

Carrying out the same type of analysis as in the free molecular case, we

obtain

N

o N(t)

(73)

In this case the time dependence is the same as for a compact structure. The

enhancement arising from the increased collision radius is cancelled out by

the decrease in the diffusion coefficient. The simulation yields a slightly

stronger time dependence, t ^ , than that predicted by coagulation theory.

This may be an artifact of the high initial concentration, for which the

coagulation equation is not valid, or that the simulation was not carried to

long enough time.

The analysis of the coagulation equation given above is intuitive. To be

rigorous one must verify that the similarity solution exists. Van Dongen ind

Ernst [32] have derived a generalized scaling form for the size distribution

function. They find

-40 -
z (74)

N(t)

where Z =

The quantity A is determined by the scaling characteristics of the coagulation

coefficient

The results of Van Dongen and Ernst are identical with our results for K given

by eqs. (66) and (72).

Forrest and Witten [33] have obtained fractal dimensions in the range of

1.5 to 1.7 for inorganic smokes, including Fe and Zn. The smokes were

produced by rapidly heating a tungsten wire plated with metal. The fractal

dimension was determined from the power law dependence of the density autocorrelation function.

Feder et al . [34] have inferred from Medalia and Heckmans’ measurements

on carbon black [35] a fractal dimension of 2.3. The difference from the

inorganic smoke may arise from the ultrasonic dispersion of the collected

carbon black particles in a liquid. The dispersion process apparently breaks

weak bonds in the agglomerates resulting in smaller, perhaps, more compact

structures. The present computer simulations are in fair agreement with the

inorganic smoke but are apparently less compact than carbon black.

k (av, av') = a^ k (v, v’) (75)

6. COMPARISON WITH EXPERIMENTS

-41 -
There are very little data on the kinetics of agglomeration in the gas

phase. Howard et al . [36] have measured the coagulation rate of carbon

particles in a low pressure, premixed flame. The primary particle size is

about 10 nm and the pressure about 20 Torr. This corresponds to free

molecular behavior with 8 t ® 0.002. The experimentally observed coagulation

rate was found by Howard _et_ al. to exceed by a factor of 10 the free molecular

coagulation rate approximately adjusted for electrostatic forces and Van der

Waals attraction. The number concentration decreased by a factor of six with

increasing height in the flame. The results of the free molecular simulation

shown in Fig. 5 suggest that as the total number concentration decreases by a

factor of six the agglomeration kinetics will approximately double the growth

rate compared to the solution of the coagulation equation. The electron

micrographs clearly show the formation of agglomerates. However, agglomeration by itself will not account for the entire factor of 10 difference between

the theoretical and experimental result. It is also possible that the

measured absolute number concentration was systematically underestimated by

the electron microscopy. Even under the best of circumstances, a measurement

of the absolute number concentration of an aerosol is difficult.

7. ACKNOWLEDGMENT

The authors acknowledge support by Defense Nuclear Agency in carrying out

this research. The authors also acknowledge Francis Sullivan for his

assistance in developing an efficient computer code.

-42 -
8. REFERENCES

[1] Phalen, R.F., Aerosol Sci J5_, 465 (1974).

[2] Sutherland, D.N., J. Colloid and Interface Sci. 2_5, 373 (1967).

[3] Meakin, P., Phys. Rev. A. 29_, 997 (1984).

[4] Kolb, M. , Botet, R. , and Jullien, R. , Phys. Rev. Lett. _51_, 1123 (1983)

[5] Void, M.J., J. Colloid Sci. 18_, 684 (1963).

[6] Witten, T.A. , Jr. and Sander, L.M., Phys. Rev. Lett. 4_7_, 1400 (1981);

Witten, T.A. and Sander, L.M. , Phys. Rev. B27 , 5686 (1983).

[7] Richter, R. , Sander, L.M., and Cheng, Z. , J. Colloid and Interface Sci

100 , 203 (1984).

[8] Sullivan, F. and Mountain, R.D., submitted to Computer Physics

Communication.

[9] Meakin, P., Chen, Z.Y., and Deutch, J.M., J. Chem. Phys. 82 , 3786

(1985).

[10] Ermack, D.L. and Buckholtz, H. , J. Comp. Phys. _35_, 169 (1980).

[11] Chandrasekhar, S., Rev. Mod. Phys. 1

, 1 (1943).

[12] Mountain, R.D., Phys. Rev. A26 , 2859 (1982).

[13] Sullivan, F., Mountain, R. D. and O'Connell, J., J. Coraput. Phys. (in

press )

[14] Meakin, P. , J. Chem. Phys. _81_, 4637 (1984).

[15] Epstein, P.S., Phys. Rev. _23_, 710 (1924).

[16] Proudman, I. and Pearson, J.R.A., J. Fluid Mech. 2_, 237 (1957).

[17] Van Dyke, M. , Perturbation Methods in Fluid Mechanics , Academic Press,

New York, 153 (1964).

[18] Felderhof, B.U., Physica 80A , 63 (1975).

[19] Lamb, H. , Hydrodynamics , Dover, New York, p. 610 (1932).

[20] Felderhof, B.U. and Deutch, J.M., J. Chem. Phys. 62, 2391 (1975).

[21] Overbeck, H. A., Journal fur die reine und angewandte Mathematik 81, 62

(1876). .

[22] Horvath, H. , Staub Reinhaltung der Luft 34, 197 (1974).

-4 3-
[23J Kirkwood, J.G. and Riseman, J. , J. Chem. Phys . _16_, 565 (1948).

[24] Stanley, H.E., J. Phys. A 10, L 211 (1977).

[25] Meakin, P., Phys. Rev. B29 , 2930 (1984).

[26] Lai, F.S., Friedlander, S.K. , Pich, J. , and Hidy, G.M.J., Colloid

Interface Sci. 39^, 395 (1972).

[27] Smoluchowski, M. , Z. Physik Chem. 92 , 129 (1917a).

[28] Hidy, G.M., J. Colloid and Interface Sci. 20, 123 (1965).

[29] Ziff, R.M., McGrady, E.D. , and Meakin, P. , J. Chem. Phys. 82_, 5269

(1985).

[30] Mulholland, G.W. and Mountain, R.D. , accepted for publication in

J. Chem. Phys.

[31] Friedlander, J.K. , Smoke, Dust and Haze , Wiley and Sons, New York

(1977).

[32] Van Dongen, P.G.J. and Ernst, M.H., Phys. Rev. Lett. _54_, 1369 (1985).

[33] Forrest, S.R. and Witten, T.A. , Jr., J. Phys. A. _12_, L109 (1979).

[34] Feder, J. , Jossang, T. and Rosengvist, E. , Phys. Rev. Let. _53_, 1403

(1984).

[35] Medalia, A. I. and Heckman, F.A. , J. Colloid and Interface Sci. _3_6, 173

(1971).

[36] Howard, J.B., Wersborg, B.L. , and Williams, G.C., Faraday Symposia of

the Chemical Society No. 7, Fogs and Smokes, 109 (1973).

- 44 -
NBS-114A (rev. 2 -aci U.S. DEPT. OP COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO.

NBSIR 86-3342

2. Performing Organ. Report No. 3. Publication Date

March 1986

4. TITLE AND SUBTITLE

Simulation of Aerosol Agglomeration in the Free Molecular and Continuum Flow Regimes

5. AUTHOR(S)

George W. Mulholland, Raymond D. Mountain, and Howard Baum

6. PERFORMING ORGANIZATION (If joint or other than NBS. see instructions)

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

7. Con tract'' Grant No.

8 . Type of Report & Period Coverea

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City, State, ZIP)

Defense Nuclear Agency

..Washington, DC 20305

10. SUPPLEMENTARY NOTES

| | Document describes a computer program; SF-185, FIPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual summary of most si gn i fi cant information. If document includes a si gn i fi cant

bibliography or literature survey, mention it here)

The formation of high temperature aerosol agglomerates is simulated by following

the Langevin trajectory of each particle with the boundary condition that the

particles stick upon collision. Both the free molecular and continuum flow are

treated. A new derivation of the friction force of an agglomerate in the continuum

limit is developed based on the evaluation of the surface momentum flux at the Oseen

flow limit. The agglomerates can be described as fractal, at least in regard to

power law relationship between mass and size, with a dimensionality of 1. 7-1.9

independent of the flow regime. The particle growth is shown to be much more rapid

in the free molecular regime than in the continuum. The global kinetics are shown

to be consistent with a similarity analysis of the coagulation equation with a modified coagulation coefficient. Comparison between the simulation and coagulation

theory at small time suggests a slight fluctuation enhancement in the free molecule

case and a small-time enhancement of the coagulation rate at high concentration for

the continuum case.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolon si aerosol agglomeration; agglomerates; Erownian motion; Darcy's law; fractal dimensions

free molecular regime; Smoluchowski equation; soot

13. availability

I y | Unlimited

| | For Official Distribution. Do Not Release to NTIS

f I Order From Superintendent of Documents, U.S. Government Printing Office Washington. D C

20402.

Order From National Technical Information Service (NTIS). Springfield, VA. 22161

14. NO. OF

PRINTED PAGES

15. Pr

$9.95

USCOMM-OC 4 0 43-P ,0