7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30- June 4, 2010


Lattice Boltzmann Simulation of Turbulent Flow Laden with Finite-Size 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, particle-laden flow, turbulent flow, finite-size effects




Abstract

The paper describes a particle-resolved simulation method for turbulent flow seeded with finite size particles. The
method is based on the multiple-relaxation-time lattice Boltzmann equation. The no-slip boundary condition on the
moving particle boundaries is handled by a second-order interpolated bounce-back scheme. The populations at a
newly converted fluid node are constructed by the equilibrium distribution with non-equilibrium 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 particle-laden cases are simulated and compared to
those of particle-free 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 particle-laden 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
point-particle 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 point-particle 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 point-particle model
is no longer a valid description and the finite-size effect
of the dispersed phase must be resolved togetherwith the
carrier fluid turbulence. Furthermore, in concentrated
suspensions such as fluidized beds, particle-particle 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 finite-size 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 body-fitted 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 no-slip 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 no-slip 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 Navier-Stokes 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 no-slip boundary
condition can be imposed by using a simple interpolated
bounce-back scheme. To reduce force oscillations on
the particles, the immersed-boundary-lattice-Boltzmann
method (IB-LBM) (Feng & Michaelides 2004, 2005)
has also been developed by replacing the conventional
bounce-back 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 particle-resolved sim-
ulation of turbulent particulate suspensions. Ten Cate
et al. (i2 l4) conducted a particle-resolved simulation
of forced turbulent flows seeded with solid particles via
LBM. The carrier-fluid turbulence is maintained at Tay-
lor microscale Reynolds number Re = 61 using a
spectral forcing scheme. No-slip 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
particle-fluid 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 particle-laden turbulence. A periodic
cubic domain of 643 contains 100 spherical particles
with radius equal to 39% of the initial Taylor microscale.
Particle-fluid density ratio was set to 1.02. Gravity was
neglected and particle-particle 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 point-particle model, the finite-size
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 Navier-Stokes equations were solved by a
fractional-step 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 finite-size 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 finite-size 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 finite-size 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 multiple-relaxation-time (MRT)
lattice Boltzmann equation, is applied to simulate de-
caying homogeneous isotropic turbulence seeded with
finite size spherical particles in a three-dimensional pe-
riodic domain. In the MRT-LBE, 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 M-1 m. (2)

