7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Evaluation of collision models applied to varying dense particle jets flows
W. Gumpricht, M. Chriguit, M. Braun* and A. Sadikit
ANSYS Germany GmbH, Birkenweg 14A, 64295 Darmstadt, Germany
t Institute for Energy and Powerplant Technology, TU Darmstadt, Petersenstr. 30, 64287 Darmstadt, Germany
gumprich@ekt.tudarmstadt.de and markus.braun@ansys.com
Keywords: particle collision, stochastic model, lagrange, dispersion
Abstract
Standard stochastic interparticle collision models are derived from the kinetic theory of gases. These models
differ from each other in their algorithmic approach and their modification of the collision frequency, accounting for
additional effects such as turbulence. In comparison to deterministic models, the major advantage of stochastic models
regards their lower computational cost. However, while stochastic models allow the simulation of interparticle
collisions for cases in which a deterministic model would be impracticable, the computational cost is still a critical
issue. The NanbuBabovsky model achieves a computational expense linearly proportional to the number of parcels
and conserves energy, being advantageous over other stochastic models. The model, implemented into the commercial
CFD software ANSYS FLUENT, is validated based on results of Discrete Element Method (DEM) simulations.
Two cases with different particle volume fractions are considered. The NanbuBabovsky model is also evaluated in
comparison to the O'Rourke model. These models are evaluated based on their accuracy, statistics and computational
cost.
Introduction
Stochastic collision models are straightforward ap
proaches for the simulation of interparticle collisions
in dense dispersed twophase flows. They can be eas
ily coupled to standard EulerLagrange simulations and
the physical phenomena taking place during a collision
process can be properly modeled. These phenomena
depend on the material properties of the particles. In
case of solidparticles, for instance, effects such as cen
tral collision, friction, rotation or rupture may occur,
while the collision of liquiddroplets may result in co
alescence, bouncing or breakup.
Some well known models, such as O'Rourke (1981),
Oesterle & Petitjean (1993) and Sommerfeld (1996),
are implemented in commercial CFD software and
widely used in the simulation of dispersed twophase
flows. In comparison to deterministic models, the ma
jor advantage of stochastic models is their lower com
putational cost. They allow the representation of several
real particles by a single simulated parcel, reducing the
number of particles to be computed, while the transient
timestep can be considerably larger.
However, while stochastic models allow the simula
tion of interparticle collisions for cases in which a de
terministic model would be impracticable, the compu
tational cost is still a critical issue. In order to obtain
reasonably accurate results, the required amount of com
puted parcels is, in many cases, of the order of hundreds
of thousands. O'Rourke (1981), for instance, computes
collisions for each particle pair in a computational grid
cell, resulting in a computational cost proportional to
the square of the number of parcels. Oesterl6 & Petit
jean (1993), extended by Sommerfeld (1996), achieve a
computational expense linearly proportional to the num
ber of parcels by considering a collision partner to be
virtual, with averaged properties of particles in a cell,
and computing collisions for each particle. However,
the conservation of energy cannot be guaranteed if the
collision partner is virtual.
The NanbuBabovsky stochastic interparticle colli
sion model is derived from the original Nanbu (1980)
model, which was developed to predict molecule colli
sions in rarefied gas dynamics and to approximate the
solution of the Boltzmann equation. Babovsky (1989)
has modified the original model in order to reduce its
computational cost and to guarantee the conservation of
energy during collisions. The model, which calculates
collisions for each particle in a cell, resulting in a com
putational cost linearly proportional to the number of
parcels, is implemented into the commercial CFD soft
ware ANSYS FLUENT.
The implementation of the stochastic collision model
considers the outcome of solidparticle collisions to be
a central collision with randomly determined collision
point.
In this paper, results of NanbuBabovsky, O'Rourke
and Discrete Element Method (DEM) simulations are
compared. Two cases with different particle volume
fraction of about 0.8% and 3%, in the region of higher
concentration, are considered. The aim is to validate the
stochastic models based on results of the deterministic
model and evaluate the limitations of the stochastic mod
els for denser regimes where the occurence of multiple
particle collisions is relevant. Stochastic models are lim
ited to the calculation of binary particle collisions.
DEM simulations were performed with an internal
ANSYS DEM code. They were configured to only ac
count for recoil effects on particle collisions, neglecting
other interparticle interaction mechanisms. These re
coil effects are described by spring forces resulting from
particle overlapping, based on the Hooke's Law of elas
ticity, and are equivalent to a central particle collision.
Another aim of this work is to evaluate the Nanbu
Babovsky and the O'Rourke models regarding their
computational cost and capability of meeting the ex
pected statistics for different particle and domain dis
cretizations.
Nomenclature
Roman symbols
D spring constant (N/m)
Fspring spring force (N)
k restitution coefficient
mi particle mass (kg)
np number of particles in parcel
N number of parcels in a cell
No011 expected number of collisions in a cell
Pi probability of collision
ri particle radius (m)
Ate collision timestep (s)
vi particle velocity (m/s)
vij  particles' relative velocity (m/s)
Vc cell volume (m3)
Greek symbols
vij collision frequency (1/s)
V) uniformly distributed random number
particle volume fraction in a jet
6 length of particle overlapping (m)
Subscripts
i, j particle index
n direction of collision
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Stochastic InterParticle Collision Models
The standard stochastic models are derived from the ki
netic theory of gases. They differ from each other in
their algorithmic approach and in their modification of
the collision frequency, given by
Vij (r + r)2 Vij (1)
Vc
for particles i and j. The collision frequency describes
the probability per unit time that a particle in parcel i
collides with a particle in parcel j. It assumes that par
ticles are moving in vacuum and that particle velocities
and positions are uncorrelated.
Based on equation 1, the expected number of colli
sions to occur in a cell over a timestep Ate is computed
with
1 N N
Noi = n jV At
i=1 j=1
where the 1/2 term is a result of symmetry.
While different models modify the collision fre
quency to account for the influence of the flow on inter
particle collisions, their base theory is the same and
modifications of the collision frequency can be applied
to any of these models. Sommerfeld (1996), for in
stance, takes turbulence effects into account by consider
ing the fluctuating velocity in the relative velocity vj .
O'Rourke (1981) introduced a so called collision effi
ciency to correct vij, accounting for the interaction of
the droplets with the surrounding gas flow.
Moreover, these models are limited to the occurrence
of binary particle collisions and only particles located in
the same computational grid cell are allowed to collide
with each other.
O'Rourke
The O'Rourke model calculates collisions for each
particle pair in a grid cell. It describes the probability
that a particle in parcel i collides with n particles in par
cel j with a Poisson distribution. Hence, the probability
of no collision between i and j is given by
Pi,o = exp(VijAtc)
A collision occurs when y > Pij,o, where E [0, 1]
is a uniformly distributed random number.
In case of solidparticles, however, a particle in i
may collide with at most one particle in j. Therefore,
vijAt, < 1 is required in order to achieve the expected
statistics.
NanbuBabovsky
The NanbuBabovsky model calculates collisions for
each particle i in a grid cell. The probability that i col
lides with any other particle in the cell is given by
N
P, = EPi (4)
j1
where Pij is the 1p' b.ibilii \ ih.l i collides with, defined
as
1
Pij = 2 At (5)
A collision occurs when
> Pij (6)
where y E [0, 1] is a uniformly distributed random num
ber. On average, this is equivalent to y < Pi. The colli
sion partner j is determined randomly
j = [Nj + 1 (7)
where [...J describes the integer part of the argument.
Hence, the sum in equation 4 is not solved and the com
putational cost of the model is linearly proportional to
the number of parcels.
There is a timestep limitation associated with the
NanbuBabovsky model. This limitation results from the
condition that NPj < 1. The fulfillment of this condi
tion can be guaranteed by limiting the collision timestep
2
Atc <  (8)
VijN
Hence, the maximum allowed collision timestep
At.,max is determined by estimating the maximum col
lision frequency for all particles in a cell
2
A x max(vij)N
Collision Outcome
The collision outcome is a central collision with ran
domly determined collision point. The collision point is
the point where the two colliding particles touch each
other. It results from the intersection of the straight line
connecting the centroid of both particles and their sur
face. The frame of reference of particle i is consid
ered and the collision direction is described by the ba
sis vector n, which points in the direction of the straight
line connecting the centroid of i and the collision point.
Hence, derived from the laws of momentum and energy
conservation, the postcollisional velocity components
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
of the particles in n direction is computed with follow
ing equations:
mnvi,n + ". kmj(vi, j,) (10)
mi + mj
new ivi,n + ..* .. kmi(vj,n vi,n)
n (11)
mi + mnj
where m, and my describe the mass of each particle.
Note that only the velocities in the direction of n are
changed. k is the restitution coefficient. k 1 in ideal
elastic collisions, while in ideal plastic collisions, k = 0.
Solution Strategy
Stochastic interparticle collision models are executed
after particle trajectories are calculated with the La
grangian solver. Since only particles located in the same
computational grid cell are allowed to collide with each
other, the collision routine loops over all cells and calcu
lates collisions for the particles located in each one.
The occurrence of collisions is determined through
a BernoulliExperiment. The NanbuBabovsky imple
mentation estimates the collision timestep before cal
culating collisions for the particles in a cell. If the es
timated timestep is smaller than the Lagrangian time
step, the BernoulliExperiment is repeated more than
once, until the collision time is equal to the Lagrangian
timestep. The O'Rourke implementation considers the
Lagrangian timestep to be the collision timestep.
In case of O'Rourke, the routine loops over all par
ticle pairs in the cell and tests each pair for collision.
The NanbuBabovsky routine loops over all particles in
a cell, testing each particle for collision. Once a colli
sion occurs, the direction of collision is randomly deter
mined and the postcollisional velocities are calculated
with equations 10 and 11.
Discrete Element Method
DEM is a deterministic model and particle collisions
occur only if particles overlap each other. Simulations
were configured to only account for recoil effects on par
ticle collisions, neglecting other interparticle interac
tion mechanisms. These recoil effects are described by
spring forces resulting from particle overlapping, given
by
Fspring D6
based on the Hooke's Law of elasticity, where D is the
spring constant and 6 is the length of overlapping. The
spring force acts in the direction of the straight line con
necting the centroid of both particles.
This requires a sufficiently small computational time
step At to accurately predict the outcome of collisions,
SI
Figure 1: Problem geometry.
depending on the given spring constant. A larger At
results in a larger particle overlapping and, depending
on D, the recoil force may be strongly overpredicted.
Hence, even though momentum is conserved in all coor
dinate directions, the method does not conserve energy
and the error of the postcollisional kinetic energy de
pends on the computational timestep.
For the simulations presented in this work, a spring
constant D = le+05 N/m and a At = 2e08 s were con
sidered. Investigation has shown that, with this configu
ration, the expected error is about 1.5% and 4.5% for
the cases with lower and higher particle density, respec
tivelly. This error is comparable to experimental mea
surement devices and is sufficient to assess the quality
of the stochastic models.
Configuration
The problem geometry is illustrated in figure 1. Two
particle jets with 0.015 m of diameter and 45 degrees
of inclination relative to the boundary are injected into a
threedimensional domain. As both jets cross each other,
particle collisions occur causing the dispersion of parti
cles over the computational domain.
The injected particles have 5e04 m of diameter and
15 m/s of initial velocity in the direction of the jet axis.
The problem is configured in such a way that parti
cle motion is affected only by interparticle collisions.
Hence, particles are considered to be moving in vac
uum and drag forces as well as gravitation are set to 0.
Interparticle collisions are considered to be ideal elastic
(k 1), while other material properties are irrelevant.
Moreover, the boundary conditions are set to escape and
particles leave the domain after reaching its borders.
Overall, two cases with particle volume fraction of
about 0. and 1.5% in each jet were computed, while
in the region near the crossing point of both jets the par
ticle volume fraction may be up to twice as much. The
particle volume fraction in each jet is described with
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
.In order to have enough parcels for a statistical
application of the model, the case with '. 0. .
was computed with one particle per parcel (n, 1).
The other case, with '. 1.5%, was computed with
np 5. In addition, all stochastic simulations were
computed ti.iiiMnicill with a computational timestep At
= le04 s, and a grid with 16 x 16 x 40 cells was used.
Results and Discussion
Particle Volume Fraction of 0.4% in each Jet
Since the NanbuBabovsky and O'Rourke models
share the same collision frequency and the amount of
simulated parcels is sufficient for a statistical treatment
of the problem, both models delivered practically the
same results.
Qualitatively, the stochastic models predict the same
as DEM. One can see in figure 2 (bottom) that, as both
jets cross each other, some particles collide with each
other as are dispersed over the computational domain,
while other particles that do not collide, or collide only
slightly, keep a path similar to their original one. This is
the same process predicted by DEM, illustrated in figure
2 (top).
In figure 3, this process can be recognized in the di
mensionless particle concentration profiles over the y
coordinate, in sampling planes A, B and C. The dimen
sionless particle concentration is defined as the ratio of
the number of particles sampled in each measurement
interval to the mean number of particles. The mean is
the total number of sampled particles in each plane di
vided by the number of intervals. The profiles are plotted
separately in order to better analyze the behavior of each
A BC
I .
l
I I I DEM
injection1
injection0 ..
NanbuBabovsky
Zx
Figure 2: Particle dispersion colored by particle veloc
ity magnitude for the '. 0. case.
006
007
006
005
0.04
03
0.03
DEM
*NanbuBabovsky
**O'Rourke
a,.
k ~
''
a 01 Plane A
0 1 2 3 4 5 6 7 8 9 10
o as
008
0 07 k
006 ",
005 ,,
00 as
A.
0.04
003
0.02
Plane B 001
1 2 3 4 5 6 7 8 9 10
007 
006 
005 
0.04
003
0 02
oo01 Plane C
01 234 56789
N/Nm
Figure 3: Dimensionless
0.06
0.05
0.04
003 
0.02
0 1 2 3 4 5 6 7 8 9 10
0.08 
0.07 
0.06 .,
0.05
06 1 2 3 4 5 6 7 8 9 10
0.08
0.06
0.05
0.04 
003
0.02
001 '
1 23 4 5678
N/Nm
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
particles is considered. Again, the profiles are qualita
tively comparable, while a deviation between the results
of stochastic models and DEM is present.
The effects of interparticle collisions are also de
scribed by the change in the original velocities of the
particles. In comparison to DEM, this velocity change
is accurately predicted by the stochastic models. Figure
5 (left) shows the particle axial velocity (i.e. xvelocity)
profiles in each plane. These profiles reflect the averaged
velocities of the particles sampled in each measurement
interval. The average axial velocity is smaller in plane
A than in planes B and C, since the concentration of par
ticles near the jets crossing point is higher and several
particles with a lower axial velocity and a higher trans
verse velocity are located in this region. As particles
move away from the crossing point, the particles that re
main near the domain axis are the ones with a higher ax
ial velocity and lower transverse velocity. Consequently,
the average particle axial velocity near the center axis is
higher in planes B and C. The averaged velocity near the
9 10
particle concentration of
injection0 (left) and injection1 (right), for
the '. z 0. case.
jet, where injection0 describes the jet on the bottom and
injection1 the one on the top (see figure 2). Since many
particles do not collide, the concentration of particles
is higher within the particle jets. This is described by
the peak of each profile. These peaks are closer to the
boundaries in planes farther from the jet crossing point,
describing the movement of the jets toward the domain
borders, while the symmetry of both jets can be seen.
Even though the stochastic models qualitatively predicts
the particle concentration with good agreement, a small
deviation between the results of stochastic models and
DEM can be seen. After the crossing point, the stochas
tic models predict a higher particle concentration within
each jet. This reflects an underpredicted number of col
lisions, which may be a result of the models' grid depen
dency.
Figure 4 (left) shows the particle volume fraction pro
files. The volume fraction describes the dispersion of
particles more accurately, since the amount of sampled
* DEM
*NanbuBabovsky
*O'Rourke
Plane A
0204 06 0.8 1 1.2 14 16
.>
/
Plane B 
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
Plane C
0.2 0.4 0.6 0.8 1 1.2 1.4
DPM volume fraction [%]
0 008
00 0.5 1.5 2 2.5 3
0 08
0.07,
0.06 i
005
004
003 '
0.02
4
0.01.
0 0.5 1 1.5 2 2.5 3
0 D7
ooa6,
007
A
004
003
002
001*1
DO 0 5 1 15 2 25 3
DPM volume fraction [%]
Figure 4: Particle volume fraction for the '
(left) and 1 .'. (right) cases.
* 
008 . i i
S. DEM
007 .. NanbuBabovsky
O'Rourke
0.06
005 
I0.04 ,
0.03 
0.02
o.o01 / Plane A
0 10 15 20 25
008 i i '
0.07
0.06
0.05 
o0.04 7
0.03
002 l
0.01 PIno R
004
>. /
0.03 
002 f
o01 Plane C
0 5 l10 15 20
Particle x velocity [m/s]
0.02
0.01 *
20 15 10 5 0 5 10 15 20
008 1
007 
006 7
0.05 
004 "
0.03
0.02
0.01
20 15 10 5 0 5 10 15 20
008
t
007
oo006
005 /'
004
0.03
002 f
001
20 15 oI 5 0 5 110 15 20
Particle y velocity [m/s]
Figure 5: Axial (left) and transverse (right) velocity
profiles for the '. 0. case.
domain borders is also higher in planes B and C, in com
parison to plane A. This is explained by a higher concen
tration of particles in these regions. In plane A, few par
ticles are sampled near the borders, since only a few par
ticles collide in such a way that their ratio of transverse
to axial velocity is large enough for them to approach
the borders before crossing plane A. The velocity oscil
lations in this region is explained by the small number
of sampled particles. Also, the location of the particle
jets, which contains particles that did not collide, is rec
ognizable in the axial velocity profiles in planes B and
C. In these locations, the velocity profile is more or less
flat, describing the initial particle axial velocity.
The particle transverse velocity (i.e. yvelocity) pro
files are shown in figure 5 (right). In the domain center
axis (y 0.4), the transverse velocities are averaged to
0 due to symmetry. Closer to the borders, the flat profiles
describe the location of the particle jets, with unchanged
initial velocities. In planes A and B, the velocity oscil
lates near the domain boundaries due to the low amount
of sampled particles.
:Ir
i
injection1
iniectiono 
ZA
NatrliuBao:E. k,
Figure 6: Particle dispersion colored by particle veloc
ity magnitude for the '. 1 ".'. case.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Particle Volume Fraction of 1.5% in each Jet
In a denser regime, most particles collide and tend
to accumulate near the domain axis before being dis
persed. Qualitatively, this effect is accurately predicted
by the stochastic models, as illustrated in figure 6 for
the '. 1.5% case. In this case, the average dis
tance between particle centroids may be as low as 2.6
particle diameters in the jets crossing region. Hence, the
occurence of multiple particle collisions is considerable.
Figure 7 shows the dimensionless particle concentra
tion profiles of each particle jet. Even though the jet
symmetry is recognizable, all peaks are located near the
center axis, describing the accumulation of particles in
this region. DEM predicts a higher relative concentra
tion of particles near the center axis. However, it does
not mean that the predicted particle volume fraction is
higher, since the number of sampled particles and the
Nmean is smaller in DEM results. This is a result of a
higher particle dispersion in lateral (i.e. zcoordinate)
direction. According to figure 4 (right), DEM predicts
a lower particle volume fraction in this region. Hence,
the stochastic models predict a higher dispersion in y
coordinate direction, while DEM predicts a higher dis
persion in zcoordinate direction. This may be explained
by the high concentration of particles near the crossing
point and the occurrence of multiple particle collisions,
which is not taken into account by stochastic models. As
the particle flow coming from the top and bottom of the
domain is intense, more particles tend to disperse to the
sides. Qualitatively, however, this effect is accurately
predicted by the stochastic models, as seen in figure 8.
The axial and transverse particle velocity profiles in
each sampling plane is illustrated in figure 9. The
A BC
I 4
I A
i O E . .M'
,"1 :
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.08I'I'I
DEM
007 *NanbuBabovsky
S *.O'Rourke
0 06
005 ^
0 04
>, _''
0.03
002
0oo Plane A
0 1 2 3 4 5 6 7 8 9 10
0.08i
0074
006
005 "'
004 ">
0.03
002 T
001o Plane B
S1 2 3 4 5 6 7 8 9 10
008 
007
0.06 ,
0.05 
1004 '*
003 .'
002
0.01' Plane C
0 1 2 3 4 5 6 7 8 9 10
N/N
0.08 
0.07
0.06
0.05 3
0.02
001
0 1 2 3 4 5 6 7 8 9 10
0 1 2 3 4 5 6 7 8 9 10
0.081
0.06 
0.05 ,
0.04
003
0.02 '
0.01
0 1 2 3 4 5 6 7 8 9 10
N/N
Figure 7: Dimensionless particle concentration of
injection0 (left) and injection1 (right), for
the '. w 1.5% case.
stochastic models predict the postcollisional particle
velocities with good accuracy. Here, the tendency is
similar to results for the lower particle volume fraction
S 0. I'. However, since most particles collide and
the particle dispersion over the domain is higher, the pro
files are smoother and the presence of the particle jets is
no longer recognizable. Near the domain boundaries, the
amount of particles sampled in plane A is low, resulting
in an oscillation of averaged velocities.
Statistics Investigation
A quick investigation was conducted to verify if the
implemented NanbuBabovsky model and the O'Rourke
model are able to meet the expected statistics for differ
ent particle discretizations and grid resolutions. In order
to ensure enough particles for the application of a statis
tical treatment even for a coarser particle discretization,
a denser regime with a particle volume fraction close to
6% in each injected particle jet was considered.
I 3M601
2 S5101
2 71 01
2 53&t01
242&01
DPE.l
Figure 8: Cross section view of particle dispersion be
hind plane A for the '. 1 ".'. case.
The expected number of collisions to occur in the
computational domain, over a timestep At le04 s,
is the sum of the expected values of each grid cell, com
puted with equation 2. Figure 10 compares the number
of computed collisions for different particle discretiza
tions, ranging from 1 to 80 particles per parcel, and the
averaged expected number of collision events over time,
computed based on the case with 1 particle per parcel.
A computational grid with 16 x 16 x 40 cells was used.
The number of computed collisions with the Nanbu
Babovsky model oscillates about the expected value for
all cases, while the amplitude of oscillations increases
as particle discretization gets coarser. O'Rourke fails to
meet the expected statistics if the particle discretization
is coarse enough, underpredicting the number of colli
sions. This is due to the large vij At term in equation 3,
resulted from a large np.
As described in equation 1, the collision frequency of
stochastic models depends on the volume of the cell in
which the particles are located. Therefore, the resolu
tion of the computational grid may have a substantial
influence on the accuracy of the collision model if par
ticles are nonuniformly distributed over the computa
tional domain. This occurs due to wasted cell volume
by a coarser grid, resulting in a smaller calculated col
lision frequency. Figure 11 illustrates the dependency
of the stochastic models on the grid resolution. Simula
tions with three different computational grids were per
formed, where parcels representing one single particle
are injected. For each simulation, both the expected and
computed number of collisions were tracked. The num
ber of computed collisions decreases for coarser grids,
due to the effect of wasted volume. However, since
this dependency is related to the collision frequency, the
number of expected collisions also decreases for coarser
grids. The NanbuBabovsky model meets the expected
statistics for all grid resolutions, while O'Rourke consid
erably underpredicts interparticle collisions if the com
putational grid is fine enough. Again, this is a results of
a large ,j At, term in equation 3.
. '* *.
.,'*,*.., 
', ,* )i
131
I 3901
* rlOr X
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0 08 i i i .
SM *DEM
007, NanbuBabovsky
*. *'O'Rourke
006
005
0.04
003 .
002 
ool .* Plane A
0 5 10 1 5 20 25
008
007 
0.06
0.05
. 004
003
002 V
0.01 Plane B
00 5 10 15 20 25
0 o08
0.07 .
006 
005
0 04
0 03
0.02
o01 Plane C
0.
0 5 10 15 20 25
Particle xvelocity [m/s]
008
007
006
005
0.04
0.03
"'{
002
o01 4
020 15 10 5 0 5 10 15 20
0.08
0.07
0.06 
0.05
0.04
003 3
0.02
0.01 *
20 15 10 5 0 5 10 15 20
0.08 ,
0.07
0.07 
0.06 
0.05
0.04
0 03
0.02
001
20 15 10 5 0 5 10 15 20
Particle yvelocity [m/s]
Figure 9: Axial (left) and transverse (right) velocity
profiles for the '. z 1.5% case.
Computational Cost
The Discrete Element Method is computationally
very expensive. The serial simulation of the. ;
1.5% case, with about 13,000 particles in the domain,
for instance, required more than a week of computa
tional time in a 2.3 GHz AMD Opteron CPU, while the
stochastic models delivered results in just a few minutes.
Although the computational cost of the O'Rourke
model is substantially lower in comparison to DEM, the
model is still expensive for cases with a large amount
of simulated parcels. In order to compare the compu
tational expense of the NanbuBabovsky and O'Rourke
approaches, simulations injecting different number of
parcels into the domain were performed with each
model. A denser regime with a particle volume fraction
close to i*'. in each particle jet was considered. First,
1,000 timesteps of 2e05 s were computed to ensure a
more or less stationary state. Then, 200 additional time
steps were computed while the elapsed computing time
was measured.
012000

S10000
E
a 8000
0 6000
o 4000
 2000
E
z
O'Rourke
I ~i
0 0.01 0.02 0.03
Time [s]
Figure 10: Statistics investigation for different particle
discretizations.
Figure 12 shows the elapsed computing time of simu
lations with different number of parcels in the computa
tional domain. If the amount of parcels is small enough
(i.e. < 29,000), other processes have a higher weight
on the computing time than the collision model, and
the cost increase is not even linearly proportional to N.
Between 29,000 and 290,000 parcels, the cost increase
of NanbuBabovsky simulations is more or less linearly
proportional to N, being a bit overproportional after
92,000 parcels. This may be a result of the timestep
condition 8. The O'Rourke approach does not present
a cost increase proportional to the square of the number
of parcels. This may be explained by the model's ap
proach of computing collisions for particle pairs in each
cell. Since particles are not necessarily proportionally
distributed over the cells after refining the particle dis
cretization, the cost increase does not have to be pro
portional to the increase in the overall number of parti
cles. However, the computational cost of the O'Rourke
approach increases more rapidly than the one of the
NanbuBabovsky model. If 290,000 parcels are com
puted, for instance, the cost of the O'Rourke model is
about 167% higher.
This investigation reflects the behavior of a particular
case, however, the timestep condition 8 also depends
on the collision frequency vij. In cases with a high dis
crepancy in particle sizes and nonuniform particle dis
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
.10000
S8000
0
6000
4000
'5
S2000
E
Figure 11: Statistics investigation for different domain
discretizations.
cretization, as well as a large number of real particles,
max(vij) is high, requiring a considerably small At,.
However, it is important to mention that, if these cases
involve solidparticles, a simulation with the O'Rourke
approach also has to account for a certain timestep con
dition in order to limit the term uij Atc, whether consid
ering a Poisson distribution or not, since only one colli
sion is allowed between particles in each colliding par
cel. Otherwise, the number of collision events is under
predicted, as seen in figures 10 and 11.
Conclusions
In this work, the NanbuBabovsky stochastic inter
particle collision model, implemented into the commer
cial CFD software ANSYS FLUENT, is introduced. The
model is validated based on results of Discrete Element
Method (DEM) simulations. An internal DEM code was
used.
Even though stochastic models are advantageous be
cause of their lower computational cost in comparison
to deterministic models, the computational cost is still a
critical issue. In many cases, the amount of computed
parcels is of the order of hundreds of thousands, requir
ing a lot of hardware resources to calculate interparticle
collisions.
10000
SNanbuBabovsky
**O'Rourke
1000
a,
E
WO 100
E
NanbuBabovsky
 ., A(V
S 8x8x20 cells
... 8x8x20 cells (expected)
16x16x40 cells
16x16x40 cells (expected)
32x32x80 cells
 32x32x80 cells (expected)
0.01 0.02 0.03 0.04 0.
le+06
Figure 12: Elapsed computing time of Nanbu
Babovsky and O'Rourke simulations.
The main advantage of the NanbuBabovsky model
in comparison to other stochastic models is that it is
able to compromise the conservation of energy with a
lower computational cost. While some models consider
a particle's collision partner to be virtual in order to re
duce their computational expense, the NanbuBabovsky
model considers a real collision partner, conserving en
ergy during collision processes. For certain computa
tional timesteps, the model has a cost linearly propor
tional to the number of parcels and, although this lin
earity does not hold for all timesteps, the model's ex
pense is considerably lower in comparison to the ap
proach of the O'Rourke model, which computes colli
sions for each particle pair is a cell.
Results of both NanbuBabovsky and O'Rourke sim
ulations are compared to DEM results. If the amount
of simulated parcels is high enough for applying a sta
tistical treatment of the problem, and the collision time
step is small enough, both models deliver practically the
same results, since they are based on the same theory. In
comparison to DEM, both stochastic models are able to
predict interparticle collisions with reasonable accuracy
and low computational cost. A DEM simulation that re
quired more than a week to be concluded was performed
with both stochastic models in just a few minutes. The
accuracy of the models is investigated for particle vol
ume fractions up to :'. in regions of higher particle con
centration.
The collision statistics and, consequently, the results
of stochastic models depend on the computational grid
resolution. The NanbuBabovsky model is robust re
garding its capability of meeting the expected statis
tics. The number of computed collisions always oscil
lates about the expected values, independent on the par
ticle discretization and grid resolution. If the number
of simulated particles is high enough, the expected col
lision frequency is met with good accuracy. This is a
result of the collision timestep estimation in the im
le+04 le+05
Number of parcels
~ I T T r
plemented model. The NanbuBabovsky model con
siders the collision timestep apart from the Lagrangian
timestep, such that collisions may be computed more
than once per Lagrangian timestep. This is just a rep
etition of the BernoulliExperiment and particle trajec
tories are not calculated between consecutive collision
timesteps. The O'Rourke model, however, considers
the Lagrangian timestep and, in some cases, the num
ber of calculated collisions is underpredicted due to a
large expected number of collisions between two parti
cles over a single timestep.
However, both stochastic models have their limita
tions. In contrast to deterministic models, the stochastic
models do not take particle positions into account, being
appropriate to the simulation of free flow regimes, where
accumulation of particles does not occur. These mod
els are unable to control the concentration of particles in
cases where particles tend to accumulate in certain re
gions, such that the bulk density of the dispersed phase
may exceed the maximum allowed in some grid cells.
Also, stochastic models are limited to the occurrence
of binary particle collisions, being suitable to cases in
which the average distance between particles is much
larger than their diameters.
In order to accurately meet the expected statistics over
each computational timestep, the amount of simulated
parcels needs to be sufficiently large. Even though the
NanbuBabovsky model is robust and the number of
computed collision events always oscillates about the
expected values, the amplitude of these oscillations in
creases as the amount of computed parcels decreases,
possibly resulting in a higher simulation error.
Moreover, the implemented NanbuBabovsky model
assumes that particles are moving in a vacuum, not con
sidering the influence of the continuous phase on inter
particle collisions and not taking turbulence into ac
count, being more suitable to flow regimes with St > 1.
In order to achieve accurate results for cases in which
the flow is influential to the occurrence of interparticle
collisions, the model has to be extended. This can be
easily achieved by modifying the collision frequency.
O'Rourke, for instance, defined a collision efficiency
used to modify the collision frequency, accounting for
the interaction of the droplets with the surrounding gas
flow. In highly turbulent cases, with St < 1, an ap
proach similar to the one of Sommerfeld, considering
the fluctuating velocities, could be used.
Acknowledgements
The authors appreciate the technical support of Ralf
Kroeger in configuring the Discrete Element Method
simulations.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
References
Babovsky, H., Convergence proof for Nanbu's Boltz
mann simulation scheme, Europ. J. of Mech. B/Fluids
1 (1989), 4145.
Bicanic, N., Discrete Element Methods, In: Encyclo
pedia of Computational Mechanics, (ed. Stein, E., de
Borst, R. and Hughes, T. J. R.), Vol. 1., Wiley, 2004.
Crowe, C. T., On the relative importance of particle
particle collisions in gasparticle flows, Proc. of the
Conf. on Gas Borne Particles, Paper C78/81 (1981),
135137.
Cundall, P A. and Strack, O. D. L., A distinct element
model for granular assemblies. G6otechnique 29 (1979),
4765.
Elghobashi, S., On predicting particleladen turbulent
flows, Proc. 7h Workshop on TwoPhase Flow Predic
tions (Ed. M. Sommerfeld) (1994).
Gidaspow, D., Multiphase flow and fluidization: contin
uum and kinetic theory descriptions, Academic Press,
New York, 1994.
Gross, D., Hauger, W., Schroder, J. and Wall, W. A.,
Technische Mechanik 3, Springer Verlag, Berlin, 2008.
Kanther, W., GasFeststoffStromungen in komplexen
Geometrien, Dissertation, Universitat Dortmund, 2003.
Nanbu, K., Direct simulation scheme derived from the
Boltzmann equation. I. Monocomponent gases, J. Phys.
Soc. Japan 49 (1980), Nr. 5, 20422049.
Oesterl6, B. and Petitjean, A., Simulation of particleto
particle interactions in gassolid flows, Int. J. of Multi
phase Flow 19 (1993), Nr. 1, 199211.
O'Rourke, P. J., Collective drop effects on vaporizing
liquid sprays, Dissertation, Princeton University, 1981.
Schiller, L. and Naumann, A., Uber die grundlegende
Berechnung bei der Schwerkraftaufbereitung, Ver. Deut.
Ing. 44 (1933), 318320.
Schnell, W., Gross, D. and Hauger, W., Technische
Mechanik 2, Springer Verlag, Berlin, 2002.
Sommerfeld, M., Theoretical and Experimental Mod
eling of Particulate Flows, Lecture Series 200006,
MartinLutherUniversitat, 2000.
Sommerfeld, M., Modellierung und numerische Berech
nung von partiklebeladenen turbulenten Strimungen mit
Hilfe des Euler/LagrangeVerfahrens, Shaker Verlag,
Aachen, 1996.
