Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Particle Velocity near Vertical Boundaries A Source of Uncertainty in TwoFluid Models
Adelina P. Davis1, Efstathios E. Michaelides1 and ZhiGang Feng2
SDepartment of Mechanical Engineering University of Texas at San Antonio
San Antonio TX, 78259, USA,
2 Department of Mechanical and Energy Engineering, University of North Texas
Denton TX 76207 USA
ihm424@my.utsa.edu, stathis.michaelides@utsa.edu, feng@unt.edu
Keywords: DNS modeling, particle velocity, wall, solid boundary, boundary conditions
Abstract
One of the significant uncertainties in the use of the twofluid models is the solid phase boundary conditions. Typically the
fluid boundary condition at a solid wall is the so called noslip condition, which has become the cornerstone of singlephase
fluid dynamics. By analogy, in twofluid models the same noslip conditions are also applied as the boundary conditions of
the solid phase (or phases) at impermeable walls. However, there is not a physical reason as to why the noslip condition
should apply to fine or coarse particles and there are no detailed experiments or simulations that support such a hypothesis.
We have undertaken a study, to determine what the behavior of particles near solid impermeable walls is and what the velocity
of the particles is at the plane where their boundary conditions need to be prescribed (one particle radius far from the wall).
This will help prescribe boundary conditions for the solids phase and will remove a cause of the uncertainty in the gasparticle
modeling. The study is based on Direct Numerical Simulations (DNS) of a dense particulate mixture and the results pertain
to both dilute and dense mixtures. Since the model of particlewall interactions/collisions is important in the results, the
computations include the parameters used in the wall interaction model of DNS studies. The results show that the
solidsphase velocity very close to a vertical wall is clearly negative (in the direction of gravity) and that the absolute value of
the velocity depends on several parameters including the particlewall interaction parameters.
Introduction
Simulations of particulate flows are very important due to
its wide applications in engineering and they take places in
many natural processes (Azpitarte & Buscaglia 2003, Feng
& Michaelides vol. 52 2009, Feng & Michaelides 2005,
Uhlmann 2009, Feng & Michaelides 2004). Due to its vast
occurrence and application, the simulations by different
techniques have taken place and are presented in many
studies, some of which are described below. Some
techniques described by (Feng & Michaelides vol. 52 2009,
Gan, Chang & Howard 2003, Kim & Choi 2004, Feng &
Michaelides 2008) have taken into consideration thermal
convection that occurs during the process. Various studies
have been done to describe twophase fluid flow by
twofluid model (Azpitarte & Buscaglia 2003, Biswas,
Esmaeeli & Tryggvason 2005, Liu, Wang, Gong & Peng
2008, Chahed, Roig & Masbernat 2003, Sokolichin,
Eigenberger, Lapin & Lfibbert 1997, Wilde, Vierendeels,
Heynderickx & Marin 2005). Biswas et al. (2005)
compared the twofluid model for laminar flow with direct
numerical simulation (DNS) method of bubbly flows and
received good results between the solutions with the
necessary adjustments. Chahed et al. (2003) presented a
study on EulerianEulerian twofluid model in which
limitations of its application are discussed. Cantero,
Balachandar & Garcia (2008) introduced EulerianEulerian
model in which it took into account the inertia and particle
settling.
One of the more recent studies that was done on direct
numerical simulation (DNS) of twophase fluid flow with a
fictitious domain was implemented by (Yu & Shao 2009), in
which spherical and nonspherical particles were
investigated. Numerous studies have been done on
implementation of noslip boundary conditions (Goldstein,
Handler & Sirovich 1993, Wang, Ge & Li 2006, Zhang &
Prosperetti 2003, Feng & Michaelides vol. 38 2009). In
Goldstein et al. (1993), boundary conditions are modeled
with an external force. Wang et al. (2006) in their study
described three traditional ways in which particlewall
interaction are done: 1) Specular reflections, 2)
Bounceback reflections, and 3) Maxwellian reflection. In
numerical studies, the noslip boundary conditions are
applied for the twophase flows at the impermeable wall.
Massoudi (2007) has examined the importance of boundary
conditions in mixture theory and in Computational Fluid
Dynamics (CFD) application.
DNS method to simulate the fluidized bed reactor is applied
for 264 particles of d = 0.6 cm and for 420 particles of d =
0.4 cm and the behavior of solid particles at the
impermeable wall and near the wall is evaluated in this
study. To determine what the behavior of particles near
solid impermeable walls is and what the velocity of the
particles is at the plane where their boundary conditions
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
need to be prescribed (one particle radius far from the wall)
is given in this study. This will help prescribe boundary
conditions for the solid phase and will remove a cause of the
uncertainty in the gasparticle modeling.
DNS method and governing equations
The DNS method that is applied in this study is the
Immersed Boundary method (Feng & Michaelides vol. 52
2009, 2008), however the isothermal process has been
assumed throughout computations. Governing equations
are given below.
Governing equations in dimensionless form
Dimensionless terms in the numerical simulations given
below for the momentum equation:
u I t p
u' = =  = t' =
Uref Lref Lref Pjo Ure'
frere
Pro Uref
And for the energy equation:1)
PfoUrf
And for the energy equation:
0 TTr ,q' P= Lref ,
Tref PfocfTrefUref
ALref
PfocfTrefUref
Where u, I, t, p, f &, T q, k are velocity, length, time,
pressure, body force, dimensionless temperature,
temperature, energy source, and energy density, respectively.
Uref, Lref and Tref are the reference velocity, lengthscale and
temperature. The fluid density at the ambient temperature
is pfo and cf is the specific heat of the fluid. The convective
timescale (LreUref) is used as the reference timescale for
the problem. Tfo is the ambient fluid temperature, and in
most cases, Tref is chosen to be the initial temperature
difference between the particles and the ambient fluid.
The flow is defined by the following dimensionless
parameters:
Pr = p cr = C (3)
Pfo Cf
Where Re is the Reynolds number; Pr is the Prandtl number;
Pe is the Peclet number; and Gr is the Grashoff number and
is evaluated as follows:
Pfureflref cf
Re= rere ,Pr = ,Pe = PrRe,
1Lf kf
Gr fLref (T TO) (PfUrefLref) (
Gr = ) (4)
In which pp is the particle density and pff, vf, kf, fif are
density, dynamic viscosity, kinematic viscosity, thermal
conductivity, and thermal expansion for the fluid.
Gravitational constant is g.
The following set of dimensionless governing equations
may be obtained that describe the entire domain 0 of the
fluid and of the entire particles:
A.) For the velocity field:
+u. Vu = Vp+ 2+ Gre +, E (5)
at Re Re2 (5)
where eg is the unit vector in the direction of gravitational
acceleration.
B.) For the force density field:
f + u V + p V2 e, XE XZSi (6)
at Re Re 2 91
C.) Continuity equation:
V = 0, E n
D.) For the velocity field inside the solid particle region:
u = Up + xp x (X tp) (8)
E.) For the temperature field:
Energy equation:
0 + u VO = V20 + q + X, X E
at Pe
Energy density function:
S= 00+ ul O V20 q, E CSi
at Pe
F.) For the motion of the particles:
gref
(Pr 1)p u fdv + (Pr )p e ( 11)
(Pr 1) = (x p)X dv (12)
In two dimensions we have the following closure
expressions for the volume of particles and for the moment
of inertia, where r represents radius of cylindrical particles:
V, = 7r2; Ip = rr4 (13)
For the particle temperature the energy balance equation for
particle is used, but in this study the temperature of particles
is constant:
PrVpC, r = i 0 fids + fqpdv (14)
The uniformity of temperature in the solid region is
enforced by the condition:
0 = Op(t) (15)
The initial conditions of the above set of equations are:
Op(t = 0) = Opo; p(t = 0) = 0po; u(t = 0) =
u0; O(t = 0) = 00; Op(t = 0) = Opo (16)
The subscript 0 represents the initial value or the initial
ambient state.
All the above equations have been solved numerically by
finite difference based scheme in (Feng & Michaelides
2008).
Numerical Scheme
For the performed simulations the following dimensions of
fluidized bed reactor were chosen x = 15 cm and y = 45 cm.
Jet width is taken as 4 cm. The fluid density was taken as
p = 1.0 g/cm3 and the viscosity at t = 0.01 g/(s cm). The
number of particles in the simulation is N = 264, where the
diameter is constant and equal to d = 0.6 cm. When the
simulations for this diameter were done the simulations was
repeated for N = 420 particles with d = 0.4 cm. For this
particle diameter of d = 0.4 cm the velocity of the jet was
chosen as 80 cm/sec, and for the diameter of d = 0.4 cm jet
velocity was taken as 60 cm/sec. Dimensionless
parameters of the flow are: density pr = 2.3, specific heat
capacity cpr = 0.01, and thermal conductivity kr = 7.5. The
simulation was done for colder particles, and for this reason
Grashof number was taken as Gr = 100. Prandtl number
is Pr = 0.7. The total time of the simulation is t = 40 sec
and the time step was chosen as t = 0.02 sec (dt = 0.0001).
In order to see the behavior of particles at the impermeable
walls all the simulations were performed one radius away
from the wall and data for analyses was selected for d = 0.6
cm from 0.3 cm to 1.5 cm at the left border and 13.5 cm to
(10)
Paper No
Paper No
14.7cm at the right border, and for d = 0.4 cm from 0.2 cm
to 1.0 cm at the left border and 14.0 cm to 14.8 cm at the
right border. All the particles have been numbered and
tracked throughout the computations.
Simulations were conducted for a number of c, factors,
which is the force scale factor (Feng & Michaelides 2005).
For every case, velocity in y direction was taken and
analyzed. Velocity that was received from the simulations
were dimensionalized and all the data is presented below.
In order to dimensionalize the velocity, terminal velocity of
particles needs to be evaluated:
Knowing that drag force and gravity force are equal FD = FG,
the following equation could be written:
1 rdPf2 rfd3 (17)
Sp VTCD (pp p)9 (17)
Where CD and Re is the drag coefficient and Reynolds
number respectively, is evaluated as following:
24
CD = (1 + 0.15Re0.687) (18)
Re
Re =fVd (19)
When numerous numbers of iterations was performed the
terminal velocity was obtained. The average flow velocity
was evaluated as following:
V1 = Viet Aje (20)
Bottom
where Ajet is the jet flow cross sectional area and Abottom is
the cross sectional area of the bottom part in the fluidized
bed reactor.
The dimensionless velocity is evaluated as following:
va (21)
VT
where v, is the particle velocity obtained from the
simulations in y direction.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
particlewall interaction.
1) Collision between particles
When the gap between particles is less than a given
threshold repulsive force is applied and that force is added
to a total force that particle experiences:
(o, I1Xi x, > R, + R, + (
F C = (c IIIxx, JRR 2 x( x ,x, R,+ (24)
where c, is the force scale factor, cp is the stiffness
parameter for collisions, R, and R, are the radii of the two
particles, and is the threshold or "safe zone".
1) Collision with the wall
This method is introduced by repulsive force between the
particle and the wall is as follows:
(0, I1x, xJ, I > 2R, + <
F (,,12R,)(2,_ (25)
c ) ( 1x, x,J < 2R, + (
where e, is another stiffness parameter; x, is the position of
fictitious particle P,,, which is located symmetrically on the
other side of the wall W.
Results and Discussion
We conducted numerous simulations to analyze the behavior
of particles at the boundary.
A.) Trajectories
In order to have the visual understanding of the particulate
flow near the boundary, the trajectories for the first 30
particles in the reactor were plotted to see the motion. Its
behavior is shown for the time of 20 seconds. Fig. 1
represents the trajectories of particles during the simulation.
Collision models and force scale factor
Feng & Michaelides (2004) discussed the restoration force
that depends on the parameter c, which is the force scale
factor. The restoration force was developed to restore the
boundary points back to the reference points and is
modelled by the linear springlike relation:
Fi = kuij (22)
where k is the spring constant and Fij = xij xb is the
displacement; x, is the boundary point and x', is the
reference point. The restoration force could be rewritten in
the following form:
0, I1Fll = 0o,
F, = c I1,11 (23)
Ed Iull' > 0,
in which ed is the stiffness scale factor for particle
deformation and it is proportional to inverse of spring
constant. Force scale factor is estimated as buoyancy force,
and it is a reasonable assumption for most of particulate
flows:
Cij = Pf(a 1)Vg (24)
where V, is the particle volume and g is the gravitational
acceleration. However, in some general cases a scaling
force factor is defined by the user. The value of the scaling
factor is of the order of magnitude of the characteristic force
for the process acting on one particle.
In order to avoid particle overlap and introduce wall
reflections in the computations, the following two
techniques were introduced: particleparticle interaction and
0 0.3 0.6 0.9
x, cm
0 3
S#75
S #8
#9
~ #t 11
# 12
t 13
t 14
 15
# 16
4#t 17
#18
_* #23
a# 24
#25
#26
1.2 1.5  #28
Figure 1: Trajectories of the particles for the first 20 sec of
the simulation.
B.) Full domain snapshots
To plot the velocity as a function of position the snapshots
were taken at different time: 16 sec, 18 sec and 20 seconds.
Velocity of the particles was taken dimensionless and the
width of the reactor was also taken. Fig. 2 represents the
snapshot at t = 16 sec of all the particles taken part in the
simulation process. As can be seen, some particles are still
at rest and some are moving with high velocities.
Paper No
1.5
1
0.5
oJ
0
.2
E 0.5
1
1.5
2
x dimensionless
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0.6
y = 3.936x4 8.014x3 + 4.904x2 0.795x 0.372
0.4
0.2
0 . .1 . . 0 . . . .
0.2 0.2 0.4 0.6 0.8 1
0.4
*x di
x dimensionless
Figure 2: Full domain of the particles with d
time t=16 sec.
0.6 cm at
Fig. 3 represents the full domain of the simulator at time t =
18 sec. From comparison of Fig. 2 and Fig. 3, it can be
noticed that particles get displaced with time. Some of the
particles settle down and velocity changes as time elapses.
1.5 r
1.5 1
Figure 3:
t=18 sec.
y = 1.070x4 5.330x3 + 5.495x2 1.428x 0.308
. .*
1
x dimensionless
Full domain of the particles d = 0.6 cm at time
Figure 4: Full domain of the particles d = 0.6 cm at time
t=20 sec.
C.) Average velocity with all particles
When the simulations for different c, factors were
performed the data was extracted for the discrete even value
time period of t = 16 sec to t = 30 sec (at an interval of the
first 0.2 sec with a time step of 0.02; i.e. 16.016.2,
18.018.2, etc). The data contains the values for all
particles that are located close to the borders, and the flow at
these times is considered fully developed. The particles at
the borders were followed and their velocity versus position
is plotted accordingly. Plotting and evaluating the
trendline for different graphs helped to calculate the velocity
at the border (at 0.3 cm from the left and at 14.7 cm from
the right side). That is the reason why there are so many
values of velocity per one value of c, coefficient (Fig. 5).
2
S4
E
U
>6
8
10
12
At time of t = 20 sec, less particles have positive velocity
and more particles are settled down. From time to time the
particles change their position in the domain and their
behavior at the walls of the bed reactor could be evaluated at
different times.
5,001,000 B0,000,000 1 00,000
A
A x a
* U
1,000,000
M 3,000,000
A A 5,000,000
A X 8,000,000
x X 10,000,000
0 14,000,000
cij
Figure 5: Velocity for all particles d = 0.6 cm versus c.
Dimensionalizing the velocity of particles, as described
earlier, we obtained Fig. 6.
y = 11.67x4 25.82x3 + 17.60x2 3.533x 0.273
S*.
A
Paper No
5,000,000 10,000,000
1,000,000
3,000,000
A 5,000,000
X 8,000,000
X 10,000,000
14,000,000
15,000,000
0.1 +
0
S0.2
E
a?
bf0.3
0.4
f
y= 5.2967E10x 3.8253E01
t~ 9
ip 4
0.5 
>o5 ; i A i
0.5 I f
x
0.6 cij
Figure 6: Dimensionless velocity for all particles d = 0.6
cm versus co.
The average value of velocity was evaluated and is plotted
versus c. in Fig.7. Plotting the trendline and its equation
helps us understand how velocity varies with cy. The
sample standard deviation was also evaluated and its value
is s = 0.0089.
0 0.1
c
0
50.2
E
S0.3
C0.

cu
>0.4
5,000,000 10,000,000 15,000,000
I
y = 1.34727E09x 4.09610E01
0.5 L
Figure 7: Average dimensionless velocity for all particles d
= 0.6 cm versus c .
The simulations also were done for the particles with a
smaller diameter d = 0.4 cm for which similar analysis were
performed. In Fig. 8 is shown the dependence of average
velocity of all particles that took place in the simulation and
the c. factor.
Figure 8: Average dimensionless velocity for all particles d
= 0.4 cm versus c.
In both cases could be noticed that velocity does not depend
on the c. factor and it lies in the range of 0.38 0.41.
D.) average velocity with suspended particles
As seen from the Figs. 24, some particles are settled
throughout the simulations and some particles are
suspended. Fig. 9 represents the values of velocity that
was taken only for suspended particles versus c. factor. As
in previously described procedure the trendline for different
times was built and the value for the 0.3 cm and for 14.7 cm
of velocity was evaluated.
S0.2 5,000,000 10,000,000 15,000,000
0.2 
0.4 I I
0.6 I x o
1,000,000
0.8 x 3,000ooo,ooo000
A A 5,000,000
1 
x 8,000,000
1.2 A X 10,000,000
S14,000,000
1.4 cij
Figure 9: Dimensionless velocity for suspended particles d
= 0.6 cm versus co.
Taking the average of the velocity and plotting it versus the
value of c. factor at which the simulation was taken place
the following Fig. 10 was received. The sample of
standard deviation in this case is s = 0.0263.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0 0 I 0. I . 0 00
5,000,000 10,000,000 15,000,000
0.1
S0.2
0 
a0. 3
"0.3
E
0.4
Paper No
0
0.1
0.2
0.3
0.4
0.5
5,000,000 10,000,000 15,000,000
y = 2.94977E09x 4.74440E01
 *
0.6 L
Figure 10: Average dimensionless velocity for suspended
particles versus c. for d = 0.6 cm.
Conducting the same analysis for the d = 0.4 cm and taking
the velocity only of suspended particles in consideration the
Fig. 11 was obtained. Connecting all obtained values for
average velocity with a straight line could be noticed that
velocity for suspended particles varies with c, factor
insignificantly.
0.1
0.2
0.3
0.4
0.5
0.6
S 5,000,000 10,000,000 15,000,000
y = 4.1557E09x 4.6162E01
 ^ ~
0.7 L cij
Figure 11: Average dimensionless velocity for suspended
particles versus c. for d = 0.4 cm.
From Figs. 9 and 10, could be noticed that the average value
of the velocity for suspended particles lies in the range of
0.46 0.47.
Conclusions
In this study behavior of particles for the diameter of d = 0.6
cm and d = 0.4 cm is evaluated and the simulations of the
fluidized bed reactor were performed. Conducting
numerous simulations for the different force scale factors c,
shows that the velocity at the boundary of solid particles
exists and it is clearly negative and is directed towards the
gravity. For different c, factors the velocity at the
impermeable walls varies, however, straight correlation can
be made. Assuming there is a noslip boundary condition
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
for the solid particles at the wall brings an uncertainty in the
gasparticle modeling.
As seen from obtained results for these two cases velocity
does not depend on the c, factor and it lies in the range of
0.38 0.41, if the velocity of all the particles is taken in
consideration. But the average value of the velocity for
suspended particles lies in the range of 0.46 0.47.
The results of this study can help prescribe boundary
conditions for the solid particles at the walls.
Acknowledgements
This work was partly supported by a grant from the DOE,
through the National Energy Technology Laboratory,
(DENT0008064) Mr. Steven Seachman project manager;
and by a grant from NSF (HRD0932339) Drs. Demetris
Kazakos and Richard Smith, project managers.
References
Azpitarte, O.E. & Buscaglia, G.C. Analytical and numerical
evaluation of twofluid model solutions for laminar fully
developed bubbly twophase flows, Chem. Eng. Science 58
(2003)37653776.
Biswas, S., Esmaeeli, A. & Tryggvason, G. Comparison of
results from DNS of bubbly flows with a twofluid model
for twodimensional laminar flows, Int. J. Multiphase Flow
31(2005)10361048.
Cantero, M.I., Balachandar, S. & Garcia, M.H. An
EulerianEulerian model for gravity currents driven by
inertial particles, Int. J. Multiphase Flow 34 (2008)
484501.
Chahed, J., Roig, V & Masbernat, L. EulerianEulerian
twofluid model for turbulent gasliquid bubbly flows, Int. J.
Multiphase Flow 29 (2003) 2349.
Feng, Z.G. & Michaelides, E.E. The immersed
boundarylattice Boltzmann method for solving
fluidparticles interaction problems, J. Comput. Phys. 195
(2004) 602628.
Feng, Z.G. & Michaelides, E.E. Proteus: a direct forcing
method in the simulations of particulate flows, J. Comput.
Phys. 202 (2005) 2051.
Feng, Z.G. & Michaelides, E.E. Inclusion of heat transfer
computations for particle laden flows, Phys. Fluids 20
(2008) 110.
Feng, Z.G. & Michaelides, E.E. Heat transfer in particulate
flows with Direct Numerical Simulation (DNS), Int. J. Heat
and Mass Transfer 52 (2009) 777786.
Feng, Z.G. & Michaelides, E.E. Robust treatment of noslip
boundary condition and velocity updating for the
latticeBoltzmann simulation of particulate flows,
Computers & Fluids 38 (2009) 370381.
Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Gan, H., Chang, J.Z. & Howard, H.H. Direct numerical
simulation of the sedimentation of solid particles with
thermal convection, J. Fluid Mech. 481 (2003) 385411.
Goldstein, D., Handler, R. & Sirovich, L. Modeling a
noslip boundary with an external force field, J. Comput.
Phys. 105 (1993) 354366.
Kim, J. & Choi, H. An immersedboundary finitevolume
method for simulations of heat transfer in complex
geometries, Korean Soc. Mech. Eng. Int. J. 18 (2004)
10261035.
Massoudi, M. Boundary conditions in mixture theory and in
CFD applications of higher order models, Computers and
Mathematics Applications 53 (2007) 156167.
Liu, S., Wang, Z., Gong, Z. & Peng, Q. Simulation of
atmospheric binary mixtures based on twofluid model,
Graphical Models 70 (2008) 117124.
Sokolichin, A., Eigenberger, G., Lapin, A. & Libbert, A.
Dynamic numerical simulation of gasliquid twophase
flows Euler/Euler versus Euler/Lagrange, Chem. Eng.
Science 52 (1997) 61126.
Uhlmann, M. An immersed boundary method with direct
forcing for the simulation of particulate flows, J. Comput.
Phys. 2009 (2005) 448476.
Wang, L., Ge, W. & Li, J. A new wall boundary condition in
particle methods, Computer Phys. Communications 174
(2006) 386390.
Wilde, J.D., Vierendeels, J., Heynderickx, G.J. & Marin, G.B.
Simultaneous solution algorithms for EulerianEulerian
gassolid flow models: Stability analysis and convergence
behaviour of a point and a plane solver, J. Comput. Phys.
207(2005)309353.
Yu, Z. & Shao, X. Direct numerical simulation of particulate
flows with a fictitious domain method, Int. J. Multiphase
Flow (2009), doi:10.1016/j.ijmultiphaseflow.2009.10.001.
Zhang, Z. & Prosperetti, A. A method for particle simulation,
J. Appl. Mech. 70 (2003) 6474.
