ICMF 2010, Tampa, FL, May 30 June 4, 2010
Mesoscale simulations of particulate flows with parallel distributed Lagrange
multiplier technique
Y. Kanarska*
Lawrence Livermore National Laboratory, Livermore, CA 94550, USA
kanarskal @llnl.gov
Keywords: particleladen flows, fictitious domain method
Abstract
Fluid particulate flows are common phenomena in nature and industry. Modeling of such flows at micro and
macro levels as well establishing relationships between these approaches are needed to understand properties of
the particulate matter. We propose a computational technique based on the direct numerical simulation of the
particulate flows. The numerical method is based on the distributed Lagrange multiplier technique following the
ideas of Glowinski et al. (1999). Each particle is explicitly resolved on an Eulerian grid as a separate domain, using
solid volume fractions. The fluid equations are solved through the entire computational domain, however, Lagrange
multiplier constrains are applied inside the particle domain such that the fluid within any volume associated with a
solid particle moves as an incompressible rigid body. Mutual forces for the fluidparticle interactions are internal to
the system. Particles interact with the fluid via fluid dynamic equations, resulting in implicit fluidrigidbody coupling
relations that produce realistic fluid flow around the particles (i.e., noslip boundary conditions). The particleparticle
interactions are implemented using explicit forcedisplacement interactions for frictional inelastic particles similar to
the DEM method of Cundall et al. (1979) with some modifications using a volume of an overlapping region as an
input to the contact forces. The method is flexible enough to handle arbitrary particle shapes and size distributions. A
parallel implementation of the method is based on the SAMRAI (Structured Adaptive Mesh Refinement Application
Infrastructure) library, which allows handling of large amounts of rigid particles and enables local grid refinement.
Accuracy and convergence of the presented method has been tested against known solutions for a falling sphere as
well as by examining fluid flows through stationary particle beds (periodic and cubic packing). To evaluate code
performance and validate particle contact physics algorithm, we performed simulations of a representative experiment
conducted at the University of California at Berkley for pebble flow through a narrow opening.
Introduction
Particulate flows occur in a wide range of industrial ap
plications and in nature. Clearly, it's difficult to have
one single simulation method that can cover all length
and time scales. Currently there is a hierarchy of meth
ods that can cover different length and time scales with
different levels of details (Zhu et al. 2007).
When the computational grid size is much larger
then particle size usually a twofluid (or multiphase)
approach is used. The computational fluid dynamics
twofluid approach is often associated with methods de
scribed in Gidaspow (1994) and various implementa
tions are available, such as the commercial code FLU
ENT (FLUENT 2001) and the DOE/NFTL code MFIX
(Syamlal et al. 1993). These twofluid codes usually uti
lize a granularkinetictheory based constitutive model
to represent the "fluid" comprised of fluidized particles,
and empirical twoway coupling relations between fluid
and particles. An advantage of the multifluid model is
that in principle it can be used to compute any multi
phase flow regime. However the effective use of these
models strongly depends on the constitutive or closure
relations for the solid phase and momentum exchange
between phases. For twophase systems comprised of
billions of particles (like, for example, most fluidized
beds or pneumatictransport systems for fine particu
lates) such continuum models are the only computation
ally viable simulation methods available. In fact, devel
opment of a general theory to correctly represent granu
lar flow with fluid as a continuum is still a challenging
research area.
If the cell size is larger than the particle size, a com
bined CFDDEM coupling approach is used. In this ap
proach, the motion of individual particles is obtained by
solving Newton's equations of motion, while the flow of
continuum gas is determined by the CFD on a compu
tational cell scale. A variety of continuum fluid codes,
coupled with discreteelementmethod (DEM) particles,
are utilized by researchers around the world in com
mercial codes as well as codes developed by univer
sity researchers. An overview of these methods can
be found in Zhu et al. (2007). The algorithm relies on
the parameterizations of drag terms similar to the Er
gun equation (Ergun 1952) for static packed beds, or
the WenYu equation for moving beds (Wen and Yu
1966). While some questions regarding parameteriza
tions of drag terms are still remaining, this approach
might be suitable for simulations of intermediatescale
system (about million of particles) and promises to be
a very powerful tool. However, there are some restric
tions of the algorithm and the CFDDEM coupling al
gorithm since it assumes that the cell size in the CFD
model should be larger than the particle size. This may
result in using a fairly coarse mesh in the areas (nozzles,
openings) where only few particles across an important
geometric feature are considered as well as placing se
vere restrictions on the maximum particle size that can
be included in simulations.
And, finally, if the cell size is smaller than the parti
cle diameter the direct numerical simulations, resolved
in the particle and fluid domains, can be applied. In
this case there are no assumed drag terms. Drag ef
fects and the flow around each particle are explicitly re
solved. And particleparticle interactions are modeled,
again, using a DEMlike approach. The method has
no assumptions for drag terms and can be used to im
prove fluidcoupling terms and derive closure parameter
izations that can represent the effective coarsegrained
interactions in the larger scale models just mentioned
above. Moreover, rapid granular flows demonstrate lack
of scale separations which render local closure laws in
applicable in many applications and will require multi
scale approach if it is not feasible to use a high resolution
grid everywhere.
A variety of fixedgrid and meshfree methods have
been utilized for simulation of fluidparticle systems.
We have selected the Lagrange multiplier technique fol
lowing the ideas of Glowinski et al. (1999) and Patankar
et al. (2001) to model fluidparticle systems. The advan
tage of the method is that its finite element formulation
permits the use of a fixed structured grid. This elim
inates the need for remeshing the domain, a necessity
in the unstructured gridbased methods. The objective
of this paper is to present an efficient approach based
on the advantages of utilization of a stationary Eulerian
grid with adaptive mesh refinement (AMR). A parallel
implementation of the method is based on the SAMRAI
library, which allows handling of large amounts of rigid
ICMF 2010, Tampa, FL, May 30 June 4, 2010
particles and enables local grid refinement in the areas
where higher resolutions are needed.
Numerical Model
The method we used is based on the Distributed La
grange Multiplier (DLM) technique of Glowinski et al.
(1999) and Patankar et al. (2001) that was originally de
veloped to study particulate suspension flows. The code
uses a stationary Eulerian grid. Particle positions are
treated as Lagrangian variables. The particle domain is
explicitly resolved on the Eulerian grid using solid vol
ume of fractions. The idea of the method is to solve
the fluid equation in the entire domain, and then correct
the flow inside the rigid domain using Lagrange multi
pliers. The original works are based on elastic collision
forces that prevent particles from overlap. We extended
this approach by incorporating DEM methods for inelas
tic, frictional contact forces. The governing equations
are solved using a fractionalstep scheme for time dis
cretization. The fluid equations are solved in the entire
computational domain at the first stage. It results in a
provisional divergencefree intermediate velocity field.
At the next stage, the constraint of rigid motion (in the
form of Lagrange multiplier) is applied in the solid do
main. A rigid body motion is imposed by constraining
the deformationrate tensor within the particle domain
to be zero. The code is parallelized using the SAM
RAI framework developed at LLNL Hornung and Kohn
(2002). This framework allows to track individual parti
cle position on multiple cpus and refining the resolution
on a structured grid in areas of interest (e.g. solidfluid
interfaces, maximum vorticity zones etc.).
Colissionless governing equations
The idea of the Lagrange multiplier algorithm is based
on the formulations of Glowinski et al. (1999),Patankar
et al. (2001) and Sharma and Patankar (2005). The parti
cle domain is denoted as P(t), where OP is the interface
between the particle and the fluid. F is the fluid domain
that is not shared with the particles. The entire compu
tational domain that includes both the fluid and the par
ticles is denoted by F U P. The governing equations in
the fluid domain can be written as:
(Ou
p at
(u V)u)
Vo
pg f, inFUP (1)
Vu 0, inFU P
D[u] = 0, in P(t)
ICMF 2010, Tampa, FL, May 30 June 4, 2010
D[u] n 0, on OP(t) (4)
u u, in P(t) (5)
where u is the fluid velocity, p is the pressure, p is
the fluid viscosity, p is the density that is equal to pf in
the fluid domain and equal to ps in the particle domain,
g is a body force, and n is a unit normal on the particle
surface. The rigid body velocity inside the particle us is
represented as
us U +Qxr, (6)
where U and Q are the translational and angular veloc
ities of the particle, respectively, and r is the position
vector of the point with respect to the particle centroid.
Force f that appears in the right hand side of the momen
tum equation is nonzero only in the particle domain and
arises as a result of the rigid body motion constraint in
the particle domain. D [u] is the deformationrate tensor
defined as
D[u] = [V + VuT] (7)
2
The stress tensor is then given by:
r =pI T, (8)
where I is the identity tensor, p is the pressure and T is
the viscous stress tensor given by
/pD[u], (9)
for a Newtonian fluid, where p is the viscosity of the
fluid. The particle domain in this formulation is treated
as a fluid with an additional constraint (Lagrange mul
tiplier) to impose the rigid body motion in such a way
that the deformationrate tensor D [u] within the particle
domain is zero. It should be noted that the representa
tion (6) for solid velocity is the sufficient and necessary
for condition (7) to be true. Stresses can develop in the
rigid fluid domains, but the only displacements allowed
are strainfree translations and rotations. The numerical
algorithm to solve the fluidparticle equations of motion
is presented in the next section.
The numerical algorithm
The integration of the governing equations is done us
ing a fractional time stepping approach in the time in
terval j t+1]. In the present algorithm, the velocity
and pressure are cellcentered quantities. The velocity
is defined at integer multiplies of At, whereas the pres
sure is defined at halftimesteps. The system of equa
tions (1),(3) with boundary conditions (4)(5) is solved
numerically using an operatorsplitting technique that
combines the incompressibility condition and advection
diffusion as a first step; and the constraint of rigidbody
motion in the particle domain and the related distributed
Lagrange multiplier technique as a next step.
First step: Solve NavierStokes equations in the entire
domain. At the first step we solve the fluid equations of
motion in the entire computational region and satisfy the
provisional divergencefree intermediate velocity field fi
using an implicit pressure projection technique.
u u
PA
At
V. 1/2 Rn+1/2
V i 0,
(10)
(11)
where R"+1/2 represents all terms on the rightsides of
the momentum equations except the pressure gradient
terms. An unsplit secondorder Godunov procedure is
used to approximate the nonlinear advection term that
appears in the momentum equations using both veloci
ties defined at the centers of the Cartesian grid as well
as velocities defined at the cell faces. The MAC pro
jection method of Bell et al. (1991) that corrects diver
gence of advection velocity along with MUSCL advec
tion scheme of Colella et al. (1985) are used to advect
the fluid. The CrankNicholson scheme is used to com
pute diffusion terms. A divergence constraint is satis
fied using the approximate pressure projection method
of Almgren et al. (1996). The density is computed as
p Ps p p P(1l 4), where 4 is the solid volume frac
tion which is equal to 1 in the solid domain and equal to
0 in the fluid domain.
Second step: Rigid body projection in the particle do
main. At the next stage, the constraint of rigid body mo
tion (in the form of Lagrange multiplier) is applied in the
solid region. The particle velocity u. in a given cell is
split into translational, U, and rotational, Q, parts as
us = U + x r,
(12)
where r is a vector which connects the particle's cen
troid and a center of the considered grid cell. Particle
velocities U and Q are computed by integration of the
provisional velocity field ui in the solid domain as
LMU= f psidV,
P
IPf j r x pfidV,
(13)
(14)
where M is the mass of a particle and Ip is the moment
of inertia of a particle.
The velocity field u"+l is then updated in the solid,
fluid and mixed (where both solid and fluid are pre
ICMF 2010, Tampa, FL, May 30 June 4, 2010
sented) domains as
un"+1 = U, in solid cells
u+1 = ii in fluid cells (15)
un+1 = Ous + (0 1)iu, in mixed cells
where 4 is the solid volume of fraction, computed based
on the particle position at time n. This step could be
considered as finding force f that modifies provisional
velocity field u such that the final velocity field u"+1
satisfies the rigid body motion constraints. Therefore
the force f can be defined as
Un+1 fAt
u U 
Q"+1 = n"+l [r x Fn+1/2]At, (19)
where U"+1 and !n"+1 are particle velocities com
puted at the previous stage.
Fourth step: Update particle position: Explicitly up
date particle position X+ 1 by the following procedure:
X"+l X" UAt,
2
(20)
where X is the position of the particle centroid. If a par
ticle or a solid body is represented as a polyhedron with
vertex coordinates Xv+l, their positions are updated
as
Xvn+ = Xn + [Q+1 x r],
(21)
Also as it was mentioned before the final velocity should
provide zero deformationrate tensor as
D[u"+1] = 0.
Since in the method considered here the final velocity
in the solid domain is known at this step and it is equal
to us, condition (17) is satisfied automatically and the
final velocity is calculated from equations (15) directly.
This approach is similar to one implemented in Patankar
et al. (2001). In other versions of the Lagrange mul
tiplier method (Veeramani et al. 2006; Glowinski et al.
1999) the governing equations (1)(3) include equations
of motion of each individual particle to find U and fQ.
Therefore the resulting equations are solved iteratively.
Third step: Apply particle collision forces. If the
particle Reynolds number and solid volume fraction are
low, particles do not interact and Un+1 and Qf++1 de
fined at the previous step are the final velocities that de
scribe the velocity field inside the particle. When the
concentration of particles is high enough they begin to
interact with each other. If the fluid viscosity is high it
would prevent particles from colliding with each other.
But if fluid forces are insufficient to prevent particle con
tacts, the separately applied collision forces F"n+ sim
ilar to DEM models (Cundall et al. 1979) are used in
our code. Our formulation is based on a viscoelastic
soft collision model. Since we resolve the particle do
main on a stationary Eulerian grid, we compute collision
force in each cell of the overlapping region and then in
tegrate them over the entire overlapping volume. This
is the main difference with the classical DEM approach
Cundall et al. (1979). We will talk more about the calcu
lation of collisional force Fn+1/2 in a separate section.
Once this force is calculated, the final particle veloci
ties translationall and rotational) are updated using the
velocity field U and from the previous time step as:
Un+1 Un+1 Fn+1/2At, (18)
Once the vertex positions are defined the volume frac
tions for solid bodies are computed and used for compu
tations at the next time step.
Treatment of collisions
Different collision models have been developed for the
coupled solvers for the fluid and solid systems. These
collision models aim to capture the collision process of
solid particles by introducing shortrange forces as addi
tional body forces acting on the particles. Treatment of
collisions includes a contact detection algorithm and ap
plying a collision forces to prevent particles from over
lapping. In the present work the contact is detected in the
region where two particles overlap based on the volume
of fraction function 0. The condition 0 > 1 means that
more than one particle exists in the given grid cell, there
fore additional collision force is applied in each grid cell
where particle overlap. While this approach would not
be optimal for normal DEM algorithms, it is a natural fit
for fluidparticle problems. Moreover, it easily extends
to nonspherical objects which are hard to accommodate
in a typical DEM approach. The collision force acts
along the normal (Fncomponent) as well as the tangen
tial direction (Ftcomponent) at the point of contact be
tween two particles. For spherical particles, the normal
and tangential directions are defined by the line joining
the centers of two colliding particles and two lines per
pendicular to it, respectively. The components of the
normal collisional force are computed in each cell where
particles overlap and then integrated over the overlapped
region as
F,= (k,,Vjk
dV /mp,
(22)
where Vijk is the volume of cell ijk in the overlapped
region, v, is the relative velocity between the two inter
acting particles in the normal direction, k, is the normal
ICMF 2010, Tampa, FL, May 30 June 4, 2010
spring constant (or stiffness), d, is the damping coef
ficient in the normal direction, and mn is the particle
mass.
The tangential component is computed according to a
Coulomb friction law as
F ijEk
I IF, I t,
,Ftl < lf F I
IFtl' > lf
where p is the friction coefficient, and dt is the
damping coefficient in the tangential direction.
An important advantage of the present collision
method is that it is flexible enough to handle arbitrary
particle shapes and size distributions and doesn't re
quire extra parametrizations for fluidparticle interac
tions. The shear friction forces discussed above can al
low only slow movement in the tangential direction, but
cannot stop or reverse tangential motion. This formu
lation is inadequate for applications that require truly
static friction, such as heap formation or angles of re
pose. In such situations, there is a threshold force be
low which the grains do not move at all, opposed by
static friction. However, implementation of even a sim
ple historydependent threshold rule is algorithmically
complicated and will be a subject of future work. Thus
the applications in this paper are limited by the visco
elastic collisional model described above.
According to the Hertzian contact theory Hertz
(1882), the relation between the normal force F, and
displacement 6 is given by
F, KH /2, (24)
where KH is the Hertzian stiffness. In the case of two
spheres of the same size with radius r, KH is expressed
by
2rE
KH 3(1 2' (25)
where E is Young's modulus and a is the Poisson ratio
of the particles. Since in (22) the normal contact force
is proportional to the volume of the overlapping region,
our computational stiffness can be related to the Hertzian
stiffness as
KH = k,,r61/2 (26)
Parallelization
The code has been built on top of the SAMRAI (Struc
tured Adaptive Mesh Refinement) library developed at
LLNL Hornung and Kohn (2002). SAMRAI is a general
objectoriented software infrastructure for implement
ing parallel scientific applications that employ structured
adaptive mesh refinement. The method uses a hierarchi
cal structured grid approach first developed by Berger
and Oliger (1984). In particular, AMR is based on a se
quence of nested grids with successively finer spacing in
both time and space. Increasingly finer grids are recur
sively embedded in coarse grids until the solution is suf
ficiently resolved. An error estimation procedure evalu
ates where additional refinement is needed and grid gen
eration procedures dynamically create or remove rectan
gular fine grid patches as resolution requirements change
. Automatic regridding in time is based on Richardson
extrapolation and in space on detection of gradients (ve
locity, scalar etc) in the solution. SAMRAI provides the
backbone of our implementation, managing the locally
refined Cartesian grid patch hierarchy with both the Eu
lerian and Lagrangian data points defined on the hierar
chy. It also provides facilities for performing adaptive
regridding, load balancing, and parallel data communi
cation. To store and manage the Lagrangian data points
a version of the SAMRAI IndexData patch data type is
used. For a generalpurpose solver library, we have cho
sen PETSc (Balay et al. 2004), developed at Argonne
National Laboratory. This suite solves largescale lin
ear and nonlinear equations. We used preconditioned
Krylov methods provided by this library. A parallel data
managing and implementation is done similar to the al
gorithm described in Griffith et al. (2001).
Validation against empirical data and
experiments
Fixed particle beds
Flow behavior through packed beds of spheres or other
porousmedialike structures are of crucial importance
in industry and nature. The determination of pressure
drop through a packed bed as a function of fluid flow
rate, geometrical constrains of the bed and physical
properties of bed material is very critical, for example, in
hydraulic and pneumatic devices. The wellknown em
pirically derived equation used for that purpose has been
proposed by Ergun (Ergun 1952) based on experimental
measurements:
Ap (1 E)2 v (1
L E3 D +
 E) P2
3 DP
(27)
where Ap is the pressure drop through the packed bed,
L is the bed length or height, E vS is the bed porosity,
where Vf is the voids volume, V is the total volume,
u is the superficial velocity at the exit of bed, v is the
fluid dynamic viscosity, and Dp is the particle diameter.
Ergun equation (27) incorporates momentum loses due
to viscous effects (first term in (27) which is important in
the laminar regime) and inertia effects (the second term
ICMF 2010, Tampa, FL, May 30 June 4, 2010
in (27) which dominates in the turbulent regime). The
Ergun equation is often used in a more general form by
introducing a friction coefficient A which is defined as
c Ap D3 2
A A + BRecP (28)
L (1 c)2VU'
where Reynolds number is defined as Re J= Dp
v(l )
A standard form of the Ergun equation uses the values
of empirical constants A 150, B 1.75, C 1.
However many other studies have been performed and
published to check the appropriate choice of the empir
ical constants A,B and C in (27). Some authors pro
posed values in the range 150 200 for A and 1.7 4.0
for B as well as functional forms for these coefficients
that depend on both porosity and Reynolds number, see
overview in Plessis (2001). Vortwek and Brunn (1994)
proposed to use values A 181, B 2.01 and C
0.96 to estimate the pressure drop in randomly arranged
packed beds. Franzen (1979) found that for A 164.97,
B 1.976 and C 0.9, equation (28) gives good esti
mate of the pressure drop for regularly arranged spheres
in the cubic packing. The reason for variation in the con
stants was determined as the variations in particle geom
etry and orientation, as well as macroscopic properties
of the packing. It should be noted that the Ergun equa
tion was derived for densely packed beds, and is not ex
pected to be valid for high void ratios. For that range
(when porosity is smaller than 40 %), normally the Wen
and Yu equation (Wen and Yu 1966) is used. However
this equation significantly underpredicts the drag force at
higher Reynolds numbers (Beetstra et al. 2007). Also the
transition from the Ergun to the Wen and Yu equation is
not a smooth function. For moving particles the situation
becomes even more complicated since particles begin to
interact with each other and may dissipate additional en
ergy that affects the pressure drop. Therefore detailed
mesoscale simulations that resolve the flow around each
particle are needed to predict flow characteristics such
as the pressure drop and volumetric flow rates in these
systems.
As a first step, we validate our code against experi
mental and empirical data for both randomly and reg
ularly arranged packing in different flow regimes and
Reynolds numbers. First we consider a configuration
that consists of monodisperse spheres of radius R =
0.002m arranged in a cubic packing. We consider two
flow regimes with Reynolds numbers Re = 18 and
Re = 900. A constant inflow velocity field of u
0.5 m/s was prescribed at the left boundary and open
boundary conditions are set up at the right boundary.
Other boundaries are chosen to be periodic. The porosity
is = 47.64%. Figure 2 shows measured values of the
friction coefficient A in different data sets that include
Ergun (1952); Franzen (1979); Hovekamp (2002); Mar
tin et al. (1951); Vortwek and Brunn (1994). For small
Re numbers, difference in the data sets for both random
and cubic packing is small (Figure 2a). For high Re
numbers, the friction coefficient is found to deviate by
more than half an order magnitude in the existing empir
ical and experimental data (Figure 2b). Figure 3 shows
flow patterns at different resolutions for Re 18. We
use Richardson extrapolation f(h) fexact + Ch to
estimate the convergence order using computed values
of f(h). Here f is a calculated parameter (pressure drop
in our case), h is some measure of grid spacing, C is
a constant, and p is the order of convergence. Based
on the results of simulations we found that the conver
gence rate is close to a second order: 1.86 for Re 18
and 1.8 for Re 900 (Figure 4 and Figure 6, respec
tively). However the low Reynolds number flow requires
less resolution to achieve a converged solution then the
high Reynolds number flow. About 48 cells per particle
are needed to describe flow behavior for Re 18. In
the case of Re 900, the resolution needs to be twice
higher to get the same order of error. This is mainly be
cause of the flow separation and turbulent boundary ef
fects that require finer grid resolutions to describe them
adequately (Figure 5). In these simulations we do not
use any turbulent model. However in future studies we
may need to incorporate turbulent effects and additional
drag terms to be able to simulate high Re number flow
regimes with relatively modest resolution. Another op
tion to improve the convergence rate is to revise the in
terpolation scheme for the velocity field in the mixed
cells in (15). The grid resolution needed for a con
verged solution depends on the particular configuration
and porosity. For densely packed beds, the flow char
acteristics are constrained by the maximum resolution
available in the void space between particles. Therefore
fineresolution simulations are required to describe flow
effects through these small void spaces for high Re num
bers. The overall agreement with available data is very
good (Figure 2).
The next example is flow in a periodic domain of ran
domly packed spheres of radius R = 0.11283m. Re
number in this configuration is about 900 and the solid
packing ratio is 60%. A superficial velocity is specified
as u = 0.008 m/s. The periodic random distribution of
spheres is generated using Donev et al. (2005) algorithm.
The pressure drop through the packed bed is illustrated
in Figure 7. We perform a numerical analysis to estimate
the convergence rate of the numerical solution. The error
in the pressure drop decreases with approximately sec
ond order when refining the mesh (Figure 8). The dif
ference between our converged solution and the Ergun
equation (27) for coefficients A 150 and B 1.75
is found to be 7%. This is quite reasonable agreement
ICMF 2010, Tampa, FL, May 30 June 4, 2010
since the considered configuration with given porosity
of II'; is a good representation of the cases investigated
in the experiments by Ergun (1952). It should be noted
that for loosely packed beds with porosity larger than
!II' different random configurations may produce dif
ferent results and deviations from Ergun values, this is
discussed in Freund et al. (2003). In the next section we
show how local inhomogeneities in the porosity distribu
tion of randomly packed beds with the same global pa
rameters (porosity, particle size) may influence the flow
and transport properties.
Falling sphere
We further validate our code by investigating the sedi
mentation of a cylindrical particle in a Newtonian fluid.
We consider sedimentation of a single twodimensional
disk with radius 1cm in a channel with dimension of
0.2mxO.4m. The fluid density is 2000kg/m3 and the
particle density is 2500kg/rm3. The various Reynolds
number dependent flow regimes are obtained by varying
the fluid dynamic viscosity v. We consider three cases
here with v 3, 0.1 and 0.002 Pas. Gravity accelera
tion 9.8m/s2 acts in a negative y direction. The simu
lation is started at t Os by dropping a particle at the
center of the channel at 0.33cm depth. Different grid
resolutions are considered ranging from 4 up to 256 grid
cells per particle diameter.
Figure 9 shows vertical velocity snapshots and vortic
ity contours for Re 0.2, Re 100 and Re 2000
correspondingly. Vorticity contours are plotted as an il
lustration of the different flow regimes which depend on
the Re number. For the small Re number the flow is
laminar and the disk reaches its terminal velocity rel
atively quickly (Figure 10a). With increasing the Re
number a vortex forms behind the particle. For high
Re numbers vortex shedding occurs that lead to hori
zontal deviation of the trajectory of the particle (Fig
ure 9c). We found that the vortex shedding begins at
Re 400, which agrees well with experimental and
theoretical findings (Achenbach 1974).
We also investigated the convergence rate for differ
ent Re numbers. The results agree well with the previ
ous finding for the fixed array of spheres considered in
the previous section. Whereas for small Re numbers,
only few cells per particle diameter (around 6) are re
quired to get a converged solution, for high Re num
bers, about 32 grid cells are needed to represent the flow
around the particle and the drag effects correctly (Figure
10, the right panel).
Particle flow through opening
The examples considered in the previous sections fo
cused mainly on the flow through stationary packed
beds. When the particles are not so constrained, the mo
tion of the fluid leads to particle motion, and the mov
ing particles interact with each other as well as with
the fluid leading to complex flow patterns that depend
significantly on the packing density and particle con
tact physics. For highly viscous fluids, viscous forces
could prevent contact between particles, but for most ap
plications, these forces are insufficient to prevent par
ticle contacts. As it was mentioned in () to account
for particleparticle interactions, separately applied col
lision forces similar to those used in the distinct element
method (DEM) (Cundall et al. 1979) are implemented
in the code. To validate the code performance and the
particle contact physics algorithm, we performed simu
lations of a representative experiment conducted at UC
Berkeley. The experimental configuration is shown in
Figure 15. As shown, 2500 polypropylene pebbles, ini
tially suspended in a vertical water column, were dis
charged through a narrow opening. The pebble reservoir
was 40.6cm in diameter, while the diameter of the nar
row opening has a diameter of 10.2cm. A conical re
gion, with a 45 degree cone angle, connects the pebble
reservoir to the narrow opening. The surface of the cone
has been perforated to allow back flow of water into the
reservoir as the pebbles were evacuated. The pebbles
were all the same size, with a diameter of 2.5cm and
a density of II I ;./cm3. The bottom boundary of the
container is closed, and a physical barrier is placed in
the narrow opening to inhibit pebble motion such that,
at the beginning of the experiment, both the pebbles and
water are at rest.
The experiment is initiated by removing the physi
cal barrier, thus allowing the pebbles to move upward
through the water column. As the pebbles flow out of
the lower reservoir, the evacuated volume is filled with
water flowing into the reservoir through the porous con
ical section. In the simulations, the holes in the conical
section were not modeled explicitly. Instead, the effect
of the holes was taken into account by allowing the con
ical surface to physically constrain pebble motion while
at the same time be transparent to fluid flow.
The comparisons between simulation results and ex
perimental data are presented in Figures 11 and 12. The
overall agreement is very good, though there are some
discrepancies. Initially, simulation show that the pebbles
move faster in the central region, and slower near the
wall with some of the pebbles near the wall moving in
a downward direction (see Figure 1 la). This can be ex
plained by the fact that the counter fluid flow is not dis
tributed uniformly through the cone section. There could
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 1: Geometry of AMR patches for simulation of fluid flow through the cubic array of spheres
0
. O
e _.o
<random> <CUoIC>
Figure 2: Comparison of the computed friction coefficient with available data for radnom and cubic packing at
Re = 18 (left panel) and Re = 900 (right panel). Data from Hovekamp (2002),Vortwek and Brunn (1994), Ergun
(1952), Franzen (1979) and Martin et al. (1951) are shown.
ICMF 2010, Tampa, FL, May 30 June 4, 2010
..........12 cells/particle diameter
U
24 cells/particle diameter
EL r
a~. 
*6 
es&.L ..~W~~W W~j
48 cells/particle diameter
mu
a iii Iiii ii
XtajieieZ
S a a a
Figure 3: Horizontal velocities at the central slice for the cubic packing of spheres at Re = 18.
e
LU
0 0 1
001
0001
Numerical Error 0
Richardson Extrapolation, slope=1 86
2nd order convergence slope ..........
.. ..
Cells/particle diameter
Figure 4: Convergence rate at low Reynolds numbers, Re = 18. Numerical error versus resolution is shown.
;"
E.  I

I
ICMF 2010, Tampa, FL, May 30 June 4, 2010
.: 12 cells/particle diameter 24 cells/particle diameter
_
...... .....
S 48 cells/particle diameter 96 cels/parice diameter
i .
i 
: a .
Figure 5: Horizontal velocities at the central slice for the cubic package of spheres at Re 900
100
Numerical Error 0
Richardson Extrapolation, slope=1 8
2nd order convergence slope .........
1 . ...0
i...............
...............
S 01
001
0001
10 100
Cells/particle diameter
Figure 6: Convergence rate at low Reynolds numbers, Re = 900. Numerical error versus resolution is shown.
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 7: Pressure drop through the randomly packed bed.
Numerical Error 0
Richardson Extrapolation, slope=1 8
2nd order convergence slope ...........
...........
.............
. . .. ..
Cells/particle diameter
Figure 8: Convergence rate. Numerical error versus resolution is shown.
[Re=100]
I [R
[Re=2000]
Figure 9: Numerical simulation of falling particle in a Newtonian fluid. A vertical velocity in the fluid is shown for
different Reynold numbers. Vorticity contours are plotted as an illustration of the different flow regimes
[Re=0.2]
ICMF 2010, Tampa, FL, May 30 June 4, 2010
8 cells partcle diameter 
16 cells/part cle di ameter
32 cells/particle diameter 
64 cells/particle diameter 
0 05 1 15 2 25 3 35 4
Time, s
[Re=100]
0 02 04 06 08 1 12 14
Time, s
[Re=2000]
0 005 01 015 02 025
Time, s
Figure 10: Settling velocity and convergence rate for Re =
01
001
0001
00001
1e05
1e06
Numerical Error O
Richardson Extrapolation, slope=2
2nd order convergence slope ...........
.. .. .. .
3 100
Cells/particle diameter
Numerical Error 0
''. Richardson Extrapolation, slope=1 97
2nd order convergence slope .........
""***..
0 100
Cells/particle diameter
Numerical Error 0
... Richardson Extrapolation, slope=1 46
2nd order convergence slope ..........
.".O..
"\^.
00001
10 100
Cells/particle diameter
0.2 (a), Re = 100 (b) and Re = 2000 (c)
[Re=0.2]
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 11: Comparison of the pebble flow in experiment and simulations at time 9.5 s (left panel) and 28.5 s (right
pannel)
0.05
experiment 0
0 simulations 0
0.05 
0.1 C
E 0.15 
0.2 
0.25
G 0.3 
0.35 
0.4
0.45 C
E 0.15 
0.5
5 0 5 10 15 20 25 30 35
Time, s
Figure 12: Comparison between experiment and simulations for the bed position below the chute outlet
Figure 12: Comparison between experiment and simulations for the bed position below the chute outlet
ICMF 2010, Tampa, FL, May 30 June 4, 2010
be several reasons for this behavior. One reason may be
that the holes in the conical section were not modeled
explicitly. Another, more likely reason, is that our ini
tial packing was not ideally conformed to the boundary
of the container, thus leading to higher porosity near the
skin of the container than in the interior of the pebble
bed, thus leading to higher fluid flow rates near the outer
boundary of the cylinder. However, it appears that those
local fluctuations do not affect the integral characteris
tics of the flow field. A comparison of the simulated
evolution of the bed bottom position with that observed
during the experiment is shown in Figure 12. The bed
bottom position was taken as the lowest pebble position
in the central cylindrical slice of radius 5.1cm.
Conclusions
We have presented a distributed Lagrange multiplier al
gorithm for particulate flows. The method is following
ideas of Patankar et al. (2001); Sharma and Patankar
(2005); Glowinski et al. (1999). The idea of the method
is to use operatorsplitting technique to solve the fluid
equations in the entire domain first and then correct the
flow inside the rigid domain using Lagrange multipliers.
The parallel implementation of the algorithm is done us
ing SAMRAI library. Previous work is did not include
collisions or base them on elastic forces that prevent par
ticles from overlap. We extend this approach by incor
porating DEM methods for inelastic, frictional contact
forces. Following Patankar et al. (2001), the proposed
numerical method does not include additional equations
of motion for the particle translational and angular ve
locities. Code validation was done by comparing numer
ical results with known experimental and empirical data
for a falling sphere in a Newtonian fluid, flow through
the stationary packed beds, and pebble release through
a narrow opening. We performed numereous conver
gence tests for different applications and found that for
the low and moderate Re numbers the convergence of
the method is close to the secondorder. For the high Re
numbers the convergence becomes slightly slower and
could be explained by deficiencies in the flow represen
tation near the rigid walls, in particular,the velocity in
terpolation in the mixed cells where both fluid and solid
coexist and necessity to resolve turbulent structures in
the fluid domain. The overall performance and accu
racy of the code is very good and promises to be a valu
able tool both for simulations of flow involving up to a
few hundred thousand particles as well as for calibrating
the phasecoupling relationships in twofluid continuum
simulation models.
Acknowledgements
This work performed under the auspices of the U.S. De
partment of Energy by Lawrence Livermore National
Laboratory under Contract DEAC5207NA27344.
References
Achenbach E., 1974: Vortex shedding from spheres. J.
Fluid Mech., 62(2), 209221.
Almgren A.S., J.B. Bell, W. G. Szymczak, 1996: A
numerical method for the incompressible NavierStokes
equations based on an approximate projection. SIAM J.
Sci. Comput, 17(2), 358369.
Balay S., K. Buschelman, V. Eijkhout, W.D. Gropp,
D. Kaushik, M.G. Knepley, L.C. McInnes, B.F. Smith,
and H. Zhang, 2004: PETSc Users Manual. Technical
Report ABNL95/11Revision 2.1.5, Argonne National
Laboratory, 2004.
Beetstra R., M.A. van der Hoef, J.A.M. Kuipers,
2007: Numerical study of segregation using a new drag
force correlation for polydisperse systems derived from
latticeBoltzman simulations. Chemical Eng. Science,
62, 246255.
Bell J.B., P. Colella, L.H. Howell, 1991: An efficient
secondorder projection method for viscouse incom
pressible flow. In 10th AIAA Computational Fluid Dy
namics Conference, Honolulu, June 2427, 1991.
Berger J., J. Oliger, 1984: Adaptive mesh refinement
for hyperbolic partial differential equations. J. Comp.
Phys., 53,484494
Colella P., 1985: A direct Eulerian MUSCL scheme for
gas dynamics, SIAM J. Comput, 6, 104117.
Colella P., 1990: Multidimensional upwind method for
hyperbolic conservation laws.J. Comp. Phys.,87, 171
200.
Cundall PA. and O.D.L. Strack, 1979: A discrete nu
merical model for granular assemblies. Geotechnique,
29, 4765.
Donev A., F. H. Stillinger, and S. Torquato, 2005:Neigh
bor List CollisionDriven Molecular Dynamics Simula
tion for Nonspherical Particles. I. Algorithmic Details.
J. Comp. Phys., 202(2), 737764.
Ergun, S., 1952: Fluid flow through packed columns.
Chemical Engineering Progress 48, 8994.
Franzen P., 1979: Zum EinfluB der Porengeometrie auf
den Druckverlust bei der Durchstr"omung von Poren
systemen, I. Versuche an Modellkan"alen mir variablem
Querschnitt. Rheol. Acta, 18:392 423.
Freund, H., T. Zeiser, F Huber, E. Klemm, G. Bren
ner, F Durst, G. Emig, 2003: Numerical simulations of
single phase reacting flows in randomly packed fixed
bed reactors and experimental validation. Chemical En
gineering Science, 58, 903910.
FLUENT 6.0 User's Guide Book, 2001.
Gidaspow, D., 1994: Multiphase flow and fluidization.
Continuum and Kinetic theory descriptions. Academic
Press, NY.
Glowinski, R., T.W. Pan, T.I. Hesla, D.D. Joseph, 1999:
A distributed Lagrange multiplier fictitiouss domain
method for particulate flow. Int. J. Multiphase Flow, 25,
755794.
Griffith B.E., R.D. Homung, D.M. McQueen, C.S. Pe
skin, 2001: Parallel and adaptive simulation of cardiac
fluid dynamics.
Freund H., T. Zeiser, F. Huber, E. Klemm, G. Brenner, F
Durst, G. Emig, 2003: Numerical simulations of single
phase reacting flows in randomly packed fixedbed reac
tors and experimental validation. Chemical Engineering
Science, 58, 903910.
Hertz H., 1882: Uber die Beriihrung fester elastische
Kirper, Journal of Reine und Angewandte Mathematik
92, 156171.
Hornung, R.D., A.M. Wissink, S.R. Kohn, 2002: Man
aging complex data and geometry in parallel structured
AMR applications. Engineering with Computers, 22(3
4), 181195.
Hovekamp, T.B., 2002: Experimental and numerical
investigation of porous media flow with regard to the
emulsion process. PhD thesis, 1101.
Marshall J.S., 2009: Discreteelement modeling of par
ticulate aerosol flows. J. Comp. Phys., 228, 15411561.
Martin, J.J., W.L. Mccabe, C.C. Monrad, 1951: Pres
sure drop through stacked spheres: Effect of orienta
tion. Chemical Engineering Progress, 47(2): 9194. En
gineering with Computers, 22(34), 181195.
Patankar N.A., 2001: A formulation for fast computa
tions of rigid particulate flows. Center Turbul. Res., Ann.
Res. Briefs, 185196.
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Plessis J., 2001: Analytical Quantification of Coeffi
cients in the Ergun Equation for Fluid Friction in a
Packed Bed.!!YEAR AND JOURNAL!!!
Sharma N. and N.A.Patankar, 2005: A fast computation
technique for the direct numerical simulation of rigid
particle flows. J. Comp. Phys., 205,439457.
Syamlal M., W. Rogers and T.J.O. Brien, MFIX Doc
umentation: Theory Guide U.S. Department of Energy
(DOE), Morgantown Energy Technology Center, Mor
gantown, West Virginia (1993).
Tsuiji Y., T. Tanaka and T. Ishida, 1992: Lagrangian nu
merical simulation of plug flow of cohesionless particles
in a horizontal pipe. Powder Technology, 71, 239250.
Tsuji Y., T. Kawaguchi and T. Tanaka, 1993: Discrete
particle simulation of twodimensional fluidized bed.
Powder Technology, 77, 7987.
Xu B.H., A.B. Yu, 1997: Numerical simulation of the
gassolid flow in a fluidized bed by combining dis
crete particle method with computational fluid dynam
ics. Chemical Eng. Science, 52, 27852809.
Veeramani C., P.D. Minev, K. Nandakumar, 2006: A
fictitious domain formulation for flows with rigid par
ticles: A nonLagrange multiplier version. J. Comp.
Phys., 224,867879.
Vorwerk J., P.O. Brunn, 1994: Shearing effects for
the flow of surfactant and polymer solutions through a
packedbed of spheres. Journal of NonNewtonian Fluid
Mechanics, 51(1), 7995.
Wen, C.Y and Y.H. Yu, 1966: Mechanics of fluidization,
Chem. Eng. Prog. Symp. Ser 62, 100111.
Zhu H.P, Z.Y. Zhou, R.Y Yang, A.B. Yu, 2007: Discrete
particle simulation of particulate systems: Theoretical
developments. Chemical Eng. Science, 62, 33783396.
