PARTICLE COLLISION RATE AND SMALL-SCALE STRUCTURE OF
PARTICLE CONCENTRATION IN TURBULENT FLOWS
by
Kevin Cunkuo Hu
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF
THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF
THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF
PHILOSOPHY
UNIVERSITY OF FLORIDA
1998
ACKNOWLEDGMENTS
I would like to take this opportunity to express my sincere thanks to my advisor,
Professor Renwei Mei, for his patience, encouragement, and help. This work would be
impossible without his guidance and insights. I appreciate my advisory committee
members, Professor James F. Klausner, Professor Corin Segal, Professor Ulrich H.
Kurzweg, and Professor Roger Tran-Son-Tay, for their suggestions and help. I also
express gratitude to Professor Michael W. Reeks (in England), Professor Ibrahim K.
Ebcioglu, Professor Otis Walton and Professor Nicolae D. Cristescu for their help.
Thanks to Dr. Guobao Guo and my fellow graduate students Xueliang Zhang,
Hong Shang, and others for their helpful discussions and friendship.
I acknowledge the financial support of the Engineering Research Center (ERC)
for Particle Science & Technology at the University of Florida, the National Science
Foundation (EEC-9402989), industrial partners of the ERC, and the ALCOA Foundation
Award.
I am deeply indebted to my wife, Rui, for her love, patience, and encouragement
during my thesis work, which were more than I could ever have asked.
TABLE OF CONTENTS
Page
ACKNOW LEDGM ENTS ........................................................................................ ii
LIST OF TABLES...................................................................................................vi
LIST OF FIGURES ................................................................................................ vii
ABSTRACT............................................................................................................ xii
CHAPTER 1 INTRODUCTION ........................................................................... 1
1.1 Smoluchowski's Prediction in Lam inar Shear Flow .......................................... 4
1.2 Various Theories for Collision Rate in Turbulence........................................... 5
1.3 Recent Development in Computer Simulation on Particle Collision Rate..... 7
1.4 Recent Development in Particle Concentration Non-uniformity ..................... 8
1.5 Scope of Present Study ........................................................................................ 11
CHAPTER 2 TURBULENCE AND PARTICLE COLLISION
DETECTION ......................................................................................................... 13
2.1 Turbulence Representation ................................................................................. 13
2.1.1 Isotropic, Gaussian Turbulence........................................................................ 13
2.1.2 Rapidly Sheared Homogeneous Turbulence .................................................... 16
2.1.3 Isotropic Turbulence Generation in a Periodic Box......................................... 17
2.2 Numerical Integration of Particle Motion and Error Analysis ....................... 19
2.3 Particle Collision Treatm ent ............................................................................... 21
2.3.1 Collision Detection Scheme ............................................................................. 21
2.3.2 Post-Collision Treatment.................................................................................. 24
2.4 Computer Simulation........................................................................................... 25
2.4.1 Time Step ......................................................................................................... 25
2.4.2 Turbulence Realization..................................................................................... 27
2.4.3 Validation of the Collision Detection Scheme ................................................. 27
iii
CHAPTER 3 A NEW THEORETICAL FRAMEWORK FOR PREDICTING
COLLISION RATE OF SMALL PARTICLES IN GENERAL TURBULENT
F L O W S .............................................................................................. ............................... 33
3.1 Introduction........................................................................................................... 33
3.2 Analysis................................................................................................................. 37
3.3 Numerical Simulation.......................................................................................... 43
3.4 Results and Discussions ...................................................................................... 44
3.4.1 Comparison with Simulation Results in a Laminar Shear Flow ...................... 44
3.4.2 Collision Rate and Velocity in a Gaussian, Isotropic Turbulence.................... 47
3.4.3 Collision in a Rapidly Sheared Homogeneous Turbulence.............................. 52
3.5 Summary and Conclusions.................................................................................. 55
CHAPTER 4 EFFECT OF INERTIA ON THE PARTICLE COLLISION
VELOCITY AND COLLISION ANGLE IN GAUSSIAN TURBULENCE .........59
4.1 Introduction .................................................................................................. 59
4.2 Saffnan & Turner's Analysis....................................................................... 62
4.3 Present Analysis ................................................................................................... 65
4.3.1 Asymptotic Analysis......................................................................................... 65
4.3.2 PDF-Based Method..........................................................................................68
4.4 Results and Discussions................................................................................ 69
4.5 Summary and Conclusions.................................................................................. 76
CHAPTER 5 QUANTIFICATION OF PARTICLE CONCENTRATION NON-
UNIFORMITY IN TURBULENCE............................................................................. 77
5.1 Introduction .................................................................................................. 78
5.2 Quantification Method.................................................................................. 79
5.2.1 Inherent Noise of the Concentration Non-uniformity...................................... 79
5.2.2 Concentration Non-uniformity Parameter........................................................ 82
5.3 Results and Discussions................................................................................ 89
5.4 Summary and Conclusions.................................................................................. 96
CHAPTER 6 COLLISION RATE AT FINITE PARTICLE INERTIA AND
SIZE IN GAUSSIAN TURBULENCE ........................................................................ 97
6.1 Introduction .................................................................................................. 98
6.2 Correction of Smoluchowski's Theory....................................................... 100
6.3 An Analysis on the Effect of Particle Inertia on the Collision Rate............ 106
6.3.1 Saffinan & Turner's Analysis ....................................................................... 106
6.3.2 Effect of Small Inertia on the Collision Kernel.............................................. 107
6.3.3 A Proposed Expression on Collision Kernel.................................................. 110
6.4 Comparison of the Analysis with Simulation Results.................................. 111
6.4.1 Effect of Particle Inertia................................................................................. 111
6.4.2 Effect of Particle Size..................................................................................... 113
6.4.3 Effect of Particle Mean Concentration........................................................... 114
6.4.4 Effect of Particle Gravitational Settling......................................................... 114
6.5 Collision Rate in Polydisperse Suspension..................................................... 115
6.5.1 Particle Size Spectrum.................................................................................... 122
6.5.2 Effect of Particle Inertia ................................................................................. 124
6.5.3 Effect of Particle Gravitational Settling ......................................................... 125
6.6 Collision Rate in Hard Sphere Suspension and Assessment of Collision
M odels ................................................................................................................. 127
6.7 Summary and Conclusions................................................................................ 129
CHAPTER 7 SUMM ARY AND CONCLUSIONS .............................................. 137
REFERENCES...................................................................................................... 140
BIOGRAPHICAL SKETCH ................................................................................ 146
LIST OF TABLES
Table 3.1 Comparison of the collision velocity among monosized
particles between prediction and direct numerical simulation
in lam inar shear flow................................................................................. 45
Table 3.2 Comparison of the collision velocity among monosized
particles between prediction and direct numerical simulation
in Gaussian, isotropic turbulence.............................................................. 50
LIST OF FIGURES
Figure 2.1 The effect of time step on the cumulative average of the normalized
collision kernel a', for inertialess particles in turbulence. D=0.05.............................. 26
Figure 2.2 Turbulence dissipation rate e against turbulence realization. (a) Variation of
ensemble average. (b) Cumulative average .................................................................. 29
Figure 2.3 Turbulence small eddy shear rate (e / v) 1/2 against turbulence realization.
(a) Variation of ensemble average. (b) Cumulative average......................................... 30
Figure 2.4 Inverse of total number of monodisperse particles against time in laminar
shear flow F=1.0, D=0.02............................................ ............................................... 31
Figure 2.5 Multiple collisions of HARD SPHERE particles in laminar shear flow...... 32
Figure 3.1 Sketch of particle collision velocity. The separation
distance is x2-xI = Rij= (ri+rj) ................................................................................... 35
Figure 3.2 Variations of J(), (Q), 74() with i where
fl= ), f2= G( ) and f3= ............................................................. ....... 36
Figure 3.3 Local coordinates for collision in laminar shear flow. 42
Figure 3.4 Comparison of the pdfs for the collision angles 9c between
the theory and the simulation in laminar shear flow ...................................................... 47
Figure 3.5 Convergence of a*2 as Nk -*oo ....................................................................... 50
Figure 3.6 Particle collision kernel a*, from direct numerical simulation in isotropic
turbulence. (a) Variation of the ensemble average. (b) Cumulative time average....... 51
Figure 3.7 (a) Normalized particle collision kernel in a highly sheared homogeneous
turbulence without contribution from the mean shear rate.............................................. 56
Figure 3.7 (b) Collision kernel in a rapidly sheared homogeneous turbulence. 6 denotes
the dissipation rate of the isotropic state at t=0 ............................................................... 57
Figure 3.8 Collision velocity in a rapidly sheared homogeneous turbulence.
so denotes the dissipation rate of the isotropic state at t=0 ............................................ 58
Figure 4.1 Sketch of particle collision scheme ............................................................... 63
Figure 4.2 Sketch of particle collision velocity and collision angles (012 & Oc ).
The separation distance is x2 x1 = R = (r, + r2)..................................................... 63
Figure 4.3 Comparison for the increased collision velocity due to the particle inertia
between numerical simulation and theoretical prediction, D/TI=0.957 .......................... 70
Figure 4.4 Effect of particle inertia on the probability density function of the particle
collision velocity wr. (a) D/I=0.957. (b) D/rh=2.522.................................................. 71
Figure 4.5 Effect of particle inertia on the probability density function of the particle
collision angle 012. (a) D/rl=0.957. (b) D/r9=2.522...................................................... 73
Figure 4.5 (c) Effect of particle inertia on the probability density function of the particle
collision angle 0 D/ri=0.957 ...................................................................................... 74
Figure 4.6 Effect of particle inertia on the mean collision angle 012 in isotropic,
G aussian turbulence ....................................................................................................... 74
Figure 4.7 Effect of particle inertia on the standard deviation of collision angle 012 in
isotropic, Gaussian turbulence ....................................................................................... 75
Figure 4.8 Effects of particle inertia and size on the probability density function of the
particle collision angle 012, pTk=1.023 ....................................................................... 75
Figure 5.1 The normalized variance of the concentration as function of cell size for
inertialess particles in isotropic turbulence. PTtk =1.023 .............................................. 83
Figure 5.2 The expressions (5.17) & (5.18) of the concentration non-uniformity as
function of cell size for inertialess particles in uniform shear flow. /F=oo ................. 84
Figure 5.3 The expressions (5.17) & (5.18) of the concentration non-uniformity as
function of cell size for inertialess particles in isotropic turbulence. PT k = ............... 86
Figure 5.4 The concentration non-uniformity parameter AP as function of cell size in
uniform shear flow. r=1.0, 03/F=l.0 ............................................................................ 87
Figure 5.5 The concentration non-uniformity as function of cell size for small inertia
particles in isotropic turbulence ..................................................................................... 88
Figure 5.6 The concentration non-uniformity as function of cell size for various
particle number densities at a given inertia in isotropic turbulence. p3T k =1.023........ 89
Figure 5.7 The concentration non-uniformity of sizeless particles as function of cell
size for various inertia in isotropic turbulence ............................................................... 91
Figure 5.8 The concentration non-uniformity parameter AD of finite-size particles as
function of cell size for various inertia in isotropic turbulence. (a) D/rn=0.957.
(b) D /il=2.52................................................................................................................... 93
Figure 5.9 The maximum concentration non-uniformity parameter AP of finite-size
particles as function of particle inertia in isotropic turbulence ...................................... 94
Figure 5.10 The concentration non-uniformity parameter AO as function of Ax/h for
small-size and zero-size particles in isotropic turbulence. P3 k =1.023 ......................... 94
Figure 5.11 The concentration non-uniformity of particles at a given inertia as function
of Ax/D cell size for various sizes in isotropic turbulence. (a) Log scale in x-axis.
(b) in norm al scale in x-axis ........................................................................................... 95
Figure 6.1 Collision kernel a ii against time in laminar shear flow. F=1.0, D=0.04.
(a) Variation of ensemble average. (b) Cumulative average...................................... 105
Figure 6.2 The effect of particle inertia on the increased collision kernel Aa 1 for
monodisperse particles for (a) D/rl=0.957. (b) D/rj=2.52 ........................................ 116
Figure 6.3 The effect of particle size on the collision kernel aot i for monodisperse
particles at a very dilute system. Ou k =oo................................................................... 117
Figure 6.4 The effect of particle size on the collision kernel a 1 for monodisperse
particles at 'not' very dilute system. p3- k = oo............................................................. 117
Figure 6.5 Variation of the collision kernel at 1 against particle number concentration
for monodisperse particles. 13PTk = oo, D/TI=0.435...................................................... 118
Figure 6.6 The effect of gravitational settling on the collision kernel a*i for
monodisperse particles. 3rk=1.023, D/rTI=2.52.......................................................... 118
Figure 6.7 The effect of gravitational settling for monodisperse particles. D/ri=2.52,
PT k = 1.023. (a) on the collision velocity. (b) on the concentration non-uniformity... 119
Figure 6.8 Collision kernel a41 of polydisperse particles against time. Prtk = O,
D/TI=0.957. (a) Variation of the collision kernel. (b) Cumulative average................. 120
Figure 6.9 Collision kernel cta of polydisperse particles with various inertia against
turbulence realization. D/rp=0.957. (a) Variation of the collision kernel.
(b) Cumulative average of the collision kernel............................................................ 121
Figure 6.10 Inverse of total number of particles against time for polydisperse particles.
D/TI=0.957. (a) primary size (v0) particles. (b) secondary size (2 v0) particles......... 122
Figure 6.11 The evolution of particle size distribution at various instants for
polydisperse particles. D/1=2.17, r k =1.023............................................................. 123
Figure 6.12 The effect of particle inertia on the increased collision kernel Aoal for
polydisperse particles for (a) D/il=0.957. (b) D/I=2.52 ........................................... 126
Figure 6.13 The effect of gravitational settling on the collision kernel a*, for
polydisperse particles. PTk=1.023, D/r1=2.52............................................................ 127
Figure 6.14 The effect of gravitational settling for polydisperse particles. D/rl=2.52,
PT k=1.023. (a) on the collision velocity. (b) on the concentration non-uniformity.... 128
Figure 6.15 Collision kernel a 1 of hard sphere particles with various inertia against
time in laminar shear flow. D=0.04. (a) Variation of the collision kernel.
(b) Cum ulative average................................................................................................. 130
Figure 6.16 The effects of the collision schemes on the collision kernel a I1 against time
in laminar shear flow. F=1.0, D=0.04. (a) Variation of the collision kernel. (b)
Cumulative average of the collision kernel................................................................... 131
*
Figure 6.17 The effects of the inertia on the collision kernel a11 against time for hard
sphere particles in Gaussian turbulence. D/Ti=0.957. (a) Variation of the collision
kernel. (b) Cumulative average of the collision kernel............................................... 132
Figure 6.18 The effects of the collision schemes on the collision kernel a* against time
in Gaussian turbulence. D/Ti=2.174, Prk= 1.023. (a) Variation of the collision kernel.
(b) Cumulative average of the collision kernel............................................................ 133
Figure 6.19 The effects of the collision schemes on the collision kernel a*,I of hard
sphere particles in Gaussian turbulence. D/il=0.957, Pt k =0.292. (a) Variation of the
collision kernel. (b) Cumulative average of the collision kernel................................ 134
Abstract of Dissertation Presented to the Graduate School of the University
of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor
of Philosophy
PARTICLE COLLISION RATE AND SMALL-SCALE STRUCTURE OF
PARTICLE CONCENTRATION IN TURBULENT FLOWS
by
Kevin Cunkuo Hu
May, 1998
Chairman: Renwei Mei
Major Department: Aerospace Engineering, Mechanics and Engineering Science
A study of particle collision and concentration non-uniformity in turbulent flows
is carried out. A correction to the classical result of Smoluchowski for the collision rate
of monodisperse particles in a laminar shear flow has been made. A theoretical
framework has been developed to evaluate the collision rate and collision velocity of very
small particles in general turbulent flows. The present approach differs significantly
from the classical approach of Saffinan & Turner in that the ensemble average is taken
after the collision rate for a given flow realization is calculated. This avoids the
assumption of isotropy, as needed in the classical approach, and allows for the evaluation
of the collision rate in general turbulent flows.
Direct numerical simulations and analytical studies on particle collisions in an
isotropic, Gaussian turbulence have been carried out. Several different collision models
(or post-collision treatments) are implemented in the simulations and the difference
among various collision models is examined. The effects of particle inertia, size, gravity,
and mean-concentration on the collision rate, collision velocity, and collision angle are
investigated. Asymptotic analyses have been developed to predict the effects of small
inertia on the increases in collision rate and collision velocity. Asymptotic prediction for
the increased collision velocity agrees well with the results of both numerical simulation
and PDF-based method. The increased collision velocity is found to contain terms of
both O((Ot k )-1) and O( (Tk )[2 ), where P1 is the particle response time and tk is the
kolmogorov time scale. The increases in particle size, gravitational acceleration, and
mean-concentration in general reduce the normalized particle collision rate in isotropic
turbulence for monosized particles as opposed to very small particles in the very dilute
limit.
A new method has been developed to quantify the concentration non-uniformity
of particles in turbulent flows. The effects of particle size and inertia on the concentration
non-uniformity are examined and elucidated. The concentration non-uniformity can be
greatly overestimated if the effect of finite particle size is neglected.
CHAPTER 1
INTRODUCTION
Particles and turbulence are around us everywhere and greatly impact our daily
life. Dusts, sand, rain drops, fine powders, glass beads, viruses and bacteria are all
'particles.' Most of the fluid flows occurring in nature such as air current and water flows
in rivers and seas are turbulent flows. The particles in the atmosphere remain airborne for
long periods, lowering visibility. The brown haze that is often seen over large cities in
autumn, winter and spring is due mainly to the interactions between particles and the air.
The improvement of particulate-fluid systems requires a better understanding of
the fundamental physics of the particle-particle and particle-fluid interactions. The
behavior of particulate-fluid mixtures may vary widely depending upon the conditions
under which the particle-fluid system functions. Particle-particle interactions involve the
collisions among particles which depend on inter-particle forces, particle size distribution,
and the concentration distribution. Particle-fluid interactions include the hydrodynamic
forces acting on particles by fluid and the feedback to fluid from particles, etc. Particle
collision, coagulation, and aggregation are the examples of most important scenarios of
those interactions that are determined by particle properties and flow structures. The fact
that collisions do occur among the particles of size ranging from a few Angstroms to a
few millimeters explains the vast amount of scientific and industrial applications of the
collision process. For example, molecular collisions help explain the physical properties
such as viscosity, diffusivity, and conductivity in gases and liquids. Numerous industrial
applications can be cited: design of fine spray combustion nozzles, control of industrial
emissions (Presser & Gupta 1993, Tu et al. 1996) and the understanding of the collision
process among suspended particles in pneumatic transport systems or in gas cleaning
chambers (Abrahamson 1975, William & Crane 1983), and in sewage disposal devices
(Hunt 1982), where particle size distribution and concentration non-uniformity are
important characteristics that depend on particle collision and coagulation. Most of
previous studies on particulate removal have focused on the filters, with little attention to
the particle-particle interactions and particle-fluid interactions. Despite many years of
research and development in particulate removal technologies, the collection of small
particles with size less than 1 un is still a rather difficult task. It is extremely difficult to
remove fine organic active particles (<0.03 mrn) from a waste stream because the particles
are much smaller than the pores in existing filters. Typical collection efficiencies of the
filtration devices decrease rapidly with the decrease of particle size in the range 0.2-30pm
(Davies 1953). A good alternative is to aggregate fine particles which can be removed
using conventional technologies. For the coagulation and aggregation to occur, particles
have to collide with each other and adhere. It is very desirable that one can predict how
the small particles are driven to collide with each other so they aggregate or grow in size
and become large enough to be captured by the filters. Obviously, the efficiency of
particle removal mainly depends on the control of particle collisions and concentration
non-uniformity in filtration systems (Fathikalajahi et al. 1996). It is evident that particle
collision is one of the most important interactions for understanding the fundamentals of
particle-fluid systems.
Particle collision kernel and collision rate are critical quantities in the
development of collision theories. The rate of flow induced particle collisions affects
directly the particle size distribution in particle production. Examples of turbulence
induced particle collision can be seen in the formation of rain drops in clouds, pulverized
coal combustion, agglomeration of fine powders in gas flows, air filtration equipment,
sewage disposal devices, fast fluidized beds, dust and spray burners, and so on. Particle
collisions can take place through a variety of collision mechanisms. Particle collisions
may be induced by centrifugal and gravitational forces, thermal forces in Brownian
motion, and van der Waal forces. In an actual system, the scenario after every contact
between particles may be different from one to another, depending on the underlying
physics of the system. The collisions among raindrops enhance the growth of droplets in
the air, since the contact between droplets may most likely result in the coalescence of
two droplets. The growth of bubbles in a liquid flow is another example of particle
coagulation upon collisions. The contacts between two solid particles may behave in
various different ways. When two solid particles are brought together the colliding pair
may stick or separate after collision upon contact. Many theoretical efforts were made to
predict the collision kernel and rate. To validate a theoretical analysis, direct numerical
simulations were carried out to obtain statistical average of the collision rate in a particle-
fluid system. Different interpretations of the immediate consequence after the contact
between a pair of particles have resulted in different collision models or so-called "post-
collision treatment" schemes. To quantify the post-collision treatment, a statistical
parameter called collision efficiency has been introduced. The collision efficiency is one
if particles coagulate to form a larger one. Usually, it is less than one since some of the
colliding particles may bounce back. Although many theoretical and experimental studies
have been devoted to estimates of the collision efficiency for water droplets settling
through still air (Klett & Davis 1973, Lin & Lee 1975, Beard & Ochs 1983) and droplets
settling down through laminar shear flow (Jonas & Goldsmith 1972), even an
approximate estimate of the collision efficiency in a complicated situation can be very
difficult to obtain. At this stage, it is still not possible to obtain an accurate solution for
the evaluation of collision efficiency in complex systems. To better focus on the issues of
our interests, the collision rate, the collision efficiency is assumed to be either one or
zero in this study. Although a unity collision efficiency will overestimate the size growth,
however, theoretical and numerical exploration on the governing physics may still hold
its significance for a better understanding of the interactions in particle-fluid systems.
Population balance equations are often used for predicting the evolution of
particle size distribution in particle production systems or in studying particle
coagulation. They are typically given in the form of
dn1
d t -anll1n1 a12nln2 a13nln3 -
dn2
dt 2 al jlnn a12nln2 CC22n2n2- 232n3 .- (1.1)
dn3
dt = a12nln2 a13nln3 X23n2n3 c33n3n3 -
dn4 1
dt a13n1n3 + 2 22n2n2 "a14n1ln4 24n2n4 -
In the above, ni is the number concentration of the ith particle size group. For primary
particles, i=l, with a volume v1, the particle in the ith group has a volume of iv1. The
quantity aijninj, denoted as ncij in this thesis,
ncij= aijninj (1.2)
is called collision rate which represents the number of collisions among ith particle with
jth particle per unit volume and per unit time. The coefficient aj is called the collision
kernel or collision function, it depends on particle size, particle inertia and flow structure,
and must be evaluated separately. The collision rate or collision kernel is a measure of
the ability of particles to collide, coagulate, and aggregate. It is essential to determine the
instantaneous particle size distribution and momentum transfer among particles in a
particulate system. Over past 80 years, many efforts have been devoted to developing the
particle collision theories for predicting the particle collision rate in fluid flows. A brief
review is given below.
1.1 Smoluchowski's Prediction in Laminar Shear Flow
A geometrical collision occurs when two particles reach a separation distance less
than R-ri+rj, in which ri and rj are the respective radii of two particle groups and R is the
radius of the colliding surface. Theoretical analysis of geometrical collisions in a uniform
shear flow was first carried out by Smoluchowski (1917). At a given instant the particle
number concentrations per unit volume are n, and nj for each size group. Taking any
particle of the ith-group as a target particle, the volume flow rate of thejth-group particles
reaching the colliding surface of radius R of the given target particle is nj f ,, wrdS in
which w is the relative velocity between ith-group and jth-group particles and Wr is the
radial component of w. Thus the number of collisions occurring among ith-group and
jth-group particles per unit time per unit volume is given by
ncij = ninj L, owdS= ninj F(ri+rj)3 = a, nin (1.3)
where F is the shear rate and
aij= F(ri+rj)3 (1.4)
is often called the collision kernel. For like particles, the number of collisions per unit
time and per unit volume was often simply taken as
1 4 3
ncii= nini (2r,)3. (1.5)
2 3
The factor was introduced to avoid double counting of the collisions among the like
2
particles. When the above results are used in the population balance equations, for
example, the rate of the "destruction" of ith-group particles due to collisions is -2 n ci +
4 4 4
Y n = nini F(2ri)3 ninj r(ri+rj)3 since each collision among the like particles
results in the loss of two particles of that size group. Equation (1.5) has been extensively
used in the literature in calculating the particle collision rate and particle size evolution.
In addition to Smoluchowski's work, Hocking & Jonas (1970) also calculated the
collision rate for the sedimenting particles in a uniform shear flow using Stokesian
theory. They studied the collision efficiency by which the particle coagulation was
assumed if the surface distance between two neighboring particles was less than a
prescribed fraction of particle diameter.
1.2 Various Theories for Collision Rate in Turbulence
When particles collide in complex flows, such as turbulent flows, or when more
than one collision mechanism are involved in the particulate system, it is quite difficult to
derive an exact expression for collision rate. Numerous efforts have been made to
estimate the collision rate by means of evaluating the collision velocity in turbulence
during past decades. Camp and Stein (1943) employed the concept of mean velocity
gradient related to the turbulent energy dissipation rate E and kinematics viscosity v and
obtained the collision kernel for turbulence shear-induced collision rate based on
Smoluchowski work (1917) in laminar shear flow.
4 R(__1/2
a = R' 2. (1.6)
3 v
Saffman and Turner (1956) gave two predictions of particle collision kernel in an
isotropic, Gaussian turbulence. In the first prediction, it was assumed that particles are
inertialess so that they follow the fluid completely and particle size is small compared
with the Kolmogorov length scale of the turbulence. The collision kernel was given
ca = 1294 R3 (E)1/2. (1.7)
IJ V
Their second prediction was for the particles with small inertia and gravity.
ai = 2(2n)12 R2 [(1 -)2 -t)2 (-)2
PO Dt
+(I )2i (T-2)2g2 +R2 ]1/2. (1.8)
3 po 9 v
However, collision kernel for like particles or the particles with same inertia is reduced to
a, =1.67 R3 ( ) 1/2 (1.9)
v
Thus the effect of particle inertia and gravity disappeared in the monodisperse case,
which is not physically sound. In addition, the result of their second analysis was
inconsistent with the first one, as easily seen from (1.7) and (1.9).
Abrahamson (1975) extended the concept of kinetic theory of gases for the
collisions among large particles in a high Reynolds number turbulence. It was assumed
that the velocities of colliding particles are uncorrelated, normally distributed in the
absence of gravitational settling. The variance of particle velocity v was estimated by
2 "2
= (1.10)
1+ 15 Tp 8/U2
where u2 is the mean-squared fluctuating velocity of the fluid and Tp is the particle
response time. The collision kernel, without gravitational settling, was given by
2 ( 2 > 2 > 1/2.
aij = 5.0 R ( + )12 (1.11)
Balachandar (1988) performed an analysis on the collision kernel by assuming
that the probability distribution function (PDF) of the velocity gradient in turbulence
follows a log-normal PDF distribution. He found a correction to Saffminan & Turner's
prediction for collision kernel in isotropic turbulence. Williams & Crane (1983) and
Kruis & Kusters (1996) performed stochastic analysis to obtain the collision kernel with
inclusion of the combined effects of local shear and particle inertia. Zhou et al. (1997)
carried out a similar analysis and modified the Kruis & Kusters' result in turbulence.
Those analyses on the collision kernels were obtained on the basis of two-point velocity
correlations without tracking the trajectories of individual particles in the turbulent flow
(hereinafter referred as Eulerian method). However, we have found that the collection of
the relative velocities (or collision velocities) upon collisions based on tracking the
trajectories of individual particles in turbulent flow is a biased average in favor of those
regions where the collision velocity is high. Hereinafter, this method based on particle
tracking is referred as Lagrangian method. The detailed discussions on the difference
between the results of Eulerian and Lagrangian methods on collision velocity will be
presented in chapter IV.
1.3 Recent Development in Computer Simulation on Particle Collision Rate
Recent developments of direct numerical simulations (DNS) in turbulence have
greatly enhanced the understanding of the turbulence-induced particle collisions.
Balachandar (1988) also carried out DNS of particle collision in an isotropic turbulence.
It was found that the probability distribution function of the velocity gradient follows a
log-normal distribution instead of Gaussian distribution. The collision kernel was
obtained in a frozen turbulence with a maximum of 500 particles in the flow field.
Artificially large particles were used in order to increase the collision events within a
reasonable computational time. It was concluded that the particle inertia does not affect
the average collision kernel for large, monosize particles. Large discrepancy in the
collision rate for monosize inertialess particles between his numerical results and
Saffinan & Turner's prediction was reported. However it is not clear whether this
discrepancy is due to the difference in the small scale turbulence structure or not.
Chen et al. (1995) carried out DNS results for the aerosols in a turbulent channel
flow. They mainly focused on the possible effects of collisions on the deposition rate of
the droplets on the channel walls. They found that the particle concentration in the wall
region can be significant while the bulk concentration is small. The collision kernel and
rate were found to peak in the wall region, and particle depositions can be enhanced by
the formation of larger particles.
Wang et al. (1997) also carried out DNS to investigate the collision kernel in
homogeneous turbulence. They allowed particles in the system to overlap in space which
is unphysical. The effects of finite inertia on collision kernel were also reported in their
DNS results (Zhou, Wexler & Wang 1997). Zhou et al. (1997) developed an analysis
showing that the increases in particle collision velocity, the non-uniformity of particle
concentration and the compressibility of particle velocity in the continuum sense are three
possible mechanisms to affect the overall collision kernel for particles with finite inertia.
However, they did not quantify these terms separately in their work. Thus how the
collision velocity and the concentration non-uniformity play in enhancing the overall
collision kernel remains unknown at this stage.
1.4 Recent Development in Particle Concentration Non-uniformity
Recent studies have shown that the particle concentration in turbulence may be
highly nonuniform. Heavy particles with finite inertia tend to accumulate in the regions
of high strain rate and to move away from the vortex cores. The occurrence of particle
accumulations is referred as concentration non-uniformity or preferential concentration.
The coherent vortical structures are the mechanisms that cause the nonuniform
distribution of particles. The coherent vortices and concentration non-uniformity have
been observed by many numerical simulations and experimental visualizations. The
numerical simulations for demonstration of the concentration non-uniformity were carried
out in plane mixing layers (Crowe et al. 1988, Chien & Chung 1987, Wen et al. 1992,
Tang et al. 1992, Wang 1992), in wake flows (Chien & Chung 1988, Tang et al. 1992),
in jet flows (Chung & Trout 1988, Hansell 1992), in wall-bounded flows (Kallio & Reeks
1989, McLaughlin 1989, Yonemura et al. 1993, Rouson & Eaton 1994, Tanaka et al.
1996), in homogeneous flows (Maxey & Corrsin 1986, Hunt et al. 1987, Fung & Perkins
1989, Squires & Eaton 1991, Wang & Maxey 1993, Fessler et al. 1993). The non-uniform
distribution of particles has also been found by the experiments in plane mixing layers
(Kobayashi et al. 1988, Kamalu et al. 1988, Lazaro & Lasheras 1992, Ishima et al. 1993),
in wake flows (Tang et al. 1992, Yang et al. 1993, Bachalo 1993), in jet flows (Longmire
& Eaton 1992), in wall-bounded flows (Rashidi et al. 1990, Young & Hanratty 1991).
Wang & Maxey (1993) have, through direct numerical simulations (DNS) in isotropic
turbulence, shown that the preferential concentration achieves maximum at P3k =1 in
which P3 is the reciprocal particle response time and Tk is the Kolmogorov time scale of
the turbulence. One of the consequences of the preferential concentration is the
enhancement of the settling velocity of heavy particles in the Stokes drag range (Maxey
1987, Mei 1994). One of the related industrial applications is that the particulate
removal efficiency of the venturi type scrubber was affected by the nonuniform
concentration distribution of droplets (Fathikalajahi et al. 1996).
Particle trajectories, therefore individual particle positions, can be obtained
accurately only in Lagrangian frame by solving the equations of particle motions.
Particles have finite response times in general but they do not have an intrinsic diffusivity
or viscosity associated with the motion of individual particles. Turbulence particle
diffusivity is the consequence of ensemble averaging over many realizations and it is used
for predicting the variation of average concentration field on a large scale. It cannot be
used to describe the small-scale spatial variation of the particle concentration. From the
computational point of view, a DNS using 1283 grid resolution for an isotropic
turbulence becomes routine and DNS using 2563 or even 5123 have been carried out.
For a reasonable contour representation of particle concentration, it is desirable to have at
least, say, 10 particles in each cell so that the statistical noise does not become
overwhelming. This requires 1.7x108 particles with a 2563 grid resolution. Due to the
preferential concentration, there will be regions and cells that do not contain any particle
or only one or two particles at a given instant. Thus, an accurate, three-dimensional
contour representation of instantaneous particle concentration with such high resolution
will still be heavily contaminated by the inherent statistical noise. No such three-
dimensional concentration contour has been reported to date. Due to the computer power
limitation, previous numerical simulations involved anywhere from 104 to 106 particles,
which is far below the level of reaching a reasonable statistics.
It has been difficult to quantify the extent/degree of particle concentration non-
uniformity induced by the small-scale turbulence structure. Squires & Eaton (1991) and
Wang & Maxey (1993) showed in their DNS the similarity between the instantaneous
particle concentration contours with the contours of the vorticity magnitude (or
enstrophy) in the same plane at the same instant. To quantify the deviation of the particle
concentration field from a statistically uniform distribution, Wang & Maxey (1993) used
a statistical quantity (P(k) Pb(k))2 averaged over space, where P(k) is the
probability of finding k particles in a given cell and Pb (k) follows the Bernoulli
distribution for a random distribution of the particles. While this quantity shows that the
preferential concentration maximizes near P3tk =1 (Wang & Maxey 1993), it varies with
total number of particles in each of their simulations. The particle-particle interactions
are neglected in their simulation and no systematic investigation on its dependence on the
cell size was carried out. Hence, its usefulness is limited. Although numerous
experimental studies have clearly shown the non-uniform particle distribution in turbulent
flows, it is quite difficult to measure the 3D non-uniformity so as to isolate the effects of
flow shear, particle inertia and gravity.
Obviously, the local accumulation will affect the encounters among the particles
and thus the overall collision rate is expected to be modified accordingly. However, the
effect of the preferential concentration on the collision kernel has never been estimated or
separated from other factors such as the increased particle collision velocity due to small
particle inertia.
1.5 Scope of Present Study
Particle collisions are extremely difficult to observe experimentally, even with
today's sophisticated laser instrumentations. Direct numerical simulation, however, has
its advantages of providing detailed information on particle trajectories, flow structures
and is ideally suited for parametrical investigations of particle and fluid flow properties.
Numerical implementation of particle interactions in a simulation, while posing a
considerable numerical challenge, introduces no fundamental difficulties. In this work,
extensive data obtained by numerical simulations on the collision velocity, collision
angle, collision kernel and collision rate are vital to understanding of the fundamental
physics in the particle-fluid interactions. One of the major objectives of this thesis is to
develop analytical and numerical methods for predicting the particle collision rate and
concentration non-uniformity with respect to the effects of turbulent shear, particle inertia
and gravity in particle-fluid systems.
As mentioned in 1.2, previous theories on the collision rate were restricted on
assumed turbulent flow. In Chapter III, we have developed a new theoretical framework
for prediction of the collision rate and collision velocity of small particles at zero inertia
in general turbulent flows. The collision velocity is a critical variable related to the
collision rate. However, the effects of particle inertia on the collision velocity and
collision angle have never been predicted using Lagrangian tracking method. In Chapter
IV, we carried out an asymptotic analysis to predict the effects of small particle inertia on
the collision velocity and presented the effect of the inertia on collision angle in isotropic
Gaussian turbulence. As pointed out in 1.4, no reliable method exists for quantifying
the concentration non-uniformity and the effect of the preferential concentration on the
collision rate has never been estimated. In Chapter V, a new quantification method has
been developed to reliably quantify the concentration non-uniformity of particles in
turbulent flows, and the effect of the preferential concentration on the collision rate is
predicted in isotropic Gaussian turbulence. The first result of Chapter VI is our
correction to the classical result of Smoluchowski for the collision rate of monodisperse
particles in a laminar shear flow. The combined effects of both the increased collision
velocity and increased concentration non-uniformity on the collision rate due to the
inertia are then examined. The effects of particle size, gravity and mean-concentration on
the collision rate are investigated through direct numerical simulations in Chapter VI.
The effects of the collision models on the collision rate are also presented in Chapter VI.
CHAPTER 2
TURBULENCE AND PARTICLE COLLISION DETECTION
2.1 Turbulence Representation
2.1.1 Isotropic. Gaussian Turbulence
The following model for the energy spectrum developed in Mei & Adrian (1995)
is assumed for turbulence,
3 2 k' 1 22
E(k) =2 uo Vk [l+(k/ko)217/6 exp(-To k2) (2.1)
where uo is the root-mean-squared turbulent velocity, k0 is a typical wave number and the
dimensionless parameter ro ---oko is related to the turbulent Reynolds number Re,.
Large Re. correspond to small rio and vice versa. The normalizing coefficient W in E(k)
is determined from (2.1) by satisfying the total energy requirement,
0 -4
k 2 (-1k 2'd
7 /6 2 exp2(-71 dk (2.2)
0 (1+k )
For small r10, E(k)~k4 when k = k/ko l and E(k)~k53, which is the scaling law in the
inertial subrange, when 1 k /1"qo. For large 10o, the above energy spectrum recovers
that of Kraichnan (1970). The relationship between the Eulerian integral length L, ko
and ri (or Re.) is given in Mei & Adrian (1994).
A turbulent eddy loses its identity as it is convected and dissipated; this behavior
is described by the eddy-self decay function D(x) which has an integral time scale TO.
The Fourier transformation of D(T) gives the power spectrum D(o)). A composite form
for the power spectrum D(o)) was constructed in Mei & Adrian (1995) as
2
exp(-ri Io)2)
() =A +(To(2.3)
where TI, is intended for the viscous dissipation time scale if the turbulent Reynolds
number is high, and T is close to To if ii I/To 1. The coefficient A is determined as
T exp(-rl /T2)
(2.4)
nt l-erf(rTl/T)
by satisfying J D(o) do) =1. The integral time scale To is related to T as To = 7tD(0) =
-oo
7lA .
The relationship between L I and To is not known in general; it can be taken to
have the form
To = cE(Re)LI 1/UO, (2.5)
where cE(Rex) was determined approximately in Mei & Adrian (1995) based on the
experimental data of Sato & Yamamotto (1988) for the fluid dispersion and the value of
high Reynolds number turbulent Prandtl number.
Using random Fourier modes representation, an isotropic, Gaussian, pseudo
turbulence
N
Ui (X, t) E bm) cos k(m) +(,(m) t
u ( ,t = [bi (k -x + .t )
m=!
+cIm) sin(k(m) x + co( m t ) ] (2.6)
is constructed to simulate the turbulent flow with a specified energy spectrum. In the
above Nk is the number of the random Fourier modes in one fluid realization,
k(m)and o)(m) are the wavenumber and frequency of the m-th mode. The random
coefficients b(m) and c(m) are chosen as follows,
b(m) (m)(6 (m) k(m)/k (M) f(k, (nM) (2.7)
i j qJ i j
=b ((Si)-ki W k ) fk^m ca (2.7)
where b. follows a normal distribution with
= 0 & = 8". (2.8)
The factor (ij-k(M) k(m) /k (M) in Eq (2.26) ensures the incompressibility V-u=-0 for
every Fourier mode. The scale factor f(k, o) depends on the energy spectrum and the
probability density functions (pdf), pI(k) and p2(00), of k and co,
f2(k, ) = E(k) D(o) (2.9)
2~,o = 5 (2.9)
47N2kpl(k)p2((O)
Since E(k) decays as k-5/3 in the inertia subrange, a Gaussian distribution for k would lead
to a very slow statistical convergence for ui(x, t). To achieve convergence of the statistics
involving the derivatives such as &', the following algebraic pdf
axj
p- Ii(ki)= (1 + Iki)" fori=l, 2 & 3 (2.10Oa)
2
and P1(k) = P11(k1)P12(k2)Pi3(k3) (2.10Ob)
is used to sample k, for ko= 1. The frequency om) is generated with the following pdf,
1 1
P2() = t +(/uko)2 (2.11)
The scale factor is finally set to be
E(k)D(co)
f(k, (.) = E-- (2.12)
47rtNkk P l(kl)Pi2(k2)p13(k3)p2((O)
It is noted that in the present studies on particle collisions, the collision rate is dictated by
the spatial structure and the eddy-self decay has no effect on the collision rate. In the
direct numerical simulation of particle collision in isotropic turbulence, however, the
eddy self-decay can be easily incorporated. Hence, the dependence of b(m) and c(m) on (o
is still included.
2.1.2 Rapidly Sheared Homogeneous Turbulence
For an initially isotropic turbulence under rapid shear, F=a, the variation of
b(m) with time was given in Townsend (1976) based on rapid distortion theory (RDT).
Hunt & Carruthers (1990) gave an overview on the use of RDT for representing various
types of homogeneous turbulence. Denoting the total shear S as
S= Ft (2.13)
the new wavenumber vector x(m) after the shearing is given by
(M) = k(m)
I I
() = k(m) -Sk(m) (2.14)
x2 2 1
(M k(m)
x k 3 3
where k(m) is the wavenumber vector of the m-th mode before the mean shear is applied.
The amplitude of each Fourier mode becomes
b(f)(S) = b^m)(0) + alb(2)(0)
b()(S) = a2 b(2m)(O) (2.15)
b(m)(S) = b(3)(0) + a3b(m(0)
where b(m)(0) is the corresponding amplitude of the m-th mode before the rapid
distortion. Dropping the superscript m for convenience, the coefficient ai for each mode
m is given as
Sk k2-2k2+Sklk2
al= 2 2
a-k22 x
k2k2 k2 k2-Ski
2 23/2 [tan-it ()tan-] ( )]
2k23/ 2, 22
(k +k )3 k, l2i- A;- f
1 33k1k
a2 = k2/x2 (2.16)
Sklk3k2-2k2+Skjk2
a3 = k2 2 X2
1 ~3
k2k3 1 k2 k2-Sk.
,2 ,.23/2 [tan- (- ) tan(- )]
2 2 32 2k22
where k= VkT+k2+k and x= 2 +X 2+ The coefficient c(m) is similarly obtained. It
1 2 3 2 3*
is clear from the above that the amplitude of each mode is entirely determined by the total
shear S=Ft. The dissipation rate s(S) can be evaluated as
1 N 2
e(S)/v= I < N {(Xm)2 3 [b()2 (S) +C(M)2(S)] } >" (2.17)
m=l i=l
Finally, the total velocity field is
uj(x, t) = [ b(m)(S) cos((m) -x + )(m) t)
m=l
+ c( (S) sin((m) x +)(o Mt)]
+ Fx2 8ii. (2.18)
Since the shear rate is high, the term o(m)t in the above acts merely as a random phase.
2.1.3 Isotropic Turbulence Generation in a Periodic Box
For the direct numerical simulations of particle collisions in this work, the
velocity field of an isotropic and Gaussian turbulence is generated by a random Fourier
modes representation (2.6). To effectively conduct a computational study on particle
collision, a region of finite size in which a large number of collisions occur should be
used with periodic boundary conditions. To render the turbulence represented by random
Fourier modes periodic in a box of volume R3, the random wave numbers, k( m), are
rounded to the nearest even integers while all other quantities are held fixed so that the
flow is periodic within the box of volume irt3. The periodic box of volume n3 yields
virtually the same collision statistics when the volume is increased to (27t)3. Hence it3 is
used throughout this study. The corresponding statistics (turbulence rmns velocity u0,
Reynolds number Rex, Taylor micro length scale X, Kolmogorove length scale T, and
Kolmogorov time scale Tk) are computed based on the modified wavenumbers. As will
be demonstrated using the present analysis outlined in 2.1.1, the collision statistics,
when normalized using R3( /v)1/2 for anI and <-Wr>, the results are the same using the
turbulence generated in an unbounded domain and in a periodic box.
Unless specifically mentioned, most of the isotropic turbulence used in the
numerical simulations throughout this study is within the periodic box of length nt and
have the following characteristics: root-mean-squared turbulence velocity u0=l.O,
integral length scale L, 11 =--0.594, Taylor micro length scale X=0.277, Kolmogorov length
scale r=O.023, Kolmogorov time scale Tk = (v / )I/2 =0.073, and Rex =40.1. The
turbulence is imposed on a box it3 with periodic boundary conditions in three
perpendicular directions. The periodic treatment of boundary conditions allows the
simulation results by a finite volume to be extended to an infinite homogeneous system.
Initially, a given number of monosize particles are randomly introduced in the
flow field of volume it3. The periodic boundary conditions are also applied to the
motions of the particles. The particles are advanced using ui(x, t) given by (2.6) with a
time step At. A semi-implicit finite difference scheme is implemented to integrate the
equation of particle motion and to detect collision. The scheme has an accuracy of second
order. The detailed error analysis of the semi-implicit finite difference scheme will be
given in next section.
2.2 Numerical Integration of Particle Motion and Error Analysis
Major concerns in numerical simulation are the accuracy and robustness of the
numerical scheme. Since the collision detection in the present simulation is based on the
calculation of particle trajectories, the collision kernel/rate/velocity/angle are affected by
the accuracy of the finite difference scheme. In this study, both tracking particles and
detecting collision are carried out in a Lagrangian frame. A semi-implicit finite
difference scheme is implemented to integrate the equation of particle motion and to
detect collision. This scheme is unconditionally stable. With consistent treatment for
both particle velocity and turbulence, The truncation error in the semi-implicit scheme is
0((At)2) and independent of 3. The accuracy analysis of the scheme is given as follows.
The particle dynamic equation in the absence of (or weak) gravity is
dv
dt = P(u v). (2.19)
The semi-implicit finite difference scheme for the particle dynamic equation (2.19) is
Vn+i -- n 1
--- n I (U.V)n+l +(U-V)n]
At 2
= p[(un+l + u-(vn++v)]. (2.20)
2
Taylor series expansions for u"n+1 and vn+1 give
~~'=~ +AtdVn +1 2d2vn dv
vn+ vn + dAt d + (At)2d + (At)3 3vn+ (2.21)
dt 2 2 dt2 6 dt
un+ =u AdU 1 Ad2u( 1 3d3un
dt 2 dt62 d-"
Each side of equation (2.20) can be evaluated as
dvn I d2vn 1 2d3Vn
LHS= ++ -(At) + -(At + O(At)3 (2.23)
dt 2 ( t dt2 6 dt3
I Atdun dvn 1 2 d~u dvn3
RHS=p(un-vn)+-PAt( dt )+-IP(At)( d 2u- 2 )+O(At)3. (2.24)
2 dt dt 4 dt dt
Then equation (2.20) can be rewritten as
Vn+l --V 1
0 = n-- [(u-v) n+' + (U-v) ]
At 2
dv" 1 d2vn 2 d3vn
+-(At) -+-(At)+d
dt 2 dt2 6 dt3 "
- pi3(u v"
I PAt( du- dv" ) P(At)2 (d2un
2 dt dt 4 dt
Thus, equation (2.20) is equivalent to
dvn 1n d2vn
dt 2)- t dt2
- (At)2 d3vn
6 dt3
I dun dvn I d2un d2vn
+- PAt( )+ (At)2( dt2)
2 dt dt 4 dt dt
I dAt[2vn
2 dt2
dun dvnA
dt -dt
6 d3v
6At)2 dt3
1 2 d2un d2vn
+ -I(At) (-2 )
4 dt2 dt2
Taking derivative with t repeatedly, we get
d2vn dut
dOt dt
dvn
dt
1 d3vn
- At[ d n
2 dtT
d2u
-d(
d2vn)]
dt2 A
- (At)2 dt-
6 dt4
+ 4 dt3
d3vn
dt3
d2un
d2v)
dt2
I d 4v
-At[ --
2 dt4
d3un d3
1 _dO -_ __d _
1 2d n + I p 2d4vn d4vn
-6 dt- 4 dt4 dt- )t
Repeatedly substituting above equations into the RHS of (2.26) results in
dv n (A)2 d3vn 1(2 d3vn
_+ (At)2_, + O((At)3)
dt 6 d 4 dt 3
d2vT
-dt2
(2.25)
(2.26)
d 3vn
dt3)
dtO
(2.27)
(2.28)
1(At)2 d3vn (2.29)
= M)2d--+ O((At)3). (2.29)
12 dt3
It is noted that the truncation error in equation (2.29) is independent of P3. Thus
motions of particles with very large P can be simulated without loss of accuracy.
2.3 Particle Collision Treatment
2.3.1 Collision Detection Scheme
Collision detection between any pairing particles is given as follows (Chen et al.
1995). At nth time step to, these two particles are initially located at
r1o =(x0, y1O, z10) and r20=(x20, Y20, z20) After a small time advancement At,
at (n+J)th time step ti, new positions of these two particles are r11 =(x11, y11, z11)
and r21 = (x21, Y21, z21 ). The particle position vectors at time t, which is between to
and t I, can be expressed as
r1(t) = {x1o0 + v1(t-to), yio0 + Vy(t to), Zio0 + vz(t to)} (2.30)
r2(t)= {X20 +v,2(t-t0), Y20 +Vy2(t-t0), z20 + Vz2(t t0)} (2.31)
The distance between two particles, s(t), is
s(t)= I r2 (t) ri (t)|. (2.32)
These two particles will collide if the condition s(t)
relative distance s(t) varies with time t, and it is not a monotonic function of t. The
minimum relative distance occurring at time tm is determined by
9s,
ast=t =0. (2.33)
a~t ttm
Combining (2.30)-(2.33) gives
D(
tm =to -- (2.34)
D2
where
DI = (xiO X2o)(Vxl Vx2) + (Yo y2o)(Vyi Vy2)
+(ZIo -Z2o)(vi Vz2) (2.35)
and
D2 =(Vxi -v,2 )2 +(Vyi Vy2)2 + (vz -v)2.) (2.36)
Equation (2.33) does not guarantee that the moment, when the above mentioned
minimum distance occurs, has to be within the time interval to and to + At. Three cases
need to be considered for calculating the minimum distance within a time step.
(i) if tm < to, then the minimum distance is at the initial time to.
Smin = s(to) (2.37)
i.e. Smin = {(X1io -X20o)2 +(yio -Y2o)2 +(Z10 Z20o)2}1/2. (2.38)
(ii) if tm > to + At, the minimum distance must occur at the time to + At.
Smin = s(to + At) (2.39)
i.e. Smin = (s2 +Sy +s2)12 (2.40)
where
Sx= xl0 -x20 +At(vxl Vx2)
Sy = 0 Y20 +At(vyl -Vy2)
sz= zl0 -z20 +At(vzl Vz2).
(iii) if to < tm to + At, the minimum distance occurs during the time interval
to and to + At. In this case, we have
Smin= S(tm) (2.41)
2 2 2 .1/2
i.e. Smin = (sx + Sy + S )1/2 (2.42)
where
Smx = XlO X20 + (tm to)(Vxl Vx2)
Smy = YI0 Y20 + (tm to)(Vyl Vy2)
Smz ZIO Z20 + (tm to)(Vzl Vz2 ) .
23
After Smin is evaluated within the time step, the examination of collision events is
still needed because Smin can be either larger or smaller than R. Collision occurs if
Smin < R (2.43)
At the collision instant tc,
s(t) = R (2.44)
and s(tc)=I r2(t0)- r(t.)1 (2.45)
Hence the collision instant tc can be obtained solving (2.44) and (2.45)
-b- ]b2 -4(s R2)w2
tc =2to+ (2.46)
2w2
where
b=2[(xio -X20)(vx, -vx2)+(YIO -y20)(VyI -Vy2)
(ZIo z20)(vzl vz2)
2 22
SO {(x10-x20) +(yIo-y20)2 +(ZI0-Z20)2}12
W = {(Vxi Vx2)2 + (Vy! Vy2)2 + (Vzl Vz2)2}/2 .
With perfect sticking collision model, any collision can produce a new born particle
(daughter particle). The daughter particle will inherit all possessions of parent particles,
i.e., conservation of mass and momentum. Let the daughter particle be numbered 3, then
the conservation law gives the size of daughter particle as
d3 =(d + d32)/3 (2.47)
where di (i=l 1, 2, 3) are particle radii. The velocity components of daughter particle are
Vx3 d + 2d2 (2.48)
d3
v = d3+ Vd3
Vy3- y 3y2 2 (2.49)
d3
vzld3 + V z2d 3
vz3 1 2d (2.50)
3and the position components at the end of time step
and the position components at the end of time step
3 3
X3 = {[xl0 + vxl(tc t0)]d1 + [x20 + Vx2(tc to)]d3} / d3 + v3(to t + At) (2.51)
Y3 = {[Y0io +Vy(tc -t0)]d + [Y20 + Vy2(tc to)]d} / d + Vy3(to tc + At) (2.52)
Z3 = {[z 0 + Vzi(tc to)]d3 + [z20 + vz2(tc t0)]d} / d3 + vz3(to tc + At) (2.53)
In general, the collision detection takes a lot of computer time since a large
number of particles are considered. The primitive collision detection scheme searches the
collision candidates by going through all particles N, so that the computing time is
proportional to N2. The collision detection method in Chen et al. (1995) is an efficient
scheme, and it is also implemented in this thesis. By dividing the computational domain
into a number of small cells, the potential collision partners for a given particle in one cell
are then searched within this cell and neighboring twenty-six cells during one time step.
This search only involves a small number of particles. The total computing time with this
scheme can be reduced dramatically. Only binary collision is considered by assuming a
negligible probability of multiple collisions in a small time step in a dilute condition. For
any particle moving out of the computational box, it is re-introduced into the box by
invoking periodicity. The collision rate or kernel can be determined, in principle, by
counting the number of collisions per time step. Initially, a given number of particles are
randomly distributed into the computational domain. The particles are evolved with
turbulence for a period of time before collision detection turned on. The collision events
are not counted until all particles with finite inertia have reached dynamical equilibrium
in turbulence.
2.3.2 Post-Collision Treatment
The collision event is counted when two particles are brought into contact, as
stated above. The collision among particles and its subsequent treatment are expected to
disturb the original system and permanently impact the subconsequent collision
evaluation. For the droplets in gas turbulence, the collision may, most likely, result in
coalesce of two droplets and a larger droplet is then formed by the colliding droplets.
Similar phenomena may be found in a system consisting of bubbles and fluids. The
collision between two solid particles may behave in different ways. Two solid particles
may stick together upon contact or the colliding solid particles will separate after
collision. In summary, the issue dealing with immediate consequence of those contacts is
referred as post-collision treatment. Several models with regard to post-collision
treatment have been developed during past decades.
In Chen et al. (1995) every collision results in the collision partners disappearing
from their respective size groups to produce a larger particle. Mass and momentum
conservation laws were applied for the birth of the new particle. This scheme is
considered to be a physically realistic model assuming a unity collision efficiency. It is
hereinafter referred as 'KEEP ALL' model. Balachandar (1988) used a different post-
collision scheme. Starting with monosize particles in system, after each collision, the
resulting larger particle was discarded in the computation so that it does not contribute to
future collisions. This post-collision treatment is hereinafter referred as the 'THROW
AWAY' scheme. Sundaram & Collins (1997) employed a 'HARD SPHERE' collision
model in which the particles were assumed to be rigid spheres and forced to bounce back
after the contact. Momentum conservation is applied to the particles after collision. In
contrast to the 'KEEP ALL' and 'THROW AWAY' schemes, the 'HARD SPHERE'
model allows a particle to have multiple collisions with other particles in a system. Wang
et al. (1997) tested three different schemes and presented one of their own. In their
scheme, collision does not disturb the system. The colliding particles are allowed to
overlap in space (hereinafter referred as 'OVERLAPPING' scheme), but no larger
particles are formed after the collision. Two colliding particles can separate after their
trajectories satisfy Ir2 (t)- r (t)l>R. Thus, a particle can have multiple collisions with
other particles. These four different post-collision treatments will, in principle, result in
different collision rate/kernel.
a(11
1.6
-1..... Dt=0.005, simulation
Dt=0.01, simulation
1.2
1.0. . - -- -
0.8
0 .6 -------- '
0.0 0.5 1.0 1.5 2.0 2.5 3.0
T
Figure 2.1 The effect of time step on the cumulative average of the normalized
collision kernel ca I for inertialess particles in turbulence. Particle size D=0.05.
2.4 Computer Simulation
2.4.1 Time Step
Three collision models: "THROW AWAY," "KEEP ALL" and "HARD
SPHERE" are employed in this work. For the following simulations with '"THROW
AWAY" and "KEEP ALL" collision models, the time step At=0.01O is used to advance the
motions of particles. Figure 2.1 shows the effect of time step on the collision kernel a *,
using "KEEP ALL" collision model. It can be seen that there is almost no difference
between the results of At=0.01 and At=0.005. Similar results are obtained by '"THROW
AWAY" scheme. The volume concentration of suspended particles is usually quite low
in many natural systems such as spray burners and air scrubbers, which is within the
interest of this work. In such dilute systems, the multiple collisions during a small time
step may be negligible, thus binary collision mechanism is considered solely in this
study.
2.4.2 Turbulence Realization
The averages of statistical quantities are evaluated in two directions: average over
all time steps and turbulence ensemble average over all turbulence realizations. The total
time steps mainly depend on the particle response time and the number of collision pairs
during whole time period. We start counting the time steps for collision detection after
the motions of particles have been evolved for three times the particle response time.
Ideally, to form the turbulence ensemble average, the numerical simulation should be
repeated an infinite number of realizations with different initial particle distributions. An
infinite number of turbulence realizations is impractical due to limited computer capacity
and may be unnecessary from the viewpoint of the accuracy of numerical scheme.
Figures 2.2 (a) & (b) show the turbulence dissipation rate against turbulence realizations.
While the ensemble averaged turbulence dissipation rate varies with the realization by a
large amplitude shown in figure 2.2 (a), the cumulative average of turbulence dissipation
rate approaches a constant for the number of turbulence realizations greater than 30, as
seen from figure 2.2 (b). Figures 2.3 (a) & (b) show the small eddy shear rate against
turbulence realizations. It can also be seen that, although the ensemble averaged small
eddy shear rate randomly varies with the turbulence realization in figure 2.3 (a), the
cumulative average of small eddy shear rate is close to a constant by the 40th realization
seen in figure 2.3 (b). Therefore, in this study forty turbulence realizations are used for
turbulence ensemble average.
2.4.3 Validation of the Collision Detection Scheme
Starting with monosize particles in system and using the "THROW AWAY"
collision model for post-collision treatment. The evolution of the particle number n1 (t)
follows the simple population balance equation:
dn1/dt=-a, In2 with n1(t = 0)=n0. (2.54)
The solution is
no / n1(t) = l +a n0t. (2.55)
It should be mentioned that the collision rate of monosize particles (2.54) is only accurate
for large number of particles in a system. For small number of particles, however, the
collision rate (2.54) must be corrected. The detailed discussions on this issue will be
given in 6.2. Since the statistics for the number of the remaining particles n1(t) is much
more accurate, the average slope of the curve 1 / n1(t) gives the collision kernel, so it can
be used to check the accuracy of the collision kernel obtained by simulation in laminar
shear flow.
To validate the collision detection, particle collision in a uniform shear flow of
velocity gradient [= 1.0 is first considered. For large number of particles the theoretical
prediction of the collision kernel aoij of Smoluchowski is accurate and is given by
ac = 4 R ir- / 3 For a small number of monosized particles the theoretical prediction
of the collision rate of Smoluchowski has to be corrected, as mentioned above. Ten
thousand particles of radius 0.075 are introduced into a box of volume 2x2x2. The time
step is Dt=0.01 and a total of 800 time steps is used. The number of shear flow
realizations is 80. Since no periodic condition may be imposed along vertical direction in
uniform shear flow, periodic boundary conditions are employed in the streamwise and
spanwise directions only. In such a case, a boundary correction of the collision rate has
to be imposed (Wang et al. 1997). Another correction on the collision rate is also applied
due to the inaccurate expression of the collision rate (6.2) Figure 2.4 shows the
particle number evolution. The collision kernel a 11 can be obtained from the average
slope of the curve. Very good agreement can be observed between simulation result and
Smoluchowski's prediction. Hence, the collision detection scheme and the numerical
implementation are validated.
9
8
(a)
7
c 6
o
L 5
S
S4
0
U
C
1e 3 .
**
p 2
1 *. *.
1 .* n I .
** *
0 5 10 15 20 25 30 35 40
Turbulence realizations
C
2.0
(b)
1.6
i ** **.** .** .
I 2 ** *12
U
C
S
0.8
I-
0.4 --- ---- --- ---- --- ---- --- ---
0 5 10 15 20 25 30 35 40
Turbulence realizations
Figure 2.2 Turbulence dissipation rate e against turbulence realization. (a) Variation of
ensemble average. (b) Cumulative average.
(S / v)1/2
25
(a) **
20 -
*o
15 .
.*** 4 *
10 *
5
0 I I "; ...; .. .. .. .......... ; .....I. .. .
0 5 10 15 20 25 30 35 40
Turbulence realizations
(s / v)"2
14.2
*
14.0
13. (b) *
13.8 6 *
13.6 \ ^
13.4
13.2
13.0
12.8
12.6 ----
0 5 10 15 20 25 30 35 40
Turbulence realizations
Figure 2.3 Turbulence small eddy shear rate (& / v) 1/2 against turbulence realization.
(a) Variation of ensemble average. (b) Cumulative average.
no / n
1.06
1.05
Smoluchow ski
1.04 ... Present Siulation
1.03-
1.02
1.01
1.00 I --I
0 1 2 3 4 5 6 7 8
T
Figure 2.4 Inverse of total number of monodisperse particles against time in laminar
shear flow. F=1.0, D=0.015.
For "HARD SPHERE" particles, the time step for advancing the motions of
particles is much less than 0.01. It is found that At=--0.001 may generate almost same
result as At=0.0005 for large inertia particles. For small inertia or zero inertia particles,
the collision kernel strongly depends on time step due to multiple collision mechanism,
which is shown in figures 2.5 for uniform shear flow. A pair of particles can have tens or
hundreds of collisions depending on the size of time step before they completely separate.
Similar multiple collisions have been found for particles in turbulence. Due to this
technical difficulty, the "HARD SPHERE" collision model is limited to particles with
large inertia.
V2I
V2
Figure 2.5 Multiple collisions of HARD SPHERE particles in laminar shear flow.
CHAPTER 3
A NEW THEORETICAL FRAMEWORK FOR PREDICTING COLLISION
RATE OF SMALL PARTICLES IN GENERAL TURBULENT FLOWS
A new theoretical framework is developed to predict the collision rate and collision
velocity of small, inertialess particles in general turbulent flows. The present approach
evaluates the collision rate for small particles in a given instantaneous flow field based on
the local eigen-values of the rate of strain tensor. An ensemble average is applied to the
instantaneous collision rate to obtain average collision rate. The collision kernels
predicted by Smoluchowski (1917) for laminar shear flow and by Saffman & Turner
(1956) for isotropic turbulence are recovered. The collision velocities presently predicted
in both laminar shear flow and isotropic turbulence agree well with the results form direct
numerical simulation for particle collision in both flows. The present theory for
evaluating the collision rate and collision velocity is also applied to a rapidly sheared
homogeneous turbulence. Using the mean turbulent shear rate (6/v)1/2 as the
characteristic shear rate to normalize the collision rate, the effect of the turbulence
structure on the collision rate and velocity can be reasonably described. The combined
effects of the mean flow shear and the turbulence shear on the collision rate and collision
velocity are elucidated.
3.1 Introduction
Smoluchowski (1917) considered collision rate among spherical particles in a
laminar shear flow, (ux, uy, u.) = (Fy, 0, 0) with a constant gradient F. For a target
particle of radius ri centered at an arbitrary position xi, any particle with a radius rj
moving toward the target particle will cause a collision if
Ixi xj| < Rj = ri +rj (3.1)
is satisfied. The collision rate can be evaluated as
= ninj fwr
"cr ^ ^o(3.2)
where wr is the radial component of the relative fluid velocity between xi and xj, wr<0
indicates an impending collision, -nln2wr is the particle flux moving inward to the target
particle, and A is a spherical surface with radius Rj-ri +rj. Using standard spherical
coordinates (R, 0, 4)) centered at particle i (see figure 3.1), wr on the spherical surface of
R=Rij can be written as Wr = FycosO = Rij sinO cosO cos). With dA= Rj sinO dOdo, it is
easily seen that
4li ~~~F R3 (3.3)
43 i
ncij nnj 3 r Rij (3.3)
or
Fij RR. (3.4)
3 i
Saffman & Turner (1956, hereinafter referred as ST) presented a classical theory
for the collision rate of small droplets in Gaussian, isotropic turbulence. The collision
rate among particle size group i and j is
c= ninj fWr. (3-5)
The ensemble average, denoted by < >, is necessary since turbulence is random. Because
of the continuity of the fluid flow, the volume influx entering the spherical surface of
radius Rj is equal to the efflux so that
nc1= ninj < entire Iwrl dA>. (3.6)
cij 2 n n entire sphere r
ST further interchanged the integration with the ensemble average so that ncij can be
evaluated in isotropic turbulence,
ncij 2 ninj entire sphere dA. (3.7)
-Wr y
0
I,
~C
u = Fy x
/ 4
,
z
Figure 3.1 Local coordinates for collision rate calculation in laminar shear flow.
For isotropic turbulence, the above becomes
ncij = ninj 27tRR (3.8)
where x is the local coordinate parallel to the line connecting the centers of the colliding
particles as shown in figure 3.2. For particle size smaller than Kolmogorove length scale
il, the relative velocity upon the collision was evaluated using the standard results of
small-scale turbulence theory as
= Rij< 01 > = Rij i 12 (3.9)
where s is the dissipation rate of turbulence and v is fluid kinematic viscosity. This gives
3 8n s ]/2 3 6 1/2
ncij = ninjRij [ ]v =.2944nnR () (3.10)
or aij=1.2944R (S 1/2 (3.11)
or C,3 =1.2944 (3.12)
oij (3-!)
*
where aii is normalized collision kernel.
V2
x
VIX V 2x = Wx
Figure 3.2 Sketch of particle collision velocity. The separation
distance is x2-xI = Rij= (ri+rj).
It needs pointed out that from equation (3.6) to equation (3.7), the interchange
between the integration with the ensemble average is only an approximation even in a
Gaussian, isotropic turbulence. Furthermore, the interchange between these two
operations severely limits the extension of ST's result and approach to other turbulent
flows such as the near wall region where the mean flow gradient contributes to . Of
course, in the theory of Smoluchowski for laminar shear flow, this issue does not arise
since the ensemble average is not needed. In a typical industrial facility where particle
collision leads to the desirable growth of particles, the inhomogeneity of the turbulent
flow structure is common. Due to the effects such as turbopheresis, particles tend to be
driven toward wall regions to cause a further inhomogeneity in the particle concentration
distribution. It was not clear how one can predict the collision rate of small particles in
the presence of both strong mean flow shear and high turbulence intensity (Kontamaris,
1995).
In this chapter, a new theoretical framework is presented for the collision rate of
inertialess particles, whose sizes are smaller than the Kolmogorove length scale, in
general turbulent flows. The present analysis starts from equation (3.5). Using the
leading order term, (xi-xj)-Vu, in the Taylor series expansion to approximate the relative
fluid velocity, w, at the collision instant, the collision rate for an arbitrarily given Vu is
evaluated first. The ensemble average is then carried out in a Gaussian, isotropic
turbulence to obtain the collision rate and collision velocity. Direct numerical
simulations are also performed to obtain the collision rate and collision velocity in a
Gaussian, isotropic turbulence. Good agreement is obtained between the present
*
prediction and the simulation for the collision kernel cai Excellent agreement is
obtained for the average collision velocity between the prediction and the simulation.
While the present result for the collision rate agrees with the prediction of ST, the present
theory predicts an average collision velocity <-Wr> that is 1.58 times the value given by
ST. The difference is caused by the bias of the collision toward higher values of collision
velocity which was not considered in previous studies on flow induced particle collisions.
Using the rapid distortion theory (RDT) for a rapidly sheared homogeneous turbulence
and the random Fourier modes representation, the dependence of particle collision kernel
aij on the mean flow shear rate and the turbulence structure in homogeneous shear flow is
studied under the present theoretical framework. It is found that the collision kernel, after
normalized by particle volume and the turbulence mean shear rate characterized by the
dissipation rate, mainly depends on the ratio of the mean flow shear rate to the turbulence
mean shear rate.
3.2 Analysis
The relative fluid velocity, w = u1 -u2, between target particle centered at x, and
the colliding particle centered at x2 in space with Ix, x21 = Irl2 = R less than the
Kolmogorove length scale -q can be expressed using Taylor series expansion as
w = U2-U1 ; r12.Vu +... (3.13)
r12 = X2- Xl (3.14)
to the leading order. The collision rate, ncl for a given Vu at a given instant is thus
nc12 =nn2 Wr< -w-ndA (3.15)
where n is the outward normal of the surface given by Ix2 X|i = Rj = R and w= w-n.
For a spherical surface, n= r12/R so that
J 1
-w-ndA = f -wJrl2 dA
wr;
SR Ow, 12 r1- r2 2dA. (3.16)
Representing Vu as the sum of a symmetrical part (rate-of-strain tensor) and an anti-
symmetrical part (vorticity tensor), it can be easily seen that the antisymmetrical part of
Vu does not contribute to wr*. Physically, it is simply because a rigid body rotation does
1
not contribute to normal velocity. For a symmetric tensor E = I (Vu +(Vu)T) or a
symmetric matrix Eij, a linear transformation, which is a rotation of the coordinate system
from (x, y, z) to (x-, z ), can be easily found to reduce Eij to a diagonal form as
=[eR 0 0 -
E = 0 ey 0 (3.17)
0 0 e2_
where (ek ey eZ ) are the eigen-values of E and are the principal values of the rate-of-
strain tensor. Due to incompressibility, ez + ey + ez = 0 so that there are only two
independent parameters (say ez and ey ) in E in order to evaluate the integral in (3.16).
Arranging the eigenvalues in the order of
ez =-(ex + ey) _
and defining
_y
-ewith -0.5<4<1 (3.19)
the integrand in (3.16) can be expressed as
r12 -Vu-rl2 = (ez R 2 + ey y 2 + e- j 2)
=- eR [-X 2 + Cy 2 -(1+g j 2]
= ex R2[cos2O+ sin20cos2-( l+)sin20sin2o] (3.20)
where a local spherical coordinate system (x =RcosO, Y =Rsin0cos+, z =Rsin0sino) has
been introduced. Hence
2itr
J Wr0-w-ndA= -R3ez ff <0[cos20+ 4sin20cos20-(1+-)sin2Osin2O] sinO dO do
W <0 JjWj.<0
fr W
0 0
=R3e e (9). (3.21)
Again, Wr<0 in the integration limit implies that only the non-positive values inside the
square bracket are evaluated in the integration. While an analytical expression for J9() is
difficult to obtain due to the dependence of the integration limits on arbitrary 4, accurate
numerical integration for Y(Q) can be easily carried out (using Simpson's rule). The
following interpolation is then constructed for (4) based on numerical integration
results,
8
8( [0.90725 + 0.2875(0.5+4)2 + 0.333(0.5+4)4] <0
8
,j [I + 0.554 + 0.7240712 0.9687443
+ 0.739244 0.230355]. >0 (3.22)
F 8
Consider the laminar shear flow for example. In this case, 4=0, eR = 2, H(0) = 3, so that
1 33
nc2 = n1n2 FR = nln2 FR3 which is the same as given by Smoluchowski (1917).
Equation (3.22) is valid for -1_<4-0.5 as well since it is based on (3.21). The restriction
for 4>-0.5 results only from equation (3.18) when the eigenvalues are arranged in such a
manner. For -1<_< 0, JY() is symmetrical with respect to C =-0.5. For example, 4=0
corresponds to a 2-D stagnation flow in (Rx z ) plane and C=-1 corresponds to a 2-D
stagnation flow in (R Y ) plane so that the collision rates are completely equivalent. It is
also interesting to note that c=-0.5 is equivalent to 4=1 since both correspond to an
axisymmetric stagnation flow with 4=-0.5 for an axisymmetric contraction and C=l for an
axisymmetric expansion; the only difference is the scale ek With these examples in
mind, a simple analysis for the flow field quickly shows that M1(_<0) can be obtained from
1;2f0) as
1--y )= I(1-- -1)-1+y :() for 4>0. (3.23)
Finally, the collision rate for small inertialess particles in a turbulent flow is simply
nC12 = nn2R3 < ej 9(4) > (3.24)
where the ensemble average is taken over all possible values of the rate-of-strain tensor.
It is emphasized that this result is applicable to general turbulent flow. Needless to say,
the isotropy assumption is not needed in the present framework. When the detailed
knowledge of Vu is available, such as from DNS or random Fourier modes
representation (Kraichnan, 1970) of turbulence, the ensemble average in (3.24) can be
carried out.
The foregoing analysis can also be extended to evaluate the average of the normal
(or radial) component of the collision velocity <-Wr>. For a given position with arbitrary
eR and 4 at any given instant, the instantaneous average of the normal collision velocity
can be calculated by integrating over the spherical surface of radius R as
0wr
= Wr< (-wr)dA (3.25)
In the above, the superscript 4. on the left hand side denotes the area-average, (-Wr)dA in
the numerator can be interpreted as proportional to the probability of a particle striking
the target surface at a position (0, i) over an area dA and wr in the numerator is the
relative radial velocity of the collision at (0, 0) on the target surface. The denominator is
the normalizing factor for the probability. The product, wr(-wr), clearly indicates that
A is a biased average in favor of these regions in (0, 0) where (-wr) is high.
Suppose that NR realizations of turbulence be taken in a turbulent flow, the total number
of collisions will be
Nc = NRnIn2< .r< (-wr)dA>
per unit time using equation (3.5). The chance of a particle striking the target surface
over an area dA at (0, 4) relative to the axis of the principal direction R of the
instantaneous rate-of-strain tensor in a given realization is then nin2(-wr)dA|w<0/Nc.
Integrating over the target surface and summing over NR realizations, the ensemble
average, , is obtained,
Snin2(-Wr) Wr
W> =R < Wr Nc d> = (3.26)
The above can be simplified as
<2 9(e ) >
=-R < > (3.27)
r
where
2n
( = f f [cos20+ 4sin20cos20-(l+)sin2Osin2w]2r<0 sinO dO do. (3.28)
0 0
The integration can again be evaluated numerically and represented in the following
piecewise interpolations
"G =0.4+1.0777(0.5+4)2 0.64767(0.5+4)4 for 4<0
=0.62831 + 0.68495 0.0910832
+0.09030143 0.0346584 for C>O. (3.29)
2
Similarly, the mean square value of the radial component velocity < wr > can be
evaluated as
< e3 Wf) >
= R (3.30)
with
=0.17114+l1.2643(0.5+)2-0.48123(0.5+;)4 for 4<0
=0.45714 + 0.976554 + 0.396062
+ 0.072430(3 0.0291C4 for _>0. (3.31)
The functions 4), q(4)/I(), and 7{)/4( ) are shown in figure 3.3. Take a laminar
shear flow for example: ez =172 and 4=0 so that G/H= 0.62832 (=-d5) and G/H=-16/35.
It is easily seen that
F2 n8 18 7t
=-R-jj/[|]=-0 FR, (3.32)
2 2F3 168 18 4 ,
=R 835 3/[3]=-(FR)R2 (3.33)
The results of these two moments of wr can be used to validate the present theoretical
formulation by comparing with direct numerical simulation for particle collisions.
0.04-
-0.5
0.0 0.5
Figure 3.3 Variations of i(, G(), K-(Q. ) with where
fl-38 ](),f2= ) andf3= .
To apply the above results to turbulent flows, we consider the following two
cases: i) Gaussian isotropic turbulence; and ii) rapidly sheared homogeneous turbulence
under a mean flow gradient F. In both cases, the turbulence can be represented using
random Fourier modes and ensemble averages can be obtained without resorting to the
direct numerical solutions to the Navier Stokes equations for the fluid flows. Pertinent
comparison can be made with ST's prediction in the first case. Insights on the effect of
mean flow shear can be gained by examining the collision rate in the rapidly sheared
turbulence.
3.3 Numerical Simulation
The detailed representation of the isotropic turbulence and rapidly sheared
homogeneous turbulence are given in Chapter 2. In short, the isotropic turbulence is
represented in the form of
N
Ui(X, t)= Y[b m) cos(k(m) .x +tom) ( t )
m=!
+clm) sin(k(m). x + o(m) .t ) ] (2.6)
where k(m)and ((m) are the wavenumber and frequency of the m-th mode. The random
coefficients b1m) and cm) are their corresponding amplitudes. For sheared
homogeneous turbulence, the flow field is given by
N
Ui(X, t)= E[bim)(S)cos( (m) -x+to)(m)t )
m=!
+ cm)(S)sin(X(m) -x +o(m) -t)
+ Fx2 8i (2.18)
where X(m) is the new wavenumber after the distortion by the mean flow shear, S = Ft is
the total shear by the mean flow, and b7m)(S) is the amplitude of the m-th random mode
after the rapid shearing.
The scheme developed in Chen et al. (1995) is employed here to search for
particle collisions. Detailed description on the collision algorithm is given in Chapter 2.
Periodic boundary conditions are employed for both turbulent flow and particle motion. A
number of particles (Np) are randomly introduced in the flow field of volume n3. The
average collision velocity is obtained based on all the collision events for whole time
period and all turbulence realizations.
As discussed in Chapter 2, four different post-collision treatments will give
different collision rate and it is not clear which treatment is the best. In the study
throughout this chapter, the 'THROW AWAY' scheme is implemented. However, an
extremely low particle volume concentration is used. The initial particle volume
concentration is chosen carefully so that the average number of collisions in each
turbulence realization is less than 1 The post-collision treatment thus has little practical
effect on the collision statistics. A large number of turbulence realizations are used to
obtain sufficient ensemble average. A subsequent cumulative time average smoothes out
statistical noise in the collision rate.
3.4 Results and Discussions
3.4.1 Comparison with Simulation Results in a Laminar Shear Flow
In an earlier paper dealing with particle collision a laminar shear flow (Hu & Mei,
4 .
1998), comparison for the collision kernel, 1 FR3, based on Smoluchowski's prediction
'3
and the direct numerical simulation was presented; excellent agreement was obtained. To
confirm the present analyses in the laminar shear flow, we compare the collision velocity
given by equations (3.32-3.33) with that based on the direct numerical simulations. Zero
inertial particles of radius 0.009, 0.0125 and 0.02, respectively, are introduced into a
cubic box with length Lx=2, Ly=2.0, and Lx=2.0, with the restriction 0.1
uniform shear flow with a shear rate r=I.0 is imposed and the flow velocity field is given
by v=(u, v, w) = (Fy, 0, 0). Periodic boundary conditions are employed in the streamwise
and spanwise directions. Particles near y=0.1 and y=1.9 have less collisions because
there are no particles in regions of y<0.1 and y>l .9 in the simulation. A correction due to
this boundary effect (Wang et al. 1997) was needed for the collision rate (Hu & Mei,
1998). However, this boundary correction has little effect on the collision velocity since
the decrease in the collision kernel affects both the numerator and denominator in
equation (3.25). Three hundred time steps with a step size At=0.01 are used to advance
particles and a total of 40 realizations were employed to further increase the statistical
accuracy. Ten thousand particles are used for r=0.009, 0.0125 while four thousand
particles are used for larger size with r=0.02. Table 3.1 shows the comparisons for both
<-Wr> and . Since is the second order moment of Wr, excellent agreement
for these two moments indicates that the basis of the present formulation is sound. It is
worth noting that a naive, Eulerian-based estimation of the average collision velocity
4
based on the incoming particle flux, i FR3ni, and the available collision area, 2nR ,
would lead an erroneous value of FR instead of n7 FR for <-wr>.
Table 3.1 Comparison of the collision velocity among monosized particles
between prediction and direct numerical simulation in laminar
shear flow.
Particle <-Wr> <-Wr>
r r
radius simulation prediction simulation prediction
0.009 5.6421x10-3 5.6549x10-3 3.7018xl0-5 3.7029x10-5
0.0125 7.8638x10-3 7.8540x10-3 7.1529x10-5 7.1429x10-5
0.02 1.2425x10-2 1.2566x10-3 1.8012x10-4 1.8286x10-4
In a laminar shear flow, Wr=-UxCOS0 and the collision angle Oc between the
collision velocity w and the axis connecting the center of colliding particle pair,
n/2_
easily evaluated so that the pdf for 0c is readily obtained. Since (-Wr)wr< sin0dOdo is
proportional to the probability of the particle striking the target surface at (0, )) with a
solid angle do)=dOd), it is readily seen that
(-wr)dA = -FR3sin2OcosO cos4) d0d4) for (0, 4)E Q,
= 0 elsewhere (3.35)
where Q is the region defined by (2 <0_
Wr<0. Interpreting sin20cos0sino as proportional to the probability density function of
incoming particles striking the target particle near (0, )), further integrating over 0 from 0
to 27t, treating 0
normalization, one obtains the pdf for the collision angle Oc
p(Oc) = -B sin2 ccos0c for in<_O
=0 for 0<09
where B is a normalizing factor which is 3 for 0 in radian and t/60 for 0 in degree. The
average angle is =f Ocpdf(0c)dec = 2.237 radian or 128.17. The numerical
7t2
simulation gives = 128.13. Figure 3.4 compares the pdfs based on the above
prediction and the direct numerical simulation using ten thousand particles of radius
0.009 over 40 realizations of initial particle positions. Excellent agreement is observed
for the whole region it/2<0c< The present theoretical framework is thus validated in the
laminar shear flow.
0.025 ....-,------
0.020
0.015
0.010
Analytical prediction
0.005 Simulation
0.000 # -,--
90 100 110 120 130 140 150 160 170 180
Collision angle 0 c
Figure 3.4 Comparison of the pdfs for the collision angles 0c between
the theory and the simulation in a laminar shear flow.
3.4.2 Collision Rate and Velocity in a Gaussian, Isotropic Turbulence
For particle collision in isotropic turbulence which is represented using equation
(2.6), the dimensionless collision kernel
aIl
S- R3(e/v)12 (s/v)1/2 (3.37)
is evaluated using NR=500 realizations. The dissipation rate E is based on the ensemble
average of the same NR realizations. In each realization, ez W/(4) and E/v are computed in
the flow field with 654 statistically independent data. It is found that a1I varies slightly
with the number of the Fourier modes, Nk, used in the turbulence representation.
Furthermore, it weakly depends on the power g in equation (2.10) with which the random
*
wavenumber is chosen. Figure 3.5 shows the variation of an as Nk increases from 100
to 4096 for t=1.2 and 1.4 for ri 0= 0.05 (Rex=53.8). The result seems to converge after
*
Nk=2048; a, 1(P.=1.2) converges to 1.2935 while a, 1(t=1.4) converges to 1.2945 with a
difference less than 0.1%. They are both very close to 1.2944, the value predicted by ST.
It is worth commenting that due to the rapid decrease of the energy spectrum E(k) in the
high k range, not all wavenumbers that are randomly generated need to be used in
equation (2.10). For those high wavenumber modes that correspond to very small
amplitudes, they can be simply neglected without compromising the accuracy. Thus the
effective number of the random modes is smaller than the specified Nk. From the
simulations, the effective number of the modes for g=1.4 is roughly three times that for
g=1.2 for the same Nk. Hence the fact that it takes larger values of Nk for p.=1.2 than for
ji=1.4 to get a converged result is not surprising. It is also interesting to note that, using
Nk=1024 and g=1.2, a higher value of Rex (=124.4 using fi 0=0.01) gives a,1 =1.2899 as
*
opposed to al =1.2912 at Rex=53.8. This clearly shows that shape of the energy
*
spectrum E(k) in isotropic turbulence has little effect on a, With Nk=200 and .=1.2,
*
a1 =1.2756 which is only 1.45% less than the converged result. Thus in the following
direct numerical simulation for particle collision, Nk=200 and p=1.2 will be used.
To validate the present analyses, direct numerical simulations of particle collision
in isotropic turbulence are performed. The turbulence in the periodic box of length n has
the following characteristics: root-mean-squared turbulence velocity uo=l1.0, integral
length scale L1 =0.594, Taylor micro length scale X=0.277, Kolmogorov length scale
71=0.023, Kolmogorov time scale Tk=(v/)l/2=0.073, and Re,=40.1. With the modified
Fourier modes in the periodic box, the dimensionless collision kernel based on the
*
prediction is ac =1.2711 as opposed to a11 =1.2756 without the modification of the
wavenumbers. Initially, five hundred particles of radius r=0.005 (r/ril=0.217) are
randomly distributed in the box of volume n3 with a particle volume concentration of
8.44xl0O6. These inertialess particles are advanced using local fluid velocity. A total
time period of T=6 (i.e. T/Tk=82.2) is simulated with time step At=0.01. The average
number of collisions within the entire period of the simulation 0
realization is about 0.4271<1 based on our new expression of the collision rate in 6.2.
To obtain a reliable statistics, a total of 10,000 realizations is used to generate an
*
estimated collision pairs Nc=4271. Figure 3.6 (a) shows the simulation result for ctaI in
0
numerical simulation data due to the small number of collisions within each time step
despite the fact that 10,000 realizations have been used. A cumulative average is applied
to reduce the noise. Figure 3.6 (b) shows the cumulative time average of cc1 from the
simulation. The direct numerical simulation is less than the predicted value (1.2711) by
about 2.5%. This difference may be caused by a small effect of particle size, since the
theoretical predictions are expected to be accurate only in the small size limit, i.e. r/r<
while the particle size in numerical simulation is r/rl=0.217. The effect of finite particle
size on the collision rate will be explored in Chapter 6. The agreement for o 1 between
the simulation and the prediction is quite good.
Table 3.2 shows the results of the prediction and the numerical simulation for the
two moments of particle collision velocity, <-Wr> and , in isotropic turbulence. The
predicted values are based on =1.2 and Nk=200 in the periodic box with NR=400
realizations. The particle sizes are r/rI = 0.2175 and 0.4875, respectively. Satisfactory
agreements are obtained. On the other hand, ST (1956) gave an expression, equation
(3.9), for the Eulerian-based average particle radial collision velocity. Using that
<-Wr>
expression, one would obtain -j
expression, one would obtain /v)2 0.206. The present prediction for the periodic
domain with NR=400, after ensemble average, gives a value of 0.328 for this
dimensionless radial collision velocity which agrees rather well with the direct numerical
simulation show in Table 2. Since the direct numerical simulation is based on the
tracking of particle trajectories, the probability-based prediction for <-w,> is consistent
with the Lagrangian approach. The Eulerian-based approach for <-Wr> does not take into
account the bias of the collision probability toward the high values of-wr; hence it gives a
lower value for <-Wr>.
cci
1.30
1.29
1.28
1.27
1.26
1.25 --
1.OE+2
1.0E+3
1.OE+4
Figure 3.5 Convergence ofacq as Nk --+8.
Table 3.2 Comparison of the collision velocity among monosized particles
between prediction and direct numerical simulation in Gaussian,
isotropic turbulence.
Particle <-Wr> <-Wr>
r r
radius (r/rhi) simulation prediction simulation prediction
0.217 4.5430x10-2 4.5611x10-2 2.7356x10-3 2.7296x10-3
0.4875 1.0093x 10- 1 1.0040x10-1 1.3405x10-2 1.3224x10-2
oa11
2.5
(a) -- siulation
....... Prediction
2.0
1.5
1.0-
0.5
0 10 20 30 40 50 60 70 80 90
t( / v)1/2
*
a11
1.8
1.7 (b) ------- Prediction
1.6 Cumulative average of
simtlation result
1.5
1.4
1.3
1.2
0 10 20 30 40 50 60 70 80 90
t(e / v)/2
Figure 3.6 Particle collision kernel ax from direct numerical simulation in Isotropic
turbulence. (a) Variation of the ensemble average. (b) Cumulative time average.
3.4.3 Collision in a Rapidly Sheared Homogeneous Turbulence
The statistics of the turbulence, such as four non-zero components of the
Reynolds stress tensor, generated using RDT and the random Fourier modes was first
checked against the analytical values at short times (Ftnl) given by RDT (Townsend,
1976, p.84) to ensure the correct implementations. It was found through simulation that
the development of the sheared turbulence with time does not depend on the Reynolds
number, Rex, of the initial state or the energy spectrum at t=0. The turbulence part of the
flow in the RDT is determined by the total shear S=Ft. The complete flow field, sheared
turbulence plus the mean shear flow, thus depends on two parameters: F and S=Ft. Since
the collision rate depends on the spatial structure of the turbulence while the turbulence
structure becomes increasingly anisotropic as the total shear S increases, finding an
appropriate representation of the collision rate in the homogeneous turbulent shear flow is
not trivial.
To develop a better understanding on the dependence of the collision rate on the
turbulent flow, it is instructive to examine the collision rate due to the sheared turbulence
but without the contribution from the mean shear F; this part of the collision kernel will
be denoted as ai1. Another word, a jj is evaluated using the flow field given by
equation (2.18) but without the last term [x25il which is deterministic. Following the
*
definition for all in the isotropic turbulence case, we use the dissipation rate (s/v)1/2,
which increases with S, as the turbulence characteristic shear rate to normalize t I I,
-* an!
c" = R3(/V) 1/2 (without contribution from mean shear F). (3.38)
*
It is easily seen that a,1 depends on the total shear S=Ft only. Figure 3.7 (a) shows the
-*
variation of &1 on S. In the high S limit, the simulation results suggest that a11
approach a constant of 1.1676 in the form of
a1 -1.1676+0.659/S forS>10. (3.39)
In this rage of 0
turbulence intensity becomes dominated by the longitudinal component in the x1-
direction. Based on RDT prediction which neglects the nonlinear interaction among
different scales, the longitudinal component of the turbulence makes up 92.7% of the
total kinetic energy at S=Ft=20. The anisotropy invariant IlIb reaches a value of -0.25 at
S=20. At such large values of total shear, S, the RDT results may become a little
inaccurate in comparison with DNS result. However, for the present work, since the
RDT-based turbulence is strongly anisotropic, it serves our purpose to evaluate the effect
of such anisotropy in the turbulence on particle collision rate. When the collision kernel
is normalized using (s/v)"2, the maximum difference in an is only about 15% at the
*
extreme limit of S-+co. If the isotropic result is used, &a, ~-1.2944, only a maximum
error of 10% will result from the limiting value of aE at very large total shear. This
*
encouraging result will be used to develop an approximation for a1 in the presence of
strong mean shear.
Now consideration is given to the particle collision kernel with contribution from
both the turbulence part and the mean flow shear part as the flow field is given by the
*
entire right-hand-side of equation (2.18). Figure 3.7 (b) shows aI as a function of the
ratio of the mean shear rate to the turbulence shear rate, F7/(- )/2, for a range of shear
V
rates. For a given F, (s/v)1/2 increases as t or S increases so that the abscissa mainly
reflect the effect of t or S. To see the effect of mean flow shear rate, F is made
dimensionless by the turbulence mean shear rate at t=0 in the isotropic state (eo/v)1/2.
Four values of F/(eo/v)1/2 are used: 0.87, 2.18,4.36, and 8.72; they correspond to the solid
symbols in figure 3.7 (b). The open symbols specifically correspond to the cases when
S=0 (or t--=0) so that the turbulence is still isotropic. The dimensionless values of aoI
seem to collapse reasonably well. At large values of r/FI( )1/2, the curve becomes linear
V
as the collision rate is increasingly dominated by the mean shear so that aI scales
linearly with F as Smoluchowski's theory predicts. A small degree of scattering exists for
F/( )/2<2. This is caused by the 15% maximum variation of a11I over the whole range
of total shear S, as already demonstrated in figure 3.7 (a). There are two possibilities for
F/(& )1/2 1: i) both F and S=Ft are very small even comparing with (So/v)1/2 of the initial
isotropic turbulence; ii) F is finite but S=Ft is large so that (E/v)1/2 become large
comparing with (so/v)/2. In the first case, the turbulence is closer to the isotropic state
_* *
so that a11 -1.2944. In the second case, the turbulence is highly anisotropic and a&,
-+ 1.1676 for very large values of S. No attempt is made to correlate d& with the total
shear S since S is not relevant in practical applications such as in turbulent pipe flow.
_*
Using ,11 =1.2944 in the isotropic state to represent the turbulence part of the collision
kernel, a simple interpolation is proposed to represent the effect of mean shear and
turbulence on the collision rate,
2.2 22 1.3333 22 1/2.2
ax 1 R3 (/v)1/2 1 {1.2944 + [ (/v)1/2 I 2 (3.40)
This interpolation is shown by the solid line in figure 3.7 (b). A satisfactory agreement
can be observed. It is also worth commenting that when F is normalized using the eddy
2
turn over time, 3u0 /&o, of the initial isotropic turbulence, the shear-rate parameter
F* F3u/o0, (3.41)
corresponds to F*=36.98 and 92.45 for F/(Eo/v)1/2 =0.87 and 2.18. In a low Reynolds
number turbulent channel flow, the maximum value of F* is around 35 at a location of 10
wall units away from the wall (Kim et al. 1987). Hence, the parameters used in figure 3.7
(b) are relevant to practical situations. Although a 10% error exists for a11 given by
(3.40), the advantage of this simple approximation is that the total shear S=Ft is absent
and it depends only on the spatial structure of the turbulence characteristics by the scalar
quantity (e/v)"12.
Figure 3.8 shows the radial collision velocity, <-Wr>, as a function of F/( )1/2, for
(V
different values of F/(' )1/2. Similarly, the data seem to collapse approximately into one
V
curve in the range studied. The following interpolation fits the radial collision velocity
well,
<-Wr> 1.7 1 r F 17 1/17 (.1
Reiv)/ {0.32451 + 1 ./v)12 } 41)
(eiv'/ 6v~l2](.
For the root-mean squared radial collision velocity, , the behavior is very similar to
r
that of <-Wr>. While the data is not shown here, the following interpolation, which fits
the data quite well, is provided,
1/2
r____ 1 7 0.3381 F 17 1/1.7 3
R(e/v)"2 {0.3681 + (e/v) } (3.42)
The application of the above expressions to turbulent channel flow, which is
different from the homogeneous, rapidly sheared turbulence, is yet to be tested. However,
since the rapidly sheared turbulence at large total shear S considered in this work is
already drastically different from the isotropic turbulence, there are reasons to expect a
similar dependence ofoa 1, <-Wr>, and on F and (e/v)1/2 in a turbulent channel flow
which is of significant industrial importance.
3.5 Summary and Conclusions
A theoretical framework has been developed to evaluate the collision rate and
collision velocity of small particles in general turbulent flows. The classical results for
the collision kernels by Smoluchowski (1917) for laminar shear flow and by Saffman &
Turner (1956) for isotropic turbulence are recovered. The present approach differs
significantly from that of Saffmnan & Turner in that the ensemble average is taken after
the collision rate for a given flow realization is calculated. This avoids the assumption of
isotropy, as needed in Saffmnan & Turner, and allows for the evaluation of the collision
rate in general turbulent flows. The collision velocity predicted in both laminar shear
flow and isotropic turbulence agrees well with the direct numerical simulation for particle
collisions. The present theory is used to evaluate the collision rate and velocity of small
particles in the rapidly sheared homogeneous turbulence. It is found that the effect of the
turbulence structure on the collision rate and velocity can be mainly captured by using the
mean turbulent shear rate (F/v)112. The combined effects of the mean shear and the
sheared turbulence on the collisions are elucidated. Simple interpolations are developed
to evaluate the collision rate and velocity for arbitrary mean flow shear rate and
turbulence mean shear rate.
1.4
1.3 -
1.2 -
I ........ ........ I.
10 1
100
102
S= Ft
Figure 3.7 (a) Normalized particle collision kernel in a highly sheared homogeneous
turbulence without contribution from the mean shear rate.
*
0
*
0*
0*
*
*
1.1676
6-
5-
4
4 x
("" /v) 1/2
(S )
3 0
-= 0.87
elOx 2.18
2-
"X x 4.36
8.72
1 0 S=O =0or t=0
Interpolation
0 1 2 3 4
(d/v) 1/2
Figure 7 (b) Collision kernel in a rapidly sheared homogeneous turbulence.
cO denotes the dissipation rate of the isotropic state at t=0.
1.5 -***---
S1.0
~F
^ >^ r
71/2
S 0C /V~
0.5 0.87
2.18
X^5 x 4.36
8.72
0 S=0 or t=0
--- Interpolatior
0.0 ,-,,- i
0 1 2 3 4
F
(E/V) 1/2
Figure 3.8 Collision velocity in a rapidly sheared homogeneous turbulence.
c0 denotes the dissipation rate of the isotropic state at t=0.
CHAPTER 4
EFFECT OF INERTIA ON THE PARTICLE COLLISION VELOCITY AND
COLLISION ANGLE IN GAUSSIAN TURBULENCE
Numerical simulation of particle collisions in isotropic, Gaussian turbulence is carried out
to study the effects of particle inertia on the turbulence induced collision velocity and
collision angle. An asymptotic analysis is developed to predict the effect of the small
particle inertia on the collision velocity. A PDF-based method is also implemented to
predict the effect of the small particle inertia on the collision velocity. The analysis
predicts that the increase in the collision velocity is proportional to the terms (P3Tk)-I
and (0Ck )-2, in which p-1 is the particle response time and Tk is the Komogorov time
scale. Good agreement is obtained among three different methods: the numerical
simulation, asymptotic prediction, and the PDF-based method, for the collision velocity
of particles with small inertia. The effects of both particle inertia and size on the collision
angle have also been investigated. It is found that collision angle can be significantly
affected by both inertia and size. Average collision angle and the variation of collision
angle increase as the inertia increases, while the decrease in the size will also reduce both
the average collision angle and the variation of collision angle.
4.1 Introduction
Particle collisions can take place through a variety of collision mechanisms.
Particles may be brought together by inertial centrifugal and gravitational forces. The
velocity gradient is responsible for collisions among particles in laminar shear flow and is
a dominant factor for particle collisions in turbulence. Brownian motion and Van der
Waal forces can also affect the collisions. One of the common characteristics of these
collision mechanisms is the ability to generate a negative relative velocity between two
particles along the line connected with two particle centers. Of significance are the
statistics of the relative velocity of two particles at contact (hereinafter referred as
collision velocity) and the relative angle between the velocities of two particles at contact
(hereinafter referred as collision angle). For solid particles, collision causes momentum
exchange among particles and may influence the local and global particulate flow
conditions. The suspended solid particles are produced during the manufacturing process
of materials, the combustion of heavy fuel and pulverized coal. The collisions among
suspended particles will greatly affect the properties of particulate flow and thus influence
the performance of the systems. Adhesion may occur when two colliding solid particles
have small collision velocity, while large collision velocity will produce large restitutive
rebounce velocities for the colliding particles and may overcome the bonding force
determined by cohesive surface energy. The collision also affects the angular momentum
of particles. Tangential velocity components will be imposed to the colliding particles if
the collision angle 012 is not equal to zero. The particle rotation will be either enhanced
or reduced depending on the collision angle 012, which determines the orientation of the
particle tangential velocities upon collision. Gunaraj et al. (1997) performed
experiments measuring the collision angle 012 between heavy particles (silica sand
particles with 1000-1500 gm in diameter) in turbulent channel flow.
The collision velocity is also a critical quantity in determining the collision kernel
and collision rate. In uniform shear flow, Smolochowski (1917) integrated the volume
flow rate using collision velocity along the spherical surface of collision radius R
(= ri + rj, where ri and rj are the radii of collision pair respectively) and obtained
collision kernel for uniform shear flow. In turbulent flow, particle collisions are affected
by local turbulence fluctuations. In the absence of any other force creating relative
motion between particles, these local fluctuations entirely control relative velocity
distribution which determines the collision kernel. In most of the previous theories, the
collision kernel in a particle-laden system is linearly related to the collision velocity, the
first moment of the relative velocity distribution of colliding particles. Thus the scale and
form of these particle relative velocity distributions directly determine the collision rate.
The fact that the collision rate and collision kernel in previous studies are all proportional
to the collision velocity explains the critical importance of the prediction of collision
velocity (Saffinan & Turner 1956, Williams & Crane 1983, Kruis & Kusters 1996, et al).
The particle collision has been studied in two different ways. One is based on
Eulerian method, which is favorable to theoretical prediction and was followed by
previous studies on the theoretical evaluation of both collision velocity and collision
kernel. The other is to use Lagrangian tracking in direct numerical simulation of particle
collisions. For inertialess particles in both uniform laminar shear flow and isotropic
turbulence, the collection of collision velocities using Lagrangian method shows a biased
average in favor of those regions where wr is high, as shown clearly in Chapter 3.
Saffman & Turner (1956) obtained the collision velocity using average integration on the
spherical collision surface by Eulerian method without considering the influence of the
biased average. Although their formulation correctly predicts the velocity correlation
between two points for turbulent flow, it has been incorrectly interpreted as the collision
velocity upon collisions. The effect of particle inertia on the collision velocity was
recently studied also by Sundaram and Collins (1997). They used two statistical
variables: g(R), the particle radial distribution function at contact, and P(w/R), the
conditional relative velocity probability density function at contact, where w is the
relative velocity of two particles at contact. It is noted that in Sundaram and Collins'
work, in additional to the assumption of the exponential form for P(w/R), the g(R) was
evaluated based on their numerical extrapolation form for distances smaller than the first
numerically evaluated point. However, the usefullness of g(R) and P(w/R) is rather
questionable.
In this work, we study the effect of particle inertia on increasing particle collision
velocity and collision angle by carrying out numerical simulation of particle collisions in
an isotropic, Gaussian turbulence, and developing an asymptotic analysis to predict the
effect of small particle inertia on the collision velocity based on Lagrangian description.
An integration method is also implemented to predict the effect of the small particle
inertia on the collision velocity. The increase in the particle collision velocity due to
inertia predicted by both the asymptotic analysis and integration method is in agreement
with the numerical simulation results. The analysis predicts that the increase in the
collision velocity contains terms that are proportional to (PT k)- and (pxrk )-2. Good
agreement is obtained by the numerical simulation, asymptotic prediction and the
integration method for the collision velocity of particles with small inertia. The effects of
particle inertia and size on the collision angle have also been investigated. It is found that
collision angle are significantly affected by both inertia and size. Average collision angle
and the variation of collision angle increase as the inertia increases, while the decease in
the size will reduce both the average collision angle and the variation of collision angle.
4.2 Saffman & Turner's Analysis
In Saffman & Turner's (1956) prediction of the collision rate among inertialess
particles, the volume flow rate of particles towards a target particle was evaluated by an
integration, , where wr is the radial component of the relative velocity.
The integration was carried out over the entire spherical surface of radius R, as shown in
figure 4.1. Because of the continuity of fluid flow,
=
= Entire sphere <|wrI> dS. (4.1)
For isotropic turbulence, which is true at small scales for anisotropic turbulence, the
above becomes
= 27iR2 (4.2)
where x is the local coordinate parallel to the line connecting the centers of the colliding
particles as shown in figure 4.2 and < > denotes ensemble averaging.
63
A y
SCollision area
Sfor the target
F particle
x
Figure 4.1 Sketch of particle collision scheme
W = V -V2
wx = VIx V2x
Figure 4.2 Sketch of particle collision velocity and collision angles (012 & 0 J).
The separation distance is x2 xi = R = (r, + r2).
In the limit of zero particle inertia, particle velocity v follows the fluid velocity u
very well so that the radial component of relative velocity w = vi v2 upon collision is
Iwxl-lulx U2xl-R-l = R( --2 )/2 (4.3)
Ox 157t v
and =5 R2 <(U)2 > IR2E (4.4)
Sax 15 v
where wx is the collision velocity.
For particles with inertia, the volume flow rate of particles to a target particle was
expressed as Jffentire sphere Wl P(w)dw, where P(w) is the probability density function
(PDF) of w and contains the effects of the spatial variations of velocity in the fluid, the
turbulent accelerations due to the relative motion with the turbulence, and gravity. The
volume flow rate can be integrated assuming P(w) to be a normal distribution function.
For monodisperse particles,
JJentire sphere wl P(w)dw = 2 R3 2 12 (4.5)
3 v
Compared with the equation (4.2), the equivalent collision velocity based on equation
(4.5) can be evaluated
iWxl1 2R 3 ( 2rt)12/(2R2) = 2R(6_)l12 (4.6)
3 v 3 27tv
and - 1-R2- (4.7)
9 v
It is noted that equation (4.7) is independent of the inertia parameter 13. For the small
inertia limit 13--+oo, comparing equations (4.6) & (4.7) with equations (4.3) & (4.4), an
inconsistency that (4.6) & (4.7) can not recover (4.3) & (4.4) occurs, which is physically
incorrect.
4.3 Present Analysis
4.3.1 Asymptotic Analysis
To understand and predict the effect of inertia on the collision velocity wx, it is
assumed that the particle size is much less than the Taylor micro-length scale, and the
gravitational force is negligible. For simplicity it is assumed that turbulence is isotropic
and Gaussian. The particle dynamic equation in the absence of(or weak ) gravity is
dv
dt =3(u-v). (4.8)
The above can be written as
I dv
v-au (4.9)
13 dt
For large values of P, repeatedly substituting the left hand side into the right hand side
gives
I du I d2V I Du I( ) Vu+ I d2v
v-u ---+-- u---+-(u-v)UVu+---
P3dt D2 dt2 P3Dt 13 1P2 dt2
I Du I
-- +()(4.10)
f3Dt p3
The variance of w = vi V2x in which particle 1 and particle 2 are separated by a
distance of I x2 x, I= R is
1 Duix 1 Du2x
var(wx) var(Ulx U2x Du- )
Pi Dt P2 Dt
=<(Ulx U2x)2>1x2-x=R +<(I Dux 1- Du2x )2> -x=R
P13 Dt P32 Dt x2-xI=-R
1 Dux 1 Du2x
2 <(Ulx U2x)( Dt ^ Dt)>lx2-xl=R
01 Dt f02 Dt
= term I + term II + term III. (4.11)
The first term can be expanded as
term I=<(ul U2x)2 >lx2-x=R = +-2Ix2-x,=R
=2 2f(r=R) = 2u2 [1-f(R)] (4.12)
where u2 is the mean square intensity of turbulence and f(r) is the standard longitudinal
velocity correlation coefficient. If the particle size is much less than the Taylor micro-
length scale X, the above can be simplified to
term I= <(ul -U)2 >x2-x=R 2 (2 = IXR2v
U~,1~1,R 2U "2X 15 v
for R<
Term II in equation (4.11) can be similarly handled,
I u1D,~ Di _!I(Dux)2>~(X
term II-- <( -L 2 )2 =R 2 <(DUI)> + 1 <(Du2t> D
'P1 Dt P2 Dt I x 1-R p Dt / 2 vDt
-2 1- (4.14)
13P2 Dt Dt x2-x=+R44)
where
DIwhere < D U2X >1 X _X =R is the spatial correlation of ux- in the x-direction with a
Dt Dt x2 x iDt-
separation distance of R. Noting that var( Du--) = var( DU2 ) and expanding the
Dt Dt
correlation for small R, we obtain
[+
DUIDU 2-xI=R <(-x )2>[l+R2 a2fD(o) + ] (4.15)
Dt Dt 2 ax2 "'"
where fD is the correlation coefficient of Dux. Defining the micro length scale of the
Dt
fluid acceleration Du" as
Dt
I I a2fD() (4.16)
?! 2 x2(
term II can be expressed as
term 1-<(DUx)2>( l + 1 2 <(DUx)2>(1R2
Dt t F2 0P12 Dt 7 )
_<(DU)2>( I_ 1)2+ 2<(DU)2>R for R< (4.17)
Dt 1 0i ) + 0 Dt)T f D
Finally, equations (4.13) & (4.17) lead to
var(wx)- IlR2 +<(Dux)2>( 1 1 2 <(DU)2 R2
15 v Dt 01 02 P1P2 Dt TD
-2<( -u, )(- DUX )>. (4.18)
1 Dt P2 Dt
For the like particles, 31 = 02 = 1,
I 2c 2 /Dux Du2 2 (Dux 2 R2
var(w.)- 15R21 <(Ux-U2x)(- )> <(__ >
DDt Dt
(i) (ii) (iii)
(4.19)
where term (i) is the mean-squared collision velocity for inertialess particles, terms (ii) &
(iii) are contributed by the inertia. The net increase of the mean-squared collision
velocity due to the inertia can be expressed as
2 AWD U lx Du2x 2 Dux 2 ( 2
< 2 (u )-uX p<)+ R )2> (4.20)
SDt Dt
(ii) (iii)
It is noted that term (ii) represents the first-order effect of the inertia and term (iii)
contains the second-order effect of the inertia. The present prediction, equations (4.19) &
(4.20), differs from the equation (4.7) in two aspects. First, the collision velocity in the
limit p3-+oo is consistent with the equation (4.4). Second, the present analysis includes
two additional terms in the effect of inertia, as seen from equations (4.19) & (4.20) when
P1 = P2. The effect of inertia does not vanish when P1 = 02 while it does in Saffman &
Turner's analysis. It should be mentioned that in (4.20) is computed based on
Lagrangian tracking of all colliding particles.
4.3.2 PDF-Based Method
For inertialess particles, the collision velocity can be evaluated by the direct
integration (Mei & Hu 1998) in both uniform shear flow and turbulent flows.
Lwr
= Jwr
fwr
and
2 w < Jw 2(_WrdS
=
Jwr
where (-wr)dS in the numerator can be interpreted as proportional to the probability of a
particle striking the target surface over an area dS. The denominator is the normalizing
factor for the probability. The product, wr(-Wr), clearly indicates that is a biased
average in favor of these regions where (-wr) is high.
Equations (4.21) & (4.22) may be extended for particles with inertia. For small
inertia, as shown in 4.3.1, the particle velocity to the first order P is
v = u- (4.23)
The collision velocity may be expressed as
1 Du1x 1 Du2(
wr = wx = Vix V2x = Uix U2x- +-- (4.24)
P3I Dt P2 Dt
Jw<1 DuIx I Du2x )dS. (4.25)
r<0 wrdS w<0 f(U2 31 Dt 02 Dt
For monodisperse particles, PI = P3 = P3, equation (4.25) may be rewritten as
fWr
From the equations (4.21), (4.24) and (4.26), the collision velocity is then found
From the equations (4.21), (4.24) and (4.26), the collision velocity is then found
w u Dt Dt
S2 Jwr
< r= (DtlX- --- (4.28)
Jw
The integration in right hand side of (4.28) are evaluated using Eulerian method.
4.4 Results and Discussions
The direct numerical simulations of the particle collisions in Gaussian turbulence
is performed to assess present asymptotic analysis and to predict collision angle. Initially,
a given number of particles are randomly distributed in the box of volume it3. The
inertial parameter PT|k varies from 0.022 to oo. The particle sizes D/irl=0.435, 0.957,
1.739, 2.522, 4.348 are considered, and the corresponding initial number of particles are
24000, 10000, 5000, 3000 and 1000 respectively. The particle collision detection is
turned on until all particles have reached dynamic equilibrium. The overlapping particles
are eliminated to avoid the false counting of the collisions before the collision detection
turned on. The particles are advanced by time step At=0.01. The total time steps, varying
from 200 to 600 after the equilibrium state, depend on the inertia. A total of 40
realizations is used for smooth statistics.
Figure 4.3 shows the comparisons of direct numerical simulations, asymptotic
analysis and PDF-based method for D/n=0.957. Very good agreement is observed among
the three methods for small inertia pT) k >3.6, so equation (4.20) is validated for prediction
of small inertia effect on the collision velocity. However, large deviation is observed for
pT k <3.6 since the asymptotic analysis cannot be applied for particles with large inertia.
70
The reason is that the first order approximation (4.23) is no longer valid when the 1
becomes small for particles with large inertia. The PDF results of collision velocity can
be found in figures 4.4 (a) & (b). The collision velocities are all negative, because
collision occurs only if the pair of particles are approaching each other. For not very
large inertia, say PT[ k> I, due to the strong correlation between particle velocity and fluid
velocity, the PDF curves are sharp and narrowed with a small variance and peak at almost
the same location, while zero inertia curve is at highest level. The PDFs are sharper for
small particle size (D/rT=0.957) at pTk>1.023, shown in figure 4.4 (a), compared with
figure 4.4 (b). However, the PDF curves become flatter and the location of the peak
shifts to larger IwrI for Pxk <0.073 due to little correlation between particle velocity and
fluid velocity at large inertia.
1.0E+O
1.OE-I
S---0-.-- prediction, Eq (420)
-... simulation
-*- integration, Eq (4.28)
1.0E-2
1.OE-3 -.
1.O0E-4 .. ..
0 5 10 15 20 25
1rk
Figure 4.3 Comparison for the increased collision velocity due to the particle inertia
between numerical simulation and theoretical prediction, D/r7=0.957.
PDF
0.0 4-*
-2.0
-1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0
collision velocity wr
PDF
-2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0
collision velocity wr
Figure 4.4 Effect of particle inertia on the probability density function of the
particle collision velocity wr. (a) D/A=0.957. (b) D/rl=2.522.
The effect of particle inertia on collision angle among monodisperse particles can
be found in figures 4.5 (a-c), 4.6 & 4.7. It is seen that the collision angles 012 & 0c
strongly depend on the inertia. For not very large inertia, fT > 1.023, the mean collision
angle 012 is less than 20 0 at D/A=0.957, as shown in figure 4.5 (a) & 4.6, and it is less
than 30 0 at D/r9=2.522, as shown in figures 4.5 (b) & 4.6. However, for very large
inertia, P3XTk <0.073, the mean collision angle 012 is greater than 80 0, which reflects little
biased orientation of collisions due to very weak correlation between the particles and
turbulent flow. It is noted that the PDF curves for PTk> 1.023 are much sharper and
narrower than those for 13Pk__0.073 and the variation of collision angle increases as the
particle inertia increases, as shown in figure 4.7. Figure 4.5 (c) shows the effect of inertia
on 0C. It is seen that the increase in inertia results in the collision angle 0c approaching
to symmetrical distribution with the mean 45 0.
The effect of particle size on the collision angle 012 can be seen in figure 4.8,
which clearly shows that the particle size has a significant effect on the collision angle.
The average collision angle 012 increases as the size increases. For D/T<11.739, the
maximum 012 is located at somewhere between 0and 10, and the average 012 is
around 10. However, for larger particles with large inertia, shown in figure 4.8, the
collision angle 012 is around 20 0, which is within the range of 14 0 32 0 obtained from
the experiments for heavy particles in a turbulent channel flow (Gunaraj et al. 1997).
Since little information is reported on the turbulent channel flow, the comparison is only
qualitative.
PDF
0.120
*(a)
0.100
P3 k
I k
0.080 -----7.308
S---A--- 3.654
|----1.023
0.060 ---0.292
I ... 0.073
0.040 ----..- 0.022
!.. -- infinite
0.020
0.000 % ^ --. "-' -- i -, _,,..
0 20 40 60 80 100 120 140 160
Collision angle 012
PDF
0.050
0.045 ,
(b)
0.040 r
fPTk
0.035 -___ 7.308
0.030 -- --- 3.654
0 .2 -a-1.023
0.025 -\ 0.292
a.
0.020 ^ -- --0.073
I o--0-- 0.022
0015 -0- infinite
0.010 .,AAAAAWAQ%
0.005 saw~x:3^ ^
AA AAAAAON LAA
0.000 m e
0 20 40 60 80 100 120 140 160
Collision angle 012
Figure 4.5 Effect of particle inertia on the probability density function of the
particle collision angle 012. (a) D/11=0.957. (b) D/r|=2.522.
180
180
PDF
o,020
0018 (C) r"n--u -
0.016- 7 \
0.014 ,
0.012 m m
0010 -7.308
,,,',/= /'y y -- 7.308 .^^A
/
0.008 *-/L -- 3.654 "
mM} -1.023
0*006 -*-O.9/ 0.292
0.004 ---0 073
/' -a- 0022 ?
0.002 -o-- infinite
0.000 "
0 10 20 30 40 50 60 70 80 90
Collision angle 0,
Figure 4.5 (c) Effect of particle inertia on the probability density function of the particle
collision angle 0c, D/I=0.957.
612
0 1 2 3 4 5 6 7
PTk
Figure 4.6 Effect of particle inertia on the mean collision angle 012 in isotropic,
Gaussian turbulence.
45
40
35
30
25
20
15
10
5
0
1 2 3 4 5 6 7
Figure 4.7 Effect of particle inertia on the standard deviation of collision angle 012 in
isotropic, Gaussian turbulence.
PDF
0.10
0.09 -
0.08
0.07 D / T1
0.06 .-..-.- 4.348
----- 2.522
0.05 --,--1.739
0.04 ". -----0.957
0.03 ,, -0.435
0.03 -/ *';
0.02 a .
0.01 -
...- ../ .--u
0.01 '/ "" -- "*-. -" -... -----
0.00 I I I m U.
0 10 20 30 40 50 60 70 80 90 100
Collision angle 012
Figure 4.8 Effects of particle inertia and size on the probability density function of
the particle collision angle 012, P1x k =1.023.
4.5 Summary and Conclusions
Numerical and analytical studies on particle collision velocity and angle in
isotropic, Gaussian turbulence have been carried out. The effects of particle inertia and
particle size on the collision velocity and collision angle are examined.
An asymptotic analysis has been developed to predict the effect of small particle
inertia on the particle collision velocity. The asymptotic prediction for the increased
collision velocity agrees well with the results of numerical simulation and PDF-based
method. The collision velocity is found to contain terms that are proportional to (IPxk )1
and (PTk)-2.
Particle collision angle strongly depends on particle inertia and size. The average
collision angle and the variation of the collision angle increase as the inertia increases.
Both the average collision angle and the variation of the collision angle decrease as the
particle size decreases.
CHAPTER 5
QUANTIFICATION OF PARTICLE CONCENTRATION NON-UNIFORMITY
IN TURBULENCE
For heavy particles with finite inertia, small scale structure of turbulence induces particle
preferential concentration. Particles tend to collect in the regions of high strain rate and
low vorticity, as has been observed in many Lagrangian simulations. Quantification of
the extent of the non-uniform particle concentration distribution has been difficult
because the particle concentration based on the Lagrangian statistics is highly noisy. A
method is developed to calculate reliably the variant of the particle concentration, < n'2 >,
in which n is the instantaneous particle number concentration and the prime is the
fluctuation in a small cube of volume (Ax)3, from the Lagrangian simulations of the
particle motions in an isotropic turbulence. The concentration non-uniformity is
measured using a new parameter A0 based on the first two lower order moments of n. It
is observed that the non-uniformity varies with Ax, suggesting a fractal nature of the
particle concentration distribution. For particles with zero size, the non-uniformity scales
with Ax as An (Ax) -b (b around 1) for small values of Ax. However, for finite-size
particles, the non-uniformity reaches a maximum value (around order one or less) at Ax/D
= 2.0-3.5, in which D is the particle diameter; it drops to zero as Ax/D gets close to one.
It is noted that the maximum value of the /2 is attained when the ratio of
particle response time to the Kolmogorov time scale is close to one, which has been
previously observed for the zero-size case based on other physical quantities.
5.1 Introduction
In many statistical problems, the lower order moments are often of interest. For
large number of monodisperse particles, the collision rate is nc=tXn1n2. The first
moment of particle concentration is clearly the average concentration n, i.e., the number
of particles per unit volume. The second order moment can then be expressed as
=n2+=n2(1+
-2 .1
n
where < > denotes ensemble average. The first moment n does not tell how the number
of particles varies with the volumes at different locations and gives no reflection of
turbulence structure at small scale levels. However, the second order moment, ,
-2 2
consists of two terms: n and < n' >. The variance of particle concentration, < n' >,
measures the variation of particle number density in the cells, i.e. the non-uniformity of
the particle distribution. Thus it is fundamentally important to the particle collision rate.
Since nfi is a simple average of particle concentration in a turbulent flow while < n'2>
depends on the turbulence structure and particle inertia, etc., our attention is given to the
behavior of < n'2 >.
2 -2
In this study, we use the normalized variance of particle concentration, / n
to quantify the concentration non-uniformity in turbulent flows. The normalized variance
is obtained via direct numerical simulation. As will demonstrated in 5.2, the inherent
statistical noise is eliminated from the /n By this new method no assumption
regarding the particle-particle interactions is needed. The effects of particle size and
inertia on the concentration non-uniformity can be quantified separately. It is found that
the concentration non-uniformity can be greatly overestimated if the particle size is
ignored. For particles with zero size, the non-uniformity scales with Ax as
/2- (Ax) -b (b around 1) for small values of Ax for particles whose response
time is comparable with turbulence Kolmogorov time scale, where Ax is the length of a
small cube of volume (Ax)3 used in the Lagrangian simulations of the particle motions in
turbulence. However, for finite-size particles, the non-uniformity reaches a maximum
value (around order one or less) at Ax/D = 2.0-3.5, in which D is the particle diameter
and rf is Kolmogorov length scale; it drops to zero as Ax/D gets close to one. The
maximum value of the /< n > 2 is obtained when the ratio of particle response time
to the Kolmogorov time scale is around one.
5.2 Quantification Method
5.2.1 Inherent Noise of the Concentration Non-uniformity
Suppose the computational domain L3 is divided into N. equal-volume cells each
with length Ax so that N, = (L/Ax)3. Obviously n = N/Ne is the first moment of n.
For zero-size, inertialess particles, the probability of a particle residing in a
particular cell is p = 1/Nc,. With N particles in the domain, the probability of k particles
residing in a specific cell of volume (Ax)3 obeys the Bernoulli distribution,
N! N
P(k)_k(Nk! Pk(1- P)N-k. (5.2)
The second moment of n is, by definition,
2 N
= Zk'p(k). (5.3)
k=0
Thus the variance of the particle concentration for a given cell is
=-n = k2P(k)- n 2. (5.4)
k=0
For Np(1- p) >>1, P(k) given by (5.2) can be expressed as
P(k) L2N1 exp[-(k-Np)2/2Np(l- p)], (5.5)
-\27tNp(l -p)
for k in the neighborhood of Np, which is known as DeMoivre-Laplace theorem. In the
limit of large Np(l1- p), which is very close to n for small p, an asymptotic integration
gives
N -2
=j k2P(k) dk n +n(1-p). (5.6)
0
Denoting as for zero-size, inertialess particles (P-.oo), we obtain
=- n 2=n(-p). (5.7)
2 -2
Denoting r as the ratio of < n'2 > to n ,
r= -2 (5.8)
n
it follows for zero-size, inertialess particles that
r.- (5.9)
n n
For n =0(1) or smaller, the DeMoivre-Laplace theorem does not hold anymore. In DNS
the number of particles used is relatively small from the point view of smooth statistics,
as mentioned in Section I. In DNS, n may be much less than one as cell size becomes
small. With finite-size particles, the ultimate number distribution of particles is reduced
to either zero or one per cell as the cell size is smaller than 0.5773D, where D is the
particle diameter. For n <<0(1), particle number concentration is more likely to obey
Poisson distribution
,ke-k
P(k)=-- -, (5.10)
k!
where X,=Np. The second moment of n is then
SN N 2 ke-k
= YkPk2) Yk (5.11)
k=0 k=0 k
and the variance of the particle concentration for a given cell is
-2 Nk2 ,'e-k -X 2
=- n = Yk ne (5.12)
k=O k!
ForN >>l,
N 2 Xke- k 2k e- oo ek -.
Y_ k-- k2 Y -_k- 1:[k(k- 1)+k]+--
k=O k! k=O k! k=O k!
00 Xke^ Ck Xk Ck
= Yk(k-1) --+ ke
k=O k! k=0 k!
X2 O k-2 0 xk-1
=Xe E +k e- Y
k=O (k 2)! k=O (k 1)!
X e~eX+Xe~eX
= X2 e-k ek.+ke-k ek
= X2 +X. (5.13)
Since X=Np= n, the variance of n can be expressed as
2 N k2 Cke- -2 2 -
= Yk2 n 2X X-n
k=O k!
-2 -2 -
=n + n-n = n. (5.14)
It follows for zero-size, inertialess particles that
< n'2 > I
-2I
r-4--2 (5.15)
n n
2 -2
For n >>1, seen from both equations (5.9) and (5.15), r. =< n'2 > /n -+0 so that
the variance of the particle concentration distribution is very small when it is normalized
by the square of the first moment of n. However, if N is fixed while N, takes a large
value in the limit, we have n << 1 and r. /n >> 1. In a practical simulation of particle
motions in turbulence, N is always finite and fixed. Since we are interested in
quantifying the small-scale variation of the particle concentration, it is necessary to
reduce successively the size of the cell, (Ax)3, on which n and are computed. It is
noted that the relatively large value of < n', > for inertialess particles is entirely due to the
fact that there is not enough particles in each cell as the size is reduced and r, -1/n
contains no contribution from the preferential concentration. As far as the inertial effect
is concerned, r-1/n is purely the inherent statistical noise which can be improved by
using large values of n in principle. As the inertia of the particles increases, < n'2> will
contain an increasing effect of the inertia bias. The fact that r -1/n >>1 for inertialess
particles when n << 1 implies, however, that r= < n'2 >/ n is dominated completely by
the inherent statistical noise r. -1/n. The contribution from the inertial effect to
or r thus must be isolated from the statistical noise. Figure 5.1 shows the comparison
between the DNS results r, for inertialess particles and the inherent noise 1/n, where
2 -2
/ n at every instant was recorded and a total of 300 time steps and 40 turbulence
realizations are used in our DNS. The details regarding turbulence representation, the
characteristics of the turbulence and particle collision detection scheme have been given
2 -2
in 2.1.3. It clearly shows that /n is indeed dominated by 1/n which has to be
subtracted.
5.2.2 Concentration Non-uniformity Parameter
For particles with finite inertia, the variance of the particle concentration can be
computed in direct numerical simulations
1 Nc
=-- (n-n)2. (5.16)
NC j=l
1.0E+5
1.0E+4
1.OE.3
1.0E+34
1.0E+2
_.E O- D=0.01, W =12,000 "
1.0E- a zero size, N=12,000
-1 zero size, N=500 -
1.0E-1 -- "-,i,
o N=500, 1/n
1.0E-2 N=12,000, 1 / n
1.0E-3
1.0E-4 i i :-: : : : :
0.1 1 10 100
Ax/ i
Figure 5.1 The normalized variance of the concentration as function of cell size for
inertialess particles in isotropic turbulence. Pf k =1.023.
The above variance of the concentration non-uniformity is dominated by the statistical
noise, as shown in figure 5.1. For zero-size, inertialess particles, the above should
approach to = n(1- p) for n >0(1), and =n for n 0(1). Defining the
normalized difference between above and as
n'2 n'2 () n'2
<-13 >= <-__->.<-*-o>
-2 -2 -2
n n n
n' (J3) 1 -p
> p for n > 0(1), (5.17)
2
n n
and
'2 n'2
<_!_ >= <_!->.<-p-4>
-2 -2 -2
n n n
n'2 ( ) I
=<--->- -
-2
n n
-1.20 F
-1.60
for n << 0(1).
1 10 100
Ax/ r
Figure 5.2 The expressions (5.17) & (5.18) of the concentration non-uniformity as
function of cell size for inertialess particles in isotropic turbulence. PT k =00.
The effect of inertia on the variance of n can be identified. Both equations (5.17)
and (5.18) represent the contribution to from the inertial effect for n > 0(1) and
n << 0(1), respectively. Figure 5.2 shows the results of equations (5.17) and (5.18) for
isotropic turbulent flow. Equation (5.17) obviously gives very small value for large Ax/D
for n >>0(1), while equation (5.18) on other hand gives small values for small Ax/D for
2 .-2
n << 0(1). For inertialess particles, / n from both expressions has 0(1) variation
2over cell size. Equation (5.17) gives
over cell size. Equation (5.17) gives /n --I- as cell size Ax approaches to zero
(5.18)
1- Q ----a
- 3-- =0.058, Eq (5.18)
- D=0.032, Eq (5.18)
-o-- D=0.022, Eq (5.18)
zero size, Eq (5.18)
--- D=0.058, Eq (5.17)
-0-D=0.032, Eq (5.17)
-- D=0.022, Eq (5.17)
- -zerosize, Eq (5.17)
-0.40
-0.80
,2 .-_2
while equation (5.18) gives /n -*1 as Ax increases. This is physically unsound
because the preferential concentration is an inertial phenomenon and should not appear in
turbulent flows with inertialess particles. Figure 5.3 shows the DNS results of
using equations (5.17) & (5.18) in a uniform shear flow for inertialess particles. The
uniform shear flow is represented by u=(Fy, 0, 0) with the velocity gradient F=1.0.
2 .-2
Comparing with figure 5.2, very similar behavior for < n' >/n is obtained even for a
uniform shear flow. The results are unphysical because there is no mechanism to drive
inertialess particles to form non-uniform concentration in a uniform shear flow. Thus,
one should not expect to apply the equation (5.17) or (5.18) to quantify the concentration
non-uniformity. However, the above results gave us a valuable hint. That is, the results
obtained by both equations (5.17) & (5.18) can be interpreted as an additional statistical
noise for both inertialess and inertial particles. It is hence suggested that the actual
contribution of the inertia to < n'2 > for an arbitrary 13 may be expressed as
= [-n(l- p)] [-n(- p) ] 0=, (5.19)
or
= [ n] [ n ] (5.20)
Apparently, zero concentration non-uniformity is obtained for inertialess particles in both
turbulence and uniform shear flows based on equations (5.19) & (5.20). It is noted that
is a function of 13 and n in equation (5.20), while is mainly determined
by 13 and n in (5.19). The influence ofp on in (5.19) may be negligible as far as
a small p is concerned, i.e. = f(13, n ), since p appears in both first and second
terms in right side of equation (5.19). A large p, however, means the large cell size (the
box size is the limit), which usually gives approximately uniform particle distribution. It
is also noted that --0 as n -+oo (or N-+oo). Thus, it is expected that the effect of
inertia on the concentration non-uniformity can be captured using equations (5.19) &
(5.20). To quantify the extent of the concentration non-uniformity using the second
moment of n, a dimensionless quantity AP called concentration non-uniformity
parameter, is defined as
<"n2> it, ) - (-P)
A --2 -- 1 -- (5.21)
n n n n
or
2 2
1 1
An [ -2 -1- J -I -2-- J=oo.
n n n n
0.00 I
-0.50 I
(5.22)
1 .5 0 I i l i- I I I I-* I I :4: : *:
0 0.2 0.4 0.6 0.8 1
Ax
Figure 5.3 The expressions (5.17) & (5.18) of the concentration non-uniformity as
function of cell size for inertialess particles in uniform shear flow. 13/=co.
Figure 5.4 shows the concentration non-uniformity parameter for the uniform shear
flow using equations (5.21) & (5.22). It is seen that AO is very small for all cases in a
uniform shear flow. Figure 5.5 shows AO in turbulent flows for a large range of sizes for
particles with zero inertia. It is also shown that, for zero-size particles, the concentration
non-uniformity at small inertia (Prk =7.3) still increases with the decreasing cell size.
Sn
--.-9. a___.._ .._____.._____ _____ -------__
---- D=0.04, Eq (5.18)
-a- D=0.02, Eq (5.18)
zero size, Eq (5.18)
- -D=0.04, Eq (5.17)
---- D=0.02, Eq (5.17)
-a- zero size, Eq (5.17)
The detailed discussions about the size effect on the collision rate will be given in 5.3.
The above results obtained by (5.21) & (5.22) are reasonable and hence strongly suggest
that equations (5.21) & (5.22) can be used to quantify and predict the particle
concentration non-uniformity induced by particle inertia in general turbulence. It should
also be noticed that both equations (5.21) and (5.22) always yield the same shape differed
by a constant (-=1). To facilitate the subsequent discussions, equation (5.22), which does
not contain p and is thus in a simpler form, is used hereinafter.
AP
004
0.03
0.02
0.00 -
Ax
Figure 5.4 The concentration non-uniformity parameter An as function of cell size in
uniform shear flow. F=1.0, /F=1.0.
Figure 5.6 shows the effect of total particle number N on the concentration non-
uniformity parameter AP. It is seen that An remains almost on the same curve when N
increases from 25,000 to 50,000 and 100,000 for all values of N, (from 23 to 2783).
o zero size, Eq (5.22)
o D=0.04, Eq (5.22)
o D=0.02, Eq (5.22)
S.- zero size, Eq (5.21)
&- ------ D=0.04, Eq (5.21)
^ l V- D=0.022, EEq (55.221)
eo0oo 0 o- -0 -. 10.0 r-TT. .-* .. .- .-