7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A numerical study of finite size particles in homogeneous turbulent flow
T. Doychev and M. Uhlmann
Institute for Hydromechanics, Karlsruhe Institute of Technology, 76128 Karlsruhe, Germany
{todor.doychev, markus.uhlmann} @kit.edu
Keywords: particulate flow, homogeneous turbulence, resolved particles, DNS, immersed boundary method
Abstract
The present study addresses finite size and finite Reynolds number effects in spatially homogeneous turbulent flow
seeded with solid particles. The ratio between particle diameter and Kolmogorov length scale is chosen in the
interval 525; the Reynolds number based upon the particle diameter and the terminal velocity in ambient flow is of
order 0(100). Under these conditions, the pointparticle approach loses its validity, and we resort to fullyresolved
simulations, realized with the aid of an immersed boundary method. The volume fraction of the solid phase is set to
5 10 3, whereby dominant effects of interparticle collisions are avoided. Two flow configurations are investigated:
pure sedimentation in ambient fluid and sedimentation with decaying turbulent background flow, both without mean
velocity gradients and with periodicity in all spatial directions. It is found that turbulence has a strong effect upon
the settling velocity, initially causing a significantly slower average vertical particle velocity, and at later times (when
the turbulence intensity has decayed to a weaker level) a slightly more rapid settling. At the same time turbulence is
found to affect the interparticle distances, bringing particles closer together than pure sedimentation. Further effects
of settling particles upon the turbulent kinetic energy, the dissipation rate and anisotropy are discussed, as well as
probability density functions of particle velocities and forces.
Introduction
The interaction between finitesize heavy particles and a
surrounding flow field is a fundamental hydrodynamic
problem with applications to various realworld con
figurations, such as raindrop formation in clouds, flu
idized bed reactors and combustion devices. In past
studies the details of the nearfield around the particles
have often been neglected, thereby excluding a descrip
tion of wake effects (cf. Balachandar and Eaton 2010,
for a recent review). Here we consider sedimentation
problems involving particles with diameters between 10
and 25 Kolmogorov length scales and average particle
Reynolds numbers (based upon diameter and settling ve
locity) of 0(100). For this purpose we have carried out
direct numerical simulations with interface resolution.
Today these computations are still costly and limited to
a narrow parameter range. However, numerical experi
ments are able to complement laboratory measurements,
in particular by providing data at high spatial and tem
poral resolution.
The problem of particles settling in a turbulent back
ground flow is characterized by a large set of non
dimensional parameters. Let us first consider the case
of sedimentation of a set of monodisperse spheres in
an ambient fluid. Given the fluid density pf, the kine
matic fluid viscosity vf, the vector of gravitational ac
celeration g on the one hand, and the particle diame
ter D, particle density pp and solid volume fraction sb
on the other hand, dimensional analysis shows that the
problem is determined by three nondimensional param
eters. One has already been mentioned (the solid vol
ume fraction); the other two can be taken as the den
sity ratio pp/pf and the ratio between gravitational and
viscous velocities, D g\D/.v. Now, when adding a
turbulent background flow to the picture, we introduce
at least two additional reference quantities to the prob
lem, namely a fluid velocity scale ... f and a fluid flow
length scale C. Therefore, there will be five independent
nondimensional parameters, the two additional ones be
ing, e.g., a fluid flow Reynolds number ... fj/v and a
length scale ratio D/i. In the present work we con
sider fully developed homogeneousisotropic turbulence
as the background flow, implying that the scales of the
singlephase flow field can be represented e.g. by the
Taylor microscale A and the r.m.s. velocity u,,, m
2q/73 (where q2 = E(K)dK is twice the integral of
the turbulent kinetic energy, here as evaluated in spectral
space over radial wavenumbers K).
The problem of settling particles can be character
ized by various timescales: the gravitational scale T7
ID/.g, the purely viscous scale 7, D2/v and the
Stokes drag scale Tp ppD2/(pf18v). The back
ground turbulence is often characterized by the Kol
mogorov timescale, T, ,/, and the largeeddy
turnover time TL k/E.
These dimensional considerations illustrate the ampli
tude of the parameter space. It is therefore in general dif
ficult to find reference studies matching a given param
eter point. The experiment of Parthasarathy and Faeth
(1990) provides one of the most relevant datasets for the
present simulations. It provides measurements of fluid
and particle data of particles settling in a fluid which is
initially at rest, carried out for parameter values mostly
comparable to those investigated in the present study, al
beit at a much lower solid volume fraction. Their value
is approximately = 104, whereas we are targeting
flows with a somewhat higher solid volume fraction of
# 5 103.
The purpose of the present work is to investigate the
influence of fluid turbulence upon the settling velocity
of solid spherical particles. Contrary to previous studies
dealing with single particles fixed in space (e.g. Bagchi
and Balachandar 2003; Zeng et al. 2010), our particles
are mobile, they are present in large numbers, and the
flow is turbulent a priori. In order to separate the tur
bulence effect from the effects of mobility and collectiv
ity, we consider configurations with and without back
ground turbulence.
Numerical method
The present simulations have been carried out with the
aid of a variant of the immersed boundary technique (Pe
skin 1972, 2002) proposed by Uhlmann (2005a). This
method employs a direct forcing approach, where a lo
calized volume force term is added to the momentum
equations. The additional forcing term is explicitly com
puted at each time step as a function of the desired parti
cle positions and velocities, without recurring to a feed
back procedure; thereby, the stability characteristics of
the underlying NavierStokes solver are maintained in
the presence of particles, allowing for relatively large
time steps. The necessary interpolation of variable val
ues from Eulerian grid positions to particlerelated La
grangian positions (and the inverse operation of spread
ing the computed force terms back to the Eulerian grid)
are performed by means of the regularized delta function
approach of Peskin (1972, 2002). This procedure yields
a smooth temporal variation of the hydrodynamic forces
acting on individual particles while these are in arbitrary
motion with respect to the fixed grid.
Since particles are free to visit any point in the com
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
putational domain and in order to ensure that the regu
larized delta function verifies important identities (such
as the conservation of the total force and torque dur
ing interpolation and spreading), a Cartesian grid with
uniform isotropic mesh width Ax A y Az is em
ployed. For reasons of efficiency, forcing is only applied
to the surface of the spheres, leaving the flow field inside
the particles to develop freely.
The immersed boundary technique is implemented
in a standard fractionalstep method for incompress
ible flow. The temporal discretization is semiimplicit,
based on the CrankNicholson scheme for the viscous
terms and a lowstorage threestep RungeKutta proce
dure for the nonlinear part (Verzicco and Orlandi 1996).
The spatial operators are evaluated by central finite
differences on a staggered grid. The temporal and spatial
accuracy of this scheme are of second order.
The particle motion is determined by the Runge
Kuttadiscretized Newton equations for translational and
rotational rigidbody motion, which are explicitly cou
pled to the fluid equations. The hydrodynamic forces
acting upon a particle are readily obtained by summing
the additional volume forcing term over all discrete forc
ing points. Thereby, the exchange of momentum be
tween the two phases cancels out identically and no spu
rious contributions are generated. The analogue proce
dure is applied for the computation of the hydrodynamic
torque driving the angular particle motion.
During the course of a simulation, particles can ap
proach each other closely. However, very thin inter
particle films cannot be resolved by a typical grid and
therefore the correct buildup of repulsive pressure is not
captured which in turn can lead to possible partial 'over
lap' of the particle positions in the numerical computa
tion. In practice, we use the artificial repulsion potential
of Glowinski et al. (1999), relying upon a shortrange
repulsion force (with a range of 2Ax), in order to pre
vent such nonphysical situations. Essentially the same
method is used for the treatment of particles approaching
solid walls.
The present numerical method has been submitted
to exhaustive validation tests (Uhlmann 2004, 2005a,b,
2006a), as well as grid convergence studies (Uhlmann
2006b). The computational code has been applied to the
case of vertical plane channel flow (Uhlmann 2008).
Setup of the simulations
We have performed two different types of numerical ex
periments. In the first one, particles are released from
rest in ambient fluid; the cases realized with this con
figuration of "pure a'ldillicill.lliln are denoted as "SA"
and "SB". The second type of configurations features a
turbulent background flow, which is homogeneous, (ini
Table 1: Particle properties. In all cases considered
here, particles possess the same relative density with re
spect to the fluid, viz. pp/pf = 1.5, as well as a global
solid volume fraction of bs = 5 10 3. Two differ
ent values for the gravitational acceleration are consid
ered, leading to the following nondimensional parame
ters. Note that the particle Reynolds number ReDo is
based upon the "nominal" settling velocity given by a
balance between drag and immersed weight, using the
standard drag formula (Clift et al. 1978).
cases D gD/v ReDoo
A, SA 93.9125 66.53
B, SB 141.4600 145.85
Table 2: Properties of the initial homogeneousisotropic
flow field.
cases Rec kmaxrI Lint/Lbox
A, B 180.6 0.884 0.0995
SA, SB (fluid at rest)
tially) isotropic, and which decays in time. The simula
tions of this type are denoted as cases "A" and "B". All
simulations are performed with the same particle/fluid
density ratio of pp/pf 1.5. Likewise, we chose
the same value for the global solid volume fraction of
t s 5 103 in all cases, corresponding to Np 3038
particles. Two different values for the gravitational ac
celeration were chosen (cf. table 1) such that a single
particle in ambient fluid experiences a relative flow at
a diameterbased Reynolds number of 66.53 (145.85)
in cases A and SA (B and SB), according to a balance
between drag and immersed weight, using the standard
drag formula (Clift et al. 1978)
The initial turbulent flow field is generated through
spectral simulation of forced homogeneousisotropic
turbulence, using the code of A. Wray (cf. Jim6nez and
Wray 1998). When the flow has reached a statistically
stationary state (verified by monitoring skewness and
normalized dissipation rate) a flow field is spectrally in
terpolated upon the staggered finitedifference grid. The
turbulence properties of the singlephase flow are re
Table 3: Twophase flow properties: the following pa
rameters apply to both cases with turbulence background
flow (based upon quantities at the initial time t 0).
cases D/r St, Stint TL/Tp
A, B 22.5 59.22 2.66 0.84
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
10
101
S102
ko
103
100 101
t/Tp
Figure 1: Turbulent kinetic energy k (normalized with
the initial value of cases A and B, ko) as a function of
time (normalized with the particle relaxation time Tp).
*, singlephase flow; case A;  case SA;
 case B;  case SB. The straight reference
