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
finnj@engr.orst.edu and sva@engr.orst.edu
Keywords: EulerLagrange, 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
EulerLagrange approach, volumeaveraged NavierStokes equations are solved using a colocated 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.
Introduction
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 EulerianEulerian
or EulerianLagrangian. In an EulerianEulerian or
twofluid 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
subgrid 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 subgrid 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
evaluated.
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 nonuniform 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 EulerianLagrangian 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 largeeddy 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 twoway 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
tions:
j (Xb)
d
mb d Ub)
EFb
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
FAM PfVbCAM (Dub
(Dt
Dub
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. Interbubble 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
24
CD (1 + 0.15Rc 687) (9)
Reb
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
CEo
D
Eo
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.22c3/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
0
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))
a
(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)
VP+
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
bubbles.
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 rewritten in a more commonly
used form by combining the first and last terms on the
righthand 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 gasfluidized beds (Gidaspow
(1994); van der Hoef et al. (2008)). In case of Reynolds
averaged NavierStokes equations or largeeddy simula
tions, the above equations should be timeaveraged or
spatially filtered using densityweighted 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;
7T
Fo
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, oneway 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 noslip 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 timestep. 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
addedmass 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
force.
FB sin 0s + FL FpD + FAM (21)
The dynamic pressure force and added mass force can
be combined,
9
FpD + FAM (1 + CAM) Pf (22)
rs
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)
3CDU2i
4dbg [b 1]
(1+ CAM)U2
[P 1] gsin(0,) +CLUoR7
If the flow field inside the vortex core is known (Gaus
sian, TaylorGreen, 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, 1way 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
V
mean
[mis]
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
Figure 3: Settling location of bubbles in the core of the
stationary Gaussian vortex: (.) analytical, (0) DEM.
All boundaries are assigned a noslip 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 XY 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
X
(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, counterclockwise 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 counterclockwise
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
04
050 '' 50 . .500 ''T0 2,D c
Time [s.
(b) Point 'B' vertical velocity
Figure 6: Fluid velocity at the points 'A' and 'B' shown
in Figure 4a.
ThreeDimensional 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 VDb
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 fuelair mixture.
They were able to classify the miscible droplet behavior
into several subregimes 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 multimode 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
(HR) 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 subregime. The
characteristic velocity, U, is defined using the HR 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 pas
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
Density
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
tU/D
Figure 7: Blob speed nondimensionalized 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 subregime 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 nondimensional 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)
5tfi
6I
(b) ' 
(d) 1 ,'l
(M)
\" ", .
,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 nondimensional 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 rollup 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
tions.
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 highshear 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 microbubbles 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
4P
(a) Initial Deformation
(b) Bag Vortex Ring RollUp
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 sixthorder polynomial
to match the experimental conditions. To keep the size
of the computation small and allow several parametric
studies, a twodimensional domain is simulated with pe
riodic conditions in the spanwise direction giving rise to
a vortextube. 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 subregime of the multimode
breakup category.
Parameter
Pf ,vf
Domain Size
Grid Size
Jet height (hjt)
Inflow Time
Vortex Strength,
Fo(m2/s)
Bubble Size (pm)
Inflow Velocity
a6, a5, a4,a3
a2, al, ao
Value
1,000 kg/m3; 106 m/s2
1m x 0.15m x 0.005m;
800 x 121 x 4
O.lm
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
nondimensional 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
X/hj,,
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
7
6
S5
r
x
4
02
1
3 A
: *
_* _
S AA
0.1 0.2 0.3
rs [cm]
0.4 0.5
Figure 11: Comparison of the nondimensional 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 twoway 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,= pfUorcV,, (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
d(L,.)
dt
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
(30)
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
twoway 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 bubblevortex 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, bubblefluid 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 EulerianLagrangian 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
as:
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 redistribution 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
variations.
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
(33)
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
SP VP VP= 0 (34)
SV 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 twoway 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
comes:
AR l nVb + g(Pb
\ rs
Pf)
This is in agreement with the general expression ob
tained by Druzhinin & Elghobashi for their twofluid
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
80F
60
0 40
20
x
x
c
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
(41)
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/2qF
A v 2 erf( fr)
6 6
 erfe V/2Tr
( r)
Note, that the integration has been postmulitplied 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
bubblescale 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.
Conclusions
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
1UU
0 A
U.,
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 twodimensional 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.
Acknowledgements
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
sity.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
References
SV Apte, K. Mahesh, P. Moin, and JC Oefelein. Large
eddy simulation of swirling particleladen flows in a
coaxialjet combustor. International Journal of Multi
phase Flow, 29(8):13111331, 2003.
SV Apte, K. Mahesh, and T. Lundgren. Accounting
for finitesize effects in simulations of disperse particle
laden flows. International Journal of Multiphase Flow,
2007.
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 CouetteTaylor 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 EulerLagrange model using mixed domain
decomposition and a mirror domain technique: Appli
cation to dispersed gasliquid twophase flow. Journal
of Computational Physics, 220(1):216248, 2006.
NG Deen, T. Solberg, and BH Hjertager. Numerical
Simulation of the GasLiquid Flow in a Square Cross
sectioned Bubble Column. In 14th Int. Congress of
Chemical and Process Engineering. PrahaCzech Re
public. Citeseer, 2000.
E. Delnoij, FA Lammers, JAM Kuipers, and WPM
Van Swaaij. Dynamic simulation of dispersed gasliquid
twophase flow using a discrete bubble model. Chemical
engineering science, 52(9):14291458, 1997.
OA Druzhinin and S. Elghobashi. Direct numerical sim
ulations of bubbleladen turbulent flows using the two
fluid formulation. Physics of Fluids, 10:685, 1998.
A. Ferrante and S. Elghobashi. On the physical mech
anisms of twoway coupling in particleladen isotropic
turbulence. Physics of Fluids, 15:315, 2003.
A. Ferrante and S.E. Elghobashi. On the effects of mi
crobubbles on TaylorGreen vortex flow. Journal of
Fluid Mechanics, 572:145177, 2007.
D. Gidaspow. Multiphase flow and fluidization: con
tinuum and kinetic theory descriptions. Academic Pr,
1994.
G. Hu and I. Celik. EulerianLagrangian based large
eddy simulation of a partially aerated flat bubble col
umn. Chemical Engineering Science, 63(1):253271,
2008.
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 fluidparticle suspensions. Int. J.
Multiphase Flow, 16(1):3542, 1990.
TS Lundgren. Slow flow through stationary random
beds and suspensions of spheres. Journal of Fluid Me
chanics, 51(02):273299, 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. LargeEddy Simulation of Re
alistic Gas TurbineCombustors. AIAA journal, 44(4):
698708,2006.
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. Finiteparticlesize effects
in disperse twophase flows. Theoretical and Computa
tional Fluid Dynamics, 7(6):429440, 1995.
PG Saffman. The lift on a small sphere in a slow shear
flow. Journal of Fluid Mechanics, 22(02):385400.
PG Saffman. Vortex dynamics. Cambridge University
Press, 1995.
L. Schiller and A. Naumann. A drag coefficient correla
tion. Vdi Zeitung, 77:318320, 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 kE turbulence model to the dynamic simu
lation of bubble columns: Part I. Detailed numerical
simulations. Chemical Engineering Science, 54(1314):
22732284, 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:171202, 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.
ANNUAL REVIEW OF FLUID MECHANICS, 40:47,
2008.
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:439454, 2007.
JH Walther and P Koumoutsakos. Threedimensional
vortex methods for particleladen flows with twoway
coupling. Journal of Computational Physics, 167(1):
3971, 2001.
DZ Zhang and A. Prosperetti. Momentum and energy
equations for disperse twophase flows and their closure
for dilute suspensions. International Journal of Multi
phase Flow, 23(3):425453, 1997.
