7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
3D Droplet Dispersion in a Cylinder Wake Flow
Albert Lee
Department of Aerospace Engineering, University of Illinois UrbanaChampaign
104 South Wright Street, Urbana, IL 61801, USA
Eric Loth
Department of Mechanical & Aerospace Engineering, University of Virginia
122 Engineer's Way, Charlottesville, VA 22904, USA
aylee2@illinois.edu and loth@virginia.edu
Keywords: Cylinder wake, Hybrid RANS/LES, Particle dispersion
Abstract
A hybrid RANS/LES scheme with a continuous random walk (CRW) model was used to investigate droplet dispersion in the
wake of a subcritical cylinder at a Reynolds number of 8000. The fluidphase flow field was solved assuming oneway
coupling using the NicholsNelson hybrid scheme with the Menter SST turbulence model and exhibited the Karman vortex
street, but with some instabilities in the wake. Lagrangian droplet trajectories were calculated using both the instantaneous
resolved flow field and a stochastic subgrid fluctuations based on a CRW model. The effect of droplet inertia was studied
using various drop sizes which were injected downstream of the cylinder in the absence of gravity. Substantial differences
were found as a function of droplet Stokes number. Small Stokes number droplets were affected by both the largescale and
smallscale vortex structures and closely followed the fluid flow. Droplets with intermediate Stokes numbers experienced a
centrifugal effect from the vortices causing them to disperse over a greater area and cluster in the low vorticity regions on the
outer edges of the wake. Large Stokes number droplets were found to be less responsive to the vortex street and dispersed
the least. Vorticitybased concentration distributions illustrated the preferential bias of the different particle sizes similar to
results seen in homogenous isotropic turbulence, but with weaker inertial trends.
numerically simulated particle flow in low Reynolds
number cylinder wake flow due to their applications in
environmental and combustion process. Other studies
have used the discrete vortex method for particles in high
Reynolds number cylinder flow. One may expect different
particle dispersions in such cases since the flow field is
dependent on Reynolds number.
At very low Reynolds numbers, there is a steady
laminar wake. The cylinder wake begins to develop a pure
Karman vortex street when 90 < Re < 300. In this region,
direct numerical simulation (DNS) have been used to study
particleladen flows. Wu et al. (2009) examined
trajectories of particles with various Stokes and Froude
numbers in the wake of a circular cylinder at Re = 100 It
was found that particle motion was dependent on only
Stokes numbers when the Froude number was large, i.e.
gravity effects were small. Particles with small Stokes
numbers filled the vortex core but particles at intermediate
Stokes numbers concentrated at the edge of the vortex street.
When the particle Froude number was also small, particles
concentrated in the lower regions of the vortex street due to
the significant gravity force. Luo's DNS study of
particles around a circular cylinder to investigate the
dispersion of coal particles in heat exchangers was
Introduction
The transport of particles in bluff body wakes is a
phenomenon common in nature and engineering
applications such as the dispersion of aerosol particles
around buildings or the dispersion of fuel particles in
combustion processes. While these two examples involve
particles flowing around the bluff bodies, some applications
involve injecting particles into the wake. In groundbased
icing research tunnels to test aircraft icing, water droplets
are sprayed into the wind tunnel upstream of the test section
just downstream of a spray bar, which has a substantial
wake. Uniform liquid water content is desired at the test
section to simulate cloud conditions. The NASA Glenn
Icing Research Tunnel and the Arnold Engineering and
Development Center (AEDC) Engine Test Facility J2 test
cell both rely on the wakes of the spraybars to promote the
dispersion of particles (DeAngelis et al. 1997). To study
the fundamental dispersion of droplets in buffbody wake a
circular cylinder was chosen since it is well characterized in
terms of wake dynamics and velocity profiles.
Cylinder flow can be categorized based on the
Reynolds number of the flow (based on freestream velocity
and cylinder diameter). The majority of studies have
conducted at a range Re =140 260 (Luo et al. 2009).
This threedimensional study found that both the spanwise
vortex structures and the streamwise vortex tube structures
effect the dispersion of particles with intermediate and large
Stokes numbers, which were found to vacate regions of high
vorticity in favor of low vorticity regions. Jacobs et al.
(2004) studied the effects of compressibility and particle
density to fluid density ratio on particle dispersion in the
wake of a square cylinder and found similar particle
clustering in regions of low vorticity.
As Reynolds number increases further, the cylinder
enters the subcritical regime where the boundary layer is
still laminar, but instabilities in the vortex street begin to
form. The cylinder stays in this regime until the Reynolds
number is increased to about 4x105 (Schlichting & Gersten
2000) when the flow transitions to the critical regime. Jin
et al. (2009) injected particles into the wake of a splitter
plate with a Reynolds number of 6500 based on the plate
thickness using largeeddy simulations (LES). Snapshots
of the particle dispersion showed small particles
concentrating in the vortex center while intermediatesized
particles concentrated on the outer edges of the vortices,
which matched the experimental results. Large particles
gathered in between the vortices instead of the outer edges
of the wake. Particle studies at much higher Reynolds
numbers include the works of Chen and Huang, who both
used the discrete vortex method to investigate
twodimensional cylinders in particleladen flows (Chen et
al. 2009 and Huang & Wu 2006) at Reynolds numbers of
2.73x105 and 106, respectively. Huang and Wu found that
intermediatedsized particles do not enter the vortex cores,
but concentrate on the outer edges of the vortex structures.
The region around the vortices where few particles exist
becomes wider as the Stokes number was increased from
0.25 to 4.0. However, these studies do not take into
account the threedimensional characteristics of the
fullyturbulent wake, which would change the particle
dispersion compared to dispersion in low Reynolds number
cylinder wakes.
The objective of this study was to investigate
numerically the dispersion of water droplets in the wake of a
cylinder to better understand inertia effects on
threedimensional dispersion and preferential concentration.
The hybrid RANS/LES scheme in the computational fluid
dynamics (CFD) code WINDUS was used to compute the
unsteady wake flow of a cylinder at a Reynolds number of
8000. This scheme takes advantage of the accuracy and
efficiency of Reynoldsaveraged NavierStokes (RANS)
schemes in the attached boundary layer regions while using
LES techniques in the separated flow and wake regions,
where RANS techniques are inaccurate. Differentsized
droplets were then injected to the wake to analyze the
transport behaviors of the different particles. The droplet
trajectories were calculated using the instantaneous flow
field along with a CRW model for the subgrid turbulence
This study builds upon the work of Rybalko and Loth
(2008), who developed the CRW for RANS/LES
simulations and evaluated its performance against DNS by
considering 20 pLm droplets in the wake of a cylinder at a
Reynolds number of 800. While Rybalko and Loth
focused on numerical scheme development, the present
study focuses on the particle dispersion physics for a range
of differentsized particles and also investigates preferential
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
concentration. To the author's knowledge, this study is the
first to investigate dispersion characteristics of a multitude
of droplet sizes in a subcritical cylinder wake using the
hybrid RANS/LES scheme. Also, it is the first study to
quantify preferential concentration for cylinder wake flow.
Nomenclature
c, CRW model calibration constant
cA CRW model calibration constant
c. RANS turbulence model constant
d Particle diameter
D Cylinder diameter
f Stokes drag correction factor
k Turbulent kinetic energy
12 Number Density
Np Number of particles
Pst Vorticityconditioned concentration distribution
Re Reynolds number
Sr Strouhal number
St Stokes number
At Timestep
uf Fluid instantaneous velocity vector
u Mean velocity
SInstantaneous fluid velocity fluctuation
 Mean velocity fluctuation