The diagonal relaxation matrix S specifies the relaxation
rates for the non-conserved 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 1-6,
( 1, 0), (l, 0, 1), (0, 1, 1), i = 7-18.
(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 no-slip boundary con-
dition on solid particle surfaces are implemented using
a second-order interpolated bounce-back 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 one-node 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 non-equilibrium 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 non-equilibrium 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 bounce-back
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 particle-particle
interactions cannot be resolved, the following pair-wise
repulsive force between particles is added (Glowinski et
al. 2001; Feng & Michaelides 2005)


S0, (


rij > Rij +

C I (rij ) 'Ij +c


(9)
whererij =- Yi-Yj, 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 Crank-Nicolson 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 single-particle 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 Finite-Size Particles


Prior to the simulation of particle-laden turbulence, the
accuracy our LBM code was verified by a comparison
with a pseudo-spectral code for particle-free 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 plane-cut 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 10-4, yielding an initial Tay-
lor microscale Reynolds number Rexo 78, energy dis-
sipation rate co 1.07 x 10-4, and Kolmogorov length
scale 1.30 x 10 2. Similar to Peng et al. (2010),
after initialization of velocity field, a pre-relaxation 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 finite-size 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- ----------------------c-e
------- 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 particle-free
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 particle-free turbulence, the TKE decreases faster in
the particle-laden 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 particle-free
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
rigid-body 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 particle-laden cases feature less TKE
than the particle-free 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 particle-laden
cases due to large mass loading and large Stokes num-
ber.


Conclusions

In this paper, we have presented a particle-resolved sim-
ulation method for turbulent flow seeded with finite
size particles. The method was based on the multiple-
relaxation-time lattice Boltzmann equation (d'Humibres
et al. 2002). The no-slip boundary condition on the
moving particles boundaries was handled by a second-
order interpolated bounce-back 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 particle-particle
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 finite-size particles. In the
particle-free case, the flow statistics matched precisely
the results obtained from the accurate pseudo-spectral
method, consistent with the observation of Peng et al.
(2010). Three particle-seeded cases reported in Lucci
et al. (2010) were simulated, with very similar results.
The results show that particles of Taylor microscale size
introduce small-scale 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 short-range particle-particle 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 OCI-0904534, ATM-
0527140, and ATM-0730766, and National Natural Sci-
ence Foundation of China (Project No. 10628206).
Computing resources are provided by National Center
for Atmospheric Research through CISL-35751010 and
CISL-35751014.

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. 287-311,
1998

Ayala O., Grabowski W.W, and Wang L.-P., A
hybrid approach for simulating turbulent collisions
of hydrodynamically-interacting particles, J. Comput.
Phys., Vol. 225, pp. 51-73, 2007

Caiazzo A., Analysis of lattice Boltzmann nodes initial-
isation in moving boundary problems, Prog. Comput.
Fluid Dyn., Vol. 8, pp. 3-10, 2008

Climent E., and Maxey M.R., Numerical simulations of
random suspensions at finite Reynolds numbers, Int. J.
Multiphase Flow, Vol. 29, pp. 579-601, 2003

d'Humi&res D., Ginzburg I., Krafczyk M., Lallemand
P., and Luo L.-S., Multiple-relaxation-time lattice Boltz-
mann models in three-dimensions, Phil. Trans. R. Soc.
Lond. A, Vol. 360, pp. 437-451, 2002

Elghobashi S., and Truesdell G., On the two-way inter-
action between homogeneous turbulence and dispersed
solid particles. I: Turbulence modification, Phys Fluids,
Vol. 5, pp. 1790-1801, 1993











Feng Z.G. and Michaelides E.E., The immersed
boundary-lattice Boltzmann method for solving fluid
particles interaction problems, J. Comput. Phys., Vol.
195, pp. 602-628, 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. 20-51, 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. 363-426, 2001

Hu H., Direct simulation of solid-liquid mixtures, Int. J.
Multiphase Flow, Vol. 22, pp. 335-352, 1996

Hu H., Patankar N.A., and Zhu M.Y, Direct numerical
simulations of solid-liquid systems using the aribitrary
Lagrangian-Eulerian technique, J. Comput. Phys., Vol.
169, pp. 427-462, 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.
285-309, 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. 311-339,
1994a

Lallemand P. and Luo L.-S., Lattice Boltzmann method
for moving boundaries, J. Comput. Phys., Vol. 184, pp.
406-421, 2003

Lucci F, Derrante A., and Elghobashi S., Modulation
of isotropic turbulence by particles of Taylor-lengthscale
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. 1603-1626, 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. 883-889, 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 pseudo-spectral meth-
ods for decaying turbulence. Part I. Low-order statistics,
Comput. Fluids, Vol. 39, pp. 568-591, 2010

Peskin C., The immersed boundary method, Acta Nu-
merica, Vol. 11, pp. 479-517, 2002

Squires K.D. and Eaton J.K., Particle response and tur-
bulence modification in isotropic turbulence, Phys. Flu-
ids A, Vol. 2, pp. 1191-1203, 1990

Squires K.D. and Eaton J.K., Preferential concentration
of particles by turbulence, Phys. Fluids A, Vol. 3, pp.
1169-1179, 1991

Sundaram S. and Collins L.R., Collision Statistics in an
isotropic particle-laden turbulent suspension. Part 1. Di-
rect numerical simulations, J. Fluid Mech., Vol. 335, pp.
75-109, 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. 233-271, 2004

Takagi S., Oguz H., Zhang Z., and Prosperetti A.,
Physalis: A new method for particle simulation. Part II:
Two-dimensional Navier-Stokes flow around cylinders,
J. Comput. Phys., Vol. 187, pp. 371-390, 2003

Uhlmann M., An immersed boundary method with di-
rect forcing for the simulation of particulate flows, J.
Comput. Phys., Vol. 209, pp. 448-476, 2005

Uhlmann M., Interface-resolved 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: point-particle based,
hybrid simulations and beyond, Int. J. Multiphase Flow,
Vol. 35, pp. 854-867, 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.
221-233, 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 second-order method for
three-dimensional particle simulation, J. Comput. Phys.,
Vol. 210, pp. 292-324, 2005

Zhou Y, Wexler A.S., Wang L.-P., Modelling turbulent
collision of bidisperse inertial particles, J. Fluid Mech.,
Vol. 433, pp. 77-104, 2001