7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A VariableDensity FictitiousDomain Method for Fully Resolved Simulation of
HighDensity Ratio FluidParticle Systems
Sourabh V. Apte* and Justin Finn*
School of Mechanical Industrial and Manufacturing Engineering
Oregon State University, Corvallis, OR 97333, USA
sva@engr.orst.edu and finnj@engr.orst.edu
Keywords: Fully resolved simulations, fictitious domain methods, highdensity ratio,
freely moving particles, particlevortex interactions.
Abstract
A numerical scheme for fully resolved simulation of fluidparticle systems with freely moving rigid particles is
developed. The approach is based on a fictitious domain method wherein the entire fluidparticle domain is assumed
to be an incompressible fluid but with variable density. The flow inside the particle domain is constrained to be a
rigid body motion using an additional rigidity constraint in a fractional step scheme. The rigidity constraint force
is obtained based on the fast computation technique proposed by Sharma and Patankar (2005). The particle is
assumed to be made up of material points moving on a fixed background mesh where the fluid flow equations are
solved. The basic finitevolume solver is based on a colocated grid incompressible but variable density flow. The
incompressibility constraint is imposed by solving a variablecoefficient pressure equation giving rise to a stable
scheme for high density ratio fluidparticle systems. Through various verification and validation test cases on fixed
and freely moving particles it is shown that the numerical approach is accurate and stable for a wide range of
fluidparticle density ratios.
Introduction
Fully resolved simulations (FRS) of fluidparticle sys
tems, wherein all scales associated with the fluid and
particle motion are completely resolved, are of great im
portance for understanding particleturbulence interac
tions. Considerable work has been done on fully re
solved simulations of freely moving particles in fluid
flows on fixed grids. For example, distributed La
grange multiplier/fictious domain (DLM) based meth
ods (Glowinski et al. (2001)) and Immersed Boundary
method (IBM by Peskin (2003); Mittal and laccarino
(2005); Uhlmann (2005); Kim and Choi (2006); Mit
tal et al. (2008)) have been developed and shown to be
very effective in computing fluidparticle systems and
fluidstructure interaction problems. Lattice Boltzmann
method (LBM by Ladd and Verberg (2001)) has been
developed and effectively used for simulations of rigid
as well as deforming particles. Combination of the
DLM, direct forcing based IBM, and LatticeBoltzmann
methods (termed as Proteus) was recently developed
by Feng and Michaelides (2005). A secondorder ac
curate fixed grid method (PHYSALIS) was developed
by Zhang and Prosperetti (2005), which gives good so
lutions for spherical particles by using local spectral rep
resentations of the solution near a spherical boundary.
Taira and Colonius (2007) proposed a new im
plementation of the immersed boundary method to
achieve secondorder accuracy. They compared IBM
with fictitiousdomain based methods to point out sub
tle differences when the immersed objects are con
strained to undergo specified motion. In the fictiti
tious domain/DLM method (see Glowinski et al. (2001);
Patankar et al. (2000)), the entire fluidparticle domain
is assumed to be a fluid and the flow in the particle do
main is constrained to be a rigidbody motion through
rigidity constraint in terms of a stress or a force. Sharma
and Patankar (2005) proposed a fast technique to obtain
the rigidity constraint force that eliminated the need for
an iterative procedure to solve for the rigid body mo
tion. Recently, Veeramani et al. (2007) proposed a sim
ilar approach in the context of finiteelement methods.
Apte et al. (2009) extended the finitevolume based fic
titious domain approach by Sharma and Patankar (2005)
to large number of particles in complex turbulent flows
on colocated grids.
Majority of the above studies have been applied to
simulate rigid particulate flows with fluid particle den
sity ratios in the range of 0.1 10. Large density ra
tios are common in many practical applications involv
ing complex flows; for example coal particles in a oxy
coal boiler, aeolian particle transport, aerosol transport,
among others. Being able to simulate such flows with
fully resolved direct or largeeddy simulations is criti
cal. The fictitious domain approach with fast compu
tation of the rigidity constraint (Sharma and Patankar
(2005); Apte et al. (2009)) when applied to large density
ratio fluidparticle systems leads to numerical difficul
ties. Sharp gradients in density across the fluidparticle
interface can cause numerical 'ringing' of the solution.
In the present work, we extend this numerical approach
to account for fluidparticle systems over a range of den
sity ratios.
The paper is arranged as follows. A mathematical
formulation of the fictitious domain scheme is first de
scribed. Numerical issues with the original formulation
for high density ratios and potential remedies are dis
cussed next. The new approach is implemented in a
colocated grid finite volume method. The approach is
first validated for basic test cases to show good predic
tive capability. Namely, flow over a fixed sphere and
a NACA hydrofoil are first investigated. Next, freely
falling/rising spherical particles at different Reynolds
numbers are considered and compared with available ex
perimental data at low fluidparticle density ratios. Fi
nally, the approach is applied to a wide range of density
ratios 103 106 to show stable solution. Specifically,
interactions of a buoyant sphere with a stationary Gaus
sian vortex are simulated to show the capability of the
approach to study particlevortex interactions.
Mathematical Formulation
Let F be the computational domain which includes both
the fluid (F (t)) and the particle (Fp(t)) domains. Let
the fluid boundary not shared with the particle be de
noted by B and have a Dirichlet condition (generaliza
tion of boundary conditions is possible). For simplic
ity, let there be a single rigid object in the domain and
the body force be assumed constant so that there is no
net torque acting on the object. The basis of fictitious
domain based approach is to extend the NavierStokes
equations for fluid motion over the entire domain F in
clusive of immersed object. The natural choice is to as
sume that the immersed object region is filled with a
Newtonian fluid of density equal to the object density
(pp) and some fluid viscosity (fp). Both the real and
fictitious fluid regions will be assumed as incompress
ible and thus incompressibility constraint applies over
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
the entire region. In the numerical approach presented
by Sharma and Patankar (2005), the particle region is
identified by an indicator (color) function O which has
unit value inside the particle region and vanishes in the
fluid region. Owing to finite number of grid cells, the
boundary region of the particle typically is smeared with
0 < O < 1. The density field over the entire domain is
then given as,
P =pp + p (1 ).
The indicator function moves with the particle resulting
in
DO
D 0, (2)
Dt
where D/Dt is the material derivative. The fluid veloc
ity field is constrained by the conservation of mass over
the entire domain given as
Op
+ V (u) 0. (3)
at
The conservation of mass together with the indicator
function advection implies that for an incompressible
fluid, V u = 0 over the entire domain.
The momentum equation for fluid motion applicable
in the entire domain F in the conservative form is then
given by:
dpu)
" + V (puu)
Vp+
V (LF (Vu + (Vu))) + pg + f, (4)
where p is the density field, u the velocity vector, p the
pressure, tp the fluid viscosity, g the gravitational ac
celeration, and f is an additional body force that enforces
rigid body motion within the immersed object region
Fp. This can be done by following the fast procedure for
obtaining the rigidity constraint force as first proposed
by Sharma and Patankar (2005) for freely moving rigid
objects in laminar flows.
Sharp variations in density over a single grid cell can
lead to numerical instabilities when the momentum and
continuity equations are solved in the above conserva
tive form. This was shown to be a problem in volume of
fluid formulations for two immiscible fluids by Rudman
(1998). This problem is remedied by performing consis
tent flux constructions for mass and momentum fluxes at
the control volume faces as shown by Rudman (1998);
Bussmann et al. (2002); Raessi (2008). Inconsistencies
in flux calculations for mass and momentum leads to in
correct accelerations of fluids near interfaces leading to
numerical instability at high density ratios.
An alternative approach, that is commonly followed
in interface tracking schemes based on level set meth
ods (Osher and Fedkiw (2003); Osher and Sethian
(1988)), is to solve the above equations in a non
conservative form, wherein computation of density vari
ations across cell faces are not required especially for
colocated grid formulations. However, level set meth
ods suffer from the loss of mass owing to the non
conservative advection of the signeddistance function
especially for deforming interfaces. For the present
work on rigid body motion, the surface area of the in
terface between the fluid and particles remains constant
over time as interface deformation is not possible. Since
the rigid particles are advanced in a purely Lagrangian
frame their mass is conserved discretely in a numerical
formulation. When cast in the nonconservative form,
and making use of the fact that the fluid velocity is diver
gence free over the entire domain, the momentum con
servation equation becomes,
u+ V (u)
at
1V+
Vp
V (/ p p Vu+ (Vu)) + g + f.
p \p
In the present work, we solve the momentum equa
tion in the above form together with the incompressibil
ity constraint V u 0 over the entire domain. In order
to enforce that the material inside the immersed object
moves in a rigid fashion, a rigidity constraint force is re
quired that leads to a nonzero forcing function f. This
can be obtained using a fractional step scheme:
1. In this first step, the rigidity constraint force f in
equation 5 is set to zero and the equation together
with the incompressibility constraint is solved by
standard fractionalstep schemes over the entire do
main. Accordingly, a variable ceffiicient pressure
Poisson equation is derived and used to project the
velocity field onto an incompressible solution. The
obtained velocity field is denoted as u"+1 inside
the fluid domain and u inside the object.
2. The velocity field for a freely moving object is ob
tained in a second step by projecting the flow field
onto a rigid body motion. Inside the object:
(un+1
PP (u At
a i
To solve for u+ 1 inside the particle region we re
quire f. This is obtained by first finding the rigid
body motion inside the particle region. The veloc
ity field in the particle domain involves only trans
lation and angular velocities. Thus u is split into a
rigid body motion (uRBM U + (2 x r) and resid
ual nonrigid motion (u'). The translational and ro
tational components of the rigid body motion are
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
obtained by conserving the linear and angular mo
menta and are given as:
MpU
'PQ
Sppudx;
J r x ppfidx,
Pp
where Mp is the mass of the particle and Ip
fpP pp [(r r)I r r]dx is the moment of inertia
tensor. Knowing U and Q2 for each particle, the
rigid body motion inside the particle region uRBM
can be calculated.
3. The rigidity constraint force is then simply obtained
as f p(uRBM ui)/At. This sets un+1
uRBM in the particle domain. Note that the rigid
ity constraint is nonzero only inside the particle
domain and zero everywhere else. This constraint
is then imposed in a third fractional step.
Numerical Approach
The preceding mathematical formulation is imple
mented in a colocated, structured grid, three
dimensional flow solver based on a fractionalstep
scheme developed by Apte et al. (2009). Accordingly,
in the present work the fluidparticle system is solved
by a threelevel fractional step scheme. First the mo
mentum equations (without the pressure and the rigidity
constraint terms) are solved. The incompressibility con
straint is then imposed by solving a variablecuefficient
Poisson equation for pressure. Finally, the rigid body
motion is enforced by constraining the flow inside the
immersed object to translational and rotational motion.
The main steps of the numerical approach are given be
low.
Immersed Object Representation In the numer
ical implementation, we create small material volumes
of cubic shape that completely occupy the immersed ob
ject. Each material volume is assigned the properties of
the immersed object (e.g. density etc.). The shape of
the object can be reconstructed from these material vol
umes by computing an indicator or color function (with
value of unity inside the object and zero outside) on a
fixed background mesh used for flow solution. In this
work, the material volumes are forced to undergo rigid
motion, based on the translational and rotational veloci
ties of the object, resulting in no relative motion among
them. At each timestep the material volumes are ad
vanced to new locations. In the present approach, the
boundary of the object is represented in a stairstepped
fashion and it is straightforward to create the material
volumes using a boundingbox algorithm as described
in Apte et al. (2009).
(a) time
n+1 u,U,,unf,
n+1/2 x,,p,O,p
n u,,U,,u,f,
n1/2 X,p,e,p
Figure 1: Schematic of the variable storage in time
and space: (a) timestaggering, (b) three
dimensional variable storage, (c) cv and face
notation, (d) index notation for a given k
index in the z direction. The velocity fields
(ui, UN) are staggered in time with respect to
the volume fraction (e), density (p), and par
ticle position (Xi), the pressure field (p), and
the rigid body force (fi,p). All variables are
collocated in space at the centroid of a con
trol volume except the facenormal velocity
uN which is stored at the centroid of the faces
of the control volume.
Discretized Equations and Numerical Algorithm
Figure 1 shows the schematic of variable storage in time
and space. All variables are stored at the control vol
ume (cv) center with the exception of the facenormal
velocity UN, located at the face centers. The facenormal
velocity is used to enforce continuity equation. Capi
tal letters are used to denote particle fields. The time
staggering is done so that the variables are located most
conveniently for the timeadvancement scheme. We fol
low the collocated spatial arrangement for velocity and
pressure field (Kim and Choi (2000), Mahesh et al.
(2006)). Accordingly, the particle positions (Xi), den
sity (p), volume fraction (e), viscosity (p), and the pres
sure (p) are located at time level t 1/2 and "+1/ 2
whereas the velocity fields (Ui, UN, and Ui) and the rigid
body constraint force fiR are located at time level
t" and t"+1. This makes the discretization symmetric
in time, a feature important to obtain good conserva
tion properties. The semidiscretization of the governing
equations in each timestep is given below.
Step 1: Starting with a solution at tn and t 1/2, the
centroids of material volumes (Xi,M) representing im
mersed objects are first advanced explicitly.
Xn+1/2 xn1/2 x 1/2
iM i,P + Rij XjM
,n1/2 (
)3,P ) +
UMAt, (9)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where Xi,M is the position vector of the material vol
ume center, Xi,p is the position vector of the immersed
object centroid, Ui,M is the translation velocity, Q2,M is
the angular velocity, and At is the timestep. Here Rij
is the rotation matrix evaluated using particle locations
at t"1/2. The details of the particle update and the ro
tation matrix are similar to that presented by Apte et al.
(2009).
Step 2: Knowing the new positions of the mate
rial volumes and particle centroid, an indicator function
(color function) "~+1/2 is evaluated at the cvcenter
of the fixed background grid. We use a discrete delta
function developed by Roma et al. (1999) to compute
the color function. The color function is unity inside the
particle region and vanishes outside with smooth varia
tion near the boundary. This thus allows identification
of the particle on the background mesh. Details of the
interpolation between the material volume centers and
the cv center are similar to that presented by Apte et al.
(2009). The density and the viscosity are then calculated
over the entire domain as:
n+1/2
Pcv
U,+1/2
Sppon 1/2 +p
c' v + PF 1
,,'t/2 + fF (
S/2) (10)
c,+ ) (11)
where pp is the density of the immersed particle and
PF is the density of the surrounding fluid. Likewise
pp is dynamic viscosity of the fictitious fluid inside the
particle region, and pF is the dynamic viscosity of the
surrounding fluid. For particles with specified motion
(microvalves) pp is assumed equal to the fluid viscosity
(fP).
Step 3: Advance the momentum equations using the
fractional step method. First, obtain a predicted veloc
ity field over the entire domain. We advance the veloc
ity field from t" to t"+ The predicted velocity fields
may not satisfy the continuity or the rigidity constraints.
These are enforced later.
1 *n+1/2
At v facee UNAface
faces of cv
1 (1
p+1/2 V
pcv cv
n
n~ 2v
: *n1/2 T
fac f c e NjfaceAacece) + gi (12)
faces of cv /
where gi is the gravitational acceleration, Vcv is the vol
ume of the cv, Aface is the area of the face of a control
volume, Nj,face is the facenormal vector and
*n+1/2
*n+1/2
Tij,face
2 (Uface + Uface)
Vlc 2 1 ( + aj
xi )J face
face
In the above expressions, the velocities at the 'face' are
obtained by using arithmetic averages of the neighboring
cvs attached to the face. For the viscous terms, the ve
locity gradients in the direction of the momentum com
ponent are obtained implicitly using CrankNicholson
scheme. A centered discretization scheme is used for
spatial gradients. Evaluation of the pressure gradients at
the cv centers is explained below.
Step 4: Solve the variable coefficient Poisson equa
tion for pressure:
1 z
cs of c A ce
faces of cv
1 p n+1/2
Aface 12
faces of c fac
faces of cv f ace
(13)
where face is obtained using arithmetic averages of den
sity in the neighboring cvs. The facenormal velocity
u* and the facenormal pressure gradient are obtained
as:
1 (nb
uN = (inbr Ui,cv)Ni,face
r 2 r1
p 1n+/2
6N
n+1/2 n+1/2
Pnbr Pcv
Scv,nbr
where nbr represents neighboring cv associated with the
face of the cv, and IScv,nbr is the distance between
the two cvs. The variablecoefficient pressure equa
tion is solved using a BiConjugate gradient algorithm
by van der Vorst (2003).
Step 5: Reconstruct the pressure gradient at the cv
centers using density and facearea weighting first pro
posed by Ham and Young (2003)
1 p "+1/2
f+1/2 6xi
Faces of cv +/2 SN 12 i,faceAface
Pfac (14)
Spaces of cv INi,faceAfacel
Step 6: Update the cvcenter and facenormal veloci
ties to satisfy the incompressibility constraint:
p n+1/2
ii,cv = c At cxi (15)
6xi
Atn+1/2
N = At (16)
6N
The facenormal velocity field UN will satisfy the incom
pressibility constraint, however, the cvbased velocity
may not satisfy the rigidbody constraint inside the par
ticle region. Note that in the absence of any rigid body,
p = PF throughout the domain, and the algorithm re
duces to the standard fractional step scheme for single
phase, incompressible flow. The above velocity field
will then be denoted as .."+. In the presence of rigid
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
bodies, the following steps are performed to enforce the
rigidity constraint within the particle domain.
Step 7: First interpolate the velocity field i,cv from
the grid cvs to the material volume centroids to ob
tain U,,M using the kernel interpolation outlined in Apte
et al. (2009). Solve for the translational and rotational
velocity fields
N
MpUP = VMPMUM
M=1
N
IpQp = pMVM(r x UM),
M=1
where subscripts P and M denote the particle and
the material volume centroids respectively, VM is the
volume and pM the density of each material volume,
M.p = CM pMVM is the total mass of the particle,
Ip is the moment of inertia of the particle about the co
ordinate axes fixed to the particle centroid, and r is the
position vector of a point within the particle region with
respect to the particle centroid.. The moment of inertia
is given as
N
P PMVM [(r r)I r o r],
M=I
where I represents the identity matrix. The rigid body
motion is then obtained as:
UBM= U4 + Qp x (XM Xp). (20)
Step 8: Compute the rigidbody constraint force and
correct the velocity field to satisfy this constraint within
the particle region.
Fn+l
iM
UM T RBM,n+l
i Ui M
The force on the grid control volumes (f,,cv) is obtained
from Fi,M through a consistent interpolation scheme
scheme as presented by Apte et al. (2009). The velocity
field inside the particle region is then modified as:
un = i,cv + At f[c+. (22)
i,cv Uicv +
Numerical Test Cases
The above numerical algorithm is implemented in a par
allel, finite volume framework and validated for a num
ber of test cases: (i) flow over a fixed sphere and hy
drofoil, (ii) particle subjected to constant acceleration
for varying fluidparticle density ratios, and (iii) freely
falling/rising particles at low and high density ratios. Fi
nally, we simulate interactions of a buoyant sphere with
a stationary Gaussian vortex at different density ratios.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 1: Drag coefficient Cd for flow over a sphere at different Reynolds number
Study Re,
20 50 100 150 300 350
Current 2.633 1.550 1.101 0.907 0.686 0.649
Mittal (1999) 1.57 1.09 0.62
Mittal et al. (2008) 1.08 0.88 0.68 0.63
Clift et al. (1978) 2.61 1.57 1.09 0.89 0.684 0.644
Johnson and Patel (1999) 1.57 1.08 0.9 0.629
Marella et al. (2005) 1.56 1.06 0.85 0.621
Kim et al. (2001) 1.087 0.657
Figure 2: Isosurface of A 0.008 for flow over a stationary sphere at Rep
line shows the symmetry axis of the structure in this plane.
Flow over a Sphere To evaluate the accuracy of the al NACA 0008 HI
gorithm for threedimensional configurations, flow over flow over a stati
a fixed sphere at different Reynolds number is first eval to show the apl
uated and compared with published data. A sphere regular geometi
of diameter d = 1.10 mm is placed in a domain of lenging for this
15d x 15d x 15d. The sphere is located at x = 5d and foil's leading e
y = z = 7.5d. The grid used is 128 x 128 x 128. The These features
grid is nonuniform but it is refined and uniform around tics associated
the sphere forming a patch of 1.5d x 1.5d x 1.5d. There any flow solver
are approximately 26 grid points along a diameter of the
A two dimen
sphere. For comparison Mittal et al. (2008) used a do C n
Cartesian grid
main of 16d x 15d x 15d and a grid of 192 x 120 x 120 is
direction is gen
for their highest Reynolds number of 350 and Marella in the streamwi
et al. (2005) employed a 130 x 110 x 110 mesh on rese
respectively. Tt
a 16d x 15d x 15d domain. The fluid properties are wise direction
a ,wise direction
p 1 kg/m and the viscosity p 10 kg/m.s. The where c is the
x direction is slightly moved towards the inlet in order te h
the hydrofoil a
to increase the size of the domain in the wake. Also
bic with spacin
the density of gridpoints is increased in the wake of s
et al. (2008) us
the sphere in order to properly resolve the flow. Table 1
directions for a
compares the predicted drag coefficients with other pub solution. Kun
lished data showing very good predictions. grid and a two
grid and a two
points in the hy
Figure 2 shows the vortical structure represented by points in the hy
Fiur number, Re =
the eigenvalues of the velocity gradient tensor A in the number
of attack c=
xz plane for Rep = 350. Qualitatively the plots show
tions are run f
very similar structures as shown by Mittal (1999). This
Relevant comp
snapshot shows the vortex ring in the near wake of the R
table 2.
sphere. Another important feature is the symmetry axis
shown in the xz plane which has been observed experi At this Reyn
mentally. dimensional an
= 350 in the xz plane. The dashdotted
rdrofoil The computational solution of
onary NACA 0008 hydrofoil is presented
plicability of the present approach to ir
ries. The thin hydrofoil geometry is chal
method because of the small radius at the
dge and the sharp tip at the trailing edge.
lead to distinct lift and drag characteris
with the geometry, and it is important for
to properly resolve them.
isional flow is assumed, and a block type
with periodic boundaries in the spanwise
berated using 504 x 200 x 4 grid points
se, crossstream, and spanwise directions
ie domain itself is 13c long in the stream
and 8c wide in the crossstream direction,
hydrofoil chord length. The grid around
nd in the near wake is uniform and cu
g equal to c/400. By comparison, Mittal
ed 926 x 211 grid points in the x and y
cartesian grid based immersed boundary
z and Kroo (2001) used a body fitted C
dimensional solver with 256 x 64 grid
drofoil plane. The chord based Reynolds
pUoc/i, is fixed at 2000. Two angles
0 and a 40 are used. The simula
)r 15 non dimensional time units, tug
utational parameters are summarized in
lolds number, the flow is laminar, two
d steady. Figure 3 shows the contours of
Table 2: Parameters for NACA0008 hydrofoil.
a Lu LY N. Ny N,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 3: Comparison of steady state drag and lift coeffi
cients with other 2D computational results
13c 8c 504 200 4 400
vorticity for a 00 and 40 after the flow has reached
steady state. At a 0, the flow is symmetric, and the
wake is quite narrow. At a = 4, the wake is widened,
and separation is indicated on the suction side by detach
ment of the high vorticity region from the surface.
0.25 
0.25 k
0 0.5 x/ 1.5 2
(a) a = 0
.,
0.25 F
0.25 F
o U.b X/C 1
1.b5
(b) a = 40
Figure 3: Contours of vorticity (out of the page) for flow
around a NACA 0008 hydrofoil at Re = 2,000
The lift and drag coefficients for the hydrofoil are
given as:
Study
Present Results
Mittal et al. (2008)
Kunz and Kroo (2001)
a 0
CD
0.078
0.078
0.076
a = 4
CD
0.083
0.081
0.080
CL
0.270
0.273
0.272
effects are negligible, this case has a simple analytical
solution of u = uo + at; where u and uo are the in
stantaneous and initial velocity fields, respectively and t
is time. Initially, a particle of diameter 5 mm is placed
at the center of a cubic box of 3 cm. Uniform cubic
grids of 100 x 100 x 100 grid points is used. The fluid
density is set to 1 kg/m3 and the viscosity and gravi
tational acceleration are turned off. The particle starts
from rest and is subjected to uniform acceleration of
(40, 40, 40) m/s2. The particle density is varied
over a wide range 103, 104 and 106 and the distance trav
eled by the particle is compared to the analytical solution
of S luolt + 1/2at2. Figure 4 shows the temporal
evolution of relative error in the distance traveled by the
particle compared to the exact solution (I sums exact )
Sexact t. 0.01
over 1000 iterations with fixed time step of 10 ps. It
is observed that the error remains small for high density
ratios.
FD .
CD '
pU ;
FL
SpU
The steady state values of CL and CD are taken after
15 time units and are compared to the published data
in table 3. The present steady state values are in good
agreement with the previous computational studies.
Sphere subjected to uniform acceleration To test
the stability of the numerical algorithm for high density
ratio between the particle and the fluid, we consider mo
tion of a spherical particle subjected to uniform accel
eration (a) in a closed box. If viscous and gravitational
time, s
Figure 4: Temporal evolution of the relative error in dis
tance traveled by the sphere under uniform ac
celeration: A pp/pf = 103, pp/pf = 104,
E pp/pf = 106.
Sphere rising in inclined channel A lighterthan
fluid particle rising in an inclined channel is considered
next. The simulation is conducted with a fluid den
sity of pf 1115 kg/m3, a particle density of pp
2,000 00, 4"
cU: 30 25 20 15 10 5 0 5 10 15 20 25 30
ac/U 30 25 20 15 10 5 0 5 10 15 20 25 30
x
x
Figure 5: Comparison of experimental and numerical
simulations  Lomholt et al. (2002) with
the current simulation . Figure 5(a) shows
the particle trajectory inside the domain (the
 line shows the initial trajectory due only
to the effect of gravity), 5(b) the velocity of
the particle in the lateral direction and 5(c) the
velocity on the vertical direction. The particle
position is expressed in [m] and the velocities
are expressed in [m/s].
1081 kg/nm3, and a fluid viscosity of v 3.125 nmm2/s.
The Reynolds number Rest" based on the Stokes set
tling velocity W is defined as:
Restokes 2dpW 4 p
p i 9v2 pf
where g = 9.82 n/s2 is the gravitational acceleration,
and dp = 2 rmm is the diameter of the particle. The
channel is inclined at an angle of 8.23 with the vertical.
This is simulated by adding components of gravitational
forces in the horizontal and vertical directions.
The computational domain consists of a rectangular
box with dimensions 10 mm in the x direction, 80 rmm
in the y direction and 40 mm in the z direction. The
grid is Cartesian and uniform over the domain with
40 x 320 x 160 grid points, respectively in the x, y and
z directions so that A 0.25 x 10 3 m. The parti
cle is injected at : = 1.4 mm, y = 1.0 mm and
z 20.0 nmm. Simulation results for a Reynolds num
ber of Re toke = 13.6 are compared with experimen
tal and numerical data from Lomholt et al. (2002). As
illustrated in Figure 5, the numerical simulation exhibits
excellent agreement with both experimental and numer
ical results. Buoyancy forces cause the particle to rise
and travel alongside the right wall of the domain. Ulti
mately, the particle follows the right wall without touch
ing it, keeping a very thin lubrication layer between the
particle and the wall.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Freely falling sphere We consider the problem of a
single sphere falling under gravity in a closed container.
The particle density is (pp 1120 kg/m3) and the di
ameter is (15 mm). The sphere is settling in a box of
dimensions 10 x 10 x 16 cm3. The particle is released
at a height H = 12 cm from the bottom of the box. The
boundaries of the box are treated as noslip walls. The
fluid properties are varied to obtain different Reynolds
numbers based on the terminal velocity of the particle.
The simulation conditions correspond to the experimen
tal study by ten Cate et al. (2002). Table 4 provides de
tailed information about the parameters used in this test
problem. We simulate the above cases on a fine uniform
Table 4: Parameters for freely falling sphere.
PF (kg/m3) pF (10 3Ns/m2)
970 373
965 212
962 113
960 58
Re
1.5
4.1
11.6
31.9
grid of 100 x 100 x 160 points with a grid resolution of
A 1 mm. This provides around 15 grid points inside
the particle domain. The material volumes are cubical
with 5, where AM is the size of the material
Am
volume. Accordingly, there are around 75 material vol
umes along the diameter of the spherical particle in each
direction. A uniform timestep (At 0.5 ms) is used
for all cases.
Figures 6ab show the comparison of the time evolu
tion of particle settling velocity and position at different
times obtained from the numerical simulations with the
experimental data by ten Cate et al. (2002). The simu
lation predictions for both the particle velocity and the
particle position show good agreement with the experi
mental data. The slowing of the particle towards the end
of the simulation are to due to the presence of the bottom
wall.
Falling sphere with high density ratio To test the ac
curacy of the numerical algorithm for high density ratio
between the particle and the fluid, we consider the same
problem of falling sphere, except vary the density of the
particle keeping the fluid density fixed at pf 1 kg/n3.
The particle density (pp) is varied over a wide range
10, 100, 500 and 1000. To keep the terminal velocity
of the particle small the viscosity of the fluid is set at
0.06 [kg/n.s]. The diameter of the particle, the domain
size, computational grid and release point at same as that
of the above case.
Mordant and Pinton (2000) performed experiments
on freely falling spherical particles in a large water tank
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
E0.01
_0.
0
a,
>0.1;
(a) Fall Velocity
time, s"
(b) Normalized Height
Figure 6: Comparison with the experimental data of the
sphere fall velocity and the normalized height
from the bottom wall for different Reynolds
numbers: (Symbols: experiment by ten Cate
et al. (2002), lines: present simulation) Here
H = 05dp where h is the height of the
dp
sphere center from the bottom wall and dp is
the particle diameter.
02
02 04 06 08 1 12 14 16 18 2
Figure 7: Temporal evolution of a spherical particle ris
ing in a water column: U* U/Uterminal,
T* t/T95: solid line by Mordant and Pinton
(2000); pp/pf 500 (o); pp/pf 100 (0);
pp/pf 10 (>) present simulations.
for various density ratios (maximum density ratio con
sidered was pp/pf 14.6). They showed that for
small particles falling in a large tank (that is, for small
values of the ratio of particle diameter to tank width
dp/L ~ 0.005) the temporal evolution of the particle
velocity can be well predicted by the the curve:
U* = 1 exp  ,
( T95
where U* is the velocity of the particle normalized by its
terminal velocity, T95 is the time it takes for the sphere
to reach 95% of its terminal velocity, and t is time. We
compare the temporal evolution of the rising spherical
particle to this curve in Figure 7. The temporal evolu
tion of the particle velocity is qualitatively similar to the
exponential curve. The domain size in simulations is
small (dp/L 0.15) and thus wall effects become im
portant. This test case confirms the stability of the nu
merical solver when applied to largedensity ratio fluid
particle systems.
Buoyant sphere in a Gaussian vortex
We now consider the entrainment of a single buoyant
sphere in a stationary Gaussian vortex as shown in fig
ure 8(a). The buoyant sphere is released in the vicinity
of the vortex with core radius, re and initial circulation,
Fo. With sufficient vortex strength, the sphere gets en
trained into the vortex after circling around several times
and reaching a settling location with relative coordinates
rs, Os measured from the vortex center. For the Gaussian
vortex, there is no radial velocity component, and the the
tangential velocity is expressed as
Ue(r) 1 e2
27r
The vorticity and maximum tangential velocity (occurs
at r = r,) are given by
w(r) Frolt I(r/o)2
7Tr
Fo
2 27r
where ri and 92 are constants. This flow has been used
previously by Oweis et al. (2005), as a model for wingtip
vortices in their study of bubble capture and cavitation
inception. Variables relevant for the setup of this test
case are summarized in table 5. The domain size is ap
proximately 7 re x 7 re x 0.4 r,. A slip condition is
imposed at boundaries in the X and Y directions, and
the domain is periodic in the Z direction. The Carte
sian grid is refined in the area of the vortex core with
a cubic spacing of Aco, 0.1 mm and has 450 x
450 x 50 grid points in the X, Y and Z directions. The
velocity field of equation 26 is applied as an initial con
dition everywhere in the domain, creating a clockwise
vortex with initial strength Fo 0.04m2s 1. During
the first timestep, a single, 1.1 mm diameter spherical
particle is released at r r, 0 0. We assume typ
ical properties of water for the fluid phase and simulate
3 different particle densities corresponding to pf/pp =
100, 500, and 1000. In each case, the simulation is ad
vanced using a timestep of At 50 ps up to a time of
1 s. Figure 8(b) shows the entrainment trajectories of
Table 5: Parameters for the Gaussian vortex case.
re [mm] 11.45
91 1.27
9/2 0.715
L1, L ; L, [mm] 80,80,5
Acore [mm] 0.1
Ngrid 425 x 425 x 50
Fo [m2/s] 0.04
dp [m] 1.1
pf [kg/m3] 1,000
pp [kg/mn3] 1, 2, 10
if [kg/m.s] 0.001
the spheres for each density ratio considered. Because
of the high density ratios, the spheres do not follow the
fluid streamlines and instead are attracted toward the up
per right hand side of the inner core. This is due to
the lift and added mass effects effects similarly observed
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 6: Average settling coordinates of the sphere in
the vortex core as shown in figure 8(a)
Density Ratio rs/rc 0s [rad]
pf/pp 100 0.30 0.005
pf/pp = 500 0.45 0.065
pf/pp =1,000 0.57 0.106
by Mazzitelli et al. (2003); Climent et al. (2007) for bub
bles entrained in vortical structures. As density ratio and
relative buoyancy force are increased, a more direct path
is taken to the settling location by the sphere. The set
tling coordinates are tabulated in table 6 for each case.
With increased density ratio, the spheres reach equilib
rium at a larger radius and more negative angle from the
horizontal. For pf/pp equal to 100 and 500, there is little
relative motion of the sphere once it becomes entrained.
At pf/pp = 1000, a much different dynamic exists be
cause the sphere does not stay in one location. Rather it
circles a point centered at (rs/r6, 0s) = (0.57, 0.106).
This motion causes a strong and highly unsteady wake to
develop behind the sphere as is shown in figure 9. Strong
vortical structures are periodically shed from the sphere,
and advected clockwise around the vortex core back to
the sphere.
Conclusion
A numerical formulation for fully resolved simulations
of freely moving rigid particles for a wide range of fluid
particle density ratios is developed based on a fictitious
domain method. In this approach, the entire computa
tional domain is first treated as a fluid of varying density
corresponding to the fluid or particle densities in their
respective regions. The incompressibility and rigidity
constraints are applied to the fluid and particle regions,
respectively, by using a fractional step algorithm. The
momentum equations are recast to avoid computations
of density variations across the fluidparticle interface.
The resultant equations are discretized using symmetric
central differences in space and time. Use of consistent
interpolations between the particle material volumes and
the background grid and parallel implementation of the
algorithm facilitates accurate and efficient simulations
of large number of particles. Implementation of this ap
proach in finitevolume, colocated grid based numerical
solver is presented. The numerical approach is applied
to simulate flow over a fixed immersed objects (sphere
and NACA0008 hydrofoil) at different Reynolds num
bers, freely falling/rising spheres for a wide range of
density ratios and compared with published experimen
tal or numerical results to show good agreement. Fi
nally, the approach was applied to simulate freely mov
ing buoyant spheres in a stationary Gaussian vortex for
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(a) (b)
Figure 8: Entrainment of a buoyant sphere by a Gaussian vortex of core radius re: (a) Problem schematic, showing
the definition of the settling coordinates (rs, 0s), (b) Entrainment trajectories for the three cases simulated,
(c) The motion of the sphere at the settling location. pf /p = 100, A pf/pp = 500, p /pp = 1000
, Cx
'
(a) t= 0
(b) t = 0.05s
(c) t = 0.ls
r
r
I
I
I
r
I
(d) t = 0.5s
(e) t = 0.6s
(f) t = 0.7s
Figure 9: Snapshots of fluid vorticity magnitude in the vortex core for the high density ratio case (pf /p = 1000)
Contours are from 20 to 250 m2/s. Isosurface of 350 m2/s is also shown.
a range of density ratios to investigate the ability of the
solver to capture lift and addedmass effects on lighter
than fluid particles in vortical flows. For large vortex
strengths, the lighterthan fluid particles are attracted to
wards the core of the vortex and settle further away from
the center in the presence of gravity as expected.
Acknowledgements
This work was supported under US Department of En
ergy's Faculty Research Participation Program by Oak
Ridge Institute for Science and Research (ORISE) at Na
tional Energy Technology Laboratory, Albany. Simu
lations were performed on the high performance com
puting cluster at Oregon State University as well as at
TACC's Lonestar.
References
S.V. Apte, M. Martin, and N.A. Patankar. A numeri
cal method for fully resolved simulation (FRS) of rigid
particleflow interactions in complex flows. Journal of
Computational Physics, 228(8):27122738, 2009.
M. Bussmann, DB Kothe, and JM Sicilian. Modeling
high density ratio incompressible interfacial flows. In
Proc. FEDSM, volume 2, pages 200231125, 2002.
R. Clift, J.R. Grace, and M.E. Weber. Bubbles, drops,
and particles. Academic Press, New York, 1978.
E. Climent, M. Simonnet, and J. Magnaudet. Prefer
ential accumulation of bubbles in CouetteTaylor flow
patterns. Physics of Fluids, 19:083301, 2007.
Z.G. Feng and E.E. Michaelides. Proteus: a direct forc
ing method in the simulations of particulate flows. Jour
nal of Computational Physics, 202(1):2051, 2005.
R. Glowinski, TW Pan, TI Hesla, DD Joseph, and J. Pe
riaux. A fictitious domain approach to the direct numeri
cal simulation of incompressible viscous flow past mov
ing rigid bodies application to particulate flow. Journal
of Computational Physics, 169(2):363426, 2001.
F. Ham and YN. Young. A Cartesian Adaptive Level Set
Method for TwoPhase Flows. Annu. Research Briefs,
2003: Center for Turbulence Research, 2003.
TA Johnson and VC Patel. Flow past a sphere up to a
reynolds number of 300. Journal of Fluid Mechanics,
378:1970, 1999.
D. Kim and H. Choi. A secondorder timeaccurate fi
nite volume method for unsteady incompressible flow
on hybrid unstructured grids. Journal of Computational
Physics, 162(2):411428, 2000.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
D. Kim and H. Choi. Immersed boundary method for
flow around an arbitrarily moving body. Journal of
Computational Physics, 212(2):662680, 2006.
J. Kim, D. Kim, and H. Choi. An immersedboundary
finitevolume method for simulations of flow in complex
geometries. Journal of Computational Physics, 171(1):
132150,2001.
PJ Kunz and I. Kroo. Analysis and design of airfoils for
use at ultralow Reynolds numbers. Progress in Astro
nautics and Aeronautics, 195:3560, 2001.
A.J.C. Ladd and R. Verberg. LatticeBoltzmann simula
tions of particlefluid suspensions. J. Statist. Phys., 104:
11911251, 2001.
S. Lomholt, B. Stenum, and MR Maxey. Experimental
verification of the force coupling method for particulate
flows. International Journal of Multiphase Flow, 28(2):
225246, 2002.
K. Mahesh, G. Constantinescu, S. Apte, G. Iaccarino,
F. Ham, and P. Moin. LargeEddy Simulation of React
ing Turbulent Flows in Complex Geometries. Journal of
Applied Mechanics, 73:374, 2006.
S. Marella, S. Krishnan, H. Liu, and HS Udaykumar.
Sharp interface cartesian grid method i: An easily im
plemented technique for 3d moving boundary computa
tions. Journal of Computational Physics, 210(1):131,
2005.
I.M. Mazzitelli, D. Lohse, and F. Toschi. On the rele
vance of the lift force in bubbly turbulence. Journal of
Fluid Mechanics, 488:283313, 2003.
R. Mittal. A fourierchebyshev spectral collocation
method for simulating flow past spheres and spheroids.
International Journal for Numerical Methods in Fluids,
30(7):921937, 1999.
R. Mittal and G. Iaccarino. Immersed boundary meth
ods. Annual Review of Fluid Mechanics, 37(1):239261,
2005.
R. Mittal, H. Dong, M. Bozkurttas, FM Najjar, A. Var
gas, and A. Von Loebbecke. A versatile sharp interface
immersed boundary method for incompressible flows
with complex boundaries. Journal of computational
physics, 227(10):48254852, 2008.
N. Mordant and J.F. Pinton. Velocity measurement of
a settling sphere. The European Physical Journal B, 18
(2):343352, 2000.
S. Osher and R.P. Fedkiw. Level set methods and dy
namic implicit surfaces. Springer Verlag, 2003.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
S. Osher and J.A. Sethian. Fronts propagating Z. Zhang and A. Prosperetti. A secondorder method for
with curvaturedependent speed: algorithms based on threedimensional particle simulation. Journal of Com
HamiltonJacobi formulations. Journal of computa putational Physics, 210(1):292324, 2005.
tional physics, 79(1):1249, 1988.
GF Oweis, IE van der Hout, C. lyer, G. Tryggvason, and
SL Ceccio. Capture and inception of bubbles near line
vortices. Physics of Fluids, 17:022105, 2005.
NA Patankar, P. Singh, DD Joseph, R. Glowinski, and
T.W. Pan. A new formulation of the distributed La
grange multiplier/fictitious domain method for particu
late flows. International Journal of Multiphase Flow,
26(9):15091524, 2000.
C.S. Peskin. The immersed boundary method. Acta Nu
merica, 11:479517, 2003.
M. Raessi. A level set based method for calculating flux
densities in twophase flows. Annual Research Briefs,
pages 467478, 2008.
AM Roma, CS Peskin, and MJ Berger. An Adaptive
Version of the Immersed Boundary Method. Journal of
Computational Physics, 153(2):509534, 1999.
M. Rudman. A volumetracking method for incompress
ible multifluid flows with large density variations. Inter
national Journal for Numerical Methods in Fluids, 28
(2):357378, 1998.
N. Sharma and N.A. Patankar. A fast computation tech
nique for the direct numerical simulation of rigid partic
ulate flows. Journal of Computational Physics, 205(2):
439457, 2005.
K. Taira and T. Colonius. The immersed boundary
method: A projection approach. Journal of Computa
tional Physics, 225(2):21182137, 2007.
A. ten Cate, CH Nieuwstad, JJ Derksen, and HEA
Van den Akker. Particle imaging velocimetry exper
iments and latticeBoltzmann simulations on a single
sphere settling under gravity. Physics ofFluids, 14:4012,
2002.
M. Uhlmann. An immersed boundary method with di
rect forcing for the simulation of particulate flows. Jour
nal of Computational Physics, 209(2):448476, 2005.
HA van der Vorst. Iterative Krylov Methods for Large
Linear Systems. Cambridge University Press, 2003.
C. Veeramani, PD Minev, and K. Nandakumar. A ficti
tious domain formulation for flows with rigid particles:
A nonLagrange multiplier version. Journal of Compu
tational Physics, 224(2):867879, 2007.
