Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 12.5.1 - Mesoscale simulations of particulate ows with parallel distributed Lagrange multiplier technique
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00308
 Material Information
Title: 12.5.1 - Mesoscale simulations of particulate ows with parallel distributed Lagrange multiplier technique Particle-Laden Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Kanarska, Y.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: particle-laden flows
fictitious domain method
 Notes
Abstract: Fluid particulate ows are common phenomena in nature and industry. Modeling of such ows 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 ows. 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 uid equations are solved through the entire computational domain, however, Lagrange multiplier constrains are applied inside the particle domain such that the uid within any volume associated with a solid particle moves as an incompressible rigid body. Mutual forces for the uid-particle interactions are internal to the system. Particles interact with the uid via uid dynamic equations, resulting in implicit uid-rigid-body coupling relations that produce realistic uid ow around the particles (i.e., no-slip boundary conditions). The particle-particle interactions are implemented using explicit force-displacement interactions for frictional inelastic particles similar to the DEM method of Cundall et al. (1979) with some modi cations using a volume of an overlapping region as an input to the contact forces. The method is exible enough to handle arbitrary particle shapes and size distributions. A parallel implementation of the method is based on the SAMRAI (Structured Adaptive Mesh Re nement Application Infrastructure) library, which allows handling of large amounts of rigid particles and enables local grid re nement. Accuracy and convergence of the presented method has been tested against known solutions for a falling sphere as well as by examining uid ows 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 ow through a narrow opening.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00308
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 1251-Kanarska-ICMF2010.pdf

Full Text
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: particle-laden 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 fluid-particle interactions are internal to
the system. Particles interact with the fluid via fluid dynamic equations, resulting in implicit fluid-rigid-body coupling
relations that produce realistic fluid flow around the particles (i.e., no-slip boundary conditions). The particle-particle
interactions are implemented using explicit force-displacement 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 two-fluid (or multiphase)
approach is used. The computational fluid dynamics
two-fluid 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 two-fluid codes usually uti-
lize a granular-kinetic-theory based constitutive model
to represent the "fluid" comprised of fluidized particles,


