7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Dispersion of heavy particles in turbulentlike flows
Mathieu Ichard* and Jens Melheim
*University of Bergen, Institute of Physics and Technology
Postboks 7800, 5020 Bergen, Norway
mIlhLlihcu L I i:2 \con coin
Keywords: Isotropic Turbulence, Kinematic Simulation, Heavy particles motions
Abstract
The dispersion of heavy particles in turbulentlike flows is investigated. Special attention is paid to the different time scales
involved in the process of heavy particles motions. The homogeneous, stationary and isotropic turbulent flow field is simulated
by means of the kinematic approach. The statistical quantities of the turbulent flow field obtained with this approach agree
with theoretical expectations. The velocity autocorrelations of the particles and of the fluid seen by the particles along their
trajectories are investigated in details. Our results show that under the effect of inertia alone, the long time diffusivity
coefficient of particles with Stokes number less than 10 is larger than the diffusivity coefficient of fluid particles. Furthermore,
our predictions of the eddyparticle interaction time are in good agreement with the expression proposed by Oesterl &
Zaichik (2006). Finally, our kinematic simulation approach is used to simulate the heavy particle dispersion experiment of
Snyder & Lumley (1971) and the results compare well with the experimental observations.
Introduction
The dispersion of heavy particles is encountered in many
industrial or environmental applications: spray evaporation,
dust transport in the atmospheric boundary layer,
combustion in engines are a few examples. In this work a
heavy particle is defined as a small sphere with a density
two orders of magnitude larger than the fluid flow density.
Although, the motion of heavy particles; also denominated
inertial particles; has received much attention during the last
century, many open questions remained. Our objective is to
investigate the influence of the particle's inertia on the
different time scales involved in the dispersion process. For
example, for practical applications Stochastic Langevin
Models are used to predict the instantaneous velocity of the
particles and require as input the eddy particle interaction
time (Pozorski & Minier 1998).
Two main effects control the dispersion of heavy particles in
turbulent flows. The first effect is known as the inertia effect.
A particle with inertia has a response time different from
that of a fluid particle and as a consequence its trajectory
through a turbulent flow will be more or less affected by the
turbulent eddies. The second effect is due to the body forces
acting on the particle. For example, the gravity force gives
an inertial particle a drift velocity. This is the "crossing
trajectories effect" first identified by Yudine (1959) and
Csanady (1963). Most of the analytical studies on heavy
particle motions (seee Tchen 1947, Meek & Jones 1973,
Reeks 1977, Wang & Stock 1993) uses as baseline the
pioneered work of Taylor (1921) on Lagrangian dispersion.
Although his work was dedicated to the study of fluid
particle motion his ideas can be applied to the motion of
heavy particles. Taylor (1921) related the particle diffusivity
coefficient DP to the particle velocity autocorrelation and
the mean square velocity of the particle:
D (r) (v, (t)v, (t + r)) dr = v, (t))f p (r)dT
0 0
Moreover, as a direct consequence of the equation of motion
for a heavy particle, the velocity of the particle v,(t) depends
on the fluid velocity at the position of the particle at the
time t. Therefore one can expect a relation between the
particle velocity autocorrelation and the velocity
autocorrelation of the fluid seen by the particle along its
trajectory (see Wang & Stock 1993 for example). We then
conclude that the dispersion of a heavy particle depends on
the particle velocity autocorrelation, the velocity
autocorrelation of the fluid seen by the particle and the
mean square velocity of the particle.
We use a kinematic approach to simulate a homogeneous,
stationary and isotropic turbulent flow field in which heavy
particles are tracked. The kinematic method was first
proposed by Kraichnan (1970) and was improved by Fung
& al (1992). The turbulent flow field is modelled by the
superposition of a large number of randomly orientated
Fourier modes. First we present the kinematic approach that
we used to simulate the turbulent flow field. Then, we
verify that the turbulent statistics of the simulated flow field
agree with the theoretical predictions. Results of heavy
particles dispersion are presented and discussed with a focus
on the particle time scale and the eddy particle interaction
time. Finally, the paper ends with the simulation of the
Snyder and Lumley (1971) dispersion experiment.
Nomenclature
g gravitational constant (ms2)
D diffusivity coefficient (m2s 1)
v particle velocity (ms1)
u Eulerian fluid velocity (ms1)
Nk number of Fourier modes
an/bn Fourier decomposition coefficient (ms1)
kn wave number (m'1)
E Eulerain energy spectrum (m3s2)
k, cutoff wave number (m'1)
ln number of modes in the largescale eddies range
C dimensionless constant in the expression for E
ums root mean square fluid velocity (ms1)
L length scale of the flow (m)
A integral length scale (m)
Us sweeping velocity (ms1)
T time scale (s)
dp particle diameter (m)
St Stokes number for the particles
vd drift velocity (ms1)
x ratio of the rms particle velocity and u.s
Greek letters
uL dynamic viscosity (kgm s ')
T aerodynamic particle response time (s)
p velocity autocorrelation ()
co frequency (s ')
p dimensionless coeff. in the expression for u,
S frequency spectrum (m2s )
particle density (kgm 3)
Subscripts
E Eulerian
L Lagrangian
P particle
F fluid seen by the particles
Kinematic Simulation
Turbulent flow fields can be simulated using two
approaches. The first approach is the dynamical method
where the turbulent flow field is explicitly computed from
the NavierStokes equations, i.e. the dynamical equations.
Direct Numerical Simulations (DNS) or Large Eddy
Simulations (LES) are examples of the dynamical approach.
The second approach is the kinematic method which lies on
the superposition of a large number of randomly orientated
Fourier modes. The turbulent flow field is modelled in a
qualitative sense and at a rather low computational cost.
Simulations based on the kinematic method, referred to as
Kinematic Simulations (KS) in the following, use an
analytical expression for the incompressible velocity field
u(x,t) and thus are not gridbased and do not require
interpolations of the velocity field. KS have been mainly
used to compute Lagrangian statistics and turbulent
diffusion properties of fluid elements and the results have
been shown to agree well with experimental and DNS data
(see Nicolleau & Yu 2004, Malik & Vassilicos 1999). Fung
& al (2003) used a KS to study the dependence of the heavy
particles diffusivities on the structure of the flow field.
In this section we present the KS approach used to simulate
a homogeneous, isotropic and stationary turbulent flow field.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The Eulerian instantaneous velocity field at a position x and
at a time t is calculated as the sum of a finite number of
Fourier modes:
Nu
u(x, t) = a. cos(k. x+,t) + b. sin(kx + at)
n1
where Nk is the total number of Fourier modes and the
vectors a. and b, are the decomposition coefficients for
the wavevector k, having the unsteadiness frequency
c, The vectors a, and b, are carefully chosen to
ensure that the orientations of the Fourier modes are random
and uncorrelated. Following Oesterl6 & Zaichik (2006), the
components of these vectors are picked from a zero mean
Gaussian distribution with a variance equals to the velocity
variance. Moreover the velocity field must satisfy the
incompressibility condition and consequently the
orientations of a, and b, are sampled independently and
randomly in a plane normal to k :
a, k = b, k. = 0
The magnitudes of the vectors a, and b, control the
amplitude of the velocity's fluctuations and are given by:
la.J = b.2 =2E(k,)dk,
 n =1
dk= lk 1k t ne[2,Nk
2
k k
n= ANk
2
Before discussing the remaining parameters of the KS;
namely the wavevector, the Eulerian energy spectrum and
the unsteadiness frequency, we introduce the sweeping
phenomenon. Recent KS techniques can incorporate the
sweeping effect which is the advection of the small
scaleeddies by the largescale eddies (Tennekes 1975). In
the KS framework, an advanced model for the sweeping
effect has been proposed by Fung & al. (1992). A simplified
and less time consuming version of this model is used in
this work. A characteristic sweeping velocity, u is
defined and it is assumed that the smallscale eddies are
advected by this characteristic velocity.
The unit wavevectors k. = k /k, are randomly selected
with an isotropic distribution. The magnitudes of the
wavevectors are uniformly distributed from k, to k,:
n = k( n ) n< n1
n +1
and they obey a geometric distribution from k, to k, :
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
27(
k (;T
(7) A
n >n
The wave number k. is a fixed cutoff wave number
delimiting the smallscale and largescale eddies ranges, k,
defines the size of the biggest eddies in the turbulent flow
field and n. is the number of Fourier modes in the
largescale eddies range. As suggested by Fung & al (1992),
we use 6 modes for the largescale eddies. We will first
present the expression of the Eulerian energy spectrum
before discussing the parameter a in Equation (6) because
the existence of this parameter is closely related to our
choice for the Eulerian energy spectrum. We adopt the
expression used by Osborne & Vassilicos (2005) for the
Eulerian energy spectrum in the inertial range:
E(k,) = CuL(Lk,)
and set the constant of proportionality to 1/5. We use the
same approach in this work. The integral length scale is
computed by:
SE(k)k ldk
A=
4 c
i E(k)dk
0
Therefore, the value of the parameter a is adjusted so that
Equation (12) is satisfied. The definition of the cutoff wave
number agrees with the prediction of Tennekes & Lumley
(1972):
(8) ko = 1.8 =2;
3.5
where L=2;f/ko and the spectrum is conformed to the
Kolmogorov's "5/3 law" in the inertial range from k, to
kN This is a general form of the Eulerian energy spectrum.
The rate of dissipation of kinetic energy does not appear
explicitly in Equation (8) which is coherent with the fact
that in KS there is no turbulent dissipation as the turbulence
dynamics are absent. The cutoff wave number associated
with the sweeping phenomenon coincides with the lower
bound of the inertial range as in Fung & al (1992). The
energy contained in the large scale eddies is given by:
5
E(k=) Cut3L(Lk ) (9)
To sumup, the expression of the Eulerian energy spectrum
is:
0 k < k
E(k, CuL(Lk ) kkI kn k,
E(k,)=
Cus L(Lk,) I ke < kn k^
0 k, > kN
The dimensionless constant C is given by:
o 2
fE(k)dk =u
0 2
The main inputs for the KS are the largest wave number,
kN the cutoff wavenumber, k., and the fluctuating rms
velocity, u which sets the level of turbulence intensity
in the flow field. An equivalent Reynolds number can be
defined as Re ~ (kk /k) .The determination of the
cutoff wavenumber is not a trivial task. Fung & al (1992)
related this parameter to the integral length scale of the
flow:
where ko is the lower cutoff wave number for the inertial
range in Tennekes & Lumley (1972)'s analysis and f is a
characteristic length scale of the flow which can be seen as
the integral length scale. The shape of the Eulerian energy
spectrum prescribed in our KS is shown on Figure 1 and the
values of the wavenumbers are given in Table 1.
Based on dimensional analysis, the characteristic sweeping
velocity is expressed by:
IU = E(kI)k
The dimensionless parameter / is in the range [0:1]. It
means that only a fraction of the energy contained in the
cutoff wave number is used to advect the smallscale
eddies.
Finally, the unsteadiness frequency t,, which controls the
time dependence of the flow field is sampled from a
Gaussian distribution with a zero mean and a given standard
deviation as in Oesterlk & Zaichik (2006).
Eulerian Energy Spectrum
102
103
10
105
106
10' 10
k (m1)
Figure 1: Shape of the Eulerian energy spectrum used in
our KS approach. The values of the wavenumbers are
given in Table 1.
n(n, +l1)
k, = k, N
Turbulent Flow Field Properties
The KS approach described in the previous section has been
used to simulate a homogeneous, isotropic and stationary
turbulent flow field whose properties are given in Table 1.
Different tests were performed to verify that the flow field
was stationary and isotropic. We have found that:
(u(x,t)) 0 i= 1,2,3
at several positions x. The averages are performed over
time and over all flow field realizations.
As explained in the introduction, our main focus is on the
Eulerian and Lagrangian velocity autocorrelation functions.
Lagrangian statistics are computed by tracking 2700 fluid
particles in the flow field and the equation of motion is
solved via a 4th order RungeKutta method. The Eulerian
velocity autocorrelation function writes:
(u, (x, t)u, (x, t + r)
(u (x,t))
and for the Lagrangian velocity autocorrelation we have
(v, (t)v (t + )
PL (v) (,= ) v(t) = u, [x(t),t] (19)
where no summation over the indices is implied and the
average is performed over time and over all flow
realizations but also over all fluid particles for the
Lagrangian velocity autocorrelation. As the turbulent flow
field is stationary the autocorrelation functions depend only
on the difference in time, r The autocorrelations are
shown on Figure 2.
The computations have been performed with and without
sweeping velocity. The first thing to note is that the Eulerian
velocity autocorrelation is above the Lagrangian velocity
autocorrelation. It implies that the Eulerian time scale TE
is greater than the Lagrangian time scale TL. The time
scales are defined by:
TE = pE ()dr
0 (20)
T, = pL (r)dz
0
Several authors have also observed ratios TL/TE <1 (see
Fung & al 1992, Yeung 2001). However, the value of the
ratio seems to be flow dependent as a significant scatter
exists. The flow dependency is also observed in our
simulation when we vary the characteristic sweeping
velocity. The sweeping process has a strong effect on the
Eulerian velocity autocorrelation and almost no effect on the
(x, t) = 2 (X, t) = 2 (x, t)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Lagrangian and Eulerian velocity autocorrelations
1.0 60 Euler Us=0
p ...***** Euler
o0.8 A A Lagrange Us=O
SLagrange
0.6
50.4
"0.2
0.0
0.0 0.5 1.0 1.5 2.0
Tune Lag (s)
Figure 2: Eulerian and Lagrangian velocity autocorrelations
for the turbulent flow field simulated with the parameters
given in Table 1. The computations have been performed
with and without a sweeping velocity.
Lagrangian velocity autocorrelation as it can be seen on
Figure 2. The Eulerian autocorrelation decreases faster
when the sweeping process is activated. As the sweeping
velocity increases the ratio TL/TE tends toward 1. This
effect can be understood intuitively; as the strength of the
sweeping process increases more and more smallscale
eddies are advected past a fixed observer over a
characteristic time t, and as a consequence, the fixed
observer sees a velocity field which is very similar to the
velocity field seen by a moving observer along a fluid
particle trajectory. With a zero sweeping velocity we obtain
TL/TE=0.347 and with a sweeping velocity u,=u~/5
(fl=4%) we obtain TL/TE=0.714which is close to the
values Yeung (2001) has observed in his DNS simulations.
The value of the sweeping velocity u,= u, /5 is used in
the following.
The sweeping effect should also be observed on the
Eulerian and Lagrangian frequency spectrum:
E1 f E ( )er d
o (21)
1L JL f ()e 'dr
0
Tennekes (1975) has shown that the Eulerian frequency
spectrum should scale as 1 5/3 in the inertial range
due to the sweeping process. This scaling is different from
the scaling derived from the Kolmogorov hypothesis, i.e.
01 2 Figure 3 (with and without a sweeping velocity)
show the Eulerian frequency spectrum computed from the
Eulerian velocity autocorrelation.
kNk k, un Nk N p dt
(m) (m) (ms) N NR Npa (s)
1000 5 1.0 100 100 2700 0.0005
Table 1: Input parameters for the KS. NR is the number of
flow realizations and Np, is the total number of particles
tracked to compute the Lagrangian statistics.
oEu erian Fruetcy SpectunmwiL S epn Eulrian FequenySpectum wlhoul Sweepi
S..... lng w 10 .... cling w
10 10
1'in 1 1 00 10t 1Co0 i? IDo io' 10
Prequcnry (c) Fnqucncy l)
1 gl.. IF eque Sped .51r wii"h Sw I Jqin Frequency Sp eeunwil~o Swepine
"10 alngw' i. ngw'
10 2 10
I'O1 W1' 10 I ice '1013 101 i02 rln 1o0
Frequency w ) Ftquncy(w)
Figure 3: Eulerian and Lagrangian frequency spectrum with
and without sweeping velocity. The KS computations are
based on the parameters given in Table 1.
The 2 and 5/3 scaling are shown. The Eulerian
spectrum scales with c 5/3 when the sweeping process is
activated and it scales with t 2 when this process is not
activated, in agreement with the theory. As far as the
Lagrangian frequency spectrum is concerned no effects of
the sweeping process are expected. This is observed in our
computations as it can be noted on Figure 3. However,
following the Kolmogorov hypothesis in the inertial range,
the Lagrangian frequency spectrum should scale as a 2
whereas our computations exhibit a 5/3 scaling. Osborne &
Vassilicos (2005) have also observed a 5/3 scaling for the
Lagrangian frequency spectrum and they have shown that
the scaling was dependent on the expression used to define
the unsteadiness frequency. Further work is needed to fully
understand the reasons for the deviations between our KS
predictions and the theory.
Heavy Particle Dispersion
We now investigate several properties of heavy particles
dispersion in isotropic turbulence. The particles have a
diameter smaller than the Kolmogorov scale and much
larger than the fluid molecules. The particles density is two
orders of magnitude greater than the fluid density and the
loading is assumed to be dilute enough so that particles
collisions and turbulence modulation due to the presence of
the particles can be neglected. The motion of a particle in a
flow field is predicted via the BassetBoussinesqOseen
(BBO) equation. This equation is very complex but for
heavy particles it can be simplified to the following form
(Maxey & Riley 1983):
dv(t) u[x(t),t] v(t)
+g
dt (22)
dx(t)
d v(t)
dt
where a Stoke's drag law has been assumed to hold. This
assumption is only valid if the particle Reynolds number is
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
less than one. The aerodynamic response time of the particle
is given by:
r18d
a 18~u
The aerodynamic response time represents the time needed
by the particle to respond to changes in the fluid flow
velocity. The Stokes number is used to relate the
aerodynamic response time of the particle to a characteristic
timescale of the flow field. A small Stokes number, i.e.
St<
turbulent eddies and inversely a large Stokes number, i.e.
St>>l, means that the particle cannot follow the turbulent
eddies. In this work the Stokes number is based on the
integral scales of the flow, St= Tr u, /A. The drift velocity
of a particle is defined as:
Vd = ag
The main difficulty in solving Equation (22) resides in the
estimation of the fluid velocity u[x(t),t] along a particle
trajectory. The fluid velocity seen by a heavy particle
dispersing in a turbulent flow field is different from the
velocity seen by a fluid particle dispersing in the same flow
field. This is a consequence of two effects. The first effect is
the inertia effect; the aerodynamic response time of a heavy
particle to a turbulent fluctuation is different from that of a
fluid particle. The second effect is associated to the body
forces which give a finite drift velocity to the heavy particle.
This latter effect is known as "the crossingtrajectories
effect" first identified by Yudine (1959).
First, in order to investigate the inertia effect alone, particles
with different Stokes numbers and without drift velocity are
tracked in an isotropic flow field simulated by means of the
KS approach described in the previous section and the
equation of motion for a particle, Equation (22), is solved
via a 4th order RungeKutta method. The parameters of the
flow field are given in Table 1. Figure 4 shows the
Lagrangian velocity autocorrelations of the particles which
have been computed by:
() v,(t)v, (t + )
(2 (t))
Particle velocity autocorrelations
,1.0! t St=0.1
..... St=0.2
.00.8 ". St=o.5
0.6 .. " St5
.. St=10
.0.4 ~ St=20
s 1
0.2
0.0
0.0 0.5 1.0 1.5 2
Time Lag (s)
Figure 4: Particle velocity autocorrelations.
The Lagrangian timescale of the particle, Tp, which is
defined as:
TP = pP(r)dr
0
is increasing with increasing Stokes number which indicates
that the memory of the particle to its previous velocity is
increased.
However, this observation is not sufficient to conclude that
the dispersion of a particle is increasing with inertia because
as discussed in the introduction the diffusivity of a particle
is a function of both the particle velocity autocorrelation and
the mean square velocity of the particle; the latter
decreasing with increasing particle's inertia. Two interesting
questions arise: does an inertial particle disperse more than a
fluid particle and how does the diffusivity of an inertial
particle vary with the Stokes number ? We have investigated
these two questions by considering the long time,
asymptotic form, diffusivity of the particles defined as:
DP ( v T, (27)
which is only valid if t > T,. The long time diffusivity of
the inertial particles is compared to the long time diffusivity
of the fluid particles and the results are shown on Figure 5.
Our simulations show that inertial particles with a Stokes
number less than 10 disperse more than fluid particles and
that the ratio DP/DF attains its maximum for St around
1.0. This a quite interesting result which could have been
expected because a Stokes number of 1.0 means that the
aerodynamic response time of the particle is the same than
the characteristic time of the flow field. Furthermore, the
maximum value of the diffusivities ratio is relatively high
indicating that at a Stokes number of around 1.0, inertial
particles can disperse as much as 25% more than fluid
particles. This maximum is within an acceptable range of
the maximum value obtained by Fung & al (2003). However,
Fung & al (2003) have also predicted that inertial particles
with a Stokes number less than around 0.5 would disperse
less than fluid particles. This result is in clear disagreement
with our observations and the DNS results of Squires &
Eaton (1991).
Figure 5: Variations with the Stokes number of the ratio of
the diffusivity coefficient of inertial particles and the
diffusivity coefficient of fluid particles.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
It is also worth looking at the Lagrangian velocity
autocorrelation of the fluid seen by the particle along its
trajectory. We define this correlation as:
F() (u, [x(t),t]u, [x(t+ ),t+ }])
P ()=,),])
T, = J(r)dT
0
with TF being the Lagrangian timescale of the fluid seen
by the particle. This timescale, also called eddyparticle
interaction time, is one of the most important parameter to
predict the dispersion of heavy particles in a turbulent flow
field. Indeed, as its name indicates, it controls the time a
particle will interact with a turbulent eddy. This is a key
input parameter for stochastic models that are extensively
used to model particles dispersion in turbulent flow fields.
Figure 6 shows the Lagrangian velocity autocorrelation of
the fluid seen by a particle with different Stokes number.
The Eulerian velocity autocorrelation and the Lagrangian
velocity autocorrelation for fluid particles are also shown.
As observed by several authors, the velocity autocorrelation
of the fluid seen by the particle varies between two limiting
autocorrelations, namely the Eulerian velocity
autocorrelation and the Lagrangian velocity autocorrelation
for fluid particles. As St 0 we have TF , T, which
means that particles with almost zero inertia behave as fluid
particles and as St > o we have T, > T meaning that
particles with such high inertia are not able to move in the
flow field.
Different expressions have been proposed to predict the
variation of TF as a function of the Stokes number.
Pozorski & Minier (1998) proposed the following relation
between TF, TE and TL:
1 x 1x
F TL E
where x is the ratio of the particle's rms velocity to the
fluid's rms velocity seen by the particle.
Autocorrelation p'(T)
Fluid Particle
St=.1
 St=0.2
St=0.5
St=2
St=5
.. St=10
St=20
Euler
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
Time Lag (s)
Figure 6: Velocity autocorrelations of the fluid seen by a
particle along its trajectory. The Eulerian and fluid particle
velocity autocorrelations are also shown.
Oesterl6 & Zaichik (2006) developed a semiempirical
model to predict the time scale T,. The model is not only
valid for heavy particles but also for light particles, i.e. the
model should work for arbitrarydensity particles. The
expression is:
T, T, +(T,,)St 9St2 1 (30)
TF T T L 1+St (l+St)2(2 +St)
Wang & Stock (1993) have also suggested an expression
relating T, to TE and TL by curve fitting the results of
their KS:
i t 1 o 4(1+0 Is T,/T
T = TE (1+ St)o 41+0 OSt)
To conclude the review of the different expressions for the
eddy particle interaction time, we mention the recent work
ofYeo & Lee (2008) based on DNS simulations at Re =47.
They have shown that the time scale of the fluid velocity
seen by a particle does not vary monotonically between TE
and TL but it varies as an oscillatory Nshape. They found
that the function T =f(St) had a maximum at StI1 and they
proposed the following expression:
T 2/1 + +(TE T )(0.025St)l51
S= L0.245e [1s12)13] 1 + (32)
T, 1+(0.025St)15
Figure 7 shows a comparison between Equations (29), (30),
(31), (32) and the results obtained from our simulations. A
modified version of Equation (29) has also been added. In
this modified version the factor x is squared. The two
asymptotic conditions are met by all the expressions and our
results. However, when the Stokes number is in the range
0.1100, 3 different tendencies are observed. Equation (30),
the modified Equation (29) and our results follow a similar
trend whereas Equation (31) and the original form of
Equation (29) follow another trend. As one could have
expected the expression proposed by Yeo & Lee (2008) has
its own Nshape trend.
1.6 Variation of TF with the Stokes number
* Simulations
1.5 Eq (29)
SEq (29)"bis"
1.4 I Eq (30)
..... Eq(31) .
,3 Eq(32)
1.2
1.1
2 101 10 101 102
St
Figure 7: Comparison of our predictions for TF with
different analytical expressions.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The second effect known as the "crossingtrajectories
effect" is due to the gravity force acting on the particle. This
effect was first identified by Yudine (1959). Csanady (1963)
has also studied this phenomenon and he derived
expressions for the particles diffusivity in the directions
parallel and perpendicular to the drift velocity. These
expressions show that the dispersion of the particles
decreases for increasing drift velocity both in the direction
perpendicular as well as parallel to this velocity. Moreover,
in the limit of very large drift velocity, the dispersion in the
perpendicular direction to Vd is onehalf of the dispersion
in the direction parallel to vd The "crossing trajectories
effect" is investigated in the next section where the
experiment of Snyder & Lumley (1971) has been simulated.
Comparison with experiment
Snyder and Lumley (1971) measured particle dispersion and
particle velocities in gridgenerated turbulence. These
experimental data are widely recognized as the most
accurate data for particles dispersion in a turbulent flow and
provide a good test case against which our KS approach can
be assessed. Gridgenerated turbulence is similar to
homogeneous decaying turbulence (Sawford & Guest 1991).
In the experiment turbulence was generated by a grid with a
mesh dimension of M=2.54cm which gave a Reynolds
number of 104. The flow was orientated vertically upward
and the mean velocity was 6.55m/s. The test section of the
wind tunnel where the properties of the particles were
measured had dimensions 0.41mx0.41mx4.88m The
particles were released on the centerline of the wind tunnel
at x/M=20 (the grid being located at x=0m ) with a
velocity equals to the mean velocity of the flow field. Four
different types of particles with different aerodynamic
response times were tested. The properties of the particles
are given in Table 2.
Particles data were obtained using a photographic system
and the measurements were performed at ten stations spaced
logarithmically between x/M=68.4 and x/M=168. The
lateral particle dispersion data are presented as a function of
the distance downstream and relative to their displacement
at the station 1. However, as the turbulence is decaying,
flow properties such as the longitudinal integral length scale
and longitudinal turbulence intensity vary along the
downstream distance. This variation cannot be accounted
for in our KS approach and we have to choose characteristic
values for the turbulent properties.
Particles Characteristics
Hollow Glass
Copper (C) Corn Pollen (CP) (HG)
Vt (ms') Ta (s) Vt (ms') Ta (s) Vt (ms') Ta (s)
0.00
0.483 0.049 0.198 0.02 0.0167
Kinematic Simulation Parameters
k k (m1) u s (ms1) Nk NR Npart
14612 84 0.130.11 100 100 2700
Table 2: Parameters for the Snyder and Lumley experiment:
vt is the terminal velocity of the particle and the timestep is
taken as the minimum between Ta and 0.5/(kNkrms).
We follow Fung (1990) and we take the values of the
turbulent properties at the first station of measurement.
Snyder & Lumley (1971) measured the longitudinal integral
length scale which is roughly half the lateral integral length
scale. Hence, in our simulations the integral length scale is
A=0.015m and consequently the cutoff wave number is
k,=84m Snyder & Lumley (1971) reported a
Kolmogorov length scale of 0,00043m at the first
measurement station which gives kN =14612m 1. The last
input is the rms velocity, the value for this parameter is
discussed further down. The initial velocity of the particles
is set to the mean velocity of the flow field. Statistics are
computed over an ensemble of 2700 particle trajectories,
each particle dispersing in a different realization of the flow
field.
Figure 8 shows the particle dispersion curves. In both the
experiment and the simulation, the smaller the aerodynamic
response time of the particle the greater the dispersion. This
observation can also be formulated in terms of the "crossing
trajectories effect": the larger the drift velocity the lower the
dispersion. Results for the copper and corn pollen particles
show nearly perfect agreement with the experimental data.
This is because we have adjusted the value for the rms
velocity for the two types of particle. We have first run a test
with u, =0.13m/s, and we found perfect agreement for the
copper particles. For the corn pollen particles the simulated
dispersion curve had the same shape than the experimental
one but it was shifted 20% upward. We then tried a value of
0.11 m/s (still above the value measured at the second
measurement station) and we obtained perfect agreement for
the corn pollen particles. However, with this new rms
velocity value the curve for the copper particle was shifted
20% downward. Therefore, results shown on Figure 8 are
obtained with u,=0.13m/s for the copper particles and
u_=0.11m/s for the corn pollen particles. A large
discrepancy is observed for the hollow glass particles
(whatever rms velocity we take). The hollow glass particles
disperse way more in the simulation compared to the
experimental observations. We have also added the
dispersion curve for a fluid particle and it can be noted that
the simulation predicts that the hollow glass particles
disperse as the fluid particles. This seems reasonable due to
the very small aerodynamic response time of the hollow
glass particle. This observation has also been made by
Sawford & Guest (1991). The reason for the
underestimation of the hollow glass particles dispersion in
the experiment could be that "as much as 40% of the energy
of the hollow glass particles could have been lost due to the
low sampling rate", as reported by Snyder & Lumley (1971)
and mentioned by Reynolds & Iacono (21" '4).
The velocity autocorrelations of the particles have also been
measured in the experiment. Figure 9 depicts a comparison
of the measured autocorrelations with the predictions of the
KS approach. The experimental data show that the
autocorrelations of the heavy particles decrease faster than
the lighter particles. This feature is nicely reproduced in our
simulations and it is due to the "crossingtrajectories effect".
The predicted correlation for the hollow glass particles is
above the measured one which is consistent with the
observation for the dispersion curve. The predicted and
measured heavy particles correlations are very close to each
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
other and within the experimental error range expect for the
negative loops. Our results show a larger magnitude for the
negative loops. The presence of negative loops for such a
particle aerodynamic response time and inertia is also
predicted by the model of Wang & Stock (1993).
Disersion curves in Snyder & Lumley experiment
Sinu C
.... sim uCP
5 SimuHG
ExpC
4 Exp 'G
u uEdPpHG rtide
S6* in FhiduPatide
U0 0.1 0.2 0.3 0.4
Tune (s)
Figure 8: Comparison of the simulated
dispersion curves
p"(r) in Snyder & Lumley exp.
0.5
and observed
0.05 0.10 0.15 0.20 0.25
S(s)
Figure 9: Comparison of the simulated
particles velocity autocorrelations.
and observed
Conclusions and Future Work
The behaviour of heavy particles dispersing in a
homogeneous, isotropic turbulent flow field has been
investigated. A kinematic approach has been presented and
used to simulate the turbulent flow. The shape of the
Eulerian velocity autocorrelation is dependent on the
magnitude of the sweeping of the smallscale eddies by the
largescale eddies. The Eulerian and Lagrangian frequency
spectrum have been computed. The Eulerain frequency
spectrum scales as predicted by the theory of Tennekes
(1975) whereas the Lagrangian frequency spectrum does not
exhibit the ) 2 Kolmogorov scaling as we would have
expected.
The results of our simulations show that heavy particles
with a Stokes number less than 10 disperse more than fluid
particles. The maximum dispersion is obtained for particles
whose Stokes number is around 1. As far as the eddy
particle interaction time is concerned, our results agree well
with the expression proposed by Oesterl6 & Zaichik (2006).
Finally, the whole approach used in this paper is able to
reproduce with good accuracy the experimental
observations of Snyder & Lumley (1971). This kinematic
approach can therefore be used with some confidence to test
stochastic models, assess their performances and understand
their weaknesses. This would be the next step of this
ongoing work.
Acknowledgements
This work has been sponsored by a "Naering PhD" grant of
the Norwegian Research Council.
References
Csanady, GT Turbulent Diffusion of heavy particles in the
atmosphere. J. Atmos. Sci, Vol. 20, 201208 (1963)
Fung, J.C.H. Kinematic simulation of turbulent flow and
particle motions. PhD Thesis, University of Cambridge, UK
(1991)
Fung, J.C.H., Hunt, J.C.R., Malik, N.A. & Perkins, R.J.
Kinematic simulation of homogeneous turbulent flows
generated by unsteady random Fourier modes. J. Fluid
Mech., Vol. 236, 281317 (1992)
Fung, J.C.H., Hunt, J.C.R. & Perkins, R.J. Diffusivities and
velocity spectra of small inertial particles in turbulentlike
flows. Proc. R. Soc. Lond. A, Vol. 459, 445493 (2003)
Kraichnan, R.H. Diffusion by a random velocity field. Phys.
Fluids, Vol. 13, 2231 (1970)
Malik, N.A. & Vassilicos, J.C. A Lagrangian model for
turbulent dispersion with turbulentlike flow structure:
comparison with direct numerical simulation for
twoparticle statistics. Phys. Fluids, Vol. 11, 1572 (1999)
Maxey, M.R. & Riley, J.J. The equation of motion for small
rigid particle sphere in a nonuniform flow. Phys. Fluids, Vol.
26, 883889 (1983)
Meek, C.C. & Jones, B.G Studies of the behaviour of heavy
particles in a turbulent fluid flow. J. Atmos. Sci., Vol. 30,
239244 (1973)
Nicolleau, N. & Yu, G Twoparticle diffusion and the
locality assumption. Phys. Fluids, Vol. 16, 2309 (2' l4)
Oesterl6, B. & Zaichik, L.I. Time scales for predicting
dispersion of arbitrarydensity particles in isotropic
turbulence. Int. J. Multiphase Flow, Vol. 32, 838849 (2006)
Osborne, D.R. & Vassilicos, J.C. Oneparticle twotime
diffusion in threedimensional homogeneous isotropic
turbulence. Phys. Fluids, Vol. 17, 5104 (2005)
Pozorski, J. & Minier, J.P On the Lagrangian turbulent
dispersion models based on the Langevin equation. Int. J.
Multiphase Flow, Vol. 24, 913945 (1998)
Reeks, M.W. On the dispersion of small particles suspended
in an isotropic turbulent fluid. J. Fluid. Mech., Vol. 83,
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
529546 (1977)
Reynolds, A.M. & Iacono, GL. On the simulation of
particle trajectories in turbulent flows. Phys. Fluids, Vol. 16,
No. 12, 4353 (2. 114)
Sawford, B.L. & Guest, F.M. Lagrangian statistical
simulation of the turbulent motion of heavy particles.
BoundaryLayer Meteo., Vol. 54, 147166 (1991)
Squires, K.D. & Eaton, J.K. Measurements of particle
dispersion obtained from direct numerical simulations of
isotropic turbulence, J. Fluid Mech., Vol. 226, 135 (1991)
Snyder, W.H. & Lumley, J.L. Some measurements of
particle velocity autocorrelation functions in a turbulent
flow. J. Fluid Mech., Vol. 48, 4771 (1971)
Taylor, GI. Diffusion by continuous movements. Proc. Lond.
Math. Soc., Vol. 20, 196211 (1921)
Tchen, C.M. Mean value and correlation problems
connected with the motion of small particles suspended in a
turbulent fluid. PhD Thesis, University of Delft, The Hague,
Netherlands.
Tennekes, H. & Lumley, J.L. A first course in turbulence.
The MIT Press, ISBN 0262200198 (1972)
Tennekes, H. Eulerian and Lagrangian time microscales in
isotropic turbulence. J. Fluid Mech., Vol. 67, 561567
(1975)
Wang, L.P & Stock, D.E. Dispersion of heavy particles by
turbulent motion. J. Atm. Sci., Vol. 50, No. 13, 18971913
(1993)
Yeo, J.J.K. & Lee, C. Behavior of heavy particles in
isotropic turbulence. Phys. Review E, Vol. 77 (2008)
Yeung, PK. Lagrangian characteristics of turbulence and
scalar transport in direct numerical simulations. J. Fluid
Mech., Vol. 427, 241274 (2001)
Yudine, M.I. Physical considerations on heavyparticle
dispersion. Advances in Geophysics, Vol. 6, 185191 (1959)
