Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 12.1.1 - Simulation of Bubble and Particle Interactions with Vortical Flows Using a Discrete Element Model
Full Citation
Permanent Link:
 Material Information
Title: 12.1.1 - Simulation of Bubble and Particle Interactions with Vortical Flows Using a Discrete Element Model Bubbly Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Finn, J.
Apte, S.V.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
Subject: Euler-Lagrange
bubbly flows
discrete element model
two way coupling
Abstract: Discrete Element Model (DEM) is used to simulate bubble and particle interactions with vortical flows. In this Euler-Lagrange approach, volume-averaged Navier-Stokes equations are solved using a co-located grid finite volume method. The dispersed phase is assumed spherical and is modeled by tracking the centroids of the particles/bubbles in a Lagrangian frame. In addition to interphase momentum exchange, the volume occupied by the disperse phase is accounted for through fluid void fractions in the momentum and continuity equations. Interactions of the dispersed phase with vortical flows are investigated for both dilute and dense loadings. First, flow generated by the disperse phase motion are simulated to reproduce experimentally observed vortical flow features: (i) the breakup of a viscous blob of particles falling under the influence of gravity and (ii) the transient migration of a buoyant bubble plume. Finally, bubble entrainment and interaction with traveling vortex tube under dilute loadings are simulated in a twodimensional approximation of the experiments by Sridhar and Katz (1999). It is shown that under some conditions, the entrainment of eight small bubbles, less than 1:1 mm in diameter, results in significant levels of vortex distortion, comparable to the experimental observations. We find that the bubble induced distortion is due almost entirely to volumetric displacement effects present in the DEM model. A relative reaction force, defined as the ratio of net bubble to fluid reaction to the local driving force of the vortex, is used to analyze the vortex decay rate. It is shown that the global increases in vortex decay rate are directly proportional to the magnitude of this highly local relative reaction.
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: VID00297
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 1211-Finn-ICMF2010.pdf

Full Text

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

Simulation of Bubble and Particle Interactions with Vortical Flows Using a Discrete
Element Model

Justin Finn* and Sourabh V. Apte*

School of Mechanical Industrial and Manufacturing Engineering, Oregon State University, Corvallis, OR 97333, USA and
Keywords: Euler-Lagrange, Bubbly Flows, Discrete Element Model, Two Way Coupling

Discrete Element Model (DEM) is used to simulate bubble and particle interactions with vortical flows. In this
Euler-Lagrange approach, volume-averaged Navier-Stokes equations are solved using a co-located grid finite volume
method. The dispersed phase is assumed spherical and is modeled by tracking the centroids of the particles/bubbles
in a Lagrangian frame. In addition to interphase momentum exchange, the volume occupied by the disperse phase is
accounted for through fluid void fractions in the momentum and continuity equations. Interactions of the dispersed
phase with vortical flows are investigated for both dilute and dense loadings. First, flow generated by the disperse
phase motion are simulated to reproduce experimentally observed vortical flow features: (i) the breakup of a viscous
blob of particles falling under the influence of gravity and (ii) the transient migration of a buoyant bubble plume.
Finally, bubble entrainment and interaction with traveling vortex tube under dilute loadings are simulated in a two-
dimensional approximation of the experiments by Sridhar and Katz (1999). It is shown that under some conditions,
the entrainment of eight small bubbles, less than 1.1 mm in diameter, results in ,igniliik.m levels of vortex distortion,
comparable to the experimental observations. We find that the bubble induced distortion is due almost entirely to
volumetric displacement effects present in the DEM model. A relative reaction force, defined as the ratio of net bubble
to fluid reaction to the local driving force of the vortex, is used to analyze the vortex decay rate. It is shown that
the global increases in vortex decay rate are directly proportional to the magnitude of this highly local relative reaction.

Bubbles, droplets, or particle laden turbulent flows are
of interest in a broad array of important disciplines with
applications to cavitation, sprays and spray combustion,
coal combustion and fluidized beds. To date, our under-
standing of the phase coupling mechanisms which gov-
ern multiphase flow development is limited to low com-
plexity problems in basic canonical geometries. Com-
putational limits continue to restrict the use of complete,
fully resolved direct numerical computations to engi-
neering scale configurations for quite some time. It is
therefore important to develop models capable of han-
dling diverse flows containing large numbers of bubbles
or particles that can accurately capture the interphase dy-
namics defining the characteristics of multiphase flows.
A recent review by van der Hoef et al. (2008) dis-
cusses present modeling strategies used for numerical
simulation of gas/solid fluidized beds. The subset of
methods currently capable of handling both large scale
geometries and a large number of bubbles or particles
can typically be classified as either Eulerian-Eulerian
or Eulerian-Lagrangian. In an Eulerian-Eulerian or

two-fluid model (for example, Ferrante and Elghobashi
(2007)), both phases are treated as a continuum with
unique fluid properties and two sets of conservation
equations are solved. Because the idea of individ-
ual particles is not supported, closure models for inter-
phase momentum transfer are required. In an Eulerian-
Lagrangian simulation such as the ones presented in this
paper, each element of the dispersed phase is treated as a
sub-grid scale, spherical particle. Bubble or particle po-
sition is advanced in time in accordance with Newton's
second law, with problem specific closures required for
the various forces acting on the dispersed phase. Be-
cause the bubbles are assumed to be sub-grid in size their
surface shape cannot be explicitly represented, and no
interphase boundary condition is enforced.

Typical treatment of the dispersed phase as 'point par-
ticles' allows the bubbles to transfer momentum as point
sources but does not account for the volume of fluid
which they displace. To further advance these mod-
els, the effects of fluid volume displaced by the dis-
persed phase, termed in this work as volumetric cou-
pling, is considered using a variable density form of

the governing equations, as in Prosperetti and Zhang
(1995), Joseph et al. (1990), van der Hoef et al. (2008).
Specifically, the importance of volumetric coupling for
dilute as well as dense loading of the disperse phase is
Test cases dealing with dispersed phase interactions
with single, distinct vortical structures can be good
building blocks for understanding more complex prob-
lems, and have been studied by several groups. Srid-
har and Katz (1995) and Van Nierop et al. (2007) used
experimental techniques to develop models for bubble
motion in non-uniform flow, with particular emphasis
placed on the determination of the bubble lift coeffi-
cient. Oweis et al. (2005) used fully resolved, front
tracking methods to solve the flow field around deform-
ing and cavitating bubbles during entrainment by a sta-
tionary Gaussian vortex. They compared bubble capture
time predicted by a passive point particle method with
the front tracking results as well as experimental data
under similar conditions. They found their point particle
method was able to accurately capture the trajectory of
bubble entrainment, up to the point of bubble cavitation
and volume growth in the vortex core. Sridhar and Katz
(1999) observed the effects that bubbles, with diameter
400pm < db < 1, 100pm, had on the structure of pis-
ton generated vortex rings. Their experimental results
show that for a small number of entrained bubbles, at
low overall volume fraction, igmnilik. ii distortion of the
ring structure was possible under certain conditions. In
significantly distorted vortices, the presence of the bub-
bles resulted in a fragmented core, with multiple regions
of higher vorticity. They supplied good analytic ratio-
nalization of their results, including the bubble settling
locations, and the observed vortex distortion at low bub-
ble volume fraction.
In more general turbulent flows, bubbles and particles
are known to accumulate in certain regions. Bubbles are
attracted to the vortex cores, while particles will collect
in the the high shear, low vorticity regions. This prefer-
ential accumulation has been studied by Ferrante and El-
ghobashi (2003); Climent et al. (2007); Mazzitelli et al.
(2002) and can be responsible for, igiilk.iil flow mod-
ulation. More complex engineering scale problems have
also been simulated with the Eulerian-Lagrangian ap-
proach. There are growing bodies of literature dedicated
to bubble column reactors (Darmana et al. (2005); Hu
and Celik (2008)) and combustion chambers (Apte et al.
(2003); Moin and Apte (2006)) which show promising
results for the technique.
In this work it is shown that volumetric displacement
effects can be of importance even under dilute loadings,
especially for bubbly flows by simulating interactions
of small number of small bubbles with traveling vortex
tube. First, the numerical model and its implementation

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

