Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 7.4.1 - A Variable-Density Fictitious-Domain Method for Fully Resolved Simulation of High-Density Ratio Fluid-Particle Systems
Full Citation
Permanent Link:
 Material Information
Title: 7.4.1 - A Variable-Density Fictitious-Domain Method for Fully Resolved Simulation of High-Density Ratio Fluid-Particle Systems Computational Techniques for Multiphase Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Apte, S.V.
Finn, J.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
Subject: fully resolved simulations
fictitious domain methods
high-density ratio
freely moving particles
particle-vortex interactions
Abstract: A numerical scheme for fully resolved simulation of fluid-particle systems with freely moving rigid particles is developed. The approach is based on a fictitious domain method wherein the entire fluid-particle 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 finite-volume solver is based on a co-located grid incompressible but variable density flow. The incompressibility constraint is imposed by solving a variable-coefficient pressure equation giving rise to a stable scheme for high density ratio fluid-particle 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 fluid-particle density ratios.
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: VID00182
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 741-Apte-ICMF2010.pdf

Full Text

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

A Variable-Density Fictitious-Domain Method for Fully Resolved Simulation of
High-Density Ratio Fluid-Particle Systems

Sourabh V. Apte* and Justin Finn*

School of Mechanical Industrial and Manufacturing Engineering
Oregon State University, Corvallis, OR 97333, USA and
Keywords: Fully resolved simulations, fictitious domain methods, high-density ratio,
freely moving particles, particle-vortex interactions.

A numerical scheme for fully resolved simulation of fluid-particle systems with freely moving rigid particles is
developed. The approach is based on a fictitious domain method wherein the entire fluid-particle 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 finite-volume solver is based on a co-located grid incompressible but variable density flow. The
incompressibility constraint is imposed by solving a variable-coefficient pressure equation giving rise to a stable
scheme for high density ratio fluid-particle 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
fluid-particle density ratios.

Fully resolved simulations (FRS) of fluid-particle sys-
tems, wherein all scales associated with the fluid and
particle motion are completely resolved, are of great im-
portance for understanding particle-turbulence 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 fluid-particle systems and
fluid-structure 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 Lattice-Boltzmann
methods (termed as Proteus) was recently developed
by Feng and Michaelides (2005). A second-order 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 second-order accuracy. They compared IBM
with fictitious-domain 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 fluid-particle domain
is assumed to be a fluid and the flow in the particle do-
main is constrained to be a rigid-body 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 finite-element methods.
Apte et al. (2009) extended the finite-volume based fic-
titious domain approach by Sharma and Patankar (2005)
to large number of particles in complex turbulent flows

on co-located 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 large-eddy 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 fluid-particle systems leads to numerical difficul-
ties. Sharp gradients in density across the fluid-particle
interface can cause numerical 'ringing' of the solution.
In the present work, we extend this numerical approach
to account for fluid-particle 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
co-located 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 fluid-particle 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 particle-vortex 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 Navier-Stokes
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
D- 0, (2)
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

+ V (u) 0. (3)
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:

" + V (puu)


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
co-located grid formulations. However, level set meth-
ods suffer from the loss of mass owing to the non-
conservative advection of the signed-distance 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 non-conservative 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)


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 non-zero 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 fractional-step 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:

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 non-rigid 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:




J r x ppfidx,

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 non-zero 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 co-located, structured grid, three-
dimensional flow solver based on a fractional-step
scheme developed by Apte et al. (2009). Accordingly,
in the present work the fluid-particle system is solved
by a three-level 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 variable-cuefficient
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-

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 time-step the material volumes are ad-
vanced to new locations. In the present approach, the
boundary of the object is represented in a stair-stepped
fashion and it is straightforward to create the material
volumes using a bounding-box 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,

n-1/2 X,p,e,p