and empirical two-way coupling relations between fluid
and particles. An advantage of the multi-fluid 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 two-phase systems comprised of
billions of particles (like, for example, most fluidized-
beds or pneumatic-transport 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 CFD-DEM 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 discrete-element-method (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 Wen-Yu 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 intermediate-scale
system (about million of particles) and promises to be
a very powerful tool. However, there are some restric-
tions of the algorithm and the CFD-DEM 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 particle-particle interactions are modeled,
again, using a DEM-like approach. The method has
no assumptions for drag terms and can be used to im-
prove fluid-coupling terms and derive closure parameter-
izations that can represent the effective coarse-grained
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 fixed-grid and mesh-free methods have
been utilized for simulation of fluid-particle systems.
We have selected the Lagrange multiplier technique fol-
lowing the ideas of Glowinski et al. (1999) and Patankar
et al. (2001) to model fluid-particle 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 grid-based 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 fractional-step scheme for time dis-
cretization. The fluid equations are solved in the entire
computational domain at the first stage. It results in a
provisional divergence-free 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 deformation-rate 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. solid-fluid
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)


V-o


pg f, inFUP (1)


V-u 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 non-zero only in the particle domain and
arises as a result of the rigid body motion constraint in
the particle domain. D [u] is the deformation-rate 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 deformation-rate 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 strain-free translations and rotations. The numerical
algorithm to solve the fluid-particle 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 cell-centered quantities. The velocity
is defined at integer multiplies of At, whereas the pres-
sure is defined at half-timesteps. The system of equa-
tions (1),(3) with boundary conditions (4)-(5) is solved
numerically using an operator-splitting technique that


combines the incompressibility condition and advection-
diffusion as a first step; and the constraint of rigid-body
motion in the particle domain and the related distributed
Lagrange multiplier technique as a next step.
First step: Solve Navier-Stokes 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 divergence-free 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 right-sides of
the momentum equations except the pressure gradient
terms. An unsplit second-order 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 Crank-Nicholson 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 deformation-rate 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 visco-elastic
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 short-range 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 fluid-particle problems. Moreover, it easily extends
to non-spherical objects which are hard to accommodate
in a typical DEM approach. The collision force acts
along the normal (Fn-component) as well as the tangen-
tial direction (Ft-component) 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 fluid-particle 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 history-dependent 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
object-oriented 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 general-purpose solver library, we have cho-
sen PETSc (Balay et al. 2004), developed at Argonne
National Laboratory. This suite solves large-scale 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
porous-media-like 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 well-known 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
fine-resolution 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 two-dimensional
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 particle-particle 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


1e-05


1e-06


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 operator-splitting 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 second-order. 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 phase-coupling relationships in two-fluid 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 DE-AC52-07NA27344.


References

Achenbach E., 1974: Vortex shedding from spheres. J.
Fluid Mech., 62(2), 209-221.

Almgren A.S., J.B. Bell, W. G. Szymczak, 1996: A
numerical method for the incompressible Navier-Stokes
equations based on an approximate projection. SIAM J.
Sci. Comput, 17(2), 358-369.

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 ABNL-95/11-Revision 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
lattice-Boltzman simulations. Chemical Eng. Science,
62, 246-255.

Bell J.B., P. Colella, L.H. Howell, 1991: An efficient
second-order projection method for viscouse incom-
pressible flow. In 10th AIAA Computational Fluid Dy-
namics Conference, Honolulu, June 24-27, 1991.

Berger J., J. Oliger, 1984: Adaptive mesh refinement
for hyperbolic partial differential equations. J. Comp.
Phys., 53,484-494

Colella P., 1985: A direct Eulerian MUSCL scheme for
gas dynamics, SIAM J. Comput, 6, 104-117.

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, 47-65.

Donev A., F. H. Stillinger, and S. Torquato, 2005:Neigh-
bor List Collision-Driven Molecular Dynamics Simula-
tion for Nonspherical Particles. I. Algorithmic Details.
J. Comp. Phys., 202(2), 737-764.

Ergun, S., 1952: Fluid flow through packed columns.
Chemical Engineering Progress 48, 89-94.






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, 903-910.

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,
755-794.

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 fixed-bed reac-
tors and experimental validation. Chemical Engineering
Science, 58, 903-910.

Hertz H., 1882: Uber die Beriihrung fester elastische
Kirper, Journal of Reine und Angewandte Mathematik
92, 156-171.

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), 181-195.

Hovekamp, T.B., 2002: Experimental and numerical
investigation of porous media flow with regard to the
emulsion process. PhD thesis, 1-101.

Marshall J.S., 2009: Discrete-element modeling of par-
ticulate aerosol flows. J. Comp. Phys., 228, 1541-1561.

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): 91-94. En-
gineering with Computers, 22(3-4), 181-195.

Patankar N.A., 2001: A formulation for fast computa-
tions of rigid particulate flows. Center Turbul. Res., Ann.
Res. Briefs, 185-196.


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,439-457.

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, 239-250.

Tsuji Y., T. Kawaguchi and T. Tanaka, 1993: Discrete
particle simulation of two-dimensional fluidized bed.
Powder Technology, 77, 79-87.

Xu B.H., A.B. Yu, 1997: Numerical simulation of the
gas-solid flow in a fluidized bed by combining dis-
crete particle method with computational fluid dynam-
ics. Chemical Eng. Science, 52, 2785-2809.

Veeramani C., P.D. Minev, K. Nandakumar, 2006: A
fictitious domain formulation for flows with rigid par-
ticles: A non-Lagrange multiplier version. J. Comp.
Phys., 224,867-879.

Vorwerk J., P.O. Brunn, 1994: Shearing effects for
the flow of surfactant and polymer solutions through a
packed-bed of spheres. Journal of Non-Newtonian Fluid
Mechanics, 51(1), 79-95.

Wen, C.Y and Y.H. Yu, 1966: Mechanics of fluidization,
Chem. Eng. Prog. Symp. Ser 62, 100-111.

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, 3378-3396.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs