7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Lattice Boltzmann Simulation of Turbulent Flow Laden with FiniteSize Particles
H. Gao* and L.P. Wang*
Department of Mechanical Engineering, University of Delaware, Newark, DE 19716, USA
hgao@udel.edu and lwang@udel.edu
Keywords: Lattice Boltzmann simulations, particleladen flow, turbulent flow, finitesize effects
Abstract
The paper describes a particleresolved simulation method for turbulent flow seeded with finite size particles. The
method is based on the multiplerelaxationtime lattice Boltzmann equation. The noslip boundary condition on the
moving particle boundaries is handled by a secondorder interpolated bounceback scheme. The populations at a
newly converted fluid node are constructed by the equilibrium distribution with nonequilibrium corrections. The
code is parallelized with MPI and is found to be computationally efficient. The method is applied to two problems.
The first is the sedimentation of a random suspension at finite particle Reynolds numbers. Results on the hindered
settling function and relative particle vertical velocity fluctuation are presented in terms of particulate volume fraction
and particle Reynolds numbers, and are found to be in reasonable agreement with previous results based on the force
coupling method (Climent & Maxey 2003). The second problem is a decaying isotropic turbulence seeded with
particles of Taylor microscale size. Turbulence statistics for three particleladen cases are simulated and compared to
those of particlefree case. We find that, at the given volume fraction, the mass loading and particle size play a major
role in the evolution of overall kinetic energy, energy spectrum, and the flow dissipation rate. The density ratio plays
a minor role. These results agree well with those obtained by Lucci et al. (2010) based on the immersed boundary
method.
Introduction
Turbulent flows laden with solid particles are relevant
to a wide variety of natural and engineering applica
tions, such as aerosol and pollutant transport, interac
tion of cloud droplets, spay combustion, and chemical
processes. Over the past two decades, advanced compu
tational and experimental methods have been developed
and utilized to quantify the interactions between the dis
persed solid particles and the carrier fluid phase. Since
turbulent particleladen flows usually encompass a wide
range of length and time scales, it is computationally
demanding to simultaneously resolve both the carrier
phase turbulent flow and the disturbance flows due to
the particles. Earlier studies typically employed the
pointparticle model (Squires and Eaton 1990; Wang &
Maxey 1993; Elghobashi & Truesdell 1993; Sundaram
& Collins 1997; Zhou et al. 2001; Ayala et al. 2007),
namely, the particle size is small compared with Kol
mogorov scale of the carrier flow. Within the point
particle model, the motion of an individual particle is
modeled by an equation of motion (Maxey & Riley
1983) and only the scales of motion that are much larger
than the particle size are directly resolved. The point
particle model has been used by many to study preferen
tial concentration of inertial particles (Squires and Eaton
1991; Wang & Maxey 1993), turbulent modulation by
inertial particles (Squires and Eaton 1990; Elghobashi
& Truesdell 1993), turbulent collision rate of inertial
particles (Sundaram & Collins 1997; Zhou et al. 2001;
Ayala et al. 2007). The pointparticle model is compu
tationally efficient, but its validity is often questionable
when particle mass loading is significant or particles
form aggregates yielding strong multiscale couplings
between the particulate phase and the carrier phase.
In many engineering applications, however, particle
size can be of the same order as or larger than Kol
mogorov scale, the scales contained in the disturbance
flows then overlap with the scales of motion in the car
rier turbulence. In this situation, the pointparticle model
is no longer a valid description and the finitesize effect
of the dispersed phase must be resolved togetherwith the
carrier fluid turbulence. Furthermore, in concentrated
suspensions such as fluidized beds, particleparticle col
lisions can also have a significant impact on the rheology
of the system.
The motivation of this work is to understand the mo
tion and hydrodynamic interactions of finitesize inertial
particles suspended in a turbulent flow. One of the ma
jor challenges is to develop an efficient and accurate ap
proach to resolve the disturbance flows around particles
suspended in a turbulent carrier fluid.
In recent years, several computational methods have
been proposed towards this goal. Finite element meth
ods (Hu 1996; Hu et al. 2001) employ bodyfitted mesh
to implement solid particle boundary conditions. The
frequent mesh regeneration due to the geometry change
as a result of particle motion is computationally expen
sive, especially in three dimensional simulations. There
fore, methods using fixed and structured grid have re
ceived particular attention more recently. The fictitious
domain method (Patankar et al. 2000; Glowinski et al.
2001) applies a field of Lagrange multipliers to enforce
constraints on the particle so that the fluid inside parti
cle domain is forced to mimic rigid body motion. The
immersed boundary method (Peskin 2002; Uhlmann
2005, 2008) realizes the noslip boundary condition on
particle surface by imposing a localized forcing field.
Similarly, the force coupling method (Maxey & Patel
2001) uses smoothed body force field to represent the
effects of particles on the fluid phase, in terms of low
order force multiple expansion. The hybrid method
Physalis (Takagi et al. 2003; Zhang & Prosperetti 2003,
2005) handles the noslip boundary condition by cou
pling an analytical Stokes expansion valid in a narrow
but finite region near a particle surface with the numeri
cal solution outside the particle.
In contrast to the above macroscopic CFD approaches
based on the NavierStokes equations, in the mesoscopic
lattice Boltzmann method (LBM) (Aidun et al. 1998;
Ladd 1994a,b; Ten Cate et al. 2004) the fluid field is
realized through local moments of a lattice Boltzmann
equation on a uniform lattice grid. The noslip boundary
condition can be imposed by using a simple interpolated
bounceback scheme. To reduce force oscillations on
the particles, the immersedboundarylatticeBoltzmann
method (IBLBM) (Feng & Michaelides 2004, 2005)
has also been developed by replacing the conventional
bounceback scheme with a direct forcing scheme ap
plied on a set of Lagrangian boundary points represent
ing particle surfaces.
Among all these particle resolved methods, few have
been successfully applied to the particleresolved sim
ulation of turbulent particulate suspensions. Ten Cate
et al. (i2 l4) conducted a particleresolved simulation
of forced turbulent flows seeded with solid particles via
LBM. The carrierfluid turbulence is maintained at Tay
lor microscale Reynolds number Re = 61 using a
spectral forcing scheme. Noslip boundary condition on
the particle surface is implemented by applying a body
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
force field to the fluid domain, which enforces the fluid
velocity near particle surface to be the same as the solid
surface velocity. Similar to Nguyen & Ladd (2002),
a subgrid lubrication force is computed based on rela
tive location and velocities of approaching particles. The
simulation contained up to 3,868 particles, correspond
ing to volume fractions ranging from _'. to 11 '. The
density ratio between the solid and fluid phase ranges
from 1.15 to 1.73. The particle radius is set to 4 grid
spacing in the simulations. It was demonstrated that
particlefluid relative motion results in an enhancement
of kinetic energy E(k) and energy dissipation rate E(k)
at large wave numbers and a reduction of both spectra at
small wave numbers. Oscillations in E(k) were found at
large wave numbers, although they are negligible when
compared to the total kinetic energy of the system.
Zhang & Prosperetti (2005) developed the hybrid
method Physalis and demonstrated its capability in sim
ulating decaying particleladen turbulence. A periodic
cubic domain of 643 contains 100 spherical particles
with radius equal to 39% of the initial Taylor microscale.
Particlefluid density ratio was set to 1.02. Gravity was
neglected and particleparticle short range interactions
were represented by an elastic model. The Taylor mi
croscale Reynolds number changed from 29 to 14 during
the time interval of the simulation. It was reported that,
compared with the pointparticle model, the finitesize
particles decrease the kinetic energy, show less diffusion
in terms of mean particle displacement, and exhibit a
stronger tendency of clustering.
Uhlmann (2008) considered particulate suspensions
in a vertical turbulent channel flow, using the immersed
boundary method to treat the solid particles. The in
compressible NavierStokes equations were solved by a
fractionalstep scheme on a staggered grid. Up to 4,096
particles were simulated in a turbulent flow sustained at a
bulk flow Reynolds number of 2,700. The particle diam
eter was set to approximately 11 wall units, correspond
ing to a volume fraction of 0. I '. The density ratio
varied from 2.2 to 10. It was found that the presence of
particles induced large scale streaklike velocity pertur
bations, although no significant aggregation of particles
was observed.
Recently, the immersed boundary method was also
used by Lucci et al. (2010) to study modulation of de
caying turbulence by solid particles of size comparable
to the Taylor microscale. The carrier flow had an ini
tial Taylor Reynolds number of 75. Density ratio varied
from 2.56 to 10.0. Volume fraction of dispersed phase
was 11i' or less, and a maximum number of 6,400 par
ticles were considered in their simulations.
Using the force coupling method, Yeo et al. (2010)
studied turbulence modulation by finitesize particles
and bubbles. The flow was forced at large scales. They
demonstrated that the pivoting wavenumber characteriz
ing the transition from damped to enhanced energy con
tent is mainly a finitesize effect of particles and that the
transition scale is almost independent of the particle to
fluid density ratio.
The main objectives of this paper are to describe our
own implementation of LBM using interpolated bounce
back and to validate our preliminary results. The sim
ulation method is described in Section 2. In Section 3,
we first apply the method to study statistics in a random
suspension at finite particle Reynolds numbers and com
pare our results with those reported in Climent & Maxey
(2003). Then in Section 4, we consider decaying tur
bulent flows seeded with finitesize particles. We show
that our method can reproduce the results in Lucci et al.
(2010). We conclude the paper with a summary and fu
ture outlook.
Simulation Method
In this study, the mesoscopic lattice Boltzmann ap
proach, based on the multiplerelaxationtime (MRT)
lattice Boltzmann equation, is applied to simulate de
caying homogeneous isotropic turbulence seeded with
finite size spherical particles in a threedimensional pe
riodic domain. In the MRTLBE, the particle distribution
function is governed by
f(x + ei6t, t + 6t)
f(x,t)
M 1 S m m(eq) (1)
where M is an orthogonal transformation matrix con
verting the distribution function f from discrete velocity
space to the moment space m, in which the collision re
laxation is performed. The transformation between the
particle velocity space and the moment space is given as
m M f, f M1 m. (2)
The diagonal relaxation matrix S specifies the relaxation
rates for the nonconserved moments. In this work, the
D3Q19 model is utilized with the discrete velocities or
dered as
(0,0, 0), i =0,
e= ( o, 0), (o, 1, o), (0, o, i), i 16,
( 1, 0), (l, 0, 1), (0, 1, 1), i = 718.
(3)
The corresponding 19 x 19 transform matrix M, the 19
components of moment m and its equilibrium counter
part m(eq), and the relaxation matrix S are described in
detail in d'Humi&res et al. (2002).
The macroscopic variables, including density, mo
mentum, and pressure, are obtained from the moments
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
of the mesoscopic distribution function f, namely,
p = po + p, po = 1
8p = fi, pou = fiee, p = 6pca (5)
i i
where u is the macroscopic fluid velocity, and the sound
speed c, is equal to 1//3 in lattice units.
The computation domain is covered with a uniform
cubic lattice. The periodic boundary condition is ap
plied in all three directions. The noslip boundary con
dition on solid particle surfaces are implemented using
a secondorder interpolated bounceback scheme. First,
for each fluid node near a particle surface, all the links
moving into the surface are identified, and the boundary
cutting location on the link is calculated in terms of the
percentage (a) of the link outside the surface. Before the
advection step, the missing population is properly inter
polated using a and two populations lying before and
after the path of the missing population (Lallemand &
Luo 2003). For a < 0.5, the interpolation is performed
before streaming, while for a > 0.5, it is done after
streaming. For the moving boundary, this procedure is
performed for each time step. When two particles are in
close contact, two fluid lattice nodes may not be avail
able near the missing population to allow for a second
order interpolation. In this case, a linear interpolation or
a simple onenode bounce back is used instead.
It is noted that the above LBE evolution only applies
to the fluid nodes where the particle distribution func
tions are defined. No distribution function is defined for
solid lattice nodes lying within the particles. As the par
ticles are moving, a solid node may move out of the solid
region and become a fluid node with unknown distri
bution functions. In this work, the missing distribution
functions for the new fluid node are constructed by an
equilibrium distribution plus a nonequilibrium correc
tion (Caiazzo 2008), namely,
fi(x) f q)(x; U, p) + f (x + ea)
where u, is the velocity of the moving boundary, the
magnitude of which is assumed to be much smaller than
the sound speed c,, p is the fluid density averaged over
all the available fluid nodes surrounding the new fluid
node located at x. While the equilibrium part is calcu
lated with u, and p at x, the nonequilibrium part is
obtained from a neighboring node at x + ea. Here the
ea is a specified discrete velocity along which direction
the quantity ft ea takes the maximum, where n the unit
normal vector pointing outwards of the moving bound
ary at the point through which the previous solid node
moves to fluid domain. The equilibrium population is
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
calculated as
feq) (x, t) = p p u uu : (eiei c2)1
+ Cs2 2Cs4
(7)
where I is the identity matrix, and the weight W, is given
1/3,
wi 1/18,
1/36,
1 6,
7 18.
(8)
The hydrodynamic force and torque acting on a par
ticle are calculated during the interpolated bounceback
procedure by summing up momentum changes on all the
links across a particle surface, and the cross production
of momentum change and a vector from particle center
to the surface point, respectively. In case particles are
in close contact, where the short range particleparticle
interactions cannot be resolved, the following pairwise
repulsive force between particles is added (Glowinski et
al. 2001; Feng & Michaelides 2005)
S0, (
rij > Rij +
C I (rij ) 'Ij +c
(9)
whererij = YiYj, rij = Y Yj, Rij = Ri+Rj,
cij is a force scale and is set to be the buoyancy force
in this study. Y, and Yj represent the center location
of the i'th and j'th particle with radius of R, and Rj,
respectively. ( is a threshold gap distance within which
the lubrication force becomes active. In this work, ( is
set as two lattice spacing. The stiffness parameter, Ep, is
set to be small enough to avoid particle overlapping.
At each time step, with the hydrodynamic force lu
brication force, and torque acting on particles readily ob
tained, the particle translational and angular velocities,
center position and angular displacement are updated us
ing the CrankNicolson scheme,
Vt 1 f (F + Ft
2p \1 i' i
Q*t + 1 (rt + rt
o + 2 (Qt + Q
iY ( v:
e:+ (2+
t)^
where Mp and Ip are the mass and moment of inertia of
the i'th particle, and time step size At 1 .
Results: Sedimentation of a
Random Suspension
To validate our implementation of LBM, we first con
sider a random suspension at finite particle Reynolds
S1.0
0.8
0.6
0.4
0.00 0.03 0.06 0.09
volume fraction
Figure 1: Average settling velocity normalized by Vo
in 1283 simulations: comparing the present study with
Climent & Maxey (2003).
1.0
I 
S0.8 
0.6 
. 0.4
oI 
0.2
0.0
0.
I I
T I I I
.00 0.03 0.06 0.09
volume fraction
Figure 2: Fluctuation of settling velocity normalized by
its mean value, Vrms/Vmean, in 1283 simulations: com
paring the present study with Climent & Maxey (2003).
numbers and compare our results with those of Climent
& Maxey (2003). Solid spherical particles were ran
domly seeded in an initially quiescent fluid domain of
size 1283. Periodic boundary conditions were applied
in all three directions. A series of simulations were
conducted to yield five different volume concentrations,
using 1, 40, 120, 240, 480 particles of a fixed radius
of 5 grid spacing. The singleparticle case was used
as the reference case. A pressure difference was ap
plied to balance the excess weight of the particles. The
gravitational acceleration was varied to give two particle
Reynolds numbers of roughly 4.1 and 11.2, based on the
single particle sedimentation velocity and particle diam
eter. The sedimentation velocity was corrected to ensure
that the system has a zero net vertical mass flux.
S Our data Re=11.2
Our data Re=4.1
CM03 Re=10
S* CM03 Re=5
CM03 Re=l
I
I
Our data Re=11.2
Our data Re=4.1
CM03 Re=10
CM03 Re=5
CM03 Re=l
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 1: Particles parameters at release time.
Case d pp/pf Np d/r d/A 4,v m Tp/Tk
1 0 0 0
2 8.0 2.56 6400 16.1 1.2 0.10 0.28 22.5
3 8.0 5.0 6400 16.1 1.2 0.10 0.56 57.6
4 11.0 2.56 2304 22.1 1.5 0.10 0.27 42.5
An important question for a random suspension is
how the mean velocity varies with the particulate vol
ume concentration. Figure 1 shows the mean sedimenta
tion velocity normalized by Vo, the terminal velocity of a
single particle sedimenting in the same periodic domain.
This is known as the hindered settling function and it is
plotted as a function of particulate volume fraction. Also
shown are results from Climent & Maxey (2003) at par
ticle Reynolds numbers of 1, 5, and 10. Clearly, the par
ticle average settling velocity is significantly reduced as
its volume fraction is increased, and the larger the par
ticle Reynolds number the smaller the settling velocity.
The results are in good agreement with those of Climent
& Maxey (2003), with our simulations predict a some
what smaller settling velocity. The differences could be
due to two reasons. First, there are statistical uncertain
ties in both studies. Second, the force coupling method
is an approximate method in which the disturbance flow
is not fully resolved.
The relative vertical velocity fluctuations
Vrms/Vmean are shown in Fig. 2. The overall
trends are similar: the relative vertical velocity fluctua
tion increases with the volume fraction, but decreases
with the particle Reynolds number. Quantitatively, our
simulations yield smaller values. The origin for the
quantitative differences between our simulations and
those of Climent & Maxey (2003) remains to be studied.
Results: Decaying Turbulence
Seeded with FiniteSize Particles
Prior to the simulation of particleladen turbulence, the
accuracy our LBM code was verified by a comparison
with a pseudospectral code for particlefree decaying
homogeneous isotropic turbulence. Almost identical ki
netic energy and dissipation rate spectra, as well as their
time evolution were observed from the initial time to
about 2.12 eddy turnover times (2.12Te,o). This is con
sistent with the results shown in Peng et al. (2010).
Following the work of Lucci et al. (2010), the initial
velocity field at t = 0 was specified by a Gaussian field
OA
o U
oaY
0"8
on
0.82
0.61
0.55
z
0.23
0.E7
0.10
0.850
Figure 3: Vorticity contour at (a) 3000 time step
(1.27T6,o) and (b) at 5000 time step (2.12T,o) for case
4. A layer of fluid is removed to show some of the 2304
particles.
with a prescribed kinetic energy spectrum as
E(k) ( (exp
2 k} 1 kP
(14)
(b) 250
20C
(a) 250
200
150
100
50
0
(c) 250
200
150
(d) 250
200
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
vortcity
0.i0
0.74
0.61
50 100 10 20 2 7
xI
0.55
0.48
0 42
0.35
0 29
0.23
0.16
0.10
50 100 150 20 250
x
vorti
0.40
0.74
0.61
055
048
0.42
0.35
0.29
023
0.16
010
vorticity
0.80
0.74
067
0.61
0.55
048
0.42
0.35
0.29
023
0.16
0.10
Figure 4: Vorticity contour, velocity vector, and particle location at 3000 time step (1.27Te,o) on a planecut of
z 128.5 in 2563 simulation: (a) case 1, (b) case2, (c) case 3, and (d) case 4.
where k is the wave number, kp is the wave number con
taining peak energy, and uo is the initial r.m.s. velocity.
In this paper, we use a periodic domain of 2563 grid res
olution, indicating the potential range of wave number
1 < k < 128. Our simulations match the conditions
in Lucci et al. (2010): k, 4, uo 0.0503, fluid kine
matic viscosity v 1.45 x 104, yielding an initial Tay
lor microscale Reynolds number Rexo 78, energy dis
sipation rate co 1.07 x 104, and Kolmogorov length
scale 1.30 x 10 2. Similar to Peng et al. (2010),
after initialization of velocity field, a prerelaxation pro
cedure of the density field was performed to construct a
consistent pressure field based on the prescribed velocity
field; this amounts to solving the pressure Poisson equa
tion. The moment is referred to as time t = 0. The flow
field was then evolved for sometime until a converged
skewness of about 0.50 was developed. At this mo
ment, particles were released in the domain at random
locations with a gap distance of at least two grid spacing
in between. Initial translational velocities of particles are
set as the instantaneous fluid velocity at particle center,
and angular velocity is zero. A set of four simulations,
listed in Table 1, were preformed to study the turbulence
modulationby finitesize particles. These roughly match
Cases A, D, E, and G in Lucci et al. (2010). The grav
ity was set to zero. The simulations are run from t 0
to t 2.12Te,o, corresponding to 5000 time steps. The
flow fields are well resolved as shown by kmax, > 1
for the whole time interval. All the simulations were
performed on an IBM Power6 supercomputer at NCAR
vorticity
0.80
0.74
0.67
0.61
0.55
0.48
0.42
0.35
0.29
023
0.16
0.10
1.0
l9 ce
 case 4
0.8 case 3
S case 2
case 1
0.6
W 0.4
0.2
0.0
0.0 0.4 0.8 1.2 1.6 2.0
t / To
Figure 5: Temporal evolution of turbulent kinetic en
ergy normalized by its initial value.
with MPI parallelization using 32 processors. Each sim
ulation takes less than 6 hours of wall clock time.
Figure 3 displays the vorticity contour for Case 4
on the boundaries of the computational domain at 3000
time step (1.27T,0o) and 5000 time step (2.12T,,o), re
spectively. A layer of fluid to the top is removed to show
the locations of a portion of the 2304 particles. The mag
nitude of vorticity decreases with time as expected, as
can be seen from the color difference of Fig. 3(a) and
Fig. 3(b).
In Figure 4 we provide 2D visualizations of the vor
ticity magnitude, projected velocity vectors and particle
locations, on an x y plane near the center cut of the
domain (z 128.5). It is observed that the presence
of particles is often associated with high voritcity values
(the red spots), especially in Fig. 4(c), suggesting that
motion of finite size particles can indeed generate small
scale flow structures near its surface. Note that Case 3
has the largest density ratio and likely the largest slip
velocity. The enhanced activities at small scales due to
particles are evident when compared to the particlefree
case, Case 1.
Figure 5 shows the temporal evolution of turbulent ki
netic energy (TKE) scaled by its initial value Eo com
puted at time t = 0 for all the four cases. Compared with
the particlefree turbulence, the TKE decreases faster in
the particleladen cases from the moment of particle re
lease till about 1.6T6,o, and then the decaying rate tends
to be the same in all cases. In addition, the increase of
particle number results in an evident reduction of TKE,
as shown by comparing the curves from Case 2 and Case
4. The increase of particle density has a marginal effect
on decaying rate reduction, based on a comparison of
Case 2 with Case 3. These results are very similar to
those shown in Lucci et al. (2010).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
1.5  case 4
case 3
., case 2
1.2 / case 1
0.9 
0.6
0.3
0.0
0.0 0.4 0.8 1.2 1.6 2.0
t / Te,o
Figure 6: Temporal evolution of turbulent dissipation
rate normalized by its initial value.
IIF.
Figure 7: Kinetic energy spectrum at 5000 time step
(2.12T,,o).
Figure 6 shows the modification of turbulent dissipa
tion rate due to the presence of particles. At the mo
ment of about 0.26T,o, the injection of particles into
the fluid domain causes a sudden jump of the dissipation
rate, followed by a peak value at about 0.4T,o, and then
the total dissipation decreases monotonically. Particles
with the smallest Stokes number St Tp/Tk (case 2)
produce smallest dissipation rate eventually. With the
largest mass loading, case 3 introduces the highest peak
value and consequently takes the longest time to reduce
its dissipation rate comparable to that of the particlefree
case.
Finally, in Figure 7 we compare the kinetic energy
spectra at the end of simulations t 2.12Te,o. It should
be noted that we use the full velocity field including the
rigidbody velocity inside the particles to compute the
energy spectrum. The spectrum is essentially identical
for all cases for wave numbers less than 4. In range of
case 4
case 3
 case 2
 case 1
4 < k < 20, the particleladen cases feature less TKE
than the particlefree case, and the contrary trend occurs
at k > 40. Interestingly, wave like oscillations are ob
served at the tails of the spectra for turbulent particulate
suspensions. At large wavenumbers, Case 3 possesses
noticeably larger TKE than the other two particleladen
cases due to large mass loading and large Stokes num
ber.
Conclusions
In this paper, we have presented a particleresolved sim
ulation method for turbulent flow seeded with finite
size particles. The method was based on the multiple
relaxationtime lattice Boltzmann equation (d'Humibres
et al. 2002). The noslip boundary condition on the
moving particles boundaries was handled by a second
order interpolated bounceback scheme (Lallemand &
Luo 2003). The populations at a new fluid node
were constructed by equilibrium distribution with non
equilibrium correction (Caiazzo 2008). A simple lubri
cation force model was utilized to avoid particleparticle
overlap. The code was parallelized with MPI and was
found to be computationally efficient. So far, up to 6,400
particles in 3D have been considered in our simulations.
We are in the process to testing the code for larger grid
resolution and larger number of particles.
The force and torque on a particle were computed
based on summing momentum exchanges at the bound
ary links. Apart from the refill problem, the method con
serves the overall momentum of the system at the meso
scopic level, which will likely ensure the correct balance
of kinetic energy of the whole system and viscous dissi
pation.
The method was then applied to two problems. The
first is the sedimentation of a random suspension at finite
particle Reynolds numbers. Results on the hindered set
tling function and relative particle vertical velocity fluc
tuation are presented in terms of particulate volume frac
tion and particle Reynolds numbers. These results are in
reasonable agreement with the results obtained by Cli
ment & Maxey (2003) using a force coupling method.
We then applied the method to study a decaying
turbulence seeded with finitesize particles. In the
particlefree case, the flow statistics matched precisely
the results obtained from the accurate pseudospectral
method, consistent with the observation of Peng et al.
(2010). Three particleseeded cases reported in Lucci
et al. (2010) were simulated, with very similar results.
The results show that particles of Taylor microscale size
introduce smallscale features, enhance the dissipation
rate of the system, and alter the shape of energy spec
trum. At a given particulate volume fraction, the mass
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
loading and particle size play a prominent role in the
evolution of overall kinetic energy, energy spectrum, and
the flow dissipation rate. The density ratio plays a less
important role. These results are very preliminary. Fur
ther analysis and parametric study are underway.
Further analysis of the simulation results will focus on
particle phase characteristics, including fluid and parti
cle radial and transverse relative velocities, particle ra
dial distribution function, and particle collision rates.
Towards this end, the numerical method will be criti
cally examined in order to improve the representation
of shortrange particleparticle interaction. Other funda
mental numerical issues related to moving particles such
as force oscillations and the refill problem also require
further research.
Acknowledgements
This work was supported by the National Science
Foundation (NSF) under grants OCI0904534, ATM
0527140, and ATM0730766, and National Natural Sci
ence Foundation of China (Project No. 10628206).
Computing resources are provided by National Center
for Atmospheric Research through CISL35751010 and
CISL35751014.
References
Aidun C.K., Lu Y, and Ding E.J., Direct analysis of par
ticulate suspensions with inertia using the discrete Boltz
mann equation, J. Fluid Mech., Vol. 373, pp. 287311,
1998
Ayala O., Grabowski W.W, and Wang L.P., A
hybrid approach for simulating turbulent collisions
of hydrodynamicallyinteracting particles, J. Comput.
Phys., Vol. 225, pp. 5173, 2007
Caiazzo A., Analysis of lattice Boltzmann nodes initial
isation in moving boundary problems, Prog. Comput.
Fluid Dyn., Vol. 8, pp. 310, 2008
Climent E., and Maxey M.R., Numerical simulations of
random suspensions at finite Reynolds numbers, Int. J.
Multiphase Flow, Vol. 29, pp. 579601, 2003
d'Humi&res D., Ginzburg I., Krafczyk M., Lallemand
P., and Luo L.S., Multiplerelaxationtime lattice Boltz
mann models in threedimensions, Phil. Trans. R. Soc.
Lond. A, Vol. 360, pp. 437451, 2002
Elghobashi S., and Truesdell G., On the twoway inter
action between homogeneous turbulence and dispersed
solid particles. I: Turbulence modification, Phys Fluids,
Vol. 5, pp. 17901801, 1993
Feng Z.G. and Michaelides E.E., The immersed
boundarylattice Boltzmann method for solving fluid
particles interaction problems, J. Comput. Phys., Vol.
195, pp. 602628, 2004
Feng Z.G. and Michaelides E.E., Proteus: a direct forc
ing method in the simulations of particulate flows, J.
Comput. Phys., Vol. 202, pp. 2051, 2005
Glowinski R., Pan T., Hesla T, Joseph D. and P6riaux
J., A fictitious domain approach to the direct numerical
simulation of incompressible viscous flow passed mov
ing rigid bodies: application to particulate flow, J. Com
put. Phys., Vol. 169, pp. 363426, 2001
Hu H., Direct simulation of solidliquid mixtures, Int. J.
Multiphase Flow, Vol. 22, pp. 335352, 1996
Hu H., Patankar N.A., and Zhu M.Y, Direct numerical
simulations of solidliquid systems using the aribitrary
LagrangianEulerian technique, J. Comput. Phys., Vol.
169, pp. 427462, 2001
Ladd A.J.C., Numerical simulations of particulate sus
pensions via a discretized Boltzmann equation. Part 1.
Theoretical foundation, J. Fluid Mech., Vol. 271, pp.
285309, 1994a
Ladd A.J.C., Numerical simulations of particulate sus
pensions via a discretized Boltzmann equation. Part 2.
Numerical results, J. Fluid Mech., Vol. 271, pp. 311339,
1994a
Lallemand P. and Luo L.S., Lattice Boltzmann method
for moving boundaries, J. Comput. Phys., Vol. 184, pp.
406421, 2003
Lucci F, Derrante A., and Elghobashi S., Modulation
of isotropic turbulence by particles of Taylorlengthscale
size, J. Fluid Mech., in press, 2010
Maxey M.R. and Patel B.K., Localized force represen
tations for particles sedimenting in Stokes flows, Int. J.
Multiphase Flow, Vol. 27, pp. 16031626, 2001
Maxey M.R. and Riley J.J., Equation of motion for a
small rigid sphere in a nonuniform flow, Phys. Fluids,
Vol. 26, pp. 883889, 1983
Nguyen N.Q. and Ladd A.J.C., Lubrication corrections
for lattice Boltzmann simulations of particle suspen
sions, Phys. Rev. E, Vol. 66, 046708, 2002
Patankar N.A., Singh P, Joseph D., Glowinski R., and
Pan T.W, A new formulation of the distributed La
grange multiplier/fictitious domain method for particu
late flows, Int. J. Multiphase Flow, Vol. 26, pp. 1509
1524, 2000
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Peng Y, Liao W, Luo L.S., and Wang L.P, Compari
son of the lattice Boltzmann and pseudospectral meth
ods for decaying turbulence. Part I. Loworder statistics,
Comput. Fluids, Vol. 39, pp. 568591, 2010
Peskin C., The immersed boundary method, Acta Nu
merica, Vol. 11, pp. 479517, 2002
Squires K.D. and Eaton J.K., Particle response and tur
bulence modification in isotropic turbulence, Phys. Flu
ids A, Vol. 2, pp. 11911203, 1990
Squires K.D. and Eaton J.K., Preferential concentration
of particles by turbulence, Phys. Fluids A, Vol. 3, pp.
11691179, 1991
Sundaram S. and Collins L.R., Collision Statistics in an
isotropic particleladen turbulent suspension. Part 1. Di
rect numerical simulations, J. Fluid Mech., Vol. 335, pp.
75109, 1997
Ten Cate A., Derksen J.J., Portela L.M., van den
Akker H.E.A., Fully resolved simulations of colliding
monodisperse spheres in forced isotropic turbulence, J.
Fluid Mech., Vol. 519, pp. 233271, 2004
Takagi S., Oguz H., Zhang Z., and Prosperetti A.,
Physalis: A new method for particle simulation. Part II:
Twodimensional NavierStokes flow around cylinders,
J. Comput. Phys., Vol. 187, pp. 371390, 2003
Uhlmann M., An immersed boundary method with di
rect forcing for the simulation of particulate flows, J.
Comput. Phys., Vol. 209, pp. 448476, 2005
Uhlmann M., Interfaceresolved direct numerical simu
lation of vertical particulate channel flow in the turbulent
regime, Phys. Fluids, Vol. 20, 053305, 2008
Wang L.P., Maxey M.R., Settleing velocity and concen
tration distribution of heavy particles in homodeneous
isotropic turbulence, J. Fluid Mech., Vol. 256, pp. 27
68, 1993
Wang L.P., Rosa B., Gao H., He G.W, Jin G.D., Turbu
lent collision of inertial particles: pointparticle based,
hybrid simulations and beyond, Int. J. Multiphase Flow,
Vol. 35, pp. 854867, 2009
Yeo K., Dong S., Climent E., and Maxey M.R., Modula
tion of homogeneous turbulence seeded with finite size
bubbles or particles, Int. J. Multiphase Flow, Vol. 36, pp.
221233, 2010
Zhang Z., and Prosperetti A., A method for particle sim
ulation, J. Appl. Mech. Trans. ASME, Vol. 70, pp. 64
74, 2003
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Zhang Z., and Prosperetti A., A secondorder method for
threedimensional particle simulation, J. Comput. Phys.,
Vol. 210, pp. 292324, 2005
Zhou Y, Wexler A.S., Wang L.P., Modelling turbulent
collision of bidisperse inertial particles, J. Fluid Mech.,
Vol. 433, pp. 77104, 2001