SVelocity correction for NAST
Ux Freestream streamwise velocity
v Particle instantaneous velocity vector
x: Streamwise direction position
Y Lateral direction position
SSpanwise direction position
Greek letters
SRandom number vector
AGrid length scale
E Turbulent dissipation
Ai Kolmogorov length scale
A Turbulent length scale
Pu Viscosity
P Density
z Time scale
v Kinematic viscosity
n Vorticity magnitude
5 Hybrid RANS/LES clipping function
'P Decorrelation variable for discrete Markov
chain
Subscripts
D Domain
f Fluidphase
hyb Subgrid for hybrid RANS/LES
int Eddy interaction
P Dispersedphase
RANS RANS turbulence model
res Resolved
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Root mean squared
Total
Eddy traversal
Streamwise component
Lateral component
Spanwise component
Subgrid
Kolmogorov
Integral length scale
~'~I'''
7~~nc F

s~
 
(b)
Figure 1: (a) Computational grid consisting of three
zones: the small Ogrid (red), the block grid in the wake
(green), and larger Cgrid (blue) and (b) a closeup of the
cylinder and wake regions of the computational grid.
"hyb IRMNS+1 A)1 (5)
where va is the subgrid eddy viscosity and is defined as
v, = min 0.0854A ~;;l,vRAN) (6)
Equation 5 ensures that the transition between the RANS
turbulence model and the subgrid LES turbulence model is
smooth.
This hybrid RANS/LES scheme, along with a similar
computational grid, was used to predict the wake of a
cylinder at a Reynolds number of 800, which was compared
to the simulation of Rybalko, who also used the hybrid
RANS/LES scheme. All flow solver settings and
boundary conditions were the same as the current study. In
order to validate the hybrid scheme, the turbulent kinetic
energy levels in the wake were compared to DNS. The
hybrid RANS/LES and DNS solvers predicted similar k
values in zone 2 as shown in Fig. 2, providing confidence
that the hybrid scheme gives a valid cylinder wake flow
field result (Rybalko et al. 2008). The slightly different
results in the two hybrid scheme simulations may be
because a different version of the WINDUS code was used.
The version used in this study shows better performance in
the zone coupling between zones 1 and 2, which is evident
in the smaller artificial dip in k at x /D = 1.5, but also
Numerical Methodology
1. FluidPhase.l if,. , ,....
A structured 3D grid consisting of three zones was
constructed around the circular cylinder with a diameter of
6.1 cm (1D) as shown in Fig. 1. The first zone is a Ogrid
around the cylinder (161x55x43) while the second zone is a
block grid in the wake region (184x73x43). The third zone
is a larger Cgrid around zone 1 (45x93x43). The grid has
an overall diameter of 37.8 cm (40D) and a span of 12.2 cm
(2D). The initial spacing at the cylinder wall was set to
Ay = 0.0017D in order to have a significant number of grid
points to resolve the laminar boundary layer. Periodic
boundary conditions were used for the boundaries in the
spanwise direction, and the inflow conditions were set at a
Mach number of 0.1, corresponding to a velocity of 34 m/s,
and a temperature of 288.15 K. The pressure was adjusted
to 5744 Pa to obtain the correct Reynolds number. The
Roe 5thorder upwindbiased scheme was used for the
convective terms with a 2ndorder implicit time integration
scheme. The fluidphase solver was run using a timestep
of 6x106 seconds, which was found to provide a
numerically accurate and stable solution until the wake was
fullydeveloped
To model the turbulence, the NicholsNelson hybrid
RANS/LES approach with the Menter SST kco turbulence
model was chosen (Nichols & Nelson 2003). The
turbulence produced from the hybrid scheme is composed of
the resolved turbulence and the subgrid turbulence, i.e.
k = kres + khyb (1)
In RANS regions, all of the turbulent kinetic energy is
modeled using the Menter SST equations so that kres =0
and khyb = kRANS The turbulent length scale is given by
ARANilS = max 6.0 vIN~RMS, i RASk ERANS) (2)
In this expression, v/RANS is the eddy viscosity, ERANS is
the turbulent dissipation, and ORNS is the local mean
flow vorticity. In the LES regions, the subgrid kinetic
energy is defined as
khyb = kRMNSS (3)
The 5 is a clipping function defined as
II tnh4/ 23 4/3
r= 2I 1, tah4/3 _2A 4/3 (4)
where A is the grid length scale. The eddy viscosity for
the hybrid scheme is calculated using
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0.16
"I CFD
a * Left Arm
x Top Arm
Right Arm
Bottom Arm
Hybrid
  Hybrid Rybalko