line ( ) varies as 14.
ported in table 2. The actual simulation with the finite
difference code is then started from this initial field (ar
bitrarily setting the time to zero, t 0), adding par
ticles at random positions and assigning them an ini
tial translational velocity equal to the local fluid veloc
ity. Further pertinent parameter values of the simula
tions with turbulent background flow are given in ta
ble 3. The ratio between the particle diameter and the
Kolmogorov length is initially 26.7 (0.85 with respect to
the initial Taylor microscale); this ratio then decreases
due to turbulence decay, as will be seen below. Con
cerning the Stokes number, two definitions will be con
sidered. In the first, the particle response time (based on
Stokes drag) is divided by the Kolmogorov time scale,
viz. St,7 = Tp/. In the second definition, the ratio be
tween the same particle time scale and the integral fluid
time scale is formed, yielding Stit = TplTit. At the
time of particle injection, the value of St, is approxi
mately 60, whereas the value of Stint is 2.66.
The present simulations are realized with 5123 grid
points. This mesh provides a moderate resolution of
initially Ax/r = 0.28 (equivalently kmax,,r 11 I
which, however, improves quickly over the course of the
simulation. The particles are resolved with D/Ax
7.5, which is approximately the same resolution em
ployed by Lucci et al. (2010). The time step is chosen
such that the CFL number remains below 0.5.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
10 20 30 40 50
Figure 2: Anisotropy of the vertical diagonal
component of the Reynolds stress tensor, b33
(w'w')f/(2k) 1/3. Linestyles as in figure 1.
Results
Figure 1 shows the temporal evolution of the turbulent
kinetic energy (TKE) defined as k = '.'..'f/2 for all
four cases. Here and in the following the notation () f is
used for spatial averages over the volume occupied ex
clusively by the fluid. As can be observed from the fig
ure, the overall energy of the fluctuating fluid velocity is
only affected at later times by the presence of the parti
cles. In case B, the difference becomes appreciable after
approximately 10 particle time units (equivalent to 12
initial largeeddy turnover times TL), while in case A a
difference of only 10% with respect to singlephase flow
is measured at t/tp 20. As will be discussed below,
the initial turbulence intensity is very strong compared to
the fluctuations induced by the particle motion. There
fore, the twophase flow in cases A and B at early times
can be considered as nearly oneway coupled, at least as
far as kinetic energy is concerned. Conversely, at later
times the flow becomes more and more influenced by the
presence of particles, and finally the statistics of cases
A/SA and B/SB should converge. The figure shows that
the asymptotic limit of "pure s'diniimci,'l.l, i has not yet
been reached in cases A and B by the end of the present
runs.
The anisotropy of the fluctuating velocity field
can be measured by the anisotropy tensor bij
'.'..')/(2k) 6ij/3. Here only the diagonal ele
ments are nonzero, and both horizontal components
are identical by symmetry. Therefore only component
b33 is shown in figure 2 (the other two being given by
bl b22 = b33/2). In both cases without tur
bulent background flow a value of approximately 0.5
is quickly reached and remains constant. The fluctua
Figure 3: The ratio between the particle diameter and
the Kolmogorov length scale as a function of time.
I 80
100
120
140
0 10 20 30 40 50
Figure 4: The temporal evolution of the average particle
settling velocity, w,l = (wp)p (w)f, scaled by the
viscous reference velocity.
tions are therefore principally concentrated in the verti
cal direction (i.e. aligned with the settling velocity). The
experimental study of P.nili.s.n.Il.\ and Faeth (1990)
has likewise shown that sedimenting particles lead to an
anisotropic velocity field. In their cases (at low solid vol
ume fractions of = 10 4), the anisotropy was found
to measure b33 A at various particle Reynolds num
bers (including the present range). With background
turbulence, on the other hand, the anisotropy remains
initially close to zero and only slowly approaches the
strongly anisotropic value obtained in pure sedimenta
tion. The temporal growth of b33 seems to be roughly
constant, at least for values of b33 < 1/3.
The ratio between the particle diameter and the Kol
mogorov length scale, = (v3/)(1/4), in cases A, B is
shown in figure 3. While the background turbulence de
cays the Kolmogorov scale grows, leading to a decrease
I     
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
tITp
Figure 5: The temporal evolution of the difference be
tween the average particle settling velocity in turbulent
conditions, wre, and without background flow, w(,,
normalized by the value without background flow.
in time of the relative diameter D/r. On the figure
we can observe that the value is approximately halved
over the first five time units. Asymptotically, values of
7 (10) are obtained in cases A (B). This comparison of
scales shows that the relative particle size (with respect
to the spectrum of turbulent flow scales) varies in time.
It should be noted, however, that due to the developing
strong anisotropy (favoring the vertical velocity compo
nent) the scales will also be highly anisotropic. Their
analysis should be performed by means of spatial corre
lation functions, separately by horizontal/vertical coor
dinate direction.
The average particle settling velocities wrel are shown
(in viscous scaling) in figure 4. This quantity is de
fined as the difference between the mean velocities of
the two phases, viz. wrl (wp)p (w)f (wp being the
vertical component of the particle velocity, w the corre
sponding fluid velocity component, and the operator ()p
denotes averaging over the number of particles), and it
is therefore in general different from the average over
the relative velocities seen by each particle. Fundamen
tal differences between cases A, B on the one hand and
cases SA, SB on the other hand can be observed from
the figure. Firstly, the initial mean acceleration due to
gravity is much smaller in the turbulent cases. The dif
ference diminishes with time and after t/tp 7 (12) in
cases A/SA (B/SB) a crossover takes place, after which
the particles settle on average faster in the turbulent en
vironment. The relative difference in settling velocity
between cases A and SA is approximately 4% at times
larger than 20 particle time units; for cases B, SB the rel
ative difference amounts to 1% (cf. figure 5). Although
the longtime difference is relatively small, it is system
Figure 6: The temporal evolution of the relative turbu
lence intensity I = us/Wrel.
atic, and we therefore consider it as physical, especially
between cases A and SA.
Let us first discuss the different short time behavior.
One possible mechanism acting to decrease the mean
settling velocity is the wellknown nonlinear drag ef
fect (Crowe et al. 1998). Although experiments show
considerable scatter, the general view is that turbulent
fluctuations of the velocity seen by the particles together
with the nonlinearity of the drag law lead to an increase
of the mean drag coefficient, and therefore a decrease
in settling velocity. Since the turbulence intensity is ini
tially very strong, the nonlinearity of the drag can be ex
pected to play a, ig.ilik.ii role at short times. However,
the argument hinges on the knowledge of (i) a drag law
with instantaneous validity, and (ii) the precise definition
of a relative velocity (cf. discussion in Bagchi and Bal
achandar 2004). In fact, Bagchi and Balachandar (2004)
did not observe any systematic effect of turbulence upon
the mean drag.
Concerning the longtime behavior, it should be re
marked that the apparent persistence of the effect of the
initial background turbulence is somewhat surprising. It
seems that the small residual turbulent agitation is caus
ing the particles to fall faster than in the pure sedimen
tation cases. The principal mechanism known to lead
to faster settling velocities under turbulent conditions is
the preferential sweeping mechanism (Wang and Maxey
1993), which, however, has thus far only been confirmed
for very small particles. In order to test the possibility
of preferential particle trajectories in the present cases,
conditional sampling of fluid quantities along particle
trajectories needs to be performed.
Figure 6 depicts the ratio between the square root of
the turbulence intensity (ums = vk) and the instan
taneous average settling velocity, viz. I = urm,/Wrel.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(a)
102
10
10
100 101 102
Figure 7: The temporal evolution of the average dis
sipation rate normalized by the initial value 0o, shown
for case A ( ) and case B ( ). The result for
single phase flow is shown in blue.
The quantity I is a measure of the relative intensity
of particleinduced flow structures (i.e. wakes) as com
pared to velocity fluctuations attributed to the turbulent
background flow. Its value in the present cases A and
B is initially very large (in fact it has a singularity at
the origin since particles are at rest), dropping to val
ues below unity by t/t,,f 4 (2) in case A (B). The
difference between the curves corresponding to the two
cases A and B reflects the fact that the respective settling
velocities differ by roughly a factor of two. The tempo
ral decay of the parameter I in cases A and B implies
that these simulations encompass the whole spectrum of
relative turbulence intensities from very strong agitation
(early times) to the limit of purely particleinduced fluc
tuations (long times). In fact, figure 6 shows that the
asymptotic value of the parameter I in both pure sedi
mentation cases SA, SB is approximately 0.08, indepen
dent of the particle Reynolds number. It can also be seen
that at the end of the observation interval discussed here,
cases A, B have not yet reached the regime of pure sedi
mentation, with values of I still between 1.5 and 2 times
larger than the counterparts without initial background
turbulence.
The dissipation rate of turbulent kinetic energy is
shown in figure 7. Visibly, the addition of particles is
causing an increase of energy dissipation. This is ex
pected since each particle introduces additional gradi
ents near its surface, thereby dissipating kinetic energy.
In the present cases, however, the additional dissipation
is offset by additional generation of kinetic energy due
to the work the particles exert on the fluid, as evident
from figure 1 discussed above. Concerning the differ
ence in dissipation between cases A and B, an approxi
4 6 8 10 12
R/(D/2)
(b)
2 4 6 8 10 12
R/(D/2)
Figure 8: The dissipation rate averaged over a spherical
shell between an outer radius R (measured from the par
ticle center) and the particle surface (R D /2), addi
tionally averaged over all particles. Different lines corre
spond to different instants in time: t/7p 0.1;
 t/p 4.6; t/Tp = 16.6. (a) case A,
(b) case B.
mate factor of 8 separates the two curves at large times
(t/Tp > 30). When considering that the average particle
settling velocities of these two cases differ by a factor of
2.2, it appears that the average dissipation rate (when the
background turbulence has already sufficiently decayed)
is roughly proportional to the third power of the settling
velocity.
In order to quantify the location of regions where ad
ditional dissipation is acting, we have averaged the dis
sipation rate over spherical shells between a radial dis
tance p D /2 (the particle surface) and an outer radius
p = R, viz.
/R 27r fr
D, /2 J O J O
dydydp/V,,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
09
08
07
06
05
04
Figure 9: Vorticity magnitude in case A, showing a ver
tical slice at three different instants: t/Tp 0.14, 1.92,
10.35 (from top to bottom).
where Vs = J /2 J Jd dydqydp is the volume of the
shell. In figure 8 the quantity (E)s is shown as a func
tion of the outer shell radius R for different instants in
time in cases A and B. It can be seen that dissipation is
indeed much larger than average near the particle sur
faces, decreasing towards unity with distance from the
particle. The local values increase in time relative to the
boxaverage dissipation since turbulence decays in time,
i.e. the dissipation rate becomes more and more local
ized around the particles. When plotted in loglog scale
(not shown), it becomes obvious that the dissipation rate
decays as a powerlaw R" for short distances. Specifi
cally, the measured radial decay rates for the three curves
in figure 8 are n 0.4, 2.2 and 2.5 (with increasing
time).
Snapshots of the flow field of case B at different in
stants are shown in figure 9. The visualization of the
vorticity magnitude in vertical slices demonstrates how
vortical structures in the neighborhood of the particles
03I
0 10 20 30 40 50
Figure 10: The temporal evolution of the average dis
tance to the nearest neighbor, normalized by the value
for a homogeneous distribution on a regular cubical lat
tice. The straight line ( ) indicates the value for a
random distribution with the same solid volume fraction.
The lower limit of the plot corresponds to the minimum
value of the function (when all particles are in contact
with a neighboring particle).
(located above the particle location) are nearly invisible
at early times, and then progressively start to occupy the
stage. In the final state (t/Tp 10.35) the vorticity field
is visibly dominated by wakes. A further analysis of
the properties of wakes behind collectively sedimenting
particles as well as the implications for the turbulence
spectra would be desirable.
In order to describe the spatial distribution of the dis
persed phase, we determine the average distance be
tween nearest neighbors, defined as follows:
1 Np
Pdmin = in1 i),
N i J=1i
where dij Ixi) xj) is the distance between the
centers of particles i and j. In figure 10 the quantity
dmin is normalized by the value of particles distributed
on a cubical lattice d, = (VQ/N,)1/3 (where VQ
is the volume of the computational domain). It can be
seen that without background turbulence the distance to
the nearest neighbor is somewhat larger than with tur
bulence (approximately 0.65 compared to 0.58). The
turbulent cases A, B exhibit values which are initially
very close to a random distribution, similar to findings
in vertical turbulent channel flow (Uhlmann 2008). At
later times the values of dmin in cases A, B tend slowly
towards the asymptotic value reached in pure sedimenta
tion. This result implies that the initial strong turbulence
maintains an approximately random spatial distribution
dmin
o rr.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
10
10 5 0 5 10
(b)
100
101
10,
10
10 . ..
10 5 0 5 10
Figure 11: Normalized probability density function of
(a) particle velocity, and (b) hydrodynamic forces act
ing on the particles, both for case B. The histogram was
evaluated over the temporal interval t/Tp E [22.2, 26.6].
F; F;  F. For the purpose
of comparison a Gaussian is shown as a dashed line
( ).
in cases A, B, while pure sedimentation is characterized
by a structure with values closer to a homogeneous pat
tern (i.e. with particles spread further apart on average).
Probability density functions (pdfs) of particle veloc
ity and hydrodynamical force components acting on the
particles are shown in figure 11. The pdfs are normal
ized in order to study their shape independently of the
mean value and the standard deviation. The figure shows
data for case B at t/Tp z 24, averaged over a short
interval of AT/Tp 4.4. We observe that all veloc
ity components behave nearly Gaussian (both horizontal
components should be equal by symmetry, which pro
vides a measure of the quality of the statistics). On the
. . . . . . . . . . . . . . .
.. .. .. .. .. .. .. ..
.. .. .. .. .. .. .. ..
.. .. .. .. .. .. .. ..
.. .. .. .. .. .. .. ..
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . .
other hand, the pdfs of the particle forces, which under
this normalization are identical to the particle accelera
tion statistics, exhibit a clear deviation from Gaussianity.
Firstly, the horizontal components are still symmetric (as
they should be due to symmetry), but have longer tails
than the Gaussian reference curve. Secondly, the verti
cal force component is markedly asymmetric with a vis
ible positive skewness. It should be noted that the force
statistics have been filtered such that data from particles
which are instantaneously "in contact" (i.e. those which
are within the range of the artificial repulsion force) have
been removed from the dataset. Now, a positive fluctu
ation F' implies an instantaneous fluctuation which is
directed upward. A positive skewness of F' could be
generated by a nonlinear drag effect. If one assumes
that the drag is quasisteady (i.e. instantaneously given
by the standard drag law), and the characteristic relative
flow velocity has a symmetric pdf, then the nonlinearity
of the velocity/dragforce relation would indeed lead to a
positively skewed pdf for the force. A positive skewness
of the pdf for the vertical force has also been observed in
the case of vertical channel flow (unpublished data from
the simulation of Uhlmann 2008). This persistent result
merits further investigation in the future.
Conclusions
We have simulated the settling of spherical particles at
moderate Reynolds numbers and low solid volume frac
tions. The solid/fluid interfaces were fully resolved by
means of an immersed boundary method. Two spatially
homogeneous configurations were investigated: pure
sedimentation in ambient fluid and sedimentation with
decaying turbulent background flow.
The results show that turbulence initially leads to a
,ig2ilk.ili decrease in the mean settling velocity. For
large times, however, particles settle on average slightly
more rapidly in turbulent flow. Our analysis of the spa
tial distribution of particles has revealed that the parti
cles are more evenly distributed in the cases of pure sed
imentation. This effect is small, but clearly noticeable
from the average distance to the nearest neighbor.
It was found that particle velocity pdfs are close to
a Gaussian function, but hydrodynamic particle forces
possess wider tails. Interestingly, the pdf of the force
component in the vertical direction exhibits a positive
skewness.
The effect of particles upon the fluid turbulence is
only felt at later times, when the initial background
turbulence has already sufficiently decayed. The two
way coupling effect is such that turbulent kinetic en
ergy and average dissipation rate are both enhanced. The
later quantity was shown to be more and more localized
around the particles. Moreover, anisotropy is generated
such that the vertical fluctuations are strongly enhanced.
In the future the current results will first be verified at
a higher spatial resolution. It is planned to investigate
in more detail the effect of the size of the computational
domain upon the results, especially in the vertical di
rection. We will then turn to higher particle Reynolds
numbers, with the purpose of studying possible cluster
formation through wake attraction.
Acknowledgements
The authors thankfully acknowledge the computer re
sources, technical expertise and assistance provided by
the Barcelona Supercomputing Center Centro Nacional
de Supercomputaci6n as well as by the Steinbuch Centre
for Computing (SCC) at KIT. We thank A. Wray for pro
viding the spectral simulation code. A. Dejoan has pro
vided assistance in the adaptation of the spectral code.
References
P. Bagchi and S. Balachandar. Effect of turbulence on
the drag and lift of a particle. Phys. Fluids, 15(11):3496
3513, 2003.
P. Bagchi and S. Balachandar. Response of the wake
of an isolated particle to an isotropic turbulent flow. J.
Fluid Mech., 518:95123, 2004.
S. Balachandar and J.K. Eaton. Turbulent dispersed mul
tiphase flow. Ann. Rev. Fluid Mech., 42:111133, 2010.
R. Clift, J.R. Grace, and M.E. Weber. Bubbles, drops
and particles. Academic Press, 1978.
C. Crowe, M. Sommerfeld, and Y. Tsuji. Multiphase
flows with droplets and particles. CRC Press, 1998.
R. Glowinski, T.W. Pan, T.I. Hesla, and D.D. Joseph.
A distributed Lagrange multiplier/fictitious domain
method for particulate flows. Int. J. Multiphase Flow,
25:755794, 1999.
J. Jim6nez and A.A. Wray. On the characteristics of
vortex filaments in isotropic turbulence. J. Fluid Mech.,
373:255282, 1998.
F. Lucci, A. Ferrante, and S. Elghobashi. Modulation of
isotropic turbulence by particles of Taylorlengthscale
size. to appear in J. Fluid Mech., 2010.
R.N. Parthasarathy and G.M. Faeth. Turbulence mod
ulation in homogeneous dilute particleladen flows. J.
Fluid Mech., 220:485514, 1990.
C.S. Peskin. The immersed boundary method. Acta Nu
merica, 11:479517, 2002.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
C.S. Peskin. Flow patterns around heart valves: A digi
tal computer method for solving the equations of motion.
PhD thesis, Albert Einstein College of Medicine, 1972.
M. Uhlmann. An immersed boundary method with di
rect forcing for the simulation of particulate flows. J.
Comput. Phys., 209(2):448476, 2005a.
M. Uhlmann. New results on the simulation of par
ticulate flows. Technical Report No. 1038, CIEMAT,
Madrid, Spain, 2004. ISSN 11359420.
M. Uhlmann. An improved fluidsolid coupling method
for DNS of particulate flow on a fixed mesh. In M. Som
merfeld, editor, Proc. 11th Workshop TwoPhase Flow
Predictions, Merseburg, Germany, 2005b. Universitit
Halle. ISBN 3860107674.
M. Uhlmann. Direct numerical simulation of sediment
transport in a horizontal channel. Technical Report No.
1088, CIEMAT, Madrid, Spain, 2006a. ISSN 1135
9420.
M. Uhlmann. Experience with DNS of particulate flow
using a variant of the immersed boundary method. In P.
Wesseling, E. Oiate, and J. P6riaux, editors, Proc. EC
COMAS CFD 2006, Egmond aan Zee, The Netherlands,
2006b. TU Delft. ISBN 9090209700.
M. Uhlmann. Interfaceresolved direct numerical simu
lation of vertical particulate channel flow in the turbulent
regime. Phys. Fluids, 20(5):053305, 2008.
R. Verzicco and P. Orlandi. A finitedifference scheme
for threedimensional incompressible flows in cylindri
cal coordinates. J. Comput. Phys., 123:402414, 1996.
L.P. Wang and M.R. Maxey. Settling velocity and
concentration distribution of heavy particles in homoge
neous isotropic turbulence. J. Fluid Mech., 256:2768,
1993.
L. Zeng, S. Balachandar, and F.M. Najjar. Wake re
sponse of a stationary finitesized particle in a turbu
lent channel flow. Int. J. Multiphase Flow, 36:406422,
2010.