Figure 1: Schematic of the variable storage in time
and space: (a) time-staggering, (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 face-normal 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 face-normal
velocity UN, located at the face centers. The face-normal
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 time-advancement 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 semi-discretization of the governing
equations in each time-step 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 xn-1/2 x 1/2
iM i,P + Rij XjM

,n-1/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 time-step. 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.
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 cv-center
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:


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
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~ 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 face-normal vector and



2 (Uface + Uface)

Vlc 2 1 ( + aj

xi )J 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 Crank-Nicholson
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

where face is obtained using arithmetic averages of den-
sity in the neighboring cvs. The face-normal velocity
u* and the face-normal pressure gradient are obtained

1 (nb
uN = -(inbr Ui,cv)Ni,face
r 2 r1

p 1n+/2

n+1/2 n+1/2
Pnbr Pcv

where nbr represents neighboring cv associated with the
face of the cv, and IScv,nbr| is the distance between
the two cvs. The variable-coefficient pressure equa-
tion is solved using a Bi-Conjugate gradient algorithm
by van der Vorst (2003).
Step 5: Reconstruct the pressure gradient at the cv
centers using density and face-area 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 cv-center and face-normal veloci-
ties to satisfy the incompressibility constraint:

p n+1/2
ii,cv = c At cxi (15)
N = At (16)
The face-normal velocity field UN will satisfy the incom-
pressibility constraint, however, the cv-based velocity
may not satisfy the rigid-body 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

IpQp = pMVM(r x UM),

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

P PMVM [(r r)I r o r],

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 rigid-body constraint force and
correct the velocity field to satisfy this constraint within
the particle region.


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 fluid-particle 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: Iso-surface 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 three-dimensional 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 non-uniform 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 grid-points 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 dash-dotted

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


(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:


Present Results
Mittal et al. (2008)
Kunz and Kroo (2001)

a 0

a = 4


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/2|a|t2. Figure 4 shows the temporal
evolution of relative error in the distance traveled by the
particle compared to the exact solution (I sum-s 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

FD .
CD '
pU ;


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 lighter-than-
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



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 no-slip 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


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
volume. Accordingly, there are around 75 material vol-
umes along the diameter of the spherical particle in each
direction. A uniform time-step (At 0.5 ms) is used
for all cases.
Figures 6a-b 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

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



(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 = 0-5dp where h is the height of the
sphere center from the bottom wall and dp is
the particle diameter.


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 large-density 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

The vorticity and maximum tangential velocity (occurs
at r = r,) are given by

w(r) Frolt I(r/o)2

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.

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 fluid-particle 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 finite-volume, co-located 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 =



(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 added-mass effects on lighter-
than fluid particles in vortical flows. For large vortex
strengths, the lighter-than 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.

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.


S.V. Apte, M. Martin, and N.A. Patankar. A numeri-
cal method for fully resolved simulation (FRS) of rigid
particle-flow interactions in complex flows. Journal of
Computational Physics, 228(8):2712-2738, 2009.

M. Bussmann, DB Kothe, and JM Sicilian. Modeling
high density ratio incompressible interfacial flows. In
Proc. FEDSM, volume 2, pages 2002-31125, 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 Couette-Taylor 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):20-51, 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):363-426, 2001.

F. Ham and YN. Young. A Cartesian Adaptive Level Set
Method for Two-Phase 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:19-70, 1999.

D. Kim and H. Choi. A second-order time-accurate fi-
nite volume method for unsteady incompressible flow
on hybrid unstructured grids. Journal of Computational
Physics, 162(2):411-428, 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):662-680, 2006.

J. Kim, D. Kim, and H. Choi. An immersed-boundary
finite-volume method for simulations of flow in complex
geometries. Journal of Computational Physics, 171(1):

PJ Kunz and I. Kroo. Analysis and design of airfoils for
use at ultra-low Reynolds numbers. Progress in Astro-
nautics and Aeronautics, 195:35-60, 2001.

A.J.C. Ladd and R. Verberg. Lattice-Boltzmann simula-
tions of particle-fluid suspensions. J. Statist. Phys., 104:
1191-1251, 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):
225-246, 2002.

K. Mahesh, G. Constantinescu, S. Apte, G. Iaccarino,
F. Ham, and P. Moin. Large-Eddy 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):1-31,

I.M. Mazzitelli, D. Lohse, and F. Toschi. On the rele-
vance of the lift force in bubbly turbulence. Journal of
Fluid Mechanics, 488:283-313, 2003.

R. Mittal. A fourier-chebyshev spectral collocation
method for simulating flow past spheres and spheroids.
International Journal for Numerical Methods in Fluids,
30(7):921-937, 1999.

R. Mittal and G. Iaccarino. Immersed boundary meth-
ods. Annual Review of Fluid Mechanics, 37(1):239-261,

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):4825-4852, 2008.

N. Mordant and J.F. Pinton. Velocity measurement of
a settling sphere. The European Physical Journal B, 18
(2):343-352, 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 second-order method for
with curvature-dependent speed: algorithms based on three-dimensional particle simulation. Journal of Com-
Hamilton-Jacobi formulations. Journal of computa- putational Physics, 210(1):292-324, 2005.
tional physics, 79(1):12-49, 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):1509-1524, 2000.

C.S. Peskin. The immersed boundary method. Acta Nu-
merica, 11:479-517, 2003.

M. Raessi. A level set based method for calculating flux
densities in two-phase flows. Annual Research Briefs,
pages 467-478, 2008.

AM Roma, CS Peskin, and MJ Berger. An Adaptive
Version of the Immersed Boundary Method. Journal of
Computational Physics, 153(2):509-534, 1999.

M. Rudman. A volume-tracking method for incompress-
ible multifluid flows with large density variations. Inter-
national Journal for Numerical Methods in Fluids, 28
(2):357-378, 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):
439-457, 2005.

K. Taira and T. Colonius. The immersed boundary
method: A projection approach. Journal of Computa-
tional Physics, 225(2):2118-2137, 2007.

A. ten Cate, CH Nieuwstad, JJ Derksen, and HEA
Van den Akker. Particle imaging velocimetry exper-
iments and lattice-Boltzmann simulations on a single
sphere settling under gravity. Physics ofFluids, 14:4012,

M. Uhlmann. An immersed boundary method with di-
rect forcing for the simulation of particulate flows. Jour-
nal of Computational Physics, 209(2):448-476, 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 non-Lagrange multiplier version. Journal of Compu-
tational Physics, 224(2):867-879, 2007.

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