Investigation of the Effect of Radiation on the Thermophoretic
Motion of Soot Particles in Free?Molecular Regime
Except where reference is made to the work of others, the work described in this
thesis is my own or was done in collaboration with my advisory committee. This
thesis does not include proprietary or classifled information.
Vamshi Krishna Sarangapani
Certiflcate of Approval:
Jay M. Khodadadi
Professor
Mechanical Engineering
Daniel W. Mackowski, Chair
Associate Professor
Mechanical Engineering
Daniel K. Harris
Associate Professor
Mechanical Engineering
Stephen L. McFarland
Acting Dean
Graduate School
Investigation of the Effect of Radiation on the Thermophoretic
Motion of Soot Particles in Free?Molecular Regime
Vamshi Krishna Sarangapani
A Thesis
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulflllment of the
Requirements for the
Degree of
Master of Science
Auburn, Alabama
August 08, 2005
Investigation of the Effect of Radiation on the Thermophoretic
Motion of Soot Particles in Free?Molecular Regime
Vamshi Krishna Sarangapani
Permission is granted to Auburn University to make copies of this thesis at its
discretion, upon the request of individuals or institutions and at their expense.
The author reserves all publication rights.
Signature of Author
Date
Copy sent to:
Name Date
iii
Vita
Vamshi Krishna Sarangapani, son of Venkataramana Chary Sarangapani and Ra-
madevi Sarangapani, was born on August 12, 1980, in Dharmapuri, Andhra Pradesh,
India. He graduated with the degree of Bachelor of Technology in Mechanical Engi-
neering in 2003 from Indian Institute of Technology, Madras, India. He joined the
Masters program in the department of Mechanical Engineering at Auburn University
in August 2003.
iv
Thesis Abstract
Investigation of the Effect of Radiation on the Thermophoretic
Motion of Soot Particles in Free?Molecular Regime
Vamshi Krishna Sarangapani
Master of Science, August 08, 2005
(B.Tech., Indian Institute of Technology, Madras, India. June 2003)
68 Typed Pages
Directed by Daniel W. Mackowski
The coupling of radiation and thermophoresis is hypothesized to result in attrac-
tive forces among soot particles in combustion environments. The efiect of radiation
from the soot particles on their thermophoretic motion in the free?molecular regime
is studied by developing a ?synthetic? simulation model. A Monte Carlo technique is
used to carry out this study and the models are developed both for the two?particle
system and the aggregate system that mimics cluster?cluster aggregation. The trans-
fer of momentum and energy to and from the soot particles are computed via the
monte carlo method. The thermophoretic force and the coagulation ratios are calcu-
lated for the developed models. The results indicate that thermophoresis would be a
signiflcant mode of soot particle coagulation for larger particle sizes and higher gas
temperatures. The sphere aggregates are compared to single spheres with equivalent
volume, surface area, and radius of gyration and the results show that the aggregate
can be approximated to its volume equivalent sphere in a two?sphere model.
v
Acknowledgments
I would like to thank my advisor, Professor Daniel W. Mackowski for his constant
support and guidance throughout my research at Auburn University and providing me
with an opportunity to work in this exciting fleld. His methodology has greatly helped
me hone the thought process that goes into approaching and solving a problem. I
would also like to thank Professors Jay M. Khodadadi and Daniel K. Harris for serving
on my graduate committee. I thank all the professors, stafi, and students I have had
interactions with throughout this study. Finally, I thank all my friends for their help
and support when I needed them the most.
I wish to dedicate this work to my parents and my younger brother whose love
and encouragement provide me the strength to sail through life and achieve my goals.
vi
Style manual or journal used LATEX{ A Document Preparation System, Leslie
Lamport, Addison-Wesley Publishing Company, 2nd edition (1994). Bibliography
follows IEEE Transactions.
Computer software used The document preparation package TEX (speciflcally
LATEX) together with the departmental style-flle aums.sty. The plots were generated
using MATLAB?.
vii
Table of Contents
List of Figures x
List of Tables xi
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Hypothesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Thesis Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2 Literature Review 6
2.1 Experimental Approach . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.1 Soot Aggregates . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.2 Thermophoresis . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.3 Knudsen Number . . . . . . . . . . . . . . . . . . . . . . . . . 9
3 The Model 12
3.1 The velocity distribution function . . . . . . . . . . . . . . . . . . . . 12
3.2 Sampling of the distribution function . . . . . . . . . . . . . . . . . . 15
3.3 The two sphere model . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3.1 Energy transport . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.3.2 Momentum transfer . . . . . . . . . . . . . . . . . . . . . . . . 22
3.3.3 Soot radiation heat source function . . . . . . . . . . . . . . . 23
3.3.4 Thermophoretric force analysis . . . . . . . . . . . . . . . . . 27
3.3.5 Monte Carlo method . . . . . . . . . . . . . . . . . . . . . . . 29
3.4 Aggregate model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.4.1 Modifled Monte Carlo method . . . . . . . . . . . . . . . . . . 32
4 Results and Conclusions 35
4.1 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.1.1 The two sphere model . . . . . . . . . . . . . . . . . . . . . . 35
4.1.2 Aggregate model . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.1.3 Two?sphere equivalent of aggregate model . . . . . . . . . . . 41
4.1.4 Coagulation ratio . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.3 Suggestions for future research . . . . . . . . . . . . . . . . . . . . . . 44
viii
Bibliography 46
Appendices 48
A Monte Carlo code for the two-sphere model 49
B Monte Carlo code for the aggregate model 52
ix
List of Figures
1.1 TEM photograph of a soot aggregate . . . . . . . . . . . . . . . . . . 4
3.1 Two?sphere model showing the ray trace of a molecule . . . . . . . . 18
3.2 Soot aggregate generated using the power law . . . . . . . . . . . . . 33
3.3 Aggregate?model showing the ray trace of a molecule . . . . . . . . . 34
4.1 Two sphere model ? Force vs Distance . . . . . . . . . . . . . . . . . 36
4.2 Log?Log scale of Fig.4.1 . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.3 F1?cluster vs Nsamp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.4 Aggregate model ? Force vs Distance . . . . . . . . . . . . . . . . . . 39
4.5 Log?Log scale of Fig.4.4 . . . . . . . . . . . . . . . . . . . . . . . . . 40
x
List of Tables
4.1 Comparison of thermophoretic force magnitudes . . . . . . . . . . . . 41
4.2 Coagulation ratios for various a and T . . . . . . . . . . . . . . . . . 43
4.3 Coagulation ratios for various a and T . . . . . . . . . . . . . . . . . 43
xi
Chapter 1
Introduction
1.1 Motivation
Combustion plays an important role in automobile industry, industrial burners
and furnaces, gas turbines, etc., as the process involves conversion of chemical energy
in fuels to thermal energy that generates power. In addition to playing a helpful role,
combustion processes also have detrimental afiects on the human lives in the form of
harmful emissions such as NOx, and soot. Soot is formed in gas-phase combustion
at high temperature. At the microscopic level, soot forms when hydrocarbons are
heated with insu?cient air due to poor mixing. In  ames, soot can be observed in
difiusion  ames as opposed to premixed  ames, where a clean blue  ame can be seen.
Keeping in view the importance of soot, its study has been a major area of interest
for quite some time now. In boilers, soot fouling is a big concern as it brings down the
e?ciency. Soot is suspended in air and because of its extremely small size, penetrates
deep in to the lungs, thus afiecting our respiratory system. In some devices, such as
furnaces, the thermal emission from soot enhances heat transfer process via radiation.
Soot is also an important industrial product that flnds application, such as flller in
tires, toner in copiers, etc..
Soot collected from  ames consists of chain-like aggregates of spherical units.
These spherical units have a hexagonal structure similar to graphite [1], and the sizes
range between 10 nm and 80 nm. Soot formation is believed to take place in three
steps: [1] Particle inception or Nucleation: this step involves a series of homogenous
1
reactions between hydrocarbon species leading to a larger sized particles which can
be seen under an electron microscope as tiny spherical condensed particles [2]. The
rate of particle inception is extremely high, though the particles formed in this step
constitute only a minor fraction of the flnal soot mass. (2) Surface Growth: in this
step, two processes are believed to occur simultaneously. First, the small spherical
particles formed in step [1] collide and coalesce to form larger spheres. Second, the
hydrocarbon \growth species" in the gas phase react heterogeneously on the soot
surfaces [2]. (3) Aggregation: the spherical particles formed in step-(2) collide and
stick (but they no longer coalesce) to form chains. The soot formed can be controlled
only by Oxidation. It is the sole mechanism of removal of soot emission. In this
process, the soot particles formed are changed back into gas-phase species.
While soot is a major pollutant, it is an important industrial product and is a
major source of radiation from  ames. In most of the hydrocarbon fuels combustion,
the dominant part of radiation comes from the soot particles. The fourth power
dependence on temperature also makes radiation a prime mode of heat transfer in
most  ames. Radiation allows transfer of energy directly from hot product regions to
cold regions, exerting its efiects even at a distance.
Modelling of the formation, growth, and deposition of soot requires an under-
standing of the mechanisms which transport the soot in a gas. These mechanisms
include convection, difiusion, sedimentation, and thermophoresis. A signiflcant trans-
port of soot in  ames is by thermophoresis. The term ?thermophoresis? is given to
the motion of aerosol particles that move, with a constant velocity, towards lower
temperature regions, under the in uence of a temperature gradient. Thermophoresis
2
occurs because the momentum carried by the hot gas molecules that sink the par-
ticles is greater than that carried by the particles coming from cold region. Studies
have conflrmed that thermophoresis results in a larger deposition rate of aerosols than
either of sedimentation or difiusion processes [3].
1.2 Hypothesis
It can be observed, in  ames, that when a steel rod which is at room temperature
is introduced, the soot particles get deposited at a quick rate on that rod. This
behavior is observed as a result of thermophoresis. The temperature gradient existing
between the emanating soot particles and the steel rod causes the movement of the
particles towards the rod. It is submitted that a similar efiect can occur among the
soot particles themselves. As a result of continuous radiation from the soot particles,
the local temperature of the gas is afiected resulting in temperature gradients. These
temperature gradients lead to an attractive thermophoretic motion among the soot
particles. The objective of this thesis is to carry out investigations on the efiect of
radiation from soot particles on the thermophoretic motion of the soot particles, and
establish that the coupling of radiation and thermophoresis could result in attractive
forces among soot particles.
1.3 Objectives
The goal of this thesis can be achieved either by investigating the process ex-
perimentally or by building a simulation model. There may be some practicalities,
in studying the process by building an experimental set-up, such as low residence
3
Figure 1.1: TEM photograph of a soot aggregate
times. Likewise, building a simulation model of the efiect of soot radiation on the
thermophoretic motion of the particles is a complicated one because of the 3-D nature
of the problem. Nevertheless, the complexity can be reduced if we initially build a
simulation model for a simple system. Therefore, the initial step would be to model
a simple system of soot particles, essentially, a two-particle model. Extension of this
model to a more general system of particles could be carried out at a later stage. In
order to get started, even with the simplest of cases, we need to make some basic
assumptions. First among those would be the shape of the soot particles we are sim-
ulating. Transmission Electron Microscopy (TEM) measurements indicate that soot
aggregates consist of nearly spherical primary particles. Figure 1.1 is a TEM image
of a typical soot aggregate found within the annular region of a difiusion  ame [15].
The measurements show that the aggregates exhibit mass fractal behavior. The  ame
generated soot aggregates exhibit mass fractal-like behavior with a fractal dimension,
Df < 2, even when the number of primary particles in an aggregate is small [4]. The
second assumption would be the size of the particles. The soot particles we consider
would be nearly nano-sized, making them considerably smaller than the mean free
path. This assumption is supported by the gas temperatures we would be dealing
4
with (? 1500 K). At such high temperatures, the mean free path increases, thus
justifying our assumption of knudsen number, Kn >> 1 (i.e., the mean free path >>
the particle size).
1.4 Thesis Statement
This thesis consists broadly of two parts: the flrst part discusses the efiect of
radiation on the thermophoretic motion between two primary aerosol particles (as-
sumed to be spheres), and the second part discusses the efiect of radiation on the
thermophoretic motion between a cluster and a single particle.
This thesis is organized as follows: Chapter 2 gives a more detailed description
of the work done in the area of interest and the existing methodologies (literature
review). Chapter 3 presents details of the proposed methodology and the algorithm
employed. Chapter 4 presents the results of the proposed method and various obser-
vations that can be drawn from those results. Chapter 4 also includes suggestions for
future research.
5
Chapter 2
Literature Review
2.1 Experimental Approach
The main aim of the experimental approach is to visually observe the coagula-
tion of soot particles as a result of the proposed hypothesis. Studies were previously
carried to observe this efiect using co? ow difiusion burners. As a result of low resi-
dence times, it was di?cult to keep the soot particles for su?cient time in the  ame
environment. Burner assemblies that ofier longer residence times than those possible
in the co? ow difiusion burner assembly were looked in for from the available litera-
ture. A  at counter? ow difiusion burner assembly that ofiers longer residence time
was then built in the laboratory. A  at counter? ow difiusion  ame was established
using methane as the fuel. Nitrogen was the inert gas used. The study with this
assembly was also inconclusive, mainly because of two factors. The flrst factor was
again the low residence time (we obtained longer residence time than in the case of
co? ow difiusion  ame, but not good enough to observe the efiect). The order of
residence time we obtained was about 1s. The second factor is that the gelation
process was overpowering all the other efiects, making it di?cult to observe the in-
tended process. Once the aggregate size reaches a large value, the aggregation of
soot particles is largely because of the gelation efiect. The soot cloud forms a sort of
\spider-web," thereby attracting the soot particles to stick to the cloud. This is the
efiect we primarily notice if we run an experiment, thus making the study of efiect of
radiation on thermophoretic motion of soot particles inconclusive.
6
2.2 Simulation Model
Taking into consideration the practicalities in running an experiment, we opted
for a simulation model. The subsequent sections in this chapter discuss brie y on
soot aggregates and thermophoresis.
2.2.1 Soot Aggregates
Computer simulations were carried out to develop models for random cluster
formation by researchers [5,6] since early 80?s. Sorensen, in his paper discusses about
various aggregation algorithms for simulating aggregates [7] using the power law,
N = K0(Rg=a)Df (2.1)
Where N is the number of monomers in the aggregate, Rg is the radius of gyration, a is
the monomer radius, Df is the fractal dimension, and the proportionality constant K0
is the prefactor. In [7], random aggregates have been computer synthesized using both
Difiusion Limited Aggregation (DLA) and Difiusion Limited Cluster Aggregation
(DLCA).
In DLA, a monomer is chosen from a set of N monomers and placed at the center
of the sphere. The second monomer is now brought near the flrst and random-walked
until it is attached at a random angular position. Subsequent monomers are then
introduced one-by-one at a random angular position.
In DLCA, two monomers are chosen at random from a set of N; one of them is
placed at the center and the second one is introduced at a random angular position,
7
random-walked around the flrst one. This cluster is now put back into the set if the
pair of monomers are joined making the set, now, (N-1) long. The process is repeated
to obtain the desired length of cluster.
Aplotoflog(N)vslog(Rg)willproduceastraightlinewithslopeDf andy?intercept
K0.
2.2.2 Thermophoresis
Since its discovery, numerous applications have been recognized where ther-
mophoresis can play either positive or adverse roles [8]. Principle of thermophoresis
has been extensively used in the design of thermal precipitators, and aerosol sampling
methods. Thermophoresis can result in particle deposition on boiler pipes, reducing
the e?ciency of heat exchange. The efiect of thermophoresis on the transfer of the
radioactive aerosols generated in a nuclear reactor accident has been recognized as an
important factor in reactor safety assessment. The principle of thermophoresis can
be used to enhance chemical vapor deposition process which is a key in fabrication
of optical flbers. It has been used as an efiective method for micro-contamination
control in the semiconductor industry [8]. Application areas of efiect of thermophore-
sis also include: aerosol instruments and devices, microelectronics, xerography, drug
delivery and pharmaceutical, and atmospheric dispersion.
Let us now see why thermophoresis occurs. At the molecular level, the gas
molecules coming from the hot region carry more momentum than those moving
from the cold region. This imbalance results in a net force from hot region to a cold
8
region. It should also be noted that the entire imbalance in momentum transfer is only
because of the incoming  ux of molecules. The re?emission of molecules (outward
 ux) is assumed to be difiuse because of the perfect accommodation of spheres (soot
particles). Thermophoretic force can be determined by solving momentum transfer
and energy transfer equations using the distribution functions. These distribution
functions are governed by the Boltzman?s equation. In general, the solution for the
Boltzman?s equation is very di?cult unless for some limiting cases. These limiting
cases can be identifled with the help of Knudsen number. The Knudsen number is
deflned as the ratio of the gas mean free path l to a characteristic length scale of the
particle, Kn = l/a. The characteristic length scale is assumed to be the equivalent-
volume radius a of the particle.
2.2.3 Knudsen Number
The Knudsen number can be classifled into three regimes: the continuum regime
(Kn << 1), the transition regime (Kn ? 1), and the free-molecule regime (Kn >>
1). The Boltzmann equation is an integro-difierential equation and solution of this
equation is limited to few cases. The continuum based models are the Navier-Stokes
equations. Euler equations correspond to inviscid continuum limit which shows a
singular limit since the  uid is assumed to be inviscid and non-conducting. Euler
 ow corresponds to Kn = 0. The Navier-Stokes equations can be derived from the
Boltzmann equation using the Chapman-Enskog expansion. At Knudsen numbers
larger than 0.1 the Navier-Stokes equations break-down and a higher level of approx-
imation is obtained by carrying second order terms (in Kn) in the Chapman-Enskog
9
expansion. A special form of such an equation is called the Burnett equation, for
which the solution requires second-order accurate slip boundary conditions in Kn.
The Burnett equations and consistent second-order slip boundary conditions is sub-
ject to some controversy and a better way of solving high Knudsen number  ow is
through molecular based direct simulation techniques such as the Direct Simulation
Monte Carlo method (DSMC) [8].
Continuum Regime
Thethermophoresistheoryincontinuumregime(Kn<<1)canbebuiltuponthe
solution of the Navier-Stokes equations combined with the appropriate slip boundary
conditions. The boundary conditions are based on the temperature jump and thermal
creep. This accounts for the fact that the gas cannot be treated as a continuous
medium within a few mean free paths from the particle surface [8].
Transition Regime
In this range of Knudsen number (Kn?1), the theoretical solutions are the most
di?cult to obtain [8]. However, interpolation models were advocated by researchers
[9,10] for flnding results in transition regime. Brock?s near-continuum solution was
a widely used formula for interpolation until more accurate solutions to Boltzmann?s
equation became available [11]. In the limit Kn ! 1, Brock?s near continuum
solution got reduced to a form very similar (difiers only by a constant) to the free
molecular solution of Waldmann. Therefore, this formula was used for obtaining
results in the transition regime.
10
Free-molecular Regime
In the limit of Kn >> 1 the velocity distribution function of the incoming mole-
cules is not afiected by the presence of the particle. To calculate the thermophoretic
force, it is necessary to model the behavior in which the gas molecules are re ected
by the particle surface after collisions with the particle [8].
With the literature presented in the above sections, we can begin with the pro-
posed methodology. Chapter 3 will present details of the proposed methodology and
the algorithm employed. The chapter will be divided broadly into two major sections:
(1) two-particle model, and (2) aggregate model.
11
Chapter 3
The Model
3.1 The velocity distribution function
The molecules in a gas posses a velocity distribution function, which characterizes
the speed and direction of molecular travel. Typically the distribution is a function of
6 coordinates: 3 spacial coordinates (x, y, z) and three velocity coordinates (u, v,w).
Together, these 6 coordinates are referred to as phase space.
We can assume for now that the properties of the gas are uniform over space, so
that the distribution function only depends on the three velocity components. For
this case, the distribution, denoted as f, is deflned so that
f(u;v;w) du dv dw
is the number of molecules per unit volume that have velocities within du of u,
dv of v, and dw of w. In addition,
n =
Z 1
?1
du
Z 1
?1
dv
Z 1
?1
dw f(u;v;w) (3.1)
in which n is the number density (number of molecules per unit volume, units
of 1/m3). This shows that the distribution function must have units of number
density/velocity3, or s/m6. For an equilibrium gas (stationary and at a uniform tem-
perature), the distribution is given by the Maxwellian formula:
f = n?3=2C3
T
exp [?(u2 +v2 +w2)=C2T] (3.2)
12
in which CT = p2RT is the thermal speed of the molecules, with R being the speciflc
gas constant.
A polar coordinate system is often used to characterize the molecular velocity, with
du dv dw = c2 dc sin d d`
in which c is the molecular speed,
c2 = u2 + v2 + w2
and  and ` are the polar and azimuthal angles of the molecular trajectory;
u = c sin cos`
v = c sin sin`
w = c cos 
Averages (or moments) of the distribution function are obtained via
h?i = 1n
Z 1
c=0
Z ?
 =0
Z 2?
`=0
? f c2 dc sin d d` (3.3)
inwhich? issomemolecularquantity. Forexample, theaveragespeedisobtained
by setting ? = c. For the maxwellian distribution, this gives
< c > = 2p? CT =
q
8RT
?
and the mean kinetic energy, per unit mass, is obtained from setting ? = c2/2, or
< c > = 34C2T = 32RT
13
The  ux of a molecular quantity ? refers to the transport of ? across a surface
(real or imaginary) in the direction normal to the surface and per unit area of the
surface. It would have unit of ? /m2/s. Say that the normal to the surface points
in the z direction. The formula for the  ux is obtained from
j? =
Z 1
u=?1
Z 1
v=?1
Z 1
w=0
? f w du dv dw (3.4)
or, in polar coordinates,
j? =
Z 1
c=0
Z ?=2
 =0
Z 2?
`=0
? f c cos c2 dc sin d d` (3.5)
The number  ux is obtained by setting ? = 1, and for the maxwellian distribu-
tion,
jn = 12p? nCT
The momentum  ux in the z direction is obtained from ? = mw = mc cos (with m
the molecule mass), and
jmom = 14 mnC2T = 12?RT
in which ? = mn is the mass density. From the ideal gas law the momentum  ux is
equal to P/2 which represents the normal stress on the (imaginary) surface due to
molecules leaving the surface. If we calculated the  ux of momentum arriving at the
surface we would also get P/2, and the total stress would be P, as expected.
The net force acting on an object (a particle, for example) would be obtained by
integrating the net momentum  ux (incident and re ected) over the surface of the
particle.
14
3.2 Sampling of the distribution function
In our work we assume, beforehand, that we know the distribution function(s) of
the molecules incident on the particles (targets) and emitted from the particles (i.e.,
re ected). From this, we can compute the transfer of quantities (momentum, energy)
to and from the targets. However, the targets will possess a geometry which would
make di?cult an analytical evaluation of the integrals per the previous formulas.
Consequently, we will numerically compute the transport of momentum and energy
via a monte carlo method.
The MC method is conceptually straightforward. We simulate molecular tra-
jectories and observe (on a computer) how the molecules interact with the target.
By collecting averages over the simulation of a large number of molecules, we can
determine the net rate of momentum and energy transfer to the target.
A required element to implementing the MC method is the sampling of a dis-
tribution function. For example, we may know that the molecules which we are
simulating have a velocity that is described by a given velocity distribution function.
We want to randomly assign values of velocity to individual molecules so that, when
averaged over a large number of molecules, the velocities fall into a distribution that
is consistent with our modelled distribution.
The sampling approach that we follow for this thesis would be cumulative dis-
tribution approach. Say our distribution is a function of one variable, x, and that x
runs from 0 to 1 (x is arbitrary here; it does not have a physical interpretation). The
15
cumulative distribution function, denoted g(x), is deflned by
g(x) =
Rx
0 f(t) dtR
1
0 f(t) dt
in which t is a dummy variable of integration. The cumulative distribution g(x)
represents the probability of choosing a sample from f that is between 0 and x. If x
= 1 then g = 1, i.e., we are certain that x lies between 0 and 1. On the other hand,
if we set g = 1=2 we would get an equation for xm (referred to as the mode of the
distribution):
1=2 = g(xm) ! xm = g?1(1=2)
in which g?1 is the inverse of the cumulative distribution function. The interpretation
of xm is that it is equally likely (50/50 chance) that a sampled x will lie either above
or below xm. This does not imply that xm = 1/2; the value of xm depends on the
form of the distribution function.
To construct a sampling scheme, we flrst note that all values of g are equally
likely. That is, from a probabilistic point of view, g will be uniformly distributed
between 0 and 1. We can then set g equal to a uniform random number between 0
and 1 to obtain
g(x) = R
in which R is the random number between 0 and 1. The sampled x is obtained from
the inverse function;
x = g?1(R)
16
3.3 The two sphere model
A basic starting model for our analysis is as follows. A pair of identical spheres
of radii a are separated by a distance R. The spheres are contained in a gas, and the
mean free path of the gas molecules is signiflcantly larger than the largest geometrical
length. The gas is at temperature T0 and pressure P0. The spheres are at a uniform
temperature TS < T0; this temperature difierence is maintained by radiation heat
emission from the spheres. We wish to calculate the following:
1. The rate of radiation heat transfer necessary to maintain the given temperature
difierence T0 ? TS. In principle, we would know the heat transfer rate from
an analysis of radiation and from this we would calculate TS. However, the
problem will be somewhat simplifled if we assume that TS is given.
2. The net force acting between the two spheres as a result of difierences in the
molecular momentum  ux on the sphere surfaces.
Assumptions we will make are
1. The spheres are stationary: they are held in place by some invisible means.
This is a contrived problem, but we need to start somewhere.
2. The gas is stationary: no bulk motion.
3. The surfaces of the spheres have perfect momentum and energy accommodation.
All previous history of the molecular trajectory is lost upon collision with the
surface. Re-emission of the molecules is difiuse.
17
Figure 3.1: Two?sphere model showing the ray trace of a molecule
4. Free molecular limit conditions: Kn  1. There are no molecule-molecule colli-
sions. The soot particles we consider would be nearly nano-sized, making them
considerably smaller than the mean free path. This assumption is supported
by the gas temperatures we would be dealing with (? 1500 K). At such high
temperatures, the mean free path increases, thus justifying our assumption of
Kn >> 1 (i.e., the mean free path >> the particle size). Typically, the soot
particles are 0:01 ?m in size, whereas the mean free path at STP is of the order
0:065 ?m. At around a temperature of 2000 K, it is of the order of 0:3 ?m.
5. Theincomingmoleculesfromthebackgroundgasarecharacterizedbyamaxwellian
velocity distribution at temperature T0 and number density n0.
6. There-emittedmolecules(i.e., re ectedmolecules)arecharacterizedbyamaxwellian
distribution at TS and number density nS.
18
3.3.1 Energy transport
Label the spheres 1 and 2 (because of the symmetry of the problem, condi-
tions/transfer rates will be identical for both spheres, but we will distinguish between
the spheres anyway). The net rate of molecular kinetic energy transfer to sphere 1,
due to molecular collisions, is
_Qm;1 =
Z
A
Z
?
Z 1
c=0
(fin ?fS) 12m c2 ^n:cc2 dc d? dA (3.6)
in which m is the molecular mass, ^n the outward normal, c the molecular velocity
vector, A the surface area of the sphere, and d? = sin d d` is a difierential solid
angle. The velocity distribution functions fin and fS describe the incoming and
emitted (re ected) molecules, respectively.
The integrals in 3.6 cannot be trivially performed analytically because fin will be
a function of incoming direction and surface position. For positions on the hemisphere
facing the opposite sphere, part of the incoming molecules will originate from the
opposite sphere, and the remainder will originate from the background gas.
When the incoming direction points towards the opposite sphere the incoming
distribution will be fin = fS; this is because we assume that the temperatures of both
spheres are identical. Consequently, the integrand in 3.6 will vanish for directions ?
pointed towards the opposite sphere. For direction pointed towards the background
gas we will have fin = f0, and 3.6 becomes
_Qm;1 =
Z
A
Z
?0
Z 1
c=0
(f0 ?fS) 12m c2 ^n:cc2 dc d? dA (3.7)
19
in which ?0 denotes the directions which point towards the background gas; this will
be a function of position on the sphere.
The distribution functions are given by
fS = nS?3=2C3
T;S
exp (?c2=C2T;S) (3.8)
f0 = n0?3=2C3
T;0
exp (?c2=C2T;0) (3.9)
with CT;S = p2RTS and likewise for CT;0.
The integral over speed c can be performed analytically in 3.7, and the remaining
integrals over direction and position will deflne a conflguration factor F1?0 so that
_Qm;1 = m
2?1=2 (n0C
3
T;0 ?nSC
3
T;S) 4?a
2 F1?0 (3.10)
in which
4?a2 F1?0 = 1?
Z
A
Z
?0
^n:sd? dA (3.11)
with s being a unit vector which points in all directions except towards the opposite
sphere. The quantity F1?0 is identical to the radiation conflguration factor which is
used in radiation exchange problems; it represents the fraction of all emitted molecules
from 1 which travel to the background gas. By summation we have F1?0 + F1?2 = 1,
i.e., the emitted molecules either end up in the background gas or on sphere 2. F1?0
is a function only of the distance between the centers of the two spheres.
20
Equation 3.10 contains the unknown quantity nS, which can be determined by
application of mass transfer principles. The net molecular  ux at any point on the
surface must be zero, i.e., the incoming  ux must balance the outgoing  ux. This
leads to
Jn = 0 =
Z
?
Z 1
c=0
(fS ?fin) ^n:cc2 dc d? (3.12)
Again, theincomingdistributiondependsonthepositiononthesphere. Asbefore, the
integral over direction ? can be split into a fraction ?S which includes all directions
which point to the neighboring sphere, and ?0 which encompasses all directions which
point towards the background gas. For the former case the incoming distribution
function will be fS and for the latter the distribution will be f0. The integral is then
Jn = 0 =
Z
?0
Z 1
c=0
(fS ?f0) ^n:cc2 dc d? (3.13)
Since this result must hold at every point on the surface - and since ?0 is a function
of surface position - it follows that
Z 1
c=0
(fS ?f0) c c2 dc = 0 (3.14)
The integrals can be computed analytically, which gives
nSCT;S = n0CT;0 (3.15)
which provides a relation between the number densities and temperatures of the
incoming and emitted molecules.
21
By putting the above into 3.10, we get
_Qm;1 = ?0CT;0
2?1=2 (C
2
T;0 ?C
2
T;S) 4?a
2 F1?0 (3.16)
= ?0RCT;0?1=2 (T0 ?TS) 4?a2 F1?0 (3.17)
3.3.2 Momentum transfer
The force applied to the sphere due to the emitted molecules will be zero, due
to the fact that the emitted molecules are emitted isotropically (uniform in all direc-
tions). The net force is therefore due to the nonuniformity in the incoming molecules,
and will be given by
F1 =
Z
A
Z
?
Z 1
c=0
fin mc(^n:c) c2 dc d? dA (3.18)
As before, the integral over direction ? can be split into ?S and ?0 and fin will have
the corresponding value of fS or f0. We can then use
Z
A
Z
?0
Z 1
c=0
f0 mc(^n:c) c2 dc d? dA = (
Z
A
Z
?
Z 1
c=0
f0 mc(^n:c) c2 dc d? dA
?
Z
A
Z
?S
Z 1
c=0
f0 mc(^n:c) c2 dc d? dA)
= ?
Z
A
Z
?S
Z 1
c=0
f0 mc(^n:c) c2 dc d? dA
22
The integral over all directions ? in the above is zero because the momentum
 ux for this part is isotropic. We then get
F1 =
Z
A
Z
?S
Z 1
c=0
(fS ?f0) mc(^n:c) c2 dc d? dA (3.19)
with the integral over direction limited to directions which point from 1 to 2.
The integrals over speed can be performed analytically, to give
F1 = 3m8? (nSC2T;S ?n0C2T;0) 4?a2 G1?2 (3.20)
in which G1?2 is a vector which depends only on the geometry. For a pair of spheres,
components will be aligned with axis of symmetry.
4?a2 G1?2 = 1?
Z
A
Z
?S
s(^n:s) d? dA (3.21)
Using 3.15 in 3.20 gives
F1 = 3?CT;08? (CT;S ?CT;0) 4?a2 G1?2 (3.22)
3.3.3 Soot radiation heat source function
This section presents a simplifled derivation of the rate of thermal emission from
a carbonaceous soot particle. The derivation begins with correlations from [12] on
the absorption properties of a soot cloud, and backs out the rate at which a single
particle will emit radiation.
23
The rate of emission from a particle cloud, per unit volume of the cloud, will be
q000e = 4 ?p eb; W=m3 (3.23)
where ?p is the Planck mean absorption coe?cient, deflned by
?p =
Z 1
0
?? eb? d? (3.24)
where ? is wavelength and
eb? = C1?5 (exp (C
2=?T)?1)
; eb =
Z 1
0
eb? d? =  T4
where  = 5:67 ? 10?8 W=(m2:K4), are the spectral and total blackbody emissive
power functions.
The spectral absorption coe?cient ?? for carbon soot can be approximated as
?? t 7f? ; ?m?1
where f is the soot volume fraction (volume of the solid particle per unit volume of
medium; a dimensionless quantity). Replacing this into 3.24 and integrating gives
?p = 1:86 ? 103 f T; m?1
24
Using the above formulas, the emissive sink for the soot cloud will be
q000e = 4 ?p eb = Ck f T5
where Ck = 4:23?10?4 W=m3K5.
Assuming all particles are identical, the rate of emission per unit cloud volume
would be related to the cloud number density by
q000e = Qen
in which n is the number density of the particles and Qe is the rate of emission from
a single particle. Likewise, the volume fraction of the cloud would be given by
f = V n
with V being the particle volume. Using these relations results in
Qe = Ck V T5
with V given in m3. If we assume that the particles are spherical with radius a, we
get
Qe = 4?3 Ck a3 T5 (3.25)
25
From 3.16 the rate of molecular energy transfer to the particle was
_Qm;1 = ?0CT;0
2?1=2 (C
2
T;0 ?C
2
T;S) 4?a
2 F1?0 (3.26)
t ?0C
2
T;0
?1=2 (CT;0 ?CT;S) 4?a
2 F1?0 (3.27)
the last line using a flrst order approximation for small CT;0 ? CT;S . It should be
noted here that the rationale for assuming CT;0 ? CT;S is that TS has a value of
around 1499.54 K when T0 = 1500 K at a = 15nm. Hence the perturbation on the
spheres can be assumed to be negligible. By setting Qm;1 = ?Qe, we get
CT;S ?CT;0 =
p? Q
e
?0C2T;04?a2F1?0 (3.28)
The force acting between a pair of spheres was
F1 = 3?CT;08? (CT;S ?CT;0) 4?a2 G1?2
or, using the previous equation, 3.28,
F1 = 3Qe8p?C
T;0
G1?2
F1?0 (3.29)
Using the formula for Qe, 3.25, the above result shows that 1) the force will scale
with a3 (i.e., proportional to the particle volume), and 2) the force will scale with
T4:50 , and G1?2 will scale (asymptotically) as (a=r1?2)2.
26
3.3.4 Thermophoretric force analysis
The thermophoretic force analysis can be borrowed from thoroughly-investigated
phenomenon of coagulation between electrically charged particles [12] as discussed
by [13]. The force existing between the electrically charged particles follows an
inverse-square law and is an instantaneous force. However, the case of efiective
thermophoretic force is not precisely equivalent to the instantaneous force in elec-
trostatics [13]. This is because the interaction of the thermophoretic force takes place
through the carrier gas and an inverse-square relationship would hold good only when
the characteristic time of gas heat transfer propagation is considerably smaller than
the characteristic time of the particle motion [13]. In the present analysis, however,
an inverse-square relationship can be approximated taking into account that we are
dealing with ?m? sized particles, the ratio of gas to particle characteristic times for
which is of the order of 0.1 - good enough for the above mentioned approximation [13].
From the electrical analogy, the ratio of coagulation rate constants between the
two particles with thermophoresis to that due to Brownian motion alone, is expressed
[12] as
Z = yey ?1
where
y = F12 r12K
B Tg
where r12 is the separation distance at contact, KB = 1:3805?10?23 J=K is Boltz-
mann?s constant, Tg is the temperature of the gas, and F12 is the force at contact
27
between the spheres given [12] by
F12 = ?? e
2
r212 (3.30)
where ? and ? are the number of the elementary electrical charges (electron charges),
e. In the above equation
Comparing3.29and3.30, i.e., replacingelectrostaticforce(F12)withthermophoretic
force (F1), we have for y:
y = 2? 3Qe8p?C
T;0
G1?2
F1?0
r12
KB Tg (3.31)
with G1?2 and F1?0 evaluated at contact, i.e., at r12 = 2a. Using the expression
for Qe from 3.25 in 3.31 and 3.29 and substituting the values of all constants, we get
y = ?(4:5212?1018) a4 T3:5g (G1?2F
1?0
) (3.32)
and
FT = 2?F1 = (3:1208?10?5) a3 T4:5g (G1?2F
1?0
) (3.33)
respectively. F1 is multiplied by a factor 2 to get FT because F1 is the force acting
on a single sphere and FT is the total thermophoretic force between the two spheres.
The units of a and FT in the above equation are m and N respectively.
28
3.3.5 Monte Carlo method
Thecomputational studyofthisthesisisdone inMATLAB6.5version. Thebasic
outline of the algorithm we are going to implement using a monte carlo technique is
1. to release computational gas molecules from random points on the surface of
sphere 1 in random directions
2. track the molecule, determine where it ends up
3. repeat the code several times
4. compute F1?0 and G1?2 and
5. hence compute y and FT.
Sampling the position on sphere surface
To start with the monte carlo technique, we flrst sample the position on the
surface of sphere 1 from which molecules are released. The pseudo-random number
generator of MATLAB is used for all random sampling purposes in this thesis.
The polar and azimuth angles of the position on the surface are sampled from
the cumulative distribution functions represented by
R1 =
Z cos 
0
d(cos )
R2 = 12?
Z `
0
d`
29
where R is the random number generator between 0 and 1. Upon inversion, we have
cos 1 = R1 (3.34)
`1 = 2?R2 (3.35)
We have sampled cos 1 from 0 to 1, and not -1 to 1 because, the molecules
released from the hemisphere of sphere 1 not facing the second sphere do not end up
on it anyway.
Sampling the direction of molecules
The polar and azimuth angles of the direction of the molecule released from the
sphere surface are represented by
R3 = 2
Z cos 
0
cos d(cos )
R4 = 12?
Z `
0
d`
Upon inversion, we have
cos 2 = pR3 (3.36)
`2 = 2?R4 (3.37)
30
Determination of collision point
The method followed to determine the collision point (or, for that matter, if the
molecule collides with the sphere 2 or not) by [14] is employed here.
The distance of the molecule starting point (x0;y0;z0) to the center of sphere 2
(xs;ys;zs) is
rs0 = ((xs ?x0)2 + (ys ?y0)2 + (zs ?z0)2)1=2 (3.38)
If rms is the distance between the molecule, at any point in its trajectory, and
the center of sphere 2, then the minimum value of rms will occur when rms:^c = 0 [14].
) rms = rs0 sinfi (3.39)
where fi is the angle between the molecular trajectory and the position vector of the
center of sphere 2 relative to the starting point, i.e.,
cosfi = 1r
s0
[(xs ?x0) ^u + (ys ?y0) ^v + (zs ?z0) ^w] (3.40)
where (^u;^v; ^w), the trajectory of the molecule is given by [14]
^u = cos`1(sin 2cos`2cos 1 +cos 2sin 1)?sin`1sin 2sin`2 (3.41)
^v = sin`1(sin 2cos`2cos 1 +cos 2sin 1)?cos`1sin 2sin`2 (3.42)
^w = cos 2cos 1 ?sin 2cos`2sin 1 (3.43)
31
If rms ? a, then the molecule has collided with the sphere 2 [14].
The molecules are released one?by?one and the above process is repeated several
number of times. Finally, the conflguration factor, F1?2, is obtained by dividing the
number of molecules that collided with the sphere 2 by the total number of molecules
released. The subsequent calculations of y and FT are carried out by using equations
3.32 and 3.33.
3.4 Aggregate model
In this case, we consider the interaction between a sphere and a cluster of spheres
(of identical radii).The cluster is synthesized based on the power-law, 2.1. The algo-
rithm used for this synthesis is the one used in [14]. The program is run and a cluster
of 25 spheres is generated.
Sampling of the distribution function for the aggregate model is the same as
that done for the two?sphere model. The equations for energy transport, momen-
tum transfer, the soot radiation heat source function, and the thermophoretic force
analysis remain the same with the two?sphere model directly extended to multiple
sphere, i.e., with G1?2 ?! G1?i.
The only modiflcation for the aggregate model would appear in the monte carlo code
written for the two?sphere model.
3.4.1 Modifled Monte Carlo method
The algorithm of the modifled monte carlo technique would be
32
Figure 3.2: Soot aggregate generated using the power law
1. release of computational gas molecules from random points on the surface of
sphere 1 in random directions
2. track the molecule, determine if it ends up on the cluster
3. calculate the overall conflguration factor, F1?cluster by repeating the code several
times. The overall conflguration factor is calculated as F1?cluster = PNi=1F1?i.
4. compute G1?cluster and
5. hence compute ycluster and FTcluster.
Sampling the position on sphere surface and sampling the direction of molecules
are done in the same fashion as that for the two?sphere model. The determination
of the collision point (flgure 3.3) on the cluster is also done in a more advanced
33
Figure 3.3: Aggregate?model showing the ray trace of a molecule
fashion, noting the complex geometry that is involved in this case. The procedure
is as follows: When a molecule is released from sphere 1, its trajectory is followed
and the perpendicular distances from the individual spheres in the cluster to this
trajectory are determined. All the distances that are less than the sphere radius,
a are noted. Now, using the law of cosines the distance the molecule has travelled
at the collision point, sc is determined for all spheres that met the above condition.
Among the obtained values the one which has the minimum sc is identifled to be the
sphere with which the molecule has collided.
34
Chapter 4
Results and Conclusions
4.1 Results
The monte carlo code for both the two?sphere model and the aggregate model
was written in MatLab 6.5. The program was run on a PC with Intel Pentium 4
processor, 1.4 GHz and 512 MB of RAM.
4.1.1 The two sphere model
The two spheres under consideration had the same radius (= 0.01 ?m). The tem-
perature of the gas was taken to be 1500 K. The monte carlo code for the two?sphere
model comprised of release of 250,000 computational molecules from sphere 1. The
program took around 4 minutes to give the output, i.e., the force between the spheres,
the coagulation ratio, and the plot between the force and the distance between the
two spheres.
Plotted in flg. 4.1 is the thermophoretic force, in N, between the two spheres
versus the distance, in m, between them.As it can be seen from the flgure, the force
between the two spheres follows an inverse?squared variation with the distance.
The plot in flg. 4.2 represents the plot shown in flg. 4.1 on a log?log scale.
The plot has a slope of ?2. The distortion of the plot towards the larger values of
distance is the ?noise?. This should be expected because for higher values of distance
between the spheres, the conflguration factor, F1?2, identically approaches zero and
the monte carlo technique in such a case would fail.
35
0 0.5 1 1.5 2 2.5 3 3.5
x 10?7
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5x 10?16
r (m)
F (N)
Figure 4.1: Two sphere model ? Force vs Distance
The coagulation ratio, Z, for this model was observed to be around 1.0002, for
a sphere radius of 0.01 ?m and a gas temperature of 1500 K.
4.1.2 Aggregate model
The main task here would be to decide on the number of computational molecules
to be released from sphere 1 on to the cluster. For this, we flrst collect samples
comprising of, say, 100 molecules. The conflguration factor from sphere 1 to the
cluster (F1?cluster) is determined for each such sample. Every time, a group of 10
such samples is taken and the ratio of standard deviation to the mean of that group
is calculated. The flrst group would consist of samples 1 through 10, the next one
consists of samples 2 through 11, and so on. The algorithm is written in such a
36
10?8 10?7 10?610
?18
10?17
10?16
10?15
log(r)
log(F)
Figure 4.2: Log?Log scale of Fig.4.1
way that the release of molecules will be terminated when the above mentioned ratio
becomes less than or equal to 0.0005. This means that we are allowing 99:95% of
accuracy.
Fig. 4.3 shows the plot of the conflguration factor, F1?cluster versus the number of
samples of the computational molecules, Nsamp, released. In the monte carlo code, we
have set the total number of molecules to be released as 100000, but the computation
stops when the total number of molecules reach a value of around 53000 because the
desired accuracy level is reached by this number, thus saving the computation time.
A cluster consisting of 25 spheres was generated, all with same radii as that of
the sphere 1 (= 0.01 ?m). The temperature of the gas was taken to be 1500 K. The
accuracy for (F1?cluster) was set to 99:9% and 100,000 computational molecules from
37
0 100 200 300 400 500 6000.04
0.045
0.05
0.055
0.06
0.065
Nsamp
F 1?cluster
Figure 4.3: F1?cluster vs Nsamp
sphere 1 were released. The position of sphere 1 relative to the cluster is randomly
chosen around the cluster. The program took around 25 minutes to give the output,
i.e., the force between the spheres, the coagulation ratio, and the plot between the
force and the distance between the two spheres.
Fig. 4.4 shows the variation of the thermophoretic force in N between sphere 1
and the cluster with the distance in m between sphere 1 and the center of mass of
the cluster. As was the case with the two sphere model, the force in the aggregate
model bears an inverse?squared relationship with the distance. The magnitudes of
the force in the two cases are, however, difierent.
The plot in flg. 4.5 represents the plot shown in flg. 4.4 on a log?log scale. The
trend is same as that found in the case of a two?sphere model.
38
0.5 1 1.5 2 2.5 3 3.5 4
x 10?7
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5x 10?16
r (m)
F cluster
 (N)
Figure 4.4: Aggregate model ? Force vs Distance
Since y depends on r12, the separation between the cluster and the sphere 1 at
contact, we flrst deflne this quantity for the aggregate model. r12 = (the distance
between the center of mass of the aggregate and the center of the farthest sphere in
the aggregate + the radius of the farthest sphere + the radius of sphere 1). It should
be noted that sphere 1 in the simulation need not necessarily be in contact with the
farthest sphere when r12 is calculated because its position is chosen at random. The
coagulation ratio, Z, calculated for the aggregate model was observed to be around
1.0010, for a 25 sphere cluster with an individual sphere radius of 0.01 ?m and a gas
temperature of 1500 K.
39
10?7
10?16
log(r)
log(F
cluster
)
Figure 4.5: Log?Log scale of Fig.4.4
40
4.1.3 Two?sphere equivalent of aggregate model
Table 4.1: Comparison of thermophoretic force magnitudes
r F1?cluster Eq.Volume Eq.Surface area Eq.Radius
8.0487 0.3406 0.4802 1.485 3.633
10.0487 0.2397 0.2917 0.916 1.986
12.0487 0.1784 0.2070 0.610 1.321
14.0487 0.1388 0.1521 0.456 0.930
16.0487 0.0994 0.1184 0.346 0.710
18.0487 0.0929 0.0946 0.271 0.557
20.0487 0.0757 0.0753 0.214 0.450
22.0487 0.0559 0.0649 0.180 0.364
24.0487 0.0469 0.0502 0.145 0.303
26.0487 0.0442 0.0421 0.133 0.255
28.0487 0.0369 0.0397 0.118 0.220
30.0487 0.0326 0.0318 0.091 0.186
32.0487 0.0305 0.0290 0.080 0.166
34.0487 0.0272 0.0272 0.073 0.151
36.0487 0.0275 0.0250 0.068 0.142
Table 4.1 shows the comparison of the thermophoretic force,F1?cluster, in the
aggregate model to the thermophoretic force in the two?sphere model when the
sphere 2 is replaced with a sphere having (1) equivalent volume, (2) equivalent surface
area, and (3) equivalent radius (radius of gyration) as that of the cluster in the
aggregate model. The flrst column (r) of the table corresponds to the increasing
distance between the centers of the cluster and sphere 1, multiplied by a factor of
10?8m. All values shown in the flgure are multiplied by a factor of 10?15N.
41
4.1.4 Coagulation ratio
Let us flrst look at the variation of y, found in the expression for Z, with a, T,
and N (in the case of aggregates):
y / N V T3:5 r12 (G1?2F
1?0
)
(G1?2F
1?0
) / ( a
2
r212) and V / a
3
) y / (N a
5 T3:5
r12 )
where N is the number of spheres in the aggregate, V is the volume of a single
sphere, T is the temperature, and r12 is the separation between the centers at contact.
In the case of two?sphere model, N = 1 and r12 = 2a. Therefore,
y / a4 T3:5
In the case of the aggregates,
r12 / Rg / N(
1
Df )
) y / N(1?
1
Df ) a4 T3:5
Coagulation ratio, Z, signifles how prominent the thermophoretic force between
the particles is, when compared to the Brownian motion alone. It would be an
42
interesting quantity to look observe at various particle sizes and gas temperatures as
in the combustion environments deal a wide range of these parameters.
Table 4.2: Coagulation ratios for various a and T
a!
T# 0.01 0.02 0.03
1500 1.0002 1.0032 1.0163
1800 1.0004 1.00360 1.0304
2000 1.0006 1.0089 1.0461
Shown in Table 4.2 are the coagulation ratios obtained from the two?sphere
model for various values of particle sizes (?m) and gas temperatures (K).
Table 4.3: Coagulation ratios for various a and T
a!
T# 0.01 0.02 0.03
1500 1.0008 1.0116 1.0665
1800 1.0013 1.0229 1.1179
2000 1.0021 1.0327 1.1722
Shown in Table 4.3 are the coagulation ratios obtained from the aggregate model
for various values of particle sizes (?m) and gas temperatures (K).
4.2 Conclusions
A synthetic simulation model has been developed to generate soot targets, and
computational molecules and a monte carlo method has been used to calculate (1)
the thermophoretic force between cluster?sphere and sphere?sphere conflgurations,
and (2) coagulation ratios. Our results indicate that, at a = 0:01 ?m and T =
1500 K, the coagulation ratio in the two?sphere model is very negligible. However,
43
the coagulation ratios in the case of an aggregate model could be signiflcant with
increasing size of the aggregates (recalling the parameter?y dependency on N). It
could also be observed, by looking at Table 4.3, that the coagulation ratios for the
aggregate model have signiflcantly large magnitudes at larger particle sizes and higher
gas temperatures.
Noting the order of increase in magnitude of the coagulation ratios from the
two?sphere model to the aggregate model, it could be extrapolated that the efiect
of radiation on thermophoretic motion of the soot particles would be even more
prominent in the case of a cluster?cluster interaction.
The results also show that the sphere aggregates can be approximated by their
volume equivalent spheres. We can support the above statement, primarily, by noting
that the force scales as the volume of the spheres. However, the magnitude of the
force (and hence the comparison to the volume equivalent spheres) in the aggregate
model largely depends on the position of sphere 1 relative to the cluster. A lot of
work is yet to be done in this aspect to reach to any concluding remarks.
4.3 Suggestions for future research
The future research in the lines of the work presented could involve the dynamic
motion of the soot aggregates in the simulation model.
A quasi?random number generator could be used in the monte carlo method
implemented instead of the pseudo?random number generator used.
44
The interaction between a cluster and another cluster could also be an interesting
study to carry out taking into account the results presented in the current work for
the coagulation ratios of larger?sized particles.
45
Bibliography
[1] H.F. Calcote. Mechanisms of soot nucleation in  ames-a critical review. Combust.
Flame, 42:215-242, 1981.
[2] S.J. Harris and A.M. Weiner. Chemical kinetics of soot particle growth. Annu.
Rev. Phys. Chem., 36:31-52, 1985.
[3] Z. Dai, S.S. Krishnan, K.C. Lin, R. Sangras, J.S. Wu, and G.M. Faeth. Mix-
ing and radiation properties of buoyant luminous  ame environments. BFRL
Publications, 1997.
[4] K.M. Butler and G.W. Mulholland. Generation and transport of smoke compo-
nents. Fire Technology, 40(2):149-176, April 2004.
[5] T.A. Witten. Difiusion?limited aggregation, a kinetic critical phenomenon.
Phys. Rev. A, 47:1400-1403, 1981.
[6] P. Meakin. Difiusion controlled cluster formation in two, three, and four dimen-
sions. Phys. Rev. A, 27:604-607, 1983.
[7] C.M. Sorensen and G.C. Roberts. The prefactor of fractal aggregates. J. Colloid
Interface Sci., 186:447-452, 1997.
[8] F. Zheng. Thermophoresis of spherical and non-spherical particles: a review of
theories and experiments. Adv. Colloid Interface Sci., 97:255-278, 2002.
[9] J.R. Brock. The thermal force in the transition region. J. Colloid Interface Sci.,
23:448-452, 1967.
[10] J.R. Brock. Experiment and theory for the thermal force in the transition region.
J. Colloid Interface Sci., 25:392-395, 1967.
[11] L. Talbot, R.K. Cheng, R.W. Schefer, and D.R. Willis. Thermophoresis of par-
ticles in a heated boundary layer. J. Fluid Mech., 101:737-758, 1980.
[12] G. Zebel. Aerosol Science. Academic Press, NY., 1966.
[13] D.W. Mackowski, M. Tassopoulos,and D.E. Rosner. Efiect of Radiative Heat
Transfer on the Coagulation Dynamics of Combustion-Generated Particles.
Aerosol Sci. Technol., 20(1):83-95, 1994.
46
[14] D.W. Mackowski, V. Nayagam, and P.B. Sunderland. Monte?Carlo Simulation
of Hydrodynamic Drag and Thermophoresis of Fractal Aggregates of Spheres in
the Free?Molecule Flow Regime. J. Aerosol Sci., 2005.
[15] R.L. Vander Wal. Soot Precursor Material: Spacial Location via Simultane-
ous LIF?LII Imaging and Characterization via TEM. NASA Contract Report
198469, May 1996.
47
Appendices
48
Appendix A
Monte Carlo code for the two-sphere model
clearall;
l = 1?power(10;?8);
matrix = [2?l;4?l;6?l;8?l;10?l;12?l;14?l;16?l;18?l;20?l;22?l;24?l;26?
l;28?l;30?l];
for initial = 1 : 15;
a = 1?power(10;?8);
t = 1500;
c = matrix(initial);
theta = acos(2?rand?1);
phi = 2?pi?rand;
xs = 0;
ys = 0;
zs = c;
ni = 0;fx = 0;fy = 0;fz = 0;
n = 50000;
for i = 1 : n;
theta1 = acos(2?rand?1);
phi1 = 2?pi?rand;
theta2 = acos(sqrt(rand));
phi2 = 2?pi?rand;
49
x0 = a?sin(theta1)?cos(phi1);
y0 = a?sin(theta1)?sin(phi1);
z0 = a?cos(theta1);
rs0 = sqrt(abs((xs?x0)):2 +abs((ys?y0)):2 +abs((zs?z0)):2);
u = cos(phi1)?(sin(theta2)?cos(phi2)?cos(theta1) + cos(theta2)?sin(theta1))?
(sin(phi1)?sin(theta2)?sin(phi2));
v = sin(phi1)?(sin(theta2)?cos(phi2)?cos(theta1) + cos(theta2)?sin(theta1)) +
(cos(phi1)?sin(theta2)?sin(phi2));
w = cos(theta2)?cos(theta1)?(sin(theta2)?cos(phi2)?sin(theta1));
alpha = acos(((xs?x0)?u+(ys?y0)?v +(zs?z0)?w):=rs0);
rs = rs0:?sqrt(1?(cos(alpha)):2);
if cos(alpha) >= 0
if rs <= a
ni = ni+1;
fx = fx+u;
fy = fy +v;
fz = fz +w;
end
end
end
f12 = ni=n;
fxnet = fx=n;
fynet = fy=n;
50
fznet = fz=n;
f = 1?f12;
gbyf(initial) = abs(fznet=(f));
end
y = ?2?1:1303?power(10;18)?a:3 ?t:(3:5)?gbyf(1):?a
z = (y):=(exp(y)?1)
force = 2?1:5604?power(10;?5)?a:3 ?t:(4:5)?gbyf
plot(log(matrix);log(force))
figure
plot(matrix;force)
51
Appendix B
Monte Carlo code for the aggregate model
clearall;
a = 1?power(10;?8);
t = 1500;
no = 20000;
data = xlsread(0test4:xls0);
nspheres = length(data);
for n = 1 : nspheres;
xs = data(:;1);
ys = data(:;2);
zs = data(:;3);
end
max = 0;kmax = 0;
for k = 1 : nspheres;
r = sqrt(abs(xs(k)):2 +abs(ys(k)):2 +abs(zs(k)):2);
if max == 0
max = r;
kmax = k;
else
if r > max
max = r;
52
kmax = k;
end
end
end
for j = 1 : 15;
fall(j) = 0;fxall(j) = 0;fyall(j) = 0;fzall(j) = 0;
end
ctheta = 2?rand?1;
stheta = sqrt(1?ctheta:2);
phi = 2?pi?rand;
x11 = stheta?cos(phi);
y11 = stheta?sin(phi);
z11 = ctheta;
for e = 1 : 15;
rmax(e) = max+2?e;
x1 = rmax(e)?x11;
y1 = y11?rmax(e);
z1 = z11?rmax(e);
y = 1;isamp = 1;err = 1;eps = 0:0005;i = 1;
global mmin;
for m = 1 : nspheres;
f(m) = 0;fx(m) = 0;fy(m) = 0;fz(m) = 0;
end
53
i = 1;
isamp = 1;y = 1;err = 1;eps = 0:001;
while i <= no and err > eps;
ctheta1 = (rand);
stheta1 = sqrt(1?ctheta1:2);
phi1 = 2?pi?rand;
ctheta2 = (sqrt(rand));
stheta2 = sqrt(1?ctheta2:2);
phi2 = 2?pi?rand;
x0 = x1+stheta1?cos(phi1);
y0 = y1+stheta1?sin(phi1);
z0 = z1+ctheta1;
min = 0;mmin = 0;
for m = 1 : nspheres;
rs0 = sqrt(abs((xs(m)?x0)):2 +abs((ys(m)?y0)):2 +abs((zs(m)?z0)):2);
u = cos(phi1) ? (stheta2 ? cos(phi2) ? ctheta1 + ctheta2 ? stheta1) ? (sin(phi1) ?
stheta2?sin(phi2));
v = sin(phi1)?(stheta2?cos(phi2)?ctheta1+ctheta2?stheta1)+(cos(phi1)?stheta2?
sin(phi2));
w = ctheta2?ctheta1?(stheta2?cos(phi2)?stheta1);
calpha = (((xs(m)?x0):?u+(ys(m)?y0):?v +(zs(m)?z0):?w):=rs0);
salpha = sqrt(1?calpha:2);
rms = rs0?sqrt(1?(calpha):2);
54
if calpha >= 0
if rms <= 1
sc = rs0:?calpha?(1?rs0:2:?salpha:2):(1=2);
if min == 0
min = sc;
mmin = m;
umin = u;
vmin = v;
wmin = w;
else
if sc < min
min = sc;
mmin = m;
umin = u;
vmin = v;
wmin = w;
end
end
end
end
end
if (mmin = 0)
f(mmin) = f(mmin)+1;
55
fall(e) = fall(e)+1;
fx(mmin) = fx(mmin)+umin;
fy(mmin) = fy(mmin)+vmin;
fz(mmin) = fz(mmin)+wmin;
fxall(e) = fxall(e)+umin;
fyall(e) = fyall(e)+vmin;
fzall(e) = fzall(e)+wmin;
end
if (mod(i;100) == 0)
samp(y) = fall(e):=i;
y = y +1;
holdon;
end
if (floor(i=100)) >= 10 and mod(i;100) == 0;
fsamp = samp(isamp : (isamp+9));
isamp = isamp+1;
if mean(fsamp) = 0
err = std(fsamp):=mean(fsamp);
actual = i;
end
end
i = i+1;
end
56
f = f:=actual;
fall(e) = fall(e):=actual;
fxnet = fx:=actual;
fynet = fy:=actual;
fznet = fz:=actual;fxnetall(e) = fxall(e):=actual;
fynetall(e) = fyall(e):=actual;
fznetall(e) = fzall(e):=actual;
fnet = sqrt(fxnet:2 +fynet:2 +fznet:2);
fnetall(e) = sqrt(fxnetall(e):2 +fynetall(e):2 +fznetall(e):2);
gbyf = abs(fnet=(1?f));
gbyfall(e) = abs(fnetall(e)=(1?fall(e)));
c = sqrt(abs(x1?xs):2 +abs(y1?ys):2 +abs(z1?zs):2)?power(10;?8);
yc = ?2?1:1303?power(10;18)?a:3 ?t:(3:5)?gbyf:?c;
ycall = ?2?1:1303?power(10;18)?a:3 ?t:(3:5)?gbyfall(1):?rmax(1)?a;
z = (yc):=(exp(yc)?1);
zall = (ycall):=(exp(ycall)?1)
force = 2?1:5604?power(10;?5)?a:3 ?t:(4:5)?gbyf;
forceall(e) = 2?1:5604?power(10;?5)?a:3 ?t:(4:5)?gbyfall(e)
end
plot(rmax;forceall)
figure
plot(log(rmax);log(forceall))
57