DNS
0.3
g 0.2
0.1
0.0
0.12
10.08
0.04
8 12 16
Figure 2: Plot comparing the turbulent kinetic energy
values from the hybrid scheme in the cylinder wake along
the centerline at Re= 800 along with hybrid and DNS
results from Rybalk~o et al. (2008).
predicts slightly higher k downstream. However, the
turbulent kinetic energy decays at the same rate for both the
hybrid and DNS results.
To ensure that the code gives good performance at
higher Reynolds numbers (for which DNS data is not
available), the RANS/LES approach was compared to
experimental results for a cylinder at a Reynolds number of
8000 at locations in the far wake of the cylinder ranging
from x/D = 30 to x/D = 300 (Osaka et al. 1983).
Data measured at x/D=30.5 was extrapolated from
slices eight diameters away from the cylinder cross center in
the region where the flow behaves like single cylinder flow.
This data was used to evaluate the mean streamwise velocity
and the turbulent kinetic energy predicted from the hybrid
scheme, as shown in Fig. 3. The mean velocity deficit
matches the experimental results except at the centerline,
where the CFD mean deficit is less than the experiment.
The hybrid scheme results show good overall agreement in
terms of k near the centerline, but there is an over
prediction at distance of y /D of about 3. This could be
attributed to the grid resolution, which may be too coarse in
the far wake region to resolve the weakened vortex
structures in this region such that the subgrid model
contributes a significant amount of the kinetic energy. The
current study implemented a modified grid. The domain in
the current grid has an overall diameter of 40D as opposed
to 80D, which was needed to compare to the Osaka data,
and the grid resolution was made finer in the wake region
(zone 2).
2. DispersedPhase.l ifi,. , ./..
The particle equation of motion is given as
dvf 3 Pc Iv u 
(v u) C (7)
dt 4 p, d D
where if is the particle's velocity vector, u is the local fluid
velocity vector, and CD is the particle drag coefficient,
which is found using
24 4
CD & Rel/3()
0.012
0.010
0.008
S0.006
0.004
0.002
0.000
2 4
(b)
Figure 3: (a) Mean velocity deficit and (b) turbulent
kinetic energy results at x/D=30.5 for a Re=8000
cylinder wake flow compared to experimental data. The
four sets of experimental data points correspond to the four
arms of the cylinder cross where data was extracted (Osaka
et al.2008). Symmetry about y = was assumed.
and is valid for Rep <1000 (Putnam 1961). The particle
Reynolds number is defined as
Re =pv u (9)
As in Rybalko et al. (2008) and Jacobs et al. (2004), gravity
and twoway coupling were neglected herein to focus on the
effects of particle inertia. Additionally, oneway coupling
was assumed.
To take into account the contributions of the subgrid
turbulent kinetic energy, a CRW diffusion model was
implemented into the multiphase solver for the hybrid
RANS/LES scheme (Rybalko et al. 2008). The CRW
model uses a random number generator to simulate the
subgrid instantaneous fluid velocity fluctuation ..I
affecting the particle motion. A discrete Markov chain is
used to correlate u~yb with the fluctuation at the previous
time step and assumes isotropic turbulence, i.e. the rms of
the velocity fluctuations is (2kievi / 3) i
The resulting
equation for the discrete Markov chain is
...' .(t + t) = (C,'~ .+ Yuin(t)+
(1W )U/ B(t)(2kiev /3)1/2, (10)
with Y = exp (At / ,mt)
where Y is the decorrelation variable and fl is the
random number vector based on a Gaussian distribution.
The particle eddy interaction time scale rmt, is the
minimum of the eddy traversal time vr;. and the fluid
integral time scale zA giVen as
Astia 3/2
zr n. = max A, CA /~ kirvt'2 (1
trIv 1l ERANS
ERANS
The constant c~ is from the RANS turbulence model and
has a value of 0.09, and the constantS cA and c, are
calibration constants based on prior experimental data with
values of 0.78 and 0.124. The cca' term is the velocity
correction for nonhomogeneous anisotropic turbulence
(NAST) when particles move through region with gradients
in kievi and is defined as
3, = AtVkin, 1+ (12)
Results from validation performed by Rybalko in Fig. 4
show that without the CRW model, the hybrid RANS/LES
scheme underpredicts the mean dispersion of the particles
in the wake of a cylinder at Re =800 compared to DNS
results. When the CRW model is included, the results were
found to match the DNS dispersion (Rybalko et al. 2008).
Differentsized droplets with a density of 1000 kg/m
were released into the fullydeveloped wake of the cylinder.
Due to their small size, the droplets were treated as solid
spherical particles with constant density. Droplet sizes
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
were determined using the domain Stokes number defined
StD p D,/z (13)
where the particle response timer, is a function of the
particle density pp the particle diameter d, the fluid
viscosity pyf and the Stokes drag correction factor f
ppd2
Z (14)
P 18p ff
The domain time scale tD is the ratio of the cylinder
diameter and the freestream velocity
D = D/Ux (15)
Droplets categorized by their Stokes number can be
generalized into three main categories. When the Stokes
number is very small, droplets tend to follow the
macroscopic flow including smallscale turbulence
structures while droplets with large Stokes number are not
responsive to all but the largest flow structures. Droplets
with intermediate Stokes numbers ( St ~ ) experience a
delayed response to the flow and are unresponsive to the
smaller structures.
Droplet sizes released into the wake were determined
by assuming f = 1, which corresponds to Stokes flow over
the particle. To keep the notation simple, the particle
Stokes number obtained using the assumption f= 1 will
be denoted as St and will be the parameter used when
discussing the different cases in the figures and text. The
Stokes numbers ranged from 0.01 to 10, representing
droplet diameters from 2.4 pLm to 74.8 pm, in efforts to
observe the wide range of droplet behavior. Particles with
St= 0 were also released, representing fluid tracers. A
total of 250,000 droplets were released continuously at a
rate of five particles per iteration for each case,
corresponding to roughly 38 vortex shedding periods. The
particle injection location was set at .r /D= 3 relative to
the origin defined as the cylinder center. This location was
chosen to ensure the instantaneous streamwise velocity was
always positive so that none of the droplets would enter the
recirculation zone behind the cylinder. The injection
velocity was set to 22 m/s, which was found to be the mean
streamwise velocity at /D = 3.
Results and Discussion
*DNS
4 LES with CRW model
LES without CRW model
1. FluidPhase Instcuttemeous card TimeAveraged Results
The instantaneous vorticity magnitude and
streamwise velocity isosurfaces in Figs. Sa and 5b illustrate
the vortex street in the wake of the cylinder at Re = 8000 .
He The vortices shed off the cylinder at a regular frequency and
the mean Strouhal number was calculated from the lift force
4 8 12 16 time history using the equation
x/D D
Sr = (16
Figure 4: The normalized mean dispersion y@,,,, /D
for the hybrid scheme with and without the CRW model
compared to DNS in a Re= 800 cylinder wake flow
(Rybalko et al. 2008).
TU'
where T is the vortex shedding period. The Strouhal
number was found to be 0.22 which is close to the expected
value of Sr = 0.21 for a subcritical cylinder (Schlichting &
Gersten 2000).
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0.3
Subgrid
Resolved
0.2 A Total
0.1
0 4 8 12 16 20
x/D
Figure 6: The subgrid, resolved, and total turbulent
kinetic energy results along the centerline.
when comparing DNS results for cylinders at Re = 4000
and Re =10000 .
The mean streamwise velocity u, and streamwise
velocity fluctuation u;,,~,,x contours in Figs 7a and 7b were
used to determine the particle injection location. The
largest fluctuations in streamwise velocity occur just after
the boundary layer separates from the cylinder (and not
along the centerline). The peak normalized velocity
fluctuation u / ,U, behind the cylinder is roughly 0.4
located about 0.75D downstream of the cylinder trailing
edge. In comparison, the LES simulation of a cylinder at
Re=8000 by Jordan (2002) had a peak u:,~,;s,/U of
0.46 roughly 1D behind the cylinder in a contour plot
similar to Fig. 7b. Dong et al. (2006) had similar contour
plots for cylinders at Re =4000 and Re =10000 These
DNS results had peak u,1s,x, / Ux values of 0.45 and 0.50,
respectively. The peak for the Re =4000 was located
1.1D behind the cylinder while the peak for the
Re=10000 was located 0.6D behind the cylinder. This
suggests that the location of the peak velocity fluctuations is
comparable to other numerical cylinders with similar
Reynolds numbers, but the peak streamwise velocity
fluctuation is less than the other simulations. Using u,
and u,\;s,, ,the location x/D=3 along the centerline,
where us = 22 m/s, was chosen to inject the particles in
order to guarantee that the instantaneous velocity (based on
the mean and rms) was always positive so that the particles
do not go upstream and become entrained in the
recirculation zone.
2. DispersedPhase Test Conditions
Each particle case run in this simulation began at the
same fullydeveloped fluidphase solution. The
multiphase solver was run until all the droplets had left the
domain. As previously mentioned, the particle diameter
for each case was set assuming the Stoke drag correction
factor f= 1. However, this assumption is only valid for
very low Rep,,so the actual StD was calculated by using
a pathaveraged Stokes correction factor f in Eq. 8 ,
which was obtained by averaging f as each droplet
Figure 5: (a) Isosurface of the instantaneous normalized
vorticity magnitude n* = oD / U, =1 and (b) isosurface
of the instantaneous normalized streamwise velocity
us / Ux
The fluidphase solution was run until a
fullydeveloped wake was established. Then
timeaveraging was performed over a period of 60,000
iterations, which corresponds to about 10 sweeps of the flow
from the cylinder to the end of the domain. This was
found to give a sufficiently converged mean velocity field
and the resolved velocity fluctuations. The resolved
turbulent kinetic energy was calculated using the equation
k~es=[ (ux, + (U;)+ ( (17)
In addition, the subgrid and/or RANS turbulent kinetic
energy (k;,,,z) was averaged in time. The total turbulent
kinetic energy is the sum of these two components. Figure 6
shows the averaged k;yb,,zkrez, and k obtained along the
centerline downstream of the cylinder. The location
x / D = 0.5 corresponds to the trailing edge of the cylinder
where k,e, = 0. Near the wall, k;,,,z has a significant
contribution to the total kinetic energy as the hybrid scheme
transitions from RANS to LES. Significant turbulence is
produced in the near wake where the recirculation zone is
located and peaks at around x:/D=1.5. The turbulent
kinetic energy results mimic the results obtained by Rybalko
in Fig. 2, which is expected since both cylinders are in the
subcritical regime (Rybalko et al. 2008). The peak kinetic
energy occurs closer to the cylinder for the Re= 8000
case. This is consistent with the trend that the largescale
Karman vortices are formed closer to the cylinder as the
Reynolds number increases, which Dong et al. (2006) found
I I I i\
ill
(a)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Table 1: Summary of Time Scales and Stokes Numbers
for Each Simulation Case
St 0 0.01 0.1 1 10
d [pLm] 0 2.37 7.48 23.65 74.80
zp [s] 0 1.75x10 s 1.69x10 4 1.54x103 1.16x10 2
zD [s] 1.79x103 1.79x103 1.79x103 1.79x103 1.79x103
Zn [s] 2.08x10 s 2.05x10 s 2.23x10 s 2.51x10 s 2.26x10 s
StD 0 0.0098 0.095 0.86 6.5
Stn 0 0.85 7.59 61.6 513
The majority of the cases have large Stn where the
droplets are not sensitive to the Kolmogorovscale
turbulence, but the 2.47 pLm particles have Stn ~ 1 where
the smallest structures have some effect. Table 1
summarizes the time scales and Stokes numbers for each
case.
3. DispersedPhase Dispersion Results
The droplet trajectories differ greatly depending on
their size and inertia. Illustrated in Fig. 8 are instantaneous
snapshots of the particles with contours of the normalized
vorticity. Since the computational domain length was only
2D in the spanwise direction, particles that crossed the
periodic boundaries were adjusted in the figure to give an
idea of the spanwise dispersion in a larger domain for the
topdown view. For example, if a droplet crossed the
boundary z /D = +1, the periodic boundary would move
particle to the boundary z /D= 1 so that the droplet
would stay in the computational domain. Through
postprocessing, the droplet position was adjusted to have a
position z /D >+1.
The tracer particles in Fig. 8a have the same inertia as
the fluid, thus will follow the fluid. The tracer particles
immediately begin to follow the fluidphase flow after being
injected into the wake and will also become entrained into
the vortices. The droplets with St= 0.01 are also
similarly sensitive to the resolved vortices in the wake but
have less dispersion associated with subgrid kinetic energy,
i.e. the smallscale random diffusion is reduced as compared
to St=0, which is consistent with the finite Sth for this case
(Table 1) Droplets with St = 0.1 do not enter the vortices
at all. In fact, the droplets trajectories in Fig. 8c show that
the droplets prefer the outer edges of the turbulent structures.
This phenomenon is known as preferential concentration.
The droplets have high enough inertia that the vortices act
as centrifuges and disperse the droplets into the regions of
lower vorticity on the outer edges of the wake. The same
phenomenon is seen in Fig. 8d for the droplets with St = 1.
By comparing Figs. 8c and 8d, the smaller scales of
turbulence have an effect on the St = 0.1 droplets, causing
them to locally spread, but the St = keep to a single
continuous stream. This is quantified by the Kolmogorov
Stokes number where Stn = 7.63 for the St = 0.1
droplets and St~z=61.2 for the St = droplets. The
smallscale turbulence has some effect on the smaller
droplets, but at Stn = 61.2, zp is much to large for the
particle to respond to smallscale turbulence. In Fig. 8e, the
St =10 droplets are less affected by the wake flow.
0 44
0 4
0.32
0.28
0.24
0.2
0.08
0.04
0
Figure 7: (a) The normalized mean streamwise velocity
contour ux / U, and (b) the normalized mean streamwise
velocity fluctuations urms,x /U
traversed the domain. The actual StD values were found
to be less than the assumed values with a greater disparity at
the larger particles where the Rep <1 assumption begins
to break down.
In addition to the domain Stokes number, the
Kolmogorov Stokes number Stn was also calculated,
This value is a measure of the droplet sensitivity to the
smallest turbulent structures where velocity gradients are
dissipated into heat. The length scale z74 and time scale
Zn are defined as
1/2 (18)
where E was obtained by averaging the dissipation at each
particle location during the simulation. Then Stn was
found using
St, = z, / z, (19)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Although some of the stronger vortices impact the stream of
droplets, they do not have enough strength to disperse the
droplets to the outer edges of the wake, and the heavy
droplets cut through the vortex structures.
Similar particle behavior was found in other studies of
particles in cylinder wake flow. Luo's instantaneous
snapshots of particle dispersion show clustering of particles
with intermediate Stokes numbers ( St = 0.1 and St = 1) on
the outer edges of the vortices with no particles entering the
vortex structures while particles with small Stokes
(St=0.01) numbers become entrained in the vortices, O 0.40.8 1.2 1.6 2 2.42.8 3.236 4
which is consistent with the droplet dispersion seen in Fig. 8. (b)
Since Luo's particles were released upstream of the cylinder
instead of the directly in the wake, the heavy particles
experience preferential concentration in the outer edges of
the wake due to collisions with other particles and with the
cylinder (Luo et al. 2009). The result is particle dispersion
similar to intermediate Stokes number particles. When
particles are injected into the wake as in Fig. 8e, the
particles stay near the centerline, resulting in a very different r ~ 1
trajectory than if the heavy particle was released upstream
of the cylinder. Wu et al. (2008) also found that small
Stokes number particles become entrained in the vortices,
and particles with intermediate Stokes numbers experience 0 0 4 0 8 1 2 1 6 2 2 4 2 8 3 2 3 6 4
the centrifugal effect. Since gravitational effects were c
included, the heavy particles would sink as they were
transported downstream. The stream of droplets in Fig. 8e
is most similar to the snapshot of 50 pLm particles injected
into the wake of a thick splitter plane. Those particles
concentrate along the centerline but also move around the
vortex structures (Jin et al. 2009).
The effect of a higher Reynolds number flow is also
evident in Fig. 8. The vortex street seen in the cylinders of T
Jacobs et al. (2004) and Wu et al. (2009) are very clean and
structured such that the particle trajectories are similar
between each vortex pair. However, the vortex structures
seen in Fig. 8 are much more chaotic due to the more 0 0 4 0 8 1 2 1 6 2 2 4 2 8 3 2 3 6 4
unstable wake. This causes the droplets to spread more(d
irregularly, which is especially evident in the small Stokes
number cases seen in Figs. 8ac.
In order to quantify how much the droplets spread, the
mean lateral diffusion relative to the injection location was
calculated at flux planes spaced in increments of 0.5D
downstream of the inj section location using the equation.
p,rms2, _1 Ci2 (20)
p,tot i=1
0 0.4 0.8 1.2 1.6 2 2.4 2.8 3.2 3.6 4
(e)
Figure 8: Side views of the instantaneous droplet
dispersion with slices of the instantaneous normalized
vorticity magnitude 0 = oD / U at z /D = for
droplets with (a) St = 0 (b) St = 0.01, (c) St = 0.1, (d) St = 1,
rl Iand (e) St = 10.
O 0.4 0 8 1.2 1.6 2 2 4 2.8 3.2 3.6 4
0 0 4 0 8 1.2 1.6 2 2.4 2 8 3 2 3 6 4
0 0 4 0 8 1 2 1 6 2 2 4 2 8 3 2 3 6 4
(b)
O 0.4 08 1.2 1.6 2 2.4 2.8 3.2 3.6 4
I__P
0 0.4 0.8 123 1.6 2 2.4 2.8 3.2 3.6 4
(d)
0 0.4 0 8 1.2 1.6 2 2.4 2.8 3.2 3.6 4
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
St=0
St=0.01
4 1 A St=0.1
m St=1
St=10
S3 125k droplets
O 4 8 12 16 20
x/D
Figure 9: The normalized timeaveraged dispersion in the
lateral direction yp',,2,,/ D for 250,000 droplets along
with curves for 125,000 droplets to compare convergence.
where 3', is the droplet position in the lateral direction and
Np,tor is the total number of droplets passing through each
flux plane.
All of the droplets are injected with a streamwise
velocity meaning that initially, the droplet does not follow
the fluid velocity. Taking a look at the lateral diffusion
statistics at the injection location in Fig. 9 can give insight
to how long it takes to respond to the wake flow. For
example, the droplets with St = 10 do not begin to diffuse
until about 2D downstream of the injection location while
droplets with St = do not begin to diffuse until they
travel around 1D downstream. The droplets with St<1
begin to diffuse immediately, suggesting that they are very
responsive to the wake flow. Once the droplets lose the
initial injection momentum and begin to follow the wake
flow, the St=0.1 and St=1 begin dispersing at a faster
rate than the tracer particles while the St= 0.01 closely
follows the dispersion of the tracer particles because of the
centrifugal effect seen in Figs. 8c and 8d. The end result is a
larger dispersion range for the intermediate Stokes number
droplets, especially the St = droplets.
Figure 9 also shows the mean dispersion 3'p,;;;;,2 / D2
calculated from the first 125,000 droplets that crossed the
flux planes. The results do not change when the number of
droplets is doubled to 250,000 for the low Stokes number
cases. For the St = 1 and St = 10 droplets,
3'p~,;;; /D2 changed slightly in the latter portions of the
domain. The difference in yp',,,,,/ D going from
125,000 to 250,000 droplets with St = which had the
largest difference, was 0.16, or 4% change. Thus, it was
deemed that releasing additional particles would not have
significant impact on the results.
The dispersion in the spanwise direction is less than in
the lateral direction because the cylinder wake vortex
structures are a twodimensional phenomenon and smaller
scales of turbulence are the cause of spanwise dispersion as
shown by the more uniform spread in the topdown views of
Figure10: Top views of the instantaneous droplet
dispersion with a slice of the instantaneous normalized
vorticity magnitude O'=noD/U at v/D=0 for
droplets with (a) St = 0 (b) St = 0.01, (c) St = 0.1, (d) St = 1,
and (e) St = 10 .
II
M 
St=0.O
St=0.01
+ St=1
+ St=10
31 2 3 4
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0.5
0.4
0.3
N 0.2
0.1
0
0.32
0.28
0.24
0.20
* 0.16
0.12
0.08
0.04
0.00
01 23 4
y~/D
0 4 8 12 16 20
Figure 11: The normalized timeaveraged dispersion in the
spanwise direction Zp,,,.,,/ D 
the St = 0, St= 0.01, and St= 0.1 particles (Figs. 10a,
10b, and 10c) Generally, the particles stay within
x / D = 01, but occasionally a vortex structure will disperse
them further. Particles with St = do not spread much
until further downstream of the cylinder, and particles with
St= 10 have nearly zero dispersion in the spanwise
direction in Fig 10e.
The mean spanwise diffusion relative to the injection
location can be calculated using the equation
Sp,,.,,, N i4 (21)
p,ror i=1
where zi is the droplet position in the spanwise direction.
Figure 11 shows that there is much less spanwise dispersion
than the lateral dispersion. Even though there were distinct
trends in Yp,,",, there are not too many clear trends in
p,,,,, Again, the dispersion of the St= 0.01 droplets
closely follows the fluid tracer particles, and although the
dispersion rates, i.e. the slopes in Fig. 11, are different for
the different droplets, zp,,,2, is roughly the same at the
end of the computational domain. Luo, who calculated the
mean lateral and spanwise dispersion of all particles as a
function of time, found similar behavior for this range of
Stokes numbers, although the spanwise dispersion for
St =100 particles was significantly less than the other
particles (Luo et al. 2009).
4. DispersedPhase Preferential Concentration Results
Profiles of the timeaveraged particle number density
were calculated at each of the flux planes to provide insight
in how the concentration evolves through the cylinder wake.
When normalized by the total number of particles that pass
through the flux plane the expression for the normalized
number density is
N
up= (22)
Np,tor
where Np is the number of droplets inside the bin.
0.32
0.28
0.24
0.20
*C 0.16
0.19
0.08 1
0.04 1
0.00
0.32
0.28
0.24
0.20
* 0.16
0.12
0.08
0.04
0.00
St=0
St=0.01
+ St=0.1
 St=1
+ St=10
0 1 2
y/D
3 4
Figure 12: Timeaveraged normahized particle number
density at .x/D=8, (b) .x/D=13, and (c) .x/D=18
which correspond to SD, 10D, and 15D downstream of the
particle injection location. Symmetry about y = was
assumed.
sc+
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
nearly twice as likely to concentrate near 0 = 0
compared to tracer particles in isotropic turbulence which is
higher than what was found in this study, but the point
which the two curves intersect occurs around n* = 0.8
The difference in the amount of blas may be attributed to
the stronger vortex structures present in the cylinder wake.
These trends become exaggerated at higher Stokes numbers.
The St = droplets are about 2.3 times more likely to be
in the low vorticity regions of the wake, while the St = 1
and St = 10 droplets are roughly eight times more likely to
see vorticity magnitudes near zero. The curves for droplets
with St > 0.1 in Fig. 13 are very similar for n* > 0.5 in
that there are half as many particles in the high vorticity
regions compared to the fluid tracers.
Figure 14 shows the average vorticity seen by all the
droplets sampled to calculate Pn' as a function of the
pathaveraged Stokes number StD For the tracer
particles, the average normalized vorticity was 1.02. Small
StD have a similar n,, but for intermediate StD >
droplets cluster in lower vorticity regions yielding a lower
Figures 12a, 12b, and 12c show np profiles for the
flux planes located at x /D= 8 x/D = 13, and
x /D = 18 In these figures, symmetry about y = 0 was
assumed. As expected, the width of the curves gets wider
further downstream as the droplets disperse. The droplets
with St= 10 have the narrowest number density profiles
suggesting that droplets stay closer to the centerline than the
other cases. The St= 0.01 droplets and the fluid tracers
yield matching curves with a relatively constant np near
the centerline which decreases as the y /D increases.
The number density profile at x/D = 8 for droplets with
St= 0.1 and St = exhibit a dip at y = 0. This can be
attributed to the centrifugal effect of the stronger vortex
structures forcing the droplets to the outer edges of the wake
instead of near centerline where the vortices are located.
Over time, this would result in a reduction in concentration
droplets around y = 0. Further downstream, the dip in
the profiles for St = 0.1 is less pronounced because local
spreading of the droplets from the smallscale turbulence
seen in Fig. 8c causes a more uniform concentration of the
droplets. For the droplets with St = the spreading is
less, and the droplets stay in a continuous stream. This
keeps the majority of the droplets at the outer edges of the
wake. This preferential concentration phenomenon can be
seen in the timeaveraged number density contours of the
splitter plate wake of Jin (2009), but it is not explained in
depth. The curves in Fig. 12 are consistent with the mean
lateral dispersion yp~,rm2 in Fig. 9, with St=1 having a
much larger yp~,rms because np has a peak further away
from the centerline as opposed to a more Gaussianlike
profile. The number density profiles suggest that the
largescale turbulent structures help increase the spread of
the particles, but it is the smallscale structures that break up
the stream of droplets and give a more uniform
concentration.field into account, Po for each case was
normalized by Po for the fluid tracers. The normalized
vorticityconditioned particle concentration distribution was
defined as N
Poz = (23)
Np,totAQ* [Po Ist=o
where Np is the number of droplets in the bin, Np,tor is
the total number of droplets sampled, and AG* is the bin
size of the normalized vorticity magnitude
The instantaneous fluid vorticity magnitude at each
droplet sees was recorded every 1000 iterations, or roughly
every vortex shedding period, which were then used to
calculate Pn'. The results displayed in Fig. 13 show that
there is very little difference between the St = and
St = 0.01 droplets. However, there is a small increase in
Pn' when n* <0.5 and, on average, Pn' <1 when
n* > 0.5 because the St = 0.01 droplets still have a small
amount of inertia. The results for these droplets, which
have a Stn = 0.85, are similar to the Stn = 1 particles of
Sundaram and Collins (1997) which were also found to
cluster in lower vorticity regions. Those particles were
St=0
St=0.01
& St=0.1
& St=1
+ St=10
Figure 13: The normalized vorticityconditioned particle
concentration distribution Pa*, which were normalized by
the St = values.
1.2
1.0 
,0.8
0.6
0.4
0.2
0.0
1.0E03 1.0E02 1.0E01 1.0E+00 1.0E+01
StD
Figure 14: The average normalized vorticity plotted as a
function of the pathaveraged Stokes number StD, which is
plotted in logscale.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
References
Chen, B., Wang, C., Wang, Z., and Guo, L., "Investigation
of gassolid twophase flow across circular cylinders with
discrete vortex method," Applied Thermal Engineering, Vol.
29, pp. 14571466 (2009)
DeAngelis, B., Loth, E., Lankford, D. and Bartlett, C.S.
"Simulations of Turbulent Droplet Dispersion in
WindTunnel Icing Clouds," AIAA Journal of Aircraft, Vol.
34, No. 2, pp. 213219, MarchApril (1997)
Dong, S., Karniadakis, GE., Ekmekci, A., and Rockwell, D.,
..A combined direct numerical simulationparticle image
velocimetry study of the turbulent near wake," J. Fluid
Mech., Vol. 569, pp. 185207 (2006)
Huang, Y., and Wu, W., "Numerical Study of Particle
Distribution in Wake of LiquidParticle flows Past a
Circular Cylinder Using Discrete Vortex Method," Applied
Mathematics and Mechanics, Vol. 27, No. 4, pp. 535542
(2006)
Jacobs, GB., Kopriva, D.A., and Mashayek, E,
Compressible Subsonic ParticleLaden Flow over a Square
Cylinder," AIAA Journal of Propulsion and Power, Vol. 20,
No. 2, pp. 353359, MarchApril (2004)
Jin, H., Chen, S., Chen, L., and Fan J., "LES/FDF
simulation of particle dispersion in a gasparticle two phase
plane wake flow," Sci China Ser ETech Sci, Vol. 52, No. 10,
pp. 29432951 October (2009)
Luo, K., Fan, J., Li, W., and Cen, K., "Transient,
threedimensional simulation of particle dispersion in flows
around a circular cylinder (Re = 140260)," Fuel, Vol. 88,
pp. 12941301 (2009)
Jordan, S.A., "Investigation of the cylinder separated
shearlayer physics by largeeddy simulation," Int. J. Heat
and Fluid Flow, Vol. 23, pp. 112 (2002)
Nichols, R.H., and Nelson, C.C., "Application of Hybrid
RANS/LES Turbulence Models," 41st Aerospace Sciences
Meeting and Exhibit, Reno, NV, No. AIAA20030083,
January (2003)
Osaka, H., Nakamura, I., Yamada, H., Kuwata, Y.,
Kageyama, Y., "The Structure of a Turbulent Wake behind a
Cruciform Circular Cylinder," Bulletin of the JSME, Vol. 26,
No. 214, pp. 521528 (1983)
Putnam, A., "Integrable form of droplet drag coefficient,"
American Rocket Society Journal, Vol. 31, pp. 14671468,
October (1961)
Rybalk~o, M., Loth, E., Lankford, D., "LES SubGrid
Diffusion for Lagrangian Particles," ASME Fluids
Engineering Conference, Jacksonville, FL, No.
FEDSM200855207, August (2008)
0 At the larger StD case, the average vorticity
begins to increase and as the Stokes number increases
further, the average vorticity is expected to continue
increase as larger droplets cut through vortex structures due
to their large inertias.
Conclusions
Droplets were injected into the wake of a subcritical
cylinder at a Reynolds number of 8000 in order to study
how they disperse. The droplets with small Stokes
numbers tend to follow the fluid flow. Droplets with
intermediate Stokes numbers have a more significant inertia
and exhibit a centrifugal effect from the vortex street,
causing droplets to gather in the low vorticity regions in the
outer edges of the wake. This causes the particles to
spread further than the smaller particles that follow the
fluidphase. The droplets with St =10 are less
responsive to the cylinder wake and have the least
dispersion. Diffusion in the spanwise direction was
roughly the same for the different sized droplets released in
this study, although injecting even larger droplets may see
different results.
Timeaveraged particle number density profiles at the
flux planes show that the droplets with St = display the
preferential concentration effect as the majority of the
particles are located in the outer edges of the wake. The
small Stokes number droplets have more of a Gaussianlike
profile due to local spreading caused by smallerscale
turbulence while the large Stokes number droplets have a
Gaussianlike profile because the vortices are not strong
enough to overcome the droplet inertia. This type of
profile was not seen in previous studies and may help
quantitatively compare preferential concentration effects
different particles as they traverse downstream. The
vorticityconditioned particle concentration distribution also
shows that droplets with intermediate Stokes numbers prefer
the lowvorticity regions when compared to the fluid tracer
particles. This technique allows particles to be categorized
by vorticity and may allow preferential concentration
comparison of particles in different turbulent flows
To expand the understanding of droplet dispersion in
the wake of bluff bodies, futures studies should consist of an
expanded domain, in both the streamwise and spanwise
directions. Additionally, the different cylinder flow
regimes may yield different droplet dispersion statistics do
to the changes of the wake structure
Acknowledgements
This work was supported by the icing branch of the NASA
John H. Glenn Research Center for the Icing Research
Tunnel with Colin Bidwell as the technical manager.
Computer time was provided by the National Center for
Supercomputing Applications
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Schlichting, H., and Gersten, K., Boundary Layer Theory,
SpringerVerlag, 8th Ed., New York (2000)
Sundaram, S., and Collins, L., "Collision statistics in an
isotropic particleladen turbulent suspension. Part 1. Direct
numerical simulations," J. Fluid Mech., Vol. 335, pp. 75109
(1997)
Wu, Z., "Numerical study on heavy rigid particle motion in
a plane wake flow by spectral element method," Int. J.
Numer. Meth. Fluids, Vol. 61, pp. 536551 (2009)