is verified by simulating interactions of disperse phase
with vortical structures under dense loadings. The prob-
lems vary in their level of complexity, but each demon-
strate the ability of the present methods to capture bub-
ble or particle induced flow phenomenon with subgrid
modeling of the disperse phase. In the next section, the
governing equations and mathematical formulation are
outlined. Various verification test cases, chosen to show
the accuracy of the methods and the importance of vol-
umetric coupling, are discussed next. Finally, bubble
interactions with a traveling vortex tube are studied in
detail for dilute loadings. A correlation is developed for
the bubble induced vortex decay rate, based on the rela-
tive magnitude of the net bubble to fluid reaction.

Mathematical Formulation
In Discrete Element Method (DEM), the motion of the
fluid phase is computed using Direct Numerical Simu-
lation (DNS) or large-eddy simulation (LES), and La-
grangian particle tracking is used to solve the motion
of spherical subgrid scale dispersed phase. Dispersed
phase motion is calculated from Newton's second law,
with models for gravity, pressure, lift, drag, and added
mass forces. A bubble to fluid reaction is calculated
to account for two-way momentum exchange. In addi-
tion, the present formulation is derived from volume av-
eraging, and accounts for volumetric displacement of the
fluid by the dispersed phase. In this section, we use sub-
scripts b to denote the dispersed phase and f for the con-
tinuum phase. We refer to individual dispersed phase el-
ement as a 'bubble'. However, the formulation is equally
applicable to gas/solid particle systems.

Dispersed Phase Motion In the Lagrangian refer-
ence frame, the equations of motion may be written for
each bubble as a system of ordinary differential equa-

j (Xb)
mb d Ub)


where Fb is the net force acting on each bubble and has
the following contributions:

Fb FG + FP + FAM + FD + FL (3)

The gravitational force, FG, is the weight of the bubble.

FG = PbVbg

where Vb is the bubble volume and g is the gravitational
acceleration. The pressure force, Fp, is the force on the

bubble due to the total pressure gradient, including the
hydrostatic contribution.
Fp = -VbVP (5)
The added mass force, FAM, is the force which would
be exerted on the volume of fluid displaced by the pres-
ence of the bubble. It is given by


Dt )'

where uf,b represents the fluid velocity interpolated to
the bubble centroid location. For small, spherical bub-
bles and particles, it is generally accepted that CAM =
0.5. Inter-bubble collision forces are neglected in this
work. The bubble drag force, FD arises due to a differ-
ence in bubble and fluid velocities. It is given by:
1 7dc
FD -2CDPf IUb Uf,blUb f) (7)

