7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Modelling of dense particleladen pipe flow and a free jet
G. Mallouppas and B.G.M. van Wachem*
Department of Mechanical Engineering, Imperial College London, London, SW7 2AZ, United Kingdom
george.mallouppas05@ imperial.ac.uk and b.vanwachem@imperial.ac.uk
Keywords: Lagrangian simulations, particleladen gas, LES, Discrete Element Method
Abstract
In this contribution, we have performed Large Eddy Simulations (LES) of the dense particle laden turbulent flow
as described by prof. Jennifer Sinclair Curtis (University of Florida, USA) and coworkers. The fluid phase is
modelled with three different LES models, i.e. a Smagorinsky model where van Driest dampening is employed
near the walls, and the more recent WALE and Vreman models. Also, the effect of mesh spacing on the flow
prediction is researched. The particles are modelled by the Discrete Element Model (DEM), in which the particle
collisions are approximated by their local elastic and plastic deformation. The importance of collisions is shown by
comparing the simulations taking particle collisions into account to simulations neglecting this effect. The effect
of the geometry is also researched. Due to the dimensions of the box in which the turbulent pipe flow enters,
it is expected that this dimension has a significant effect on the spreading rate of the turbulent gasparticle jet.
Finally, the particle size is varied, in accordance with the available experimental data. The particles used in the simu
lations range from 25 pm to 70 pm as well as a mixture of these. All findings are compared with the experimental data.
Nomenclature
Roman symbols
c Vreman constant ()
d diameter (m)
g gravitational constant (nms1)
p pressure (Nm 2)
m solids loading ()
v velocity (ms 1)
u, friction velcity (ms 1)
CD Drag coefficient ()
CSGS Smagorinsky constant ()
CvD Smagorinsky van Driest constant ()
C, Wale constant ()
E Young's modulus (Nm 2)
G shear modulus (Nm 2)
R pipe radius (m)
Re Reynolds number ()
Sf Source term linear in velocity (kgm 3 s1)
T Torque (Nm)
T3 Additional source terms (kgm 2s 2)
Greek symbols
a local volume fraction ()
3 drag function (kgm 3s 1)
K von Karman constant ()
p viscosity (Nsm 2)
v kinematic viscosity (n2s 1)
T Stress tensor of the fluid (Nm 2)
p density (kgm 3)
a Poisson's ratio ()
c angular velocity radss 1)
A filter width (m)
Subscripts
f Gas
p Particle
Introduction
Particleladen turbulent flows can be found in various
industrial and environmental processes. Examples
of such processes are pneumatic transport, fluidised
beds; energy generation of fossil fuels; movement of
soot particles in the atmosphere; the flow of particles
in cyclones and many more. Understanding particle
particle, particlefluid and particlewall interactions is of
outmost importance as a more accurate implementation
of these processes will result in more reliable analysis
tools of such flows. Reliable numerical simulations
will therefore help the optimisation and better design
of industrial processes and also reliable prediction of
environmental processes.
There are various frameworks in which gassolid
flows can be predicted in the EulerianLagrangian
framework, e.g. Large Eddy Simulation (LES), Direct
Numerical Simulation (DNS) and the Reynolds Aver
aged NavierStokes (RANS) method. RANS methods
do not resolve the small scales, which are important
for turbulent flows, but only consider the averaged
quantities of the transport equations. The smaller scales
are modelled by appropriate closure models. DNS
methods, on the other hand, offer much higher accuracy
in solving the NavierStokes equations at the expense
of huge computational time; currently, DNS can only
solve flows of low Reynolds (Re) numbers, which are
outside of most the engineering and industrial interests.
Although the computational effort for LES is still
very high, it is considerably lower than for DNS and
therefore has become very fashionable for analysis of
flows in research and is an emerging tool in industrial
problems.
LES solves the NavierStokes equations up to a
particular lengthscale due to the application of a filter.
Lengthscales smaller than the cutoff width (A) of the
filter are modelled. This cutoff width is an indication
of the size of eddies that are retained in the computa
tions and the size of eddies that are rejected Versteeg
and Malalasekera (2007). Due to the filtering of the
NavierStokes equations, models are required to resolve
the subgridscale (SGS) stresses to account for the
effect of unresolved scales on the convective momentum
transport. In this paper, three different models are used
and compared with one another. These models are
the well known Smagorinsky model where van Driest
damping is employed near the wall; the WallAdapting
LocalEddy (WALE) viscosity model proposed by
Nicoud et al Nicoud and Durcos (1999); and finally
the Vreman model proposed by Vreman Vreman 21 11 14).
The particle interactions are modelled by the Dis
crete Element Model (DEM), first applied by Cundall
and Strack Cundall and Strack (1979). The particle
particle and particlewall collisions are approximated
by the elastic and plastic deformation occurring during
collision. Such deformation is mathematically described
by a springdashpotslider model, or softsphere model,
as proposed by Tsuji et al. Tsuji et al. (1991).
Due to the large number of particles, the compu
tational time and memory space required for these
simulations is very large. This is because computations
for each particle take place for every timestep. On
the other hand, DEM has a very high accuracy as no
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
empirical data, besides from material properties, are
required to predict the particleparticle and particlewall
interactions. The model is also fully deterministic 
making an accurate comparison of the effect of initial
conditions possible. To describe the fluidparticle inter
actions, the Wen and Yu Wen and Yu (1966) drag force
has been employed in this work. Although accounting
for the volume fraction effects on the drag force is
probably not of large importance in the testcases
simulated in the current research work, it might have an
effect on the particleclustering. The forces resulting
from the fluidparticle interactions, the particleparticle
interactions and the particlewall interactions are added
and used to track the individual particles employing the
Verlet integration scheme.
Theory
Fluidphase modelling
The equations for the fluidphase are similar to the
singlephase, except for including a local volume frac
tion and including source terms accounting for the pres
ence of particles. The governing equations are the mass
continuity and conservation of momentum, which are
given by equations (1) & (2) respectively,
a(af f) 9(a ff vy)
at 9xi
a(afpfv3) a(afPfvVp v (afT3T) 9p
at f 9x1 a9x
at + 0xi 0x2 f XJ
phases f
+Sf +Ti+ (,p) [v v]
p=1
(2)
where af is the local volume fraction of the fluid phase;
Sf the source term linear in the velocity field; T3 the
additional source terms; and 7)3 the stress tensor of the
fluid, given by
' 9 vi a vi (
Tf 8 'I xi V
2f ) a (3)
3 8x"V
where pf is the shear and Af represent the bulk viscosity
of the fluid.
Subgrid Scale Models
Large Eddy Simulation (LES) is used to solve the
masscontinuity and momentum equations, models are
therefore needed to determine the subgridscale (SGS)
stresses. Three models are used and compared with ex
perimental data. These are the Smagorinsky model with
van Driest damping near the wall; the WALE model
Nicoud and Durcos (1999); and finally the Vreman
model Vreman (21' '4). It is worthwhile mentioning that
the models are compared with one another by consid
ering two different geometries, by varying the length of
the inlet pipe. Predictions have been made for an inlet
pipe of length L = 20.0 and of L = 90.0. In addition,
different mesh spacing is researched in order to examine
the effect on the flow prediction.
Van Driest Wall damping
It is wellknown that the Smagorisnky model is not
suitable to account for the effect of walls, as the strong
velocity gradient perpendicular to the wall is present
due to the noslip boundary condition at the wall and
not a gradient driving the turbulence. Hence, near the
wall the socalled van Driest dampening is employed
which is based upon a simple mixing length model.
Prandtl's mixinglength hypothesis Wilcox (1998);
Pope (2000) states that,
dU
x7, = Ty (4)
where y is the length from the wall, in the case of pipe
flow given by,
y (R r) (5)
and pT represents the eddy viscosity, which is
t c2 dU
fiT i t dy
where l~ is the mixing length. The mixing lei
scale is approximated as
/mix = ly
Near the wall, however, ,ix is not well described
Equation (7). This is beacuse the mixing length me
'seems incapable of reproducing the y+3 variation r
the wall' Grifoll and Giralt (2000). The mixing length
therefore, modified by the van Driest damping funct
This function dampens the effect of viscosity near
wall. Hence equation (7) is modified to equation (8)
mix ts[1 c /A o]
where y+ is the dimensionless distance defined as,
and A+ is a constant, typically A+ 26. Although
the model seems very crude, it has proven to be quite
successful in many types of wall bounded turbulent flows.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Smagorinsky model
The Smagorinsky model assumes that the local SGS
stresses are proportional to the local rate of strain of the
resolved flow Versteeg and Malalasekera (2007). Thus,
the stresses are given as
Tij = 2psGsSij + 3T16Sij (10)
where Sj = + %j psGs is estimated using
equation (11)
PSGS = Pf (CsGSA)2 2 (11)
where Cscs is the Smagorinsky constant, A is the LES
filter width.
The Smagorinsky van Driest constant The
Smagorinsky model, however, does not accurately rep
resent the eddy viscosity near the wall. As shown by
Weickert et al Weickert et al. (2009), the CSGS is modi
fied as,
C,D =CsGsA(1 A) (12)
In the simulation the following constants were used;
CSGS 0.142 and A+ 25.0. Note that + = u'"
Vf
asG is,
ISGS PfC 2Sj S (13)
WALE model
Nicoud et al Nicoud and Durcos (1999) proposed
igth the walladapting local eddyviscosity (WALE) model.
With this model, the nonphysical behaviour of the eddy
viscosity as predicted by the Smagorinsky model near
(7) the wall is avoided. This is because the WALE model
Sby is of higher order compared to the Smagorinsky model,
)del which is 0(1). In this way, the eddy viscosity at the wall
ear is zero and hence, the only viscosity contribution at the
h is, wall is the molecular viscosity. The model for the eddy
ion. viscosity is given in Nicoud and Durcos (1999) and We
the ickert et al. (2009),
Id d j
(SG) s pf (C'A)2 (S (14)
(8) (S (+14
C, is related to the Smagorinsky constant via,
C. = (, G ) (S ad)3)
(SiJSi (Sdd
withl  4j Y and./. 
with^d 1 +^)1 ._9 a .
Vreman model
Vreman Vreman (2 '14) proposed an eddy viscosity
model used in LES, which is also of higher order accu
racy and produces the correct limit of the eddy viscosity
near the walls. This model assumes that the dissipation
is relatively small in transitional and nearwall regions.
The eddy viscosity is calculated as
Psos = PfC Bc (16)
with ,/3ij Ai and B /3 1/322 
12 +/311/33 3s +/322/333 /33. The model constant c
is related to the Smagorinsky constant as c z 2.5CSS.
Filtering Method
LES utilises lowpass filters so that only the large scale
eddies are resolved, whereas the smaller scale eddies are
removed. In the aforementioned models, a filter width
is used in order to define the cutoff lengthscale. In this
paper the volume of the cells used in the meshing is di
rectly related to the filter width (A), i.e. A .
As for successful LES calculations, the computational
cell aspect ratio has to remain near unity, this is a good
approximation.
Particle modelling
A discrete element model (DEM) proposed by Cun
dall and Strack Cundall and Strack (1979) is used to
model the particles. The individual trajectories of the
particles are determined in the Langrangian framework,
where particle collisions are modelled via the soft
sphere model proposed by Tsuji et al Tsuji et al. (1991),
which accounts for some nonreversible deformation.
Discrete Element Models (DEM) The interactions of
particles with other particles and walls are of dynamic
nature. This is because the particle movements are
essentially defined by the partcleparticle interactions;
particlewall interactions; particlefluid interactions
and/or body forces Kuang et al. (2008). The trajectories
of individual particles are considered (i.e. described
in a Langrangian framework) and Newton's 2"d law is
solved for each individual particle, accounting for the
fluidparticle, particleparticle, and particlewall inter
actions, and approximating the integral with the Verlet
algorithm. The flowchart shown in figure 1 illustrates
the DEM algorithm. The flowchart describes the logic
of how forces between individual particles with other
particles and/or walls of the domain are determined. In
addition, it illustrates how forces between particles and
the fluid are computed.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
START
1
READ P*.& Pmlu
DOES
NO
OVERLAP EISTWITH
PARTICLES IN
SAME/OR NEAREST /
PRTCLE CELLS
AMDIOR WAL.
YES
COMPUTE FORCE DUE TO
COMPUTE FORCE DUE
TO FLUID
r 
NO
P YES
SOLVE FLUID EQUTIOIBS
WITH LES wT PARTICLE
FORCES
I
END
Figure 1: DEM flowchart at a single timestep.
Newton's 2nd law for the particles is written out
as
MP d (v
dt ap (
vp + mpg Vp v Pf
N
+ [Fpw + Fpp]
and for the rotational momentum
L dw = T
where mp is the mass of the particle; Ip is the momen
tum of rotational inertia; Tp the torque of the particle;
wp is rotational velocity; vp is the translational velocity;
and 3 is the drag function as proposed by Wen and Yu
Wen and Yu (1966), where the reciprocal of the Eulerian
fluidparticle timescale is given by
3 Capafpf vf
P CD
4 dp
Svp 2.65
f
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 1: Terms of right hand side of equation (17).
Term Force Type
mpg
VpVPf
Fpw
Fp,
Drag force
Body force due to gravity
Force due to the pressure gradient
Particlewall contact force
Particleparticle contact force
and CD represents the coefficient of drag for an individ
ual particle and af represents the fluid volume fraction.
The right hand side terms of equation (17) are outlined
in table 1.
Implementation of particle collisions The particle
particle and particlewall interactions as taken into
account in this work are assumed to be local; i.e. no
longrange forces are included. For establishing the
nature of the interaction, each particle pair could be
interrogated. However, this would lead to a scaling of
the computational effort with N2, N being the number
of particles. Instead, a particlemesh algorithm is
adopted, in which each of the particles is assigned a cell
in the particle mesh based upon its location. Using this
particlemesh, each particle is only tested for interaction
against particles located in the same or directly neigh
bouring particle mesh cells. Although there is some
additional computational effort and a slightly more
complex algorithm required for this approach, it leads
to a scaling of N logN with the number of particles,
making it a lot more favourable for large numbers of
particles.
The particle collisions are modelled by the soft
sphere model as described by Cundall and Strack
Cundall and Strack (1979). In brief, this model uses
a springdashpotslider arrangement to describe the
particle behaviour before, during and after a collision.
The damping coefficient used, described by Tsuji et
al Tsuji (1994) & Tsuji et al. (1991) is, however, used
since it is related to the coefficient of restitution.
The normal and tangential contact forces are given
by the sum of forces due to the spring and dashpot.
From Hertzian contact theory the normal and tangential
contact forces are,
Fij = (kn6 ,i. G.n)n (20)
Ftij = kt6t G (21)
where G is the velocity vector of particle i relative to
particle j, Get is the slit velocity vector at the con
Table 2: Description of testcases.
TestCase Description
1 Twophase flow with 25 pm particles
at solids loading, m 1.0.
2 Twophase flow with 70 pm particles at solids
loading, m 1.0. Note that more than .I r.
of the range for the 70 pm
particles is 62 to 88pm Hadinoto et al. (2005).
3 Twophase flow with a binary mixture
(both 25 pm and 70 pm particles)
at solids loading, m 0.5 for each fraction.
tact point. The subscripts n and t represent the normal
and tangential components respectively and 6 is the dis
placement, or overlap, of particles i and j during colli
sion. Three parameters are required by the softsphere
model; stiffness (k), damping coefficient (rI) and the co
efficient of friction (p). p is a well known empirical
quantity, stiffness and the damping coefficient must be
estimated using equations
/n = a VMkn
*t = aV'Mkt64
where M = n and m is the mass of the particle.
Note that for wall collisions m oc, hence M = n.
The relationship between a and the coefficient of restitu
tion is well defined by Tsuji et al Tsuji et al. (1991). The
spring constants are, based upon elastic deformation,
k 4 1 4(1 1
3 E
 )
Ey
Sir rj
S2 aT 2a, 1 V r+rj) 2' I
kt 8( + (rir) 6n (25)
Sj / \ irj /
where r is the radius of the particles a the Poisson's ra
tio; E the Young's modulus, G the shear modulus (given
by G = 2(+E ) and 6, the magnitude of the deforma
tion in the normal direction. Note as above, when a par
ticle reaches a wall, r, oc, hence',+' 1.
Description of test cases
The objective of the simulations performed in this work
is to investigate the behaviour of a turbulent twophase
flow containing borosilicate glass beads with diameters
of 25 pm and 70 pm. Simulations with three different
test cases are pursued, which are summarised in table 2.
Table 3: Geometric data for the testcase.
Property Dimension Units
Pipe Diameter (D) 14.2 mm
Length of pipe (L) 1270.0 mm
D 90.0 ()
Length of chamber (Lchamber) 470.0 mm
Width of chamber chambere) 470.0 mm
Width of chamber chambere) 300.0 mm
Geometric data The main testcase involves a
particleladen gas jet issuing into a test chamber, as is
described in the experiment ?. The dimensions of the ge
ometry as described in the experimental work are sum
marised in table 3. Refer to figure 2 for the layout of
the experiment.
For the simulations, two geometries have been used
for the simulation cases, a pipe of length = 20 and a
pipe of length = 90. Both pipes exit in a chamber of
width and depth of = 6.5 and length of = 13. As
the chamber in which the jet exits is significantly smaller
than the chamber used in the experimental work, full slip
boundary conditions and a zero pressure gradient have
been used as boundary conditions on the chamber walls,
to minimize the influence of their presence.
Copper pipe, L/D = 90
Test Chamber
S0.300 m
0.470 m
0.470 m
Figure 2: Test chamber dimensions
Properties of the particles The physical properties of
the borosilicate glass beads as used in the simulations are
summarised in table 4.
'coefficient of restitution for inelastic particleparticle collisions
2coefficient of restitution for inelastic particlewall collisions
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 4: Summary of gas and particle properties used in
computations.
Property Dimenaion Units
Diameter of small particles (dp) 25.0 pm
Diameter of large particles (dp) 70.0 pm
Density of particles (pp) 2500.0 kgm 3
ep 1 0.94 ()
ep w 2 0.15 ()
Young's modulus (E) 100.0 MPa
Poisson's ratio (o) 0.25 ()
Table 5: Summary of gas and particle properties used in
computations.
Property Dimension Units
Temperature (Troom) 25.0 C
Density of air (pf) 1.204 kgm 3
Reynolds number (Re) 8,400.0 ()
Viscosity ( f) 1.81x10 5 kgm s
Kinematic viscosity (vf) 1.50x105 m2s 1
Properties of the fluid The fluid used in the simula
tions is air at standard conditions and its properties are
summarised in table ??.
Computational mesh
For each simulated geometry, two computational meshes
were created. All meshes are suitable for resolving the
wall boundary layer, and have 3 mesh points within the
y+ 5 layer. The coarsest geometry for the pipe with
length = 20 has 55,000 computational cells, and the
refined mesh for the same geometry has 260,000 cells.
The coarse mesh of the geometry with the pipe length of
S 90 has 200,000 mesh cells and the fine mesh has
around 600,000 mesh cells. Because resolving all the
individual particle trajectories takes around 85% of the
computational time, the refinement of the mesh does not
have a significant impact on the computational time.
Discretisation
The simulations performed in this paper use a finite vol
ume scheme, a second order backward Euler time dis
cretisation and a second order accurate central scheme
to approximate the advective terms in the continuity and
momentum equations. The solving procedure is made
parallel with the MPI libraries. During the simulations,
the CFL number is kept constant at 0.4, leading to a
slight variation in timestep over time.
Boundary and Initial Conditions
Experimental data for the singlephase flow are reported
at the pipe exit. This reported velocity profile is used to
determine the exact mass flow and the centerline veloc
ity to be specified at the inlet of the simulations. Due to
mass conservation, the mass flow rate at the pipe exit is
constant in the pipe. Using equations (26) & (27), two
simultaneous equations are set up to determine the von
Karman (K) and B constants.
/R6,, R
T= pf 27rrVl(r) dr + pf 27rrV (r) dr
(26)
U, U, ln(RU ) + B (27)
K Vf
Using the calculated values of K and B and equations
(28) and (30), the inlet velocity profile is produced.
VI(r) In((R
K
r) )+ B,
Vf
valid for 0 < Irl < (R 6,), where 6, is the linear
sublayer,
2
V (r) = (R r) (30)
vf
valid for (R 6,) < Ir < R. The frictionvelocity (u,)
is also required to determine the inlet velocity profile.
This is calculated using the following relationship,
U7 = (31)
P 1
where the shear wall stress is estimated using the friction
factor for a smooth pipe at Re = 8,400.
1 2 (3
Tw fpfV2 (32)
8
Note that the friction factor is estimated using Prandtl's
friction law for a fully developed turbulent flow Pope
(2000) for a smooth pipe. The values for the above quan
tities are given in table 6.
The matched velocity profile is set as an inflow
condition at the pipe inlet. In addition, in a number of
simulations a synthesized turbulent component is added
to this flow. The synthesized turbulence is generated by
sampling isotropic fluctuations, using the Fourier modes
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 6: Summary of gas and particle properties used in
computations.
Parameter Units
Friction velocity (u,) 0.565 ms 1
Wall shear stress (T7) 0.383 Pa
Average velocity (V) 8.80 ms 1
Linear sublayer (6,) 7.87x 105 m
Mass flow rate (T) 0.00167 kgs 1
von Karman constant (K) 0.2949 ()
B 0.79 ()
Friction factor (f) 0.033 ()
of the fully developed spectrum. The number of Fourier
modes which are considered (N) was 600, the turbulent
kinetic energy (qm) is 0.5, the boundary layer thickness
(Lt) is 0.0014 m and the timescale correlation of the
turbulence was estimated at T 2.48 x 10 4 s.
In the chamber a laminar flow of constant velocity
0.06 ms 1 through the top of the box is defined. In the
experiment of Hadinoto et al Hadinoto et al. (2005),
the top of the box is left open, allowing in order to
allow for air entrainment from the surrounding flow.
Unfortunately, the exact characteristics of this air flow
are not known, making it difficult to predict the correct
spreadingrate of the jet as it leaves the pipe.
The walls of the pipe are enforced with a noslip
condition for the velocity and pressure is extrapolated
onto the walls from the fluid. The pressure boundary
conditions were interpolated from downstream at the
inlet of the pipe and the inlet of the chamber walls. It is
found that interpolation from downstream improves the
results over the more common zerogradient boundary
conditions for pressure.
As the chamber in which the jet exits in the simulations
is significantly smaller than the chamber used in the
experimental work, full slip boundary conditions and
a zero pressure gradient have been used as boundary
conditions on the chamber walls, to minimize the
influence of their presence.
Particle Injection The particles are injected at the
pipe inlet (noted as L/D = 0) with a solids loading
very near of m = 1.0. The injection timing is set so as
to compensate for the difference in the solids loading at
each injection. This will produce an overall solids load
ing very near to 1.0 spread out through the geometry.
Assuming a uniform placement of particles at the inlet,
the location of each particle at the pipe inlet is known. A
random component is added to the spatial coordinates of
the particle. Also, the particles are given a random com
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
ponent on all velocity directions. Some form of initial
dispersion, caused by velocity fluctuations in the radial
and azimuthal direction, is crucial to obtain an accurate
representation of the flow behaviour of the particles.
70 pm particles A particle size distribution (PSD) is
used for the 70 pm particles (the particle diameter range
was between 62 pm to 88 pm). The number based distri
bution used was uniform (i.e. the number based proba
bility to find a particle is the same). Since the PSD must
be volume based, the cumulative density function (cdf)
of the number based distribution is raised to 1/3 in or
der to obtain the volume based cdf. Note that the solids
loading for the 70 pm particles is 1.0.
25 pm and 70 pm particles The procedure de
scribed for the 70 pm particles is repeated for the bi
nary mixture. A PSD is used with a cdf of 1.'. 25
pm and 5% 70 pm. This is because a solids loading
of 0.5 is required for each fraction. Therefore the ratio
of the number of particles for each fraction is given by,
n7 = 0.05.
25 (dp70
Results and Discussion
Up to now, simulations for testcases 2 and 3 have been
run, although not all testcases have sufficient statistics
to draw firm conclusions. The simulations are run on 8
cores in parallel.
Comparison of mesh refinement
Two different meshes were used to simulate the small
pipe geometry; the coarse mesh having approximately
55,000 mesh cells and the fine mesh having approx
imately 260,000 mesh cells, meaning a refinement of
about 1.75 in all directions. Figure 3 shows the effect of
mesh refinement on the predicted centerline velocity at
the pipe exit. It can be seen that the mesh refinement has
very little effect. A similar conclusion can be drawn for
the results presented in Figure 4, showing the magnitude
of the fluctuations at the pipe exit, and Figure 5 showing
the centerline velocity downstream from the pipe exit. It
can be seen that especially the fine mesh needs to be run
for longer, as the statistics are not yet steady.
Comparison of SGS LES models
Three different LES models for closing the SGS stresses
are employed. The results are shown in Figures 6,7, 8 &
9. In Figure 16 all models seem to overpredict the mass
flow rate at the pipe exit. However, a good correlation
between the experimental results is achieved especially
Figure 3: The simulation results of the mean axial ve
locity at the pipe exit (S 0) from the fine and coarse
mesh, compared with the single phase and multiphase
data.
Comparison between fine and coarse short meshes (rms axial gas velocity)
x experiment monodisperse
3.5 Vreman fine mesh
Vreman coarse mesh
S3.0 experiment singlephase
2.5
2.20
A
310
0.5
0.O0 0.1 0.2 0.3 0.4 0.5
r/D
Figure 4: The simulation results of the root mean square
of the velocity fluctuations at the pipe exit ( = 0)
from the fine and coarse mesh, compared with the single
phase and multiphase data.
the Vreman model. This is better illustrated in Figure 8
where the Vreman model is in good agreement with the
experimental results (particularly near the pipe exit).
The axial and radial rms velocities, however, are not pre
dicted well by the models. This is because there are not
enough statistical data and also this comparison is for
the short pipe.
Effect of pipe length
Fluid Statistics There is a strong effect of the length
of the pipe on the fluid velocity statistics. Figure 10
shows the mean axial flow velocity directly at the pipe
exit for the short pipe ( 20) and the long pipe
( 90) compared to experimental data. The data of
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Comparison between fine and coarse short meshes (centreline axial gas velocit
0 Vreman fine mesh
11 Vreman coarse mesh
S* experiment singlephase
10
S 2 4 6 8 10 12 14 16 18
x/D
Figure 5: The simulation results of the average center
line velocity after the pipe exit (from 0) from the
fine and coarse mesh, compared with the single phase
and multiphase data.
Axial velocity profiles of gasphase at pipe exit (shortcoarse mesh)
.
6
A
Sx x experiment monodisperse
Axa veoity profile ofga phs atpp ei sorors eh
Figure 6: The simulation results of the average center
line velocity at the pipe exit ( =' 0) as a result from the
different SGS closure models compared to experimental
data.
the long pipe is a very good match to the experimen
tal data. For both cases, it seems that the mass flow is
slightly overestimated at the inlet of the pipe. Figure 11
shows the root mean square of the velocity fluctuations
in the axial direction resulting from the simulation with
the short pipe compared, the simulation with the long
pipe and the comparison with experimental data. Both
predictions are reasonable and there is no clear differ
ence. Figure 12 shows the root mean square of the veloc
ity fluctuations in the radial direction resulting from the
simulation with the short pipe compared, the simulation
with the long pipe and the comparison with experimen
tal data. Both predictions underestimate the magnitude
of the fluctuations.
Figure 7: The simulation results of the root mean square
of the velocity velocity fluctuations at the pipe exit exit
( = 0) as a result from the different SGS closure mod
els compared to experimental data.
Figure 8: The simulation results of the average center
line velocity after the pipe exit (from = 0) as a result
from the different SGS closure models compared to ex
perimental data.
Particle Statistics There is also a strong effect of
the pipe length on the particle velocity statistics; the re
sults possibly show that the longest pipe, with length
S= 90 is not long enough to achieve fully developed
flow for the particles. Figure 13 shows the mean parti
cle velocity in the axial direction for the long pipe, the
short pipe, and the experimental data. There is a strong
improvement for the longer pipe, but the match is not
completely satisfactorily.
Figures 14 and 15 show the root mean square of the par
ticle velocity in the axial and radial directions, respec
tively. Although the trends in both figures seem to be
correctly predicted, there seems to be an offset in the
predictions compared to the experimental data.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Pipe Length comparison RMS axial gas velocity at pipe exit (Smagorinsky mod
x experiment monodisperse
short pipe
H  long pipe
1.5 experiment singlephase
0.5
V
B.0 0.1 0.2 0.3 0.4 0.5
Figure 9: The simulation results of the average axial
velocity component at the pipe exit ( = 0) as a re
sult from the different SGS closure models compared to
experimental data.
Pipe Length comparison Mean axial gas velocity at pipe exit (Smagorinsky model)
14
x experiment monodisperse
S short pipe
l ong pipe
S experiment singlephase
.0 0.1 0.2 0.3
r/D
Figure 10: The simulation results of the average fluid
axial velocity component at the pipe exit ( = 0) from
the short pipe ( 20) and the long pipe ( 90)
compared to experimental data.
Bimodal particle mixture
A case employing the Smagorinsky SGS model includ
ing van Driest Dampening was performed including a
binary particle mixture, with particles of dp 25pm
and dp 70pm with a 50:50 relative volumetric pro
portion with a total mass loading of 1.0.
Fluid velocity statistics Figures 16 and 17 show
the mean fluid velocity in the axial direction and the
root mean square of the fluid velocity fluctuations in the
axial direction, respectively. Both are compared to the
experimental findings. The results show a fairly good
prediction.
The prediction of the root mean square of the ve
locity fluctuations, compared to the experiments, as
Figure 11: The simulation results of the root mean
square of the fluid fluctuating axial velocity component
at the pipe exit ( = 0) from the short pipe ( = 20)
and the long pipe ( = 90) compared to experimental
data.
shown in figure 18 is not as good, probably due to a lack
of statistical values.
Particle statistics As the experiments have been
performed with phase doppler anemometry (PDA), the
velocity statistics of the small and the large particles in
the binary mixture have been determined separately. Al
though very similar, some changes have been found in
the statistics of the small particles compared to the large
ones experimentally.
Conclusions
This paper presents Large Eddy Simulation (LES)
modelling work and comparing this work to the ex
perimental data as presented by prof. J. Curtis and
coworkers. The simulations have been run without any
empirical parameters and material properties have been
chosen to accurately reflect the properties of the fluid
and the particles.
Generally, the simulations are too short to draw
strong conclusions. However, a few conclusions can
be drawn from the results presented so far. The first
conclusion is that the simulation results seem to be
fairly independent of mesh resolution. Changing the
mesh from 55,000 cells to 260,000 cells hardly affects
the results. Moreover, the results from the WALE and
Vreman LES models accurately capture the presence of
the wall, without the neccesity of applying van Driest
dampening. The application of such wall dampening
is computationally expensive in parallel simulation runs.
0.4 0.5
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 12: The simulation results of the root mean
square of the fluid fluctuating radial velocity component
at the pipe exit ( = 0) from the short pipe ( = 20)
and the long pipe ( 90) compared to experimental
data.
An important aspect of the simulations, is the size
of the geometry. The simulation with the shortest
geometry, with a pipe of 1 20, is too short for
the fluid and particles to be fully developed. There
are also concerns that the chamber in which the pipe
exits is too small to reflect an accurate spreading of the
jet. The uncertainty of the flow in the chamber in the
experimental work makes it difficult to match with the
simulations.
Future work
The most important goal for the future is to perform
longer simulations, so better statistics can be obtained
for the largest geometry. Once this is achieved, the ef
fect of varying the following properties in the simulation
will be performed:
The size of the chamber. The size of the cham
ber will affect the spreading rate of the jet. If this
will also affect the dispersion of particles will be
researched.
The van Driest dampening parameters at the wall.
This may have an important effect on the center line
velocity.
The magnitude of the SGS constants. All the SGS
models contain a constant, and although in theory
the magnitude of the constant should not effect the
flow characteristics this is likely not to be true for
this case, which has a fairly low Re number.
The properties of the particles. The fluid and par
ticle velocities are highly affected by the collisions
Figure 13: The simulation results of the average axial
particle velocity component at the pipe exit ( = 0)
from the short pipe (1 = 20 and the long pipe ( = 90)
compared to experimental data.
of the particles, which are affected by the material
properties.
References
P.A. Cundall and O.D.L. Strack. A discrete numerical
model for granular assemblies. Geotechnique, 29:47
65, 1979.
J. Grifoll and F Giralt. The near wall mixing length
formulation revisited. Technical report, 2000.
K. Hadinoto, E.N Jones, C. Yurteri, and J. S. Curtis.
Reynolds number dependence of gasphase turbulence
in gasparticle flows. Int J. Multiphase Flow, 31:416
434, 2005.
S.B. Kuang, K.W. Chu, A.B. Yu, Z.S. Zou, and Y.Q.
Feng. Computational investigation of horizontal slug
flow in pneumatic conveying. Ind. Eng. Chem. Res., 47:
470480,2008.
F Nicoud and F. Durcos. Subgridscale stress modelling
based on the square of the velocity gradient tensor. Flow,
Turbulence and Combustion, 62:183200, 1999.
S. B. Pope. Turbulent Flows. Cambridge Univ. Press.,
2000.
Y. Tsuji. Discrete particle simulation of gassolid flows.
KONA Powder and Particle, 11:5768, 1994.
Y. Tsuji, T. Tanaka, and T. Ishida. Lagrangian numeri
cal simulation of plug flow of cohesionless particles in a
horizontal pipe. Powder Technology, 71:239250, 1991.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Pipe Length comparison RMS Axial velocity profiles of solidphase at pipe exile
x x experiment monodisperse
25 short geometry
 long geometry
g 1.5
0.1 0.2 0.3 0.4 0.5
Figure 14: The simulation results of the root mean
square of the particle fluctuating axial velocity com
ponent at the pipe exit ( = 0) from the short pipe
( = 20 and the long pipe ( 90) compared to ex
perimental data.
H. Versteeg and W Malalasekera. An Introduction
to Computational Fluid Dynamics, the Finite Volume
Method. Prentice Hall, 2007.
A.W. Vreman. An eddyviscosity subgridscale model
for turbulent shear flow: Algebraic theory and applica
tions. Phys. Fluids, 16:36703681, 2004.
M. Weickert, G. Teike, O. Schmidt, and M. Sommerfeld.
Investigation of the les wale turbulence model within the
lattice boltzmann framework. Computers and Mathe
matics with Applications, in press, 2009.
C. Y. Wen and Y. H. Yu. Mechanics of fluidization.
Chemical Engineering Progress Symposium Series, 62:
100111, 1966.
D. C. Wilcox. Turbulence Modeling for CFD. DCW
Industries, 1998.
Pipe Length comparison RMS radial velocity profiles of solidphase at pipe exi
0.8
x x experiment monodisperse
07 H short geometry
7 long geometry
0.6
0. 0o1 02 03 04 05
Figure 15: The simulation results of the root mean
square of the particle fluctuating radial velocity com
ponent at the pipe exit ( = 0) from the short pipe
( = 20 and the long pipe ( = 90) compared to ex
perimental data.
Figure 16: The simulation results of the average fluid
axial velocity component at the pipe exit ( = 0) in
cluding the binary particle mixture compared to experi
mental results.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 17: The simulation results of the root mean
square of the fluid fluctuating axial velocity component
at the pipe exit ( = 0) for the dispersed flow contain
ing the binary particle mixture.
Radial RMS velocity profiles of gasphase at pipe exit (binary mixture)
x x experiment binary
0.5 H Smagorinsky
H experiment singlephase
Z 04
0 x x. x x
S0.2
I
Figure 18: The simulation results of the root mean
square of the fluid fluctuating radial velocity component
at the pipe exit (j = 0) of the flow containing the binary
mixture of particles compared to experimental data.