where Ab = is the frontal area of the bubble. The
bubble slip velocity, (ub Uf,b) is evaluated using the
local velocity field near the bubble of interest. The lift
force, FL, arises for bubbles in shear or rotating flow and
has been the subject of much discussion in the literature.
In general, it can be expressed as:
FL = -CLPfVb(Ub Uf,b) x (V x Uf,b) (8)
Where CL is the lift coefficient. At Reynolds numbers
much greater than one, CL and CD are know by empir-
ical relations only, and we choose expressions for these
values on a case by case basis. All cases except the bub-
ble column use the standard drag curve of Schiller and
Naumann (1935) for solid sphere drag
CD (1 + 0.15Rc 687) (9)
The bubble column case uses expressions suggested
by Darmana et al. (2006)
CL 0.5, C, = max [Cee, C]Eo
[ 16 0.48
C R min (1 + 0.15Re 687 (10)
I Reb b ReC


Eo + 4'

(Pf Pb '-''.

The Gaussian vortex case, and the traveling vortex tube
case use the strong lift coefficient suggested by Sridhar
and Katz (1995).

Cvz 0.22c-3/4;

IV x Uf,bldb
21uf,b Ub

The falling blob case uses the lift coefficient suggested
by Saffman

C ob 161 p (pfw) (12)

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


Figure 1: Illustration of the fluid domain, F containing
a subgrid scale dispersed phase.

Variable Density Fluid Formulation Consider a
domain F which contains a continuous fluid phase, and
a dispersed bubbly or particulate phase as shown in fig-
ure 1. Each element of the dispersed phase has a fi-
nite, characteristic diameter, db, and occupies a volume
Vb. The finite size of the disperse phase may be ac-
counted for by introducing the fluid and disperse phase
volume fractions, Of and eb, where Of 1 eb.
In doing so, the volume averaged continuity equation
becomes (Joseph et al. (1990); Zhang and Prosperetti
(1997); Jackson (1997))

(e)+V(eu)o0 (13)
Note that in this form, the fluid velocity field is not di-
vergence free, even for an incompressible fluid (Apte
et al. (2007)). Similarly, the fluid momentum equation
becomes (using the model formulation suggested by Gi-
daspow (1994),

(P Ofu) + V (POfuu)


V (effp(Vu + VUT)) + OfPfg + f (14)
The force f is not present in a single phase flow, and
arises due to interphase momentum transfer from the

f =- (Fp + FD+FL+FAM)

The interphase reaction force (f) can be split into (i) the
contribution from the pressure force and (ii) a net drag
from all other forces as described by van der Hoef et al.
(2008). The pressure force then results in the force den-
sity +bVP. This reaction term related to the pressure
gradient can be combined with the pressure gradient in
the momentum equation to obtain:

a (Pfu) + v. (PfeOuu)

(11) V (Ofe/(Vu + VT))+ f +

-VP Ofpfg

ObVP (16)
Fp Force Density

where ebVp is the Eulerian force density obtained from
the pressure force and f' is the Eulerian force den-
sity constructed from the Lagrangian force on the bub-
bles without the pressure force (equation 15 without the

pressure force, Fp). Noting that b + Of = 1, the
above equation can be re-written in a more commonly
used form by combining the first and last terms on the
right-hand side of the above equation (Gidaspow (1994);
van der Hoef et al. (2008)),

a (peOfu) + V (Pf Ouu)

-e V (P)+

V. (ff (Vu+VuT)) fpfg+f', (17)

where f' contains summation of all reaction forces in
equation 15 except the pressure force. This formula-
tion is commonly used in gas-fluidized beds (Gidaspow
(1994); van der Hoef et al. (2008)). In case of Reynolds
averaged Navier-Stokes equations or large-eddy simula-
tions, the above equations should be time-averaged or
spatially filtered using density-weighted Favre averag-
ing. Using the form in equation 17; however, gives rise
to an additional unclosed term fVP. It is therefore
advantageous to use the original form (equation 14), re-
sulting in standard variable density LES equations Moin
and Apte (2006). In this case, the reaction due to the
pressure force is treated I h.1;. hi'.
We choose to use the form of equation 14 for numeri-
cal implementation. The formulation for the fluid phase
given by equations 13, 14, and 15 represents what we
will refer to as volumetric coupling, where both bubble
size and momentum transfer are accounted. The numer-
ical solution of the fluid and bubble motion proceeds in
the framework of a parallel, fractional step, finite volume
solver (Shams (2010)).

Bubble entrainment in a Gaussian vortex
To verify the accuracy of the solver, the entrainment of
a single bubble into a stationary Gaussian vortex is con-
sidered. A schematic of the problem is shown in fig-
ure 2a. A buoyant bubble is released in the vicinity of
the vortex with core radius, rc. The Gaussian vortex is
a planar vortex with initial circulation F0 whose vortic-
ity distribution is a Gaussian function of radius. There
is no radial velocity component, and the the tangential
velocity can be expressed as

Ue(r) =2 (1 -e- T (18)
The vorticity and maximum tangential velocity are
The vorticity and maximum tangential velocity are

w(r) ro 1e 1(r )2;

2 27rr,

where rl1 and Tr2 are constants. Assuming axial sym-
metry, the hydrodynamic pressure gradient is ip/lr
/r. 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. Under

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

the right conditions, ie. the vortex is strong enough, the
bubble will become entrained in the vortex core. During
this process, it may circle the core several times before
eventually reaching a settling location with relative co-
ordinates rs, 0, measured from the vortex center.
A total of 14 individual cases of single bubble entrain-
ment are simulated using the passive, one-way coupling
approach. Variables relevant for the setup of the para-
metric study are summarized in table 1. The domain
size is 7re x 7rc x 0.4 re. The no-slip condition is im-
posed at boundaries in the X and Y directions, and the
domain is periodic in the Z direction. Uniform grid spac-
ing is used throughout the domain with A 1.25mm
and total number of grid points as 64 x 64 x 4 in the
X, Y and Z directions. The flow solver time step is
At = 0.4ms whereas the bubbles advanced with a time-
step of 2 pUs using subcycling per time-step. The veloc-
ity field of equation 18 is applied as an initial condition
everywhere in the domain, creating a clockwise vortex.
A single bubble is released during the first timestep at
r r, 0 0. In all cases, gravitational accelera-
tion is g = -9.81m/s2 in the Y direction, and typical
properties of water and air are assumed. For the pa-
rameters shown in table 1, the bubble Stokes numbers
(Stb .= 7. /36/, where w = F/Trr ) are between 0.5
and 5.7.

Force balance in vortical flow At the settling lo-
cation, there is no motion of the bubble relative to the
vortex, meaning that all forces acting on the bubble are
in balance. This is illustrated in figure 2b, where we have
included the lift, drag, pressure, added mass, and grav-
ity forces. Note that the pressure force has been split
into its dynamic and hydrostatic contributions. For the
clockwise vortex shown, the settling location will be in
first quadrant of the core where the fluid velocity is turn-
ing downward. Mazzitelli et al. (2002) showed that it is
primarily the lift force, which acts perpendicular to the
vorticity and slip velocity vectors, that is responsible for
bubble accumulation in the downward velocity side of
vortices such as this one. If the flow is steady and ax-
isymmetric then the directionality of the other forces in
figure 2b can also be deduced. The net buoyancy force
which acts upward is due the addition of the hydrostatic
pressure force, FpH and the bubble weight, FG. The
added-mass and dynamic pressure forces act in the di-
rection of negative dynamic pressure gradient. Neglect-
ing outside disturbances, this will be toward the vortex
center. The drag force acts in the direction of the slip
velocity vector which, for a stationary bubble in an ax-
isymmetric vortex, is perpendicular to the settling loca-
tion vector, rs.
Since the velocity, vorticity, and pressure gradient are
known functions of radius in the Gaussian vortex, it is

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

Table 1: Computational parameters used for the Gaussian vortex.
rc[mm] /i T/2 A[mm] Nyrid Fo[mT /s] db [m]
11.45 1.27 0.715 1.25 64x64x4 0.02, 0.03, 500, 700, 900,
0.04, 0.05 1,100, 1,300

+ FpD

FG + FpH

Figure 2: Schematic of a bubble entrained in a vortex and the forces which influence its settling location, (rs, 0s).

possible to obtain analytic expressions for the settling
coordinates (rs, 0s) of entrained bubbles using force bal-
ances in the radial and azimuthal directions. Van Nierop
et al. (2007) performed such an analysis to develop un-
coupled expressions for bubble settling coordinates in
forced, rotating flow. Here, we employ a similar ap-
proach to obtain coupled expressions for settling loca-
tion in the free, Gaussian vortex. Using figure 2b as a
guide, a force balance in the azimuthal direction reveals
that the drag force is balanced entirely by a component
of the net buoyancy force.

(FG + FpH) COs(O) FD (20)

In the radial direction, there is a more ambiguous bal-
ance between the lift force, the dynamic pressure force,
the added mass force and a component of the gravity

FB sin 0s + FL FpD + FAM (21)

The dynamic pressure force and added mass force can
be combined,

FpD + FAM (1 + CAM) Pf (22)

Inserting this and the expressions for FD, FL and FG
into equations 20 and 21, the force balances can be rear-

ranged into two coupled equations for rs and O,:

cos (Os)

4dbg [b -1]
(1+ CAM)U2
[P 1] gsin(0,) +CLUoR7

If the flow field inside the vortex core is known (Gaus-
sian, Taylor-Green, Rankine, etc.), and the drag and lift
coefficients are specified functions of Uo and cw, then
equations 23 and 24 can be solved iteratively for the
bubble settling coordinates. Figure 3 shows a compari-
son of the settling location predicted by the DEM model
with passive, one way coupling, to the values predicted
from solving equations 23 and 24 for the 14 individual
cases, showing good agreement. In the analytic predic-
tion, the value of Fo has been decreased by 12% in all
calculations to compensate for the viscous decay which
happens before the bubble comes to rest. Time averages
of the numerical, 1-way settling location are taken over
a period of 0.1 seconds, 1.8 seconds after injection.

LES of a Rising Bubble Column
For this test case, the geometry selected corresponds to
the 'Becker case' (Becker et al. (1994)), which has be-
come an almost standard case in the chemical engineer-
ing literature for dense bubbly flows. The domain is
shown in figure 4(a). Characteristic dimensions and rel-
evant simulation parameters are summarized in table 2.

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


Figure 3: Settling location of bubbles in the core of the
stationary Gaussian vortex: (.) analytical, (0) DEM.

All boundaries are assigned a no-slip condition except
for the top wall at Y = 1.5, where a slip condition
(.. = 0) is applied to approximate the experimental free
surface. Air bubbles are continuously injected into the
water filled domain from a disk located on the bottom
wall with a flow rate of 1.6 1/min and a superficial gas
velocity of 0.66m/s. The injection disk is 0.04m in di-
ameter and located 0.15m from the left hand wall. At
this gas flow rate, the fluid and bubble motion is known
experimentally to be periodic in nature. We use this test
case to see if the present methods can capture the tran-
sient behavior. A Dynamic Smagorinsky model for vari-
able density fluid flow is used to compute the subgrid
scale fluid stress.

Table 2: Parameters for the bubble column case.
L,, L,, L. 0.5m, 1.5m, 0.08m
n, nz 80, 150, 15
Gas Flow Rate 1.6 1/min
Superficial gas velocity 0.66mm/s
Injection location (x y, z) (0.15, 0, 0)
Injection Diameter di = 0.04m
Bubble Diameter 1.6mm
Bubble Density 1.2kg/m3
Liquid Density 1, O(0l0L/m3
Liquid Viscosity 10 3Ns/m2
Probe 'A' (x, y, z) 0.035m, 0.9m, 0.04m
Probe 'B' (x y, z) 0.45m, 1.05m, 0.04m

Figure 4(b) shows the streamlines corresponding to
the average velocity field in the X-Y midplane. In the
lower half of the domain, a large clockwise vortex is
generated by the bubbles rising along the left hand wall.
The upper half is highly transient, and marked by the pe-

(a) Geometry

0 0.1 0.2 0.3 0.4 0.5
(b) Avg. Velocity field

Figure 4: (a)Bubble column geometry showing the in-
stantaneous bubble locations. Points A and B corre-
spond to the location of two velocity probes. (b) Con-
tours of time averaged Y component of fluid velocity.
Streamtraces show the time averaged velocity field.

riodic migration of the bubble plume in the X direction.
This migration corresponds with the growth and collapse
of several secondary vortices in the upper regions of the
domain. The periodicity of the bubble plume is shown
in figure 5 where we have plotted 8 snapshots of instan-
taneous bubble positions, each 7 seconds apart. In fig-
ure 5a, the plume is firmly directed against the left wall
through most of its length due to the strong lower vor-
tex. A secondary, counter-clockwise vortex located in
the upper left comer of the domain pushes the top of the
plume to the right. In figures 5b through 5e, this sec-
ondary vortex strengthens and slides downward along
the left hand wall, creating a bulge in the plume. This
bulge cannot travel into the lower part of the plume due
to the size and strength of the main circulation region,
and eventually it collapses, as a new counter-clockwise
vortex is created in the upper left comer (see figures 5f
through 5h).
The periodic nature of the liquid phase is shown in
figure 6, where we plot the vertical velocity in time for
points 'A' and 'B' as shown in figure 4(a). The magni-
tude of the velocity peaks at these two points is similar
to the values measured in the experiments of Sokolichin
and Eigenberger (1999). We observe an average oscil-
lation period of about 49 seconds, which is 8 seconds
longer than observed in the experiments. Deen et al.
(2000) and Delnoij et al. (1997) also had trouble match-
ing the oscillation period in their computations. There

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

(a) (b) (c) (d) (e) (f) (g) (h)

Figure 5: Periodic migration of the bubble swarm. Each figure is a snapshot of instantaneous bubble positions, with 7
seconds between each frame. The first (a) and last (h) frame are chosen to correspond with the approximate beginning
and end of one cycle.

are numerous factors which play a role in the oscillation
period including grid resolution, bubble size, lift coeffi-
cient and drag coefficient. Also it was found that small
changes in the injection velocity and inlet conditions can
influence the oscillation periods. The LES results with
volumetric coupling, however, do show the experimen-
tally observed trends of periodic jet oscillation and large
scale vortical regions.

( i v l i

050 -'-'- 50 . .500 '-'T--0 2,D c-

Time [s.
(b) Point 'B' vertical velocity

Figure 6: Fluid velocity at the points 'A' and 'B' shown
in Figure 4a.

Three-Dimensional Falling 'Blob'
This test case is included to show the applicability of
methods to a flow driven entirely by particle motion.
Consider a three dimensional sphere composed of sev-
eral thousand tightly packed rigid particles. The fluid
region containing this particle 'blob' has a uniform ini-

tial particle volume fraction, epo. Both the local den-
sity, pf, and viscosity, pf, vary according to this volume
fraction. In the high particle number, Np, limit, the be-
havior of this particle/fluid system is analogous to a vis-
cous, droplet dissolving in a fluid with which it is misci-
ble. For two such fluids, the surface tension force at the
droplet interface vanishes, and the droplet dynamics are
governed by viscous and inertial forces. The problem
may be described in terms of the density ratio, viscosity
ratio, blob Reynolds number, and blob Froude number:

DR = Pb/pf; R= b/lf

Reb =p UDb; Frb U
[f V-Db

Mitts et al. (1996) used the miscible fluid analogy
to study the dynamics of deformation and breakup of
droplets in the supercritical regime. Supercritical con-
ditions are present in many rocket engines, where high
temperatures and pressures are required, and many ques-
tions remain about the behavior of the fuel-air mixture.
They were able to classify the miscible droplet behavior
into several sub-regimes based on the viscosity ratio and
droplet Reynolds number for density ratios near unity.
In a computational study, Walther and Koumoutsakos
(2001)(W&K) used adaptive vortex methods and La-
grangian particle tracking to model several of these sub-
regimes, including the multi-mode bag vortex ring sub-
regime, which will be the focus of this validation case.
In their work, they characterized the speed of a falling
blob of particles using the Hadamard and Rybczyfiski
(H-R) formula for the speed of a viscous, spherical drop
descending under the force of gravity :

ST(Pd p)gD /t + td
U UHR -3
12p / + -ld

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

In this test case, 97, 233 solid particles are arranged into
a spherical shape in an initially quiescent fluid. The par-
ticles are heavier than the liquid, and fall under the influ-
ence of gravity. The local viscosity of the fluid is varied
at the grid faces:

P = Pf + Op(,p yf) (26)

where Op is the particle volume fraction of the associ-
ated control volumes. The initial volume fraction is ob-
tained from the following expression derived by Lund-
gren (1972) for high volume fractions:

ptj 51op

The computational parameters were selected so that
comparisons could be made with both experimental and
numerical studies of the bag vortex ring sub-regime. The
characteristic velocity, U, is defined using the H-R for-
mula (equation 25). In doing so, the problem can be
specified completely by assigning the density ratio, vis-
cosity ratio, Reynolds number, and Froude number. For
clarity, all parameters used in the setup of this case are
summarized in table 3. Note that the particle density is
greater than the blob density due to the initial volume
fraction which is less than unity.

Table 3: Parameters used for the setup of the three di-
mensional falling 'blob.'
Parameter Value
Db 100 mm
dp 1.572 mm
Np 97,233
g 9.81 mn/s2
Opo 0.388
[b/pf 35.0
Pf 0.001 pa-s
Pb/Pf 1.6
f 1.0 kg/rn3
pp 2.544 kg/n3
Reb 330
Frb 3.332
Domain Size 10Db x 10Db x 10Db
Grid Size 64 x 64 x 64

Three different DEM models are considered and are
summarized in table 4. Model 1 is intended to be similar
to the one used by W&K. Here, the particle motion is af-
fected only by gravity, and drag forces, and the particle
to fluid momentum coupling is a result of the drag force
only. The fluid density is constant everywhere. In model
2, the added mass force, lift force, and pressure force
are also considered in the particle equation of motion,

and the momentum coupling to the fluid includes these
additional terms. In model 3, the variable density, volu-
metric coupling model is used, and the fluid density, pf,
varies according to the local particle volume fraction.

Table 4: Comparison of the 3 present DEM models
and the model used by Walther & Koumoutsakos in the
falling blob case
Model FG FD FL FAM Fp Variable
W&K Yes Yes No No No No
#1 Yes Yes No No No No
#2 Yes Yes Yes Yes Yes No
#3 Yes Yes Yes Yes Yes Yes


Figure 7: Blob speed non-dimensionalized by the H-
R velocity as a function of non dimensional time (-
)Model 1, (... )Model 2, (-- -)Model 3.

The speed of the center of gravity of the descending
blob is shown in figure 7. In the present case, the peak
blob Reynolds number based on blob velocity and ini-
tial diameter is 40 for models 1 & 2 and 46 for model
3. For the present viscosity ratio of 35, this puts the re-
sults within the bag vortex sub-regime observed by Mitts
et al. and simulated by W&K. The blob speed over the
time simulated is nearly the same for models 1 and 2,
showing that the net effect of lift, added mass and pres-
sure forces is negligible for this case. Further, this shows
that the increase in blob velocity shown by model 3 is
due entirely to the volumetric effects. The blob accel-
erates to a peak non-dimensional velocity of 0.12UHR
for models 1 & 2 and 0.14UHR for model 3 at a non-
dimensional time of about 27. The difference in speeds
of the falling blob suggests that the model for drag force
used in the volumetric coupling needs modifications to

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

(a) (a)



(b) ' -

(d) 1 ,'l

-\--"---- -"-, .
,f, ?^ ^..' 'U

(a) Re = 330, Fr = 3.33Model 1

(b) Re = 330, Fr = 3.33Model 2

(c) Re = 330, Fr = 3.33Model 3

Figure 8: Cross section of falling blobs with pb/p 1.6 and pb/p = 35 at different non dimensional times (tU/D) =
(a) 0.00, (b)7.18, (c) 10.77, (d) 14.36, (e) 17.94, (f) 21.53. [Left] Vortex method simulation of Walther and Koumout-
sakos (2001) at Re = 180 & Fr = 1.79. [Center & Right] Present results of model 1,2,& 3 at Re = 330, Fr = 3.33.

account for local changes in the fluid volume fractions.
Figure 8 compares cross sections of the blob at dif-
ferent non-dimensional times throughout this period of
initial acceleration for each model. During this time, the
blob deforms into a spherical cap shape as the particles
in the center of the blob descend faster then the ones on
the edge. This is due to viscous action at the edge of the
blob which causes the roll-up of a vortex ring around the
blob. The ring travels with the blob and causes continued
spreading in plane normal to the falling direction. Qual-
itatively similar spreading was observed by Walther and
Koumoutsakos (2001) with their vortex method predic-
The differences in blob shape during the acceleration
to peak velocity can be explained by the variable density
effects of model 3. With volumetric coupling consid-
ered, the fluid density in the blob region is lowered in
proportion to the particle volume fraction. This results
in a lower resistance to particle motion and helps to ac-
celerate the fluid with the blob. The blob deforms to a
spherical cap shape rapidly in models 1 & 2 because the
particles on the bottom of the blob do the work of push-
ing the initially quiescent fluid out of the way, while the
particles on the top accelerate. In model 3, this effect is
delayed by the lowered fluid resistance in areas of high
volume fraction, meaning the blob accelerates in a more
uniform fashion. At later times, the blob spreads and
forms a bag vortex ring which travels with the blob.Once

the ring is spread sufficiently wide (tU/D > 60), the net
volumetric effects are decreased, and the blob velocity
for model 3 approaches that of models 1 & 2.
The three dimensional shape of the blob is shown
at different stages of breakup examined by Mitts et al.
(1996) in figure 9. After the peak velocity is reached,
the blob continues to expand into an unstable ring like
shape. The particles collect in the high-shear regions
on the downward velocity side of the vortex ring, and
drive the circulation of the ring (figure 9(b). At still later
times, the ring becomes unstable, and the particles clus-
ter into four regions of high volume fraction(figure 9(c)).
These clusters accelerate as the original blob did, and
form four new vortex rings. The initial deformation, the
bag vortex ring roll up and instability, and subsequent
droplet splitting are all predicted and show good levels
of similarity with the experimental observations.

Bubble in a Traveling Vortex Tube
Interactions of micro-bubbles with a traveling vortex
corresponding to the experiments by Sridhar and Katz
(1995, 1999) are investigated. Specifically, they ob-
served that with dilute loading of eight small bubbles,
once entrained into a vortex ring, could deform the ring
significantly. The simulations are setup to investigate if
the DEM approach with volumetric effects can capture
this phenomenon.
The computational domain and the evolution of an


(a) Initial Deformation

(b) Bag Vortex Ring Roll-Up

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

undisturbed pair of vortex tubes are shown in Figure 10.
There is an inflow boundary at the left wall, an out-
flow condition at the right boundary, and walls on the
top and bottom. The total domain size is X/hjet 10,
Y/hjet 3 and is centered at Y/hjet 0 At the left in-
let boundary, ajet is pulsed for 0.27 seconds into the ini-
tially quiescent domain, which causes the roll up of two
symmetric vortex tubes. The contours in figure 10 show
the diffusion of high vorticity as the vortex tube travels
downstream. Table 5 lists the computational parame-
ters used in this study. The inflow velocity is a function
of time, and is described by a sixth-order polynomial
to match the experimental conditions. To keep the size
of the computation small and allow several parametric
studies, a two-dimensional domain is simulated with pe-
riodic conditions in the spanwise direction giving rise to
a vortex-tube. A uniform Cartesian grid is used through-
out the area below the line of symmetry with a total of
800 x 121 elements in the X and Y directions.

Table 5: Parameters for the traveling vortex tube case.

(c) Droplet Splitting

Figure 9: Comparison of the present results of model
3 with the experimental results of Mitts et al. (1996)
for the bag vortex ring sub-regime of the multi-mode
breakup category.

Pf ,vf
Domain Size
Grid Size
Jet height (hjt)
Inflow Time
Vortex Strength,
Bubble Size (pm)
Inflow Velocity
a6, a5, a4,a3
a2, al, ao

1,000 kg/m3; 10-6 m/s2
1m x 0.15m x 0.005m;
800 x 121 x 4
0.27 s
0.0159, 0.0207, 0.0254

500, 700, 900, 1,100
U(t)_ E o antn
62278; -47082; 13686; -2062
159.5; -1.289; 0.006

6 8 10

Figure 10: Computational domain and vortex tube evo-
lution: (a) Two, symmetric vortex tubes are created by
an inlet jet pulsed at X 0 Contours show vorticity
out of the plane as the undisturbed vortex travels down-
stream. Location of bubble injection is X/hjet 5.0

At a value of X/hjet 5.0, eight bubbles are in-
jected below and in front of the vortex core. Due to
buoyancy, the bubbles rise around the rear of the vor-
tex and are swept into the downward velocity region on
the forward side of the passing core. A parametric study
is performed to determine how bubble settling location
and vortex structure are affected by bubble size and vor-
tex strength. Depending on the Stokes number, the bub-
bles may circle the core multiple times before ultimately
reaching their final settling location, where their average
motion relative to the vortex center is zero. Once they
have reached this state, the settling coordinates (rs, 0s)
are averaged over all bubbles and in space over a dis-
tance of 5.2X/hjet < Xvx < 5.8X/hjet. The aver-
age settling radius for each case is plotted against the
non-dimensional parameter gd /8Fr (ratio of buoyancy
force to hydrodynamic pressure gradient) in figure 11
alongside the experimental data. Even at this small over-
all volume fraction, the local volume displacement ef-
fects are crucial in obtaining the correct bubble settling

0os5 9 0 0 0
Bubble In ion 0
-1 Bubble Injection-

O 2 4 X/h

radius, particularly as the ratio of bubble size to vortex
strength is increased.

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

then just the time derivative of the angular momentum.

Decay rate = c




3 A

: *

_* _


0.1 0.2 0.3
rs [cm]

0.4 0.5

Figure 11: Comparison of the non-dimensional set-
tling location with the experimental results of Sridhar
and Katz (1999). The parameter gdb/8F2 is the non-
dimensional ratio of the buoyancy force and the hydro-
dynamic pressure gradient experienced by the bubble.
(*) Experimental data, o two-way coupling (neglecting
volumetric effects), A volumetric coupling.

Inipl.ii.lill\ sigliliik.iii vortex distortion is predicted
by the numerical model for certain bubbles with param-
eters similar to those observed in the experiments. In
the absence of the bubbles, the vortex core is stable as
shown in figure 12a. In the presence of some entrained
bubbles, the vortex core deforms, altering the vorticity
distribution as shown in figure 12b. The contours of vor-
ticity in figure 12b are shown approximately 1 second
after bubble entrainment. The core has been fragmented
into several regions of higher vorticity. Figure 12c shows
how the average radial vorticity distribution is changed
because of the presence of the bubbles. In the undis-
trubed state, the vorticity distribution is very close to that
of the Gaussian vortex from the first test case. The vol-
umetric displacement of the fluid due to bubble motion
results in a decrease of inner core vorticity, while a band
of high vorticity has been created just outside of the bub-
ble settling radius. As a quantitative measure, we calcu-
late the change in the decay rate of angular momentum,
illustrated in figure 12d for each case with and without
bubbles. The instantaneous vortex angular momentum
is calculated by summing the momentum of all control
volumes in the core

L,= p-fUorcV,, (28)
core cv's

Where r,, is the distance from the vortex centroid to the
cv center and V, is the cv volume. The decay rate is


In all cases, we observe an increase of the decay rate
from c to i when considering volumetric coupling, al-
though the relative amount of this increase varies signif-
icantly. To normalize the amount across all cases, we
introduce the relative change in decay rate,

Relative change in decay rate (%) E = x 100
The overbar denotes the average decay rate measured
between X 5.2hj, to X = range 1% for 500/pm bubbles entrained in a vortex with
Fo = 0.0254m2s 1, to almost 150% for, 1, 100pm bub-
bles entrained in a vortex with Fo 0.0159m2s 1
These characteristics of vortex distortion were found
to be predominantly an effect of variations in the void
fraction as the bubbles travel to their settling location.
This was confirmed by computing the bubble trajectories
without considering the void fraction variations. With
two-way coupling and neglecting volumetric displace-
ment effects (i.e. Of set equal to 1), vortex distortion
was not obtained for any of the cases studied. This con-
firms the effectiveness of the present numerical model in
properly predicted bubble-vortex interactions.

A correlation for decay of angular momentum
The increase in vortex decay rate varies significantly de-
pending on bubble size and vortex strength. Increas-
ing bubble size (or decreasing the fluid volume fraction)
seems to have an effect, as does decreasing the vortex
strength. Here, we develop an analysis to understand
how this increase in the vortex decay rate scales with
the magnitude of the highly local, bubble-fluid interac-
tions. We start by calculating the net reaction force act-
ing on the fluid due to the presence of the bubbles with
volumetric coupling. A similar analysis was conducted
by Sridhar and Katz (1999); Druzhinin and Elghobashi
(1998). We will derive an expression here applicable to
the current Eulerian-Lagrangian framework, and use it
to calculate the net reaction from the entrained bubbles
at their settling location in a Gaussian vortex. We can
write the single phase, undisturbed momentum equation
C P + V + B, (31)
where C denotes the convective fluid acceleration term,
P denotes the pressure term, V the viscous term, and
B the gravitational body force term. Similarly for mul-
tiphase flow, the volume averaged momentum equation
(equation 14) can be expressed as

C P + V + B + f (32)

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

05 1 15

Time [s]
Time [s]

(a) Undisturbed vortex core

(b) Fragmented core due to finite (c) Redistribution of vorticity (d) Increase in angular momentum
size effects decay rate

Figure 12: Characteristic changes observed in distorted vortex cores due to the presence of bubbles with finite size: (a)
Undisturbed vorticity contours, (b) vorticity contours and average bubble location (black dot) with finite size effects,
(c) radial re-distribution of vorticity (- -) without void fraction variations, (-) with void fraction variations, (d)
increase in decay rate of angular momentum from c without bubbles to c with bubbles and accounting for void fraction

Here, C, P, V, and B are the same contributions to
the fluid momentum, but are representative of the dis-
turbed flow field containing the eight finite size bubbles.
If all terms are moved to the right hand side, and the sin-
gle phase terms are subtracted from the disturbed flow
terms, we can write an expression for the reaction force
per volume imposed on the fluid by the presence of the
dispersed phase:

AR- (C )+( -P)+( V)+( -B)+f
In order to simplify this and arrive at an expression ap-
plicable to the present case, we make the following as-
sumptions. First, the undisturbed flow is steady over the
timescale of bubble to fluid momentum transfer. Second,
the bubbles are not accelerating. Once they reach their
settling location, lIic\ .IIl.iil. rectilinearly with the vor-
tex core. This also implies that the lift, drag, added mass,
pressure, and gravity forces are in balance. Finally, we
assume the vortex core is well represented by the ax-
isymmetric, Gaussian profile. This can be seen in fig-
ure 12c. The advantage to this assumption is that we
now have good estimates for velocity, vorticity, and dy-
namic pressure gradient in the undisturbed vortex core.
The individual terms in equation 33 can then be simpli-
fied as follows:

Du D(Ofpfu)
C -C = Pf t Dt PfOb-
Dt Dr
S-P VP- VP= 0 (34)
S-V pV (/ eVu ) pV2 u smallif p << 1
B -B -= fPfg pfg = -(bPfg
f (FD + FL + FAM + FP) FG bPbg

The change in fluid body force can be combined with the
interphase reaction term, f, to form a net buoyancy force
experienced by the fluid. Note that the interphase reac-
tion term will be small when the bubbles are not accel-
erating, explaining why the point particle two-way cou-
pling approach causes almost no vortex distortion. The
only other ,igiilik.iiI term which arises in this simplifi-
cation is due to changes to the dynamic pressure gradi-
ent, *, /r. If all terms are combined, and we multiply
through by the fluid volume, then the total reaction force
to the fluid because of n bubbles having volume Vb be-

AR l nVb + g(Pb
\ rs


This is in agreement with the general expression ob-
tained by Druzhinin & Elghobashi for their two-fluid
model, and used by and Sridhar & Katz in their dis-
cussion of bubble induced vortex distortion. In order
to find the magnitude of the total reaction to the fluid,
we must decompose the reaction into components. In a
cylindrical coordinate system located at the vortex cen-
ter we have.

R, = nVb

g(Pb i)sinl(Os)] (36)

Re = nVbg(pb -pfcos(0')

In the present clockwise vortices, the magnitude of this
reaction will be directed up and to the left from the set-
tling location in first quadrant. The magnitude and di-
rection of the net reaction measured from the horizontal



0 40




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

(0 = 0) will be

Rnet V R1 + Rg (38)

OR 0= ,+ + arctan R,) (39)

The net reaction, Rnet, represents a local input to the
fluid. Intuitively, the local force which drives the rota-
tion of the vortex should also dictate the magnitude of
the observed effects. We therefore choose to normalize,
R,,t, by the local vortex force, F'l a1. The idea of the
vortex force was originally put forth by Prandtl (1918)
and is frequently used when studying vortex dynamics
(see Saffman (1995)). It can be thought of as the force
required to maintain steady rotation of a vortex. It is
defined for a closed region as

F, = ux dV (40)

If we insert the expressions for the Gaussian vortex ve-
locity and vorticity fields given in the first test case,
we can obtain an expression for the total vortex force
required to drive the Gaussian vortex having initial
strength F0 and core radius, re.

F,_ AzF, / jr (r/r')2 -2(r/r,)2] d
Here, the cross product in the integral has been simpli-
fied because the velocity, ua is orthogonal to the vor-
ticity, cz at all points in the undisturbed, 2D vortex
core. Using this expression, a more helpful quantity is
the vortex force required to maintain rotation of a small
band centered at the bubble settling radius, illustrated
schematically in figure 13. By changing the integration
limits in equation 41 to represent a band of width 6, cen-
tered at the bubble settling radius, r8, we obtain

cal Az/2q|F

A v 2 erf( fr-)

6 6

- erfe V/2Tr
( r)

Note, that the integration has been post-mulitplied by
the scaling quantity re/6. This is to account for the ef-
fect that choice of bandwidth has on the magnitude of
the integral. By rescaling in this way, the values of vor-
tex force become independent of the choice of 6 so long
as 6 << rc.
Using the bubble settling coordinates predicted with
volumetric coupling, the net reaction, R,,t and local
vortex force, Fioz1, are calculated for each of the 12
cases with volumetric coupling. The net reaction is then
normalized by the local vortex force to compute a rela-
tive reaction force, R,,e = Rt/F 1al. This relative

Figure 13: Schematic of differential band used to inte-
grate the local vortex force.

reaction can be thought of as a potential for vortex dis-
tortion, or momentum redistribution. The relative de-
crease in vortex decay rate, E, is plotted against R,,~
in figure 14. The figure shows that the relative reac-
tion force is capable of both collapsing the values of E
and also producing a highly linear trend across the dat-
apoints. There are two important implications of this
result. First, we see that the local vortex force is a valu-
able scaling parameter for normalizing the interphase re-
action force. Determining a local measure of the vor-
tex force in this case was straightforward because of the
good match with the Gaussian flow field. For an arbi-
trary flowfield such as homogenous turbulence, a more
general approach would need to be taken. Second, we
have obtained a correlation between a large scale vor-
tex property, E, and a property associated with the local
bubble-scale interactions, R,,1. This could be coupled
with a model for the settling coordinates such as the one
developed in the stationary vortex section, to arrive at a
lower order model for bubble induced vortex attenuation
based solely on the parameters F0, r6, and db.

The present work focused on using a discrete element
model (DEM) for simulation of bubble or particle in-
teractions with vortical flows. The volumetric displace-
ment of the fluid by the subgrid scale dispersed phase,
along with point particle momentum exchange were
used for realistic coupling of the dispersed phase to the
fluid. The numerical implementation of the DEM model
was first evaluated for cases involving vortical flows and
dense loadings where the volume displacement effects
are considerably large: (i) a rising bubble plume, and
(ii) a falling blob of particles. The results show good
agreements with published data. Next, the effectiveness
and importance of volumetric coupling even under di-
lute loading was evaluated by simulate interactions of


0 -A
50 a-- -
E(%) ~ 4.5e5 Rre,

1 2 3
Rrei (x104)

Figure 14: Linear relationship between E and R,, for
each of the 12 cases simulated with volumetric coupling.

few small bubbles with a traveling vortex tube. This cor-
responds to a two-dimensional approximation of the ex-
periments on traveling vortex rings by Sridhar and Katz
(1999). It was shown that accounting for the volumetric
displacement of fluid by the dispersed phase is critical
in predicting the correct bubble settling location as well
as capturing the experimentally observed distortion. In
our simulations, we notice a secondary effect, which is
an increase in the angular momentum decay rate due to
the finite size of bubbles. To better understand the vari-
ation in the magnitude of these changes, a method for
determining a relative reaction force was developed us-
ing the idea of a local vortex force. This potential can be
thought of as the relative ability of the bubbles to cause
change to the vortex flow structure. For the Gaussian
type vortices studied in this work, it is based entirely
on the bubble size, vortex size, the vortex strength. The
results show a linear correlation between this distortion
potential and the observed relative increases in vortex
decay rate. This important relationship shows how local
bubble scale interactions can be related to vortex scale
strength attenuations. This study has important conse-
quences in understanding bubble interactions with co-
herent vortical structures in turbulent flows.


This work was supported in part by the Office of Naval
Research (ONR grant number N000140610697). SVA
also acknowledges support from Department of En-
ergy's National Energy Technology Laboratory (NETL)
and Oak Ridge Institute for Science and Research
(ORISE). All simulations were performed on the high
performance computing cluster at Oregon State Univer-

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


SV Apte, K. Mahesh, P. Moin, and JC Oefelein. Large-
eddy simulation of swirling particle-laden flows in a
coaxial-jet combustor. International Journal of Multi-
phase Flow, 29(8):1311-1331, 2003.

SV Apte, K. Mahesh, and T. Lundgren. Accounting
for finite-size effects in simulations of disperse particle-
laden flows. International Journal of Multiphase Flow,

S. Becker, A. Sokolichin, and G. Eigenberger. Gas-
liquid flow in bubble columns and loop reactors. II:
Comparison of detailed experiments and flow simula-
tions. Chemical engineering science, 49(24 B):5747-
5762, 1994.

E. Climent, M. Simonnet, and J. Magnaudet. Prefer-
ential accumulation of bubbles in Couette-Taylor flow
patterns. Physics of Fluids, 19:083301, 2007.

D. Darmana, NG Deen, and JAM Kuipers. Detailed
modeling of hydrodynamics, mass transfer and chemi-
cal reactions in a bubble column using a discrete bubble
model. Chemical Engineering Science, 60(12):3383-
3404, 2005.

D. Darmana, NG Deen, and JAM Kuipers. Paralleliza-
tion of an Euler-Lagrange model using mixed domain
decomposition and a mirror domain technique: Appli-
cation to dispersed gas-liquid two-phase flow. Journal
of Computational Physics, 220(1):216-248, 2006.

NG Deen, T. Solberg, and BH Hjertager. Numerical
Simulation of the Gas-Liquid Flow in a Square Cross-
sectioned Bubble Column. In 14th Int. Congress of
Chemical and Process Engineering. Praha-Czech Re-
public. Citeseer, 2000.

E. Delnoij, FA Lammers, JAM Kuipers, and WPM
Van Swaaij. Dynamic simulation of dispersed gas-liquid
two-phase flow using a discrete bubble model. Chemical
engineering science, 52(9):1429-1458, 1997.

OA Druzhinin and S. Elghobashi. Direct numerical sim-
ulations of bubble-laden turbulent flows using the two-
fluid formulation. Physics of Fluids, 10:685, 1998.

A. Ferrante and S. Elghobashi. On the physical mech-
anisms of two-way coupling in particle-laden isotropic
turbulence. Physics of Fluids, 15:315, 2003.

A. Ferrante and S.E. Elghobashi. On the effects of mi-
crobubbles on Taylor-Green vortex flow. Journal of
Fluid Mechanics, 572:145-177, 2007.

D. Gidaspow. Multiphase flow and fluidization: con-
tinuum and kinetic theory descriptions. Academic Pr,

G. Hu and I. Celik. Eulerian-Lagrangian based large-
eddy simulation of a partially aerated flat bubble col-
umn. Chemical Engineering Science, 63(1):253-271,

R. Jackson. Locally averaged equations of motion for
a mixture of identical spherical particles and a Newto-
nian fluid. Chemical Engineering Science, 52(15):2457-
2469, 1997.

D.D. Joseph, T.S. Lundgren, R. Jackson, and DA Sav-
ille. Ensemble averaged and mixture theory equations
for incompressible fluid-particle suspensions. Int. J.
Multiphase Flow, 16(1):35-42, 1990.

TS Lundgren. Slow flow through stationary random
beds and suspensions of spheres. Journal of Fluid Me-
chanics, 51(02):273-299, 1972.

I.M. Mazzitelli, D. Lohse, and F. Toschi. The effect of
microbubbles on developed turbulence. Physics of Flu-
ids, 15:L5, 2002.

C. Mitts, D. Talley, and D. Poulikakos. A fundamental
study of supercritical droplet deformation and breakup
through a miscible fluid analog. In AIAA, ASME, SAE,
and ASEE, Joint Propulsion Conference and Exhibit, 32
nd, Lake Buena Vista, FL, 1996.

P. Moin and SV Apte. Large-Eddy Simulation of Re-
alistic Gas Turbine-Combustors. AIAA journal, 44(4):

GF Oweis, IE van der Hout, C. Iyer, G. Tryggvason, and
SL Ceccio. Capture and inception of bubbles near line
vortices. Physics of Fluids, 17:022105, 2005.

A. Prosperetti and DZ Zhang. Finite-particle-size effects
in disperse two-phase flows. Theoretical and Computa-
tional Fluid Dynamics, 7(6):429-440, 1995.

PG Saffman. The lift on a small sphere in a slow shear
flow. Journal of Fluid Mechanics, 22(02):385-400.

PG Saffman. Vortex dynamics. Cambridge University
Press, 1995.

L. Schiller and A. Naumann. A drag coefficient correla-
tion. Vdi Zeitung, 77:318-320, 1935.

E. Shams. Numerical simulation of cavitating bubble-
laden turbulent flows, Ph.D. Thesis, Oregon State Uni-
versity. Oregon State University, 2010.

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

A. Sokolichin and G. Eigenberger. Applicability of the
standard k-E turbulence model to the dynamic simu-
lation of bubble columns: Part I. Detailed numerical
simulations. Chemical Engineering Science, 54(13-14):
2273-2284, 1999.

G. Sridhar and J. Katz. Drag and lift forces on micro-
scopic bubbles entrained by a vortex. Physics of Fluids,
7:389, 1995.

G. Sridhar and J. Katz. Effect of entrained bubbles on the
structure of vortex rings. Journal of Fluid Mechanics,
397:171-202, 1999.

MA van der Hoef, M. van Sint Annaland, NG Deen,
and JAM Kuipers. Numerical Simulation of Dense Gas-
Solid Fluidized Beds: A Multiscale Modeling Strategy.

E.A. Van Nierop, S. Luther, J.J. Bluemink, J. Mag-
naudet, A. Prosperetti, and D. Lohse. Drag and lift
forces on bubbles in a rotating flow. Journal of Fluid
Mechanics, 571:439-454, 2007.

JH Walther and P Koumoutsakos. Three-dimensional
vortex methods for particle-laden flows with two-way
coupling. Journal of Computational Physics, 167(1):
39-71, 2001.

DZ Zhang and A. Prosperetti. Momentum and energy
equations for disperse two-phase flows and their closure
for dilute suspensions. International Journal of Multi-
phase Flow, 23(3):425-453, 1997.

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