7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Numerical Simulations of Solid Particle Motion in Rarefied Flow Using VOF
H. Strjm*, S. Sasict and B. Andersson*
Competence Centre for Catalysis / Department of Chemical Engineering,
Chalmers University of Technology, SE412 96 Giiteborg, Sweden
i Department of Applied Mechanics, Chalmers University of Technology, SE412 96 Giiteborg, Sweden
henrik.strom~chalmers.se, srdjan~chalmers.se and bengt.andersson~chalmers.se
Keywords: VOF, solid particles, rarefied gas
Abstract
A novel model is proposed for simulations of gassolid systems involving particles of size \ignlilk~.lll to that of the
bounding geometry. The size of the particles is comparable to the mean free path of the gas. The model is based on
the Volume of Fluid (VOF) method, but extended to allow for the incorporation of rarefaction effects. The procedure
suggested is used to investigate the change in the drag on a spherical particle translating along the axis of symmetry in
a cylindrical pore compared to the Stokes drag on an identical particle. It is found that the new model predicts a lower
increase of the normalized drag with increasing the blockage ratio than the slip boundary condition approach used by
Keh & Chang (2007). Our result is in agreement with the Monte Carlo simulation data of Stefanov, Barber, Ota &
Emerson (2005).
Nomenclatu re
Roman symbols
a particle radius (m)
b pore radius (m)
C, Cunningham correction factor ()
g gravitational acceleration (ms )
F drag force (N)
F, surface tension force (N)
L length (m)
n normal ()
P pressure (Pa)
r radial position (m)
t time (s)
u velocity (ms )
U relative velocity (ms )
V volume (m )
Greek symbols
volume fraction ()
/9 slip coefficient (kgm s )
/c viscosity ratio ()
Curvature ()
A mean free path (m)
Viscosity (Pas)
p density (kgm3)
o surface tension (Nm )
Dimensionless numbers
Kn, particle Knudsen number
Ma, particle Mach number
Re, particle Reynolds number
Subscripts
0 Stokes
g Gas
i Coordinate direction parallel to surface
j Coordinate direction normal to surface
p Particle
rRarefied
sSound
W Poiseuille
Introduction
Background. In many important environmental appli
cations, the interaction of small particles with neighbor
ing particles and with the flow field in a confined envi
ronment is of great interest. One such example is the fil
tration of soot particles from the exhaust gases of an in
ternal combustion engine, usually a process taking place
inside a porous wall. Figure 1 illustrates the concept:
The filter consists of a large number of paired channels
separated by a porous wall. The channels are plugged
in either end, so that the particleladen exhaust gas flow
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 1: The regimes of rarefied flow.
Flow regime Range of Knudsen number
Continuum flow Kn < 0.015
Slip flow 0.015 < Kn < 0.15
Transition flow 0.15 < Kn < 4.5
Free molecular flow 4.5 < Kn < 00
Results are shown for the welldefined case of a parti
cle translating along the centerline of a straight cylin
drical pore. These results are then compared to the
approach where the flow field is resolved using a slip
boundary condition on the particle surface.
The intention is to use the new model to derive sub
models for simulations of intricate industrial applica
tions, e.g. simulations of trapping of diesel particulate
matter in the porous walls of a diesel particulate filter.
Such submodels could then be used in simulations with
less computationally expensive methods.
Rarefied Flow. The degree of rarefaction effects on the
momentum transfer to a particle in rarefied gas flow can
be estimated from the Knudsen number, Kn p, defined
as''
Kn, = (1)
with A representing the mean free path of the gas (ap
proximately 67 nm in air at ambient conditions).
The regime where Kn, < 0.015 is known as contin
uum flow, and here the usual equations of continuum
fluid dynamics still hold. For higher values of Kny,
the particle will no longer experience gas as a contin
uum but will start to notice the effects of collisions with
individual gas molecules. Macroscopically, this effect
is characterized by the occurrence of a slip velocity be
tween the gas and the particle. For this reason, various
attempts have been made to resolve the flow field around
the particle using continuum fluid dynamics with a slip
boundary condition on the surface of the sphere (e.g.
Barber & Emerson (2000), Keh & Chang (2007)). How
ever, it seems the slip boundary condition approach per
forms unsatisfactorily for curved surfaces, particularly
as the Knudsen number increases (Einzel, Panzer & Liu
(1990), Bailey, Barber & Emerson (2004), Bailey, Bar
ber, Emerson, Lockerby & Reese (2005)). The different
regimes of rarefied flow are given in Table 1 (see e.g.
Sommerfeld, van Wachem & Oliemans (2008)).
It is of interest to clarify the regimes of rarefied flow
that are typically encountered in the applications for
which the model proposed here is intended. Taking the
motion of soot particles in the porous walls of a diesel
particulate filter at high temperature as an example, the
typical panticle Knudsen numbers range from 0.15 and
....................J !
Figure 1: Conceptual illustration of a wallflow diesel
particulate filter: two channels separated by a
porous wall.
follows approximately the streamlines indicated by the
dashdotted line. A closeup shows a particle moving in
side a generic pore. Depending on the panticle and lo
cal pore fluid and flow properties, the panticle will either
be able to pass through the wall and exit with the ex
haust gas, or it will be trapped inside the porous wall (or
on/among the previously trapped particles). The trap
ping of the particles is the fundamental first step in the
process by which particulate matter is removed from ex
haust gases. The next step is a recurrent temperature in
crease by which the oxidation of the accumulated mass
of trapped particles into less harmful substances is initi
ated.
A complication in any fluid dynamic simulation of
such a system is the fact that the flow is dilute in the
channels but confined inside the pores. In addition, the
gas is rarefied with regard to the typical spatial scales
of the particles and potentially also to the scales of the
pores. Numerical simulations of such systems is a chal
lenging problem, since a Lagrangian pointparticle ap
proach is not applicable when the particles approach the
size of the bounding geometry, and, at the same time,
molecular methods (e.g. Monte Carlo simulations) are
extremely expensive at small degrees of rarefaction. In
addition, the panticle size distribution spans several or
ders of magnitude, and the motion of larger particles will
be affected by inertial mechanisms.
In the current work, a novel model is introduced with
which more accurate simulations are made possible at a
moderate computational cost. It is intended for gassolid
systems where the gas flow is rarefied (i.e. the mean
free path of the gas is comparable to the length scale of
the bounding geometry and/or that of the particles) and
where the particles are \ignilk~.llll in size compared to
the bounding geometry.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where density and viscosity are calculated using:
P = QsPs + QpP,
p = appy9 + Q~py,. (5
Note here that in the expression for p (equation 4), we
use p, where pI, is found in a standard VOF model. This
will be discussed later.
The surface tension force F, is determined using the
model of Brackbill, Kothe & Zemach (1992):
F, = Iaca
5(Ps + Pp))
where /ce is the curvature, defined as:
V ti (7)
and the surface normal n is the gradient of az:
n = Va,,. (8)
The surface tension between the gas and the particle
is thus o.
The single set of equations, used for solving the flow
of both phases, (equations 2 and 3), is complemented by
tracking the volume fraction of each phase throughout
the domain. The interface between the particle and the
gas is thus obtained by solving the following continuity
equation for the particulate phase:
+ V (q,u) = 0, (9)
and the volume fraction of the gas phase is then given
as = 1 q,.(10)
As previously mentioned, the physical property p,
in equation 4 is exchanged for a "rarefied" density, p,,
through which the effects of rarefaction on particle mo
tion are taken into account. This rarefied density is cho
sen in such a way that the particle equation of motion
for a point particle is identical for the case of a particle
in an unbounded, creeping, rarefied flow and the case of
a particle in an unbounded, creeping, nonrarefied flow.
In other words, the effects of rarefication are lumped into
the rarefied density. This procedure ensures that the two
described particles will have the same acceleration.
The rarefied density can be readily calculated using a
suitable correction function, such as the one originally
proposed by Davies (1945). In the current work, we use
the parameters determined from Millikan's reevaluated
data as given by Allen & Raabe (1982).
ilu V(puu)
Dr
VP + V [p(Vu + VuT]
+PI'pg +F, (3)
Pr
'For a more thorough discussion on the different regimes of rarefied
flow, see e.g. Cercignani (2000).
Table 2: Approximate spatial scales in filtration of
diesel particulate matter.
Variable Value Units
Particle diameters 5 1, 000 1111
Pore diameters 2 20 pi
Mean free path 100 200 1111
up to well above unity, while common Knudsen numbers
for the bounding geometry (where
characteristic length L in equation 1) are up to 0.015.
This means that the particles are typically in the slip
flow or transition flow regimes, whereas the flow in the
bounding geometry is in the continuum or slip regime .
Table 2 lists the corresponding spatial scales for com
parison. It is interesting to note that, since the mean free
path of the gas increases with temperature, rarefaction
effects come into play at smaller geometrical scales if
the temperature is high.
Numerical method
Equations. In the current work, it is suggested that the
gassolid system can be resolved on the pore scale using
the Volume of Fluid (VOF) framework. In VOF, the two
(or more) phases present are treated as immiscible fluids.
The computational domain is discretized into computa
tional cells, each containing a fraction a; of each phase
i. The volume fractions of all the phases always sum up
to one in each computational cell (c.f. equation 10).
The herein proposed model is based on the VOF
framework, but extended as discussed in the following.
For more background on VOF, the reader is referred to
the recent reviews by Benson (2002) and Gopala & van
Wachem (2008).
Consider a solid particle of radius
ative velocity LT to that of its rarefied carrier gas. The
flow is isothermal, creeping (Rel, = 2pgaU/p, C 1)
and subsonic (May, U/U, < 0.3). In the VOF frame
work, where the particle is treated as a fluid parti
cle, the velocity field u(.r, y, z, t) and the pressure field
P(.r, y, z, t) can be determined from shared continuity
and momentum equations:
Vu = 0 (2)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 3: Summary of numerical schemes used in the
simulations.
Equation Scheme
Pressure PRESTO!
Momentum QUICK
Volume fraction CICSAM
Pressurevelocity coupling PISO
It should be noted here that for unbounded flow, Mil
likan's data is known to be in very good agreement with
the analytical formula of Beresney, Chemnyak & Fomya
gin (1990), which outperforms the slip boundary condi
tion approaches (Bailey, Barber, Emerson, Lockerby &
Reese (2005)).
The gravitational acceleration for a particle of den
sity p, should be pV,g, so the rarefied density must
be introduced also into equation 3. The effects out
side the particle are negligible, since for aerosol systems
p, > p,. An alternative approach would have been
to introduce gravitation as a source term that acts only
where up 1 In the simulations presented in the cur
rent work, gravity has been neglected altogether and this
is therefore not an issue.
For the considered fluid particle to behave like a solid
particle, the former must not distort and should have a
negligible tendency to internal circulation. This is en
forced as the surface tension, a, and the viscosity ratio,
/c p p/y,, goes to infinity. A thorough preliminary
investigation showed that particle deformation is negli
gible and that the decrease in drag due to internal circu
lation is limited to 33 ppm for the currently used values
of /c and a.
Numerical Schemes. The proposed model is imple
mented into the commercial CFD software ANSYS Flu
ent 12.1 using UserDefined Functions (UDFs). The
utilized discretization schemes are listed in Table 3.
The Compressive Interface Capturing Scheme for Ar
bitrary Meshes (CICSAM) by Ubbink (1997) is chosen
for its ability to handle high viscosity ratios and pre
serve shapes of interfaces (Waclawczyk & Koronowicz
(2008)). The unsteady solver uses a firstorder implicit
formulation, whereas for the time advancement of equa
tion 9, an explicit formulation is used. For the latter, the
Courant number is not allowed to exceed 0.25.
Since the CICSAM scheme is first order accurate, it
is necessary to use a very fine mesh resolution around
the gasparticle interface. It is however not necessary to
use the same fine resolution in the whole of the com
putational domain, and dynamic local grid refinement is
therefore a suitable means to reduce the computational
demand of the simulations. Such a scheme has been de
veloped but is outside the scope of the current work and
was not used for the simulations presented herein.
Slip particle in cylindrical pore
Simulation Setup. The case of a slip particle translating
along the centerline of a cylindrical pore, as examined
by e.g. Keh & Chang (2007), is illustrated in Figure 2.
The fully developed Poiseuille slip flow profile is given
assuming that the tangential momentum accommodation
coefficient of the walls equals unity (i.e. only diffusive
Figure 2: Slip particle of radius a nm translating along
the centerline of a cylindrical pore of radius b.
The Poiseuille slip flow profile aw at the pore
inlet has a mean flow velocity of 0.1 m/s. In
the current work, a 325 nm, L 9.6b and
b is varied between 406.25 nm and 1.625 pm.
reflection) as (Barber & Emerson (2000)):
aw (r) = 2Uw 1
I + 4
The pore walls have noslip boundary conditions for
fluid flow, whereas the particle experiences rarefaction
effects relative to the value of its slip coefficient, B,
which is related to the Cunningham correction via
OC = .(12)
Ba + 2p,
Thus, as B 00o, the conventional noslip boundary
condition is recovered, and as B 0 free slip over the
particle surface is obtained. In this work, two degrees of
particle rarefaction effects are considered:
a = 1(131
Ba
0.1.
Table 4: Normalized drag on a slip particle traveling
along the centerline of a cylindrical pore.
0.2 0.1 0.121)
0.4 0.1 0.261
0.6 0.1 0.6!)()
0.8 0.1 3.515
0.2 10.f)08
0.4 1 1.839)
0.6 14.1)18
0.8 124.727
Table 5: Summary of gas and particle properties used in
computations.
Variable Value Units
a 325 Im
py, 1000 kg/m13
Ps 1.181) kg/m13
p,1.81 10" kg/(n s)
In terms of Knudsen numbers, this corresponds to
Kni, 1 and 10, respectively.
The drag force on the slip particle in the pore, F, is
normalized with the Stokes drag force on an identical no
slip particle in unbounded flow, Fo:
Fo 6i7rypU. (15)
In order to obtain this data from Keh & Chang (2007),
the values in their Table 2 are multiplied by the factor
(,30 + 2p9) / (,30 + 3p,9).2
Keh & Chang (2007) use a Maxwellian firstorder slip
boundary condition on the surface of the particle:
p, cle
i (16)
8, ~ 3 cl~
where u; is the velocity component parallel to the sur
face of the particle and zr is the coordinate direction nor
mal to the surface of the particle. In the current model,
the slip on the particle is taken into account via the rar
efied density (c.f. equations 3 and 4). The simulations
are performed for ambient conditions and the material
properties used are listed in Table 5.
Results. The prediction of the newly proposed model is
compared to the results of Keh & Chang (2007) in Figure
3. Our results are also given in Table 4 for clarity.
2This somewhat cumbersome procedure is performed since normal
ization s .11. 11. 1. .11 .;~~1~~~. ...ll.. II...~ pl. .I ..I.I..ls.I .I.... com
pared to normalizing with the drag force on a slip particle in un
bounded flow, especially when the latter is calculated using an ap
proach known to be imprecise compared to experimental data.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
..a**~ pahyl = 0.1 Keh & Chang (2007)
25palrl = 1 Keh & Chang (2007)
***#** pahyl = 0.1 New model
20 0 8 par1= 1 New model
15
u..O
0o 0.2 0.4 0.6 0.8 1
alb
Figure 3: Normalized drag on a slip particle traveling
along the centerline of a cylindrical pore.
Discussion. For the case of a slip particle moving in a
pore, there are two effects on the drag experienced by
the particle. The first effect is a decrease in the mo
mentum transfer to the particle related to the degree of
rarefaction, i.e. the same effect that would be expected
for a slip particle in unbounded flow. The second ef
fect is an increase in drag due to the presence of the
walls. In both our results and the results of Keh & Chang
(2007), the wall effect outweights the rarefaction effect
at a/b > 0.4 for the lower degree of slip, so that the
drag force on the particle is larger than the correspond
ing Stokes drag. For a/b 0.2, the reduction of drag
due to rarefaction effects however is larger, so that the
normalized drag force becomes less than unity. Our re
sults fall below the curves of Keh & Chang (2007) for
all values of a/b, but the discrepancy between the two
methods grows as a/b 4 1. For the higher degree of
rarefaction, the slip boundary condition fails already at
low blockage ratios. Our results predict that the wall ef
fect will dominate only for the highest blockage ratio,
alb = g
In a comparison between a firstorder slip boundary
condition approach, such as the one used by Keh &
Chang (2007), and a more accurate Direct Simulation
Monte Carlo (DSMC) technique, Stefanov, Barber, Ota
& Emerson (2005) found that above some critical Knud
sen number, which is dependent on the blockage ratio
a/b, the Monte Carlo technique predicts lower normal
ized drag. The Knudsen numbers in the current simu
lations are well above the critical Knudsen numbers re
ported by Stefanov, Barber, Ota & Emerson (2005) for
the blockage ratios investigated. It is thus to be expected
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Acknowledgements
This work has been performed within the Competence
Centre for Catalysis, which is financially supported by
Chalmers University of Technology, the Swedish En
ergy Agency and the member companies: AB Volvo,
Volvo Car Corporation, Scania CV AB, GM Powertrain
Sweden AB, Haldor Topsoe A/S and The Swedish Space
Corporation.
References
Allen, M. D. & Raabe, O. G., Reevaluation of Mil
likan's oil drop data for the motion of small particles
in air, J. Aerosol Sci., Vol. 13, pp. 537547, 1982
Bailey, C. L., Barber, R. W. & Emerson, D. R., Is it
safe to use NavierStokes for gas microflows?, Euro
pean Congress on Computational Methods in Applied
Sciences and Engineering, ECCOMAS, JyviskylR, July
2004
Bailey, C. L., Barber, R. W., Emerson, D. R., Lockerby,
A. & Reese, J.M., A critical review of the drag force on
a sphere in the transition flow regime, Rarefied Gas Dy
namics: 24th International Symposium, American Insti
tute of Physics Conference Proceedings 762, Monopoli,
Bari, Italy, July 2005
Barber, R. W. & Emerson, D. R., A numerical study of
low Reynolds number slip flow in the hydrodynamic de
velopment region of circular and parallel plate ducts,
Technical Report DLTR00002, Daresbury Labora
tory, 2000
Benson, D. J., Volume of fluid interface reconstruction
methods for multimaterial problems, Appl. Mech. Rev.,
Vol. 55, pp. 151165, 2002
Beresney, S. A., Chernyak, V. G. & Fomyagin, G. A.,
Motion of a spherical particle in a rarefied gas. Part 2.
Drag and thermal polarization, J. Fluid Mech., Vol. 219,
pp. 405421, 1990
Brackbill, J. U., Kothe, D. B. & Zemach, C., A contin
uum method for modeling surface tension, J. Comput.
Phys., Vol. 100, pp. 335354, 1992
Cercignani, C., Rarefied gas dynamics: From basic con
cepts to actual calculations, Cambridge University Press,
2000
Davies, C. N., Definitive equations for the fluid resis
tance of spheres, Proc. Phys. Soc., Vol. 57, pp. 259270,
1945
that the results of Keh & Chang (2007) overpredict the
drag on the slip particle in the cylindrical pore. The cur
rent method predicts a less steep increase in drag with an
increasing blockage ratio, in line with the DSMC data.
It should be noted here that the current comparison
becomes unphysical for high values of a/b (since rar
efaction effects on the flow between the particle and the
pore wall will increase as a/b approaches unity and the
noslip boundary condition on the wall becomes inap
propriate). This anomaly may be addressed using a slip
boundary condition on the pore wall, which should be
an appropriate approach at least as long as the pore is
straight. In the DSMC simulations referred to above, the
pore walls exhibited slip, and the results are indeed sim
ilar to ours for a/b 0.2. As was mentioned in the
introduction, for our particular area of interest of appli
cation, it is to be expected that the particle are experi
encing large rarefaction effects whereas the flow in the
bounding geometry is still in the continuum regime.
Conclusions
A novel model is proposed for simulations of gassolid
systems where the particle sizes are \ignlilk~.lll to that
of the bounding geometry. It is applied to the case of a
particle traveling in a rarefied, creeping carrier gas along
the symmetry axis of a cylindrical pore.
When compared to available literature data using a
slip boundary condition on the surface of the sphere
(Keh & Chang (2007)), our model predicts a less pro
nounced increase in drag on the particle as the ratio of
the particle size to the pore size goes to one. This is
in line with previous results reported for the same sys
tem using Direct Simulation Monte Carlo simulations
by Stefanov, Barber, Ota & Emerson (2005). Since it
has also been shown that the slip boundary condition
approach severely overestimates the momentum transfer
to a particle in unbounded rarefied flow (Bailey, Barber,
Emerson, Lockerby & Reese (2005)), it is to be expected
that the increase in drag with increasing blockage ratio
should be less pronounced in a more accurate approach.
In summary, we demonstrate here that the new model
does a more accurate prediction of the momentum trans
fer to the particle but with a significantly smaller compu
tational cost compared to a Monte Carlo approach. The
new model is also easily implemented into a commercial
CFD software.
The proposed model will, after further validation, be
used in forthcoming attempts to derive new submodels
for soot particle behavior inside the porous walls of a
diesel particulate filter.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Einzel, D., Panzer, P. & Liu, M., Boundary condition for
fluid flow: curved or rough surfaces, Phys. Rev. Letters,
Vol. 64, pp. 22692272, 1990
Gopala, V. R. & van Wachem, B. G. M., Volume of fluid
methods for immisciblefluid and freesurface flows,
Chem. Eng. J., Vol. 141, pp. 204221, 2008
Keh, H. J. & Chang, Y. C., Creeping motion of a slip
spherical particle in a circular cylindrical pore, Int. J.
Multiphase Flow, Vol. 33, pp. 726741, 2007
Lockerby, D. A., Reese, J. M., Emerson, D. R. & Bar
ber, R. W., Velocity boundary condition at solid walls in
rarefied gas calculations, Phys. Rev. E, Vol. 70, 017303,
2004
Sommerfeld, M., van Wachem, B. & Oliemans, R., Best
practice guidelines for computational fluid dynamics
of dispersed multiphase flows, ERCOFTAC/SIAMUF,
Giiteborg, 2008
Stefanov, S. K., Barber, R. W., Ota, M. & Emerson,
D. R., Comparison between NavierStokes and DSMC
calculations for low Reynolds number slip flow past a
confined microsphere, Rarefied Gas Dynamics: 24th In
ternational Symposium, American Institute of Physics
Conference Proceedings 762, Monopoli, Bari, Italy, July
2005
Ubbink, O., Numerical prediction of two fluid systems
with sharp interfaces, PhD thesis, Imperial College of
Science, Technology & Medicine, London, 1997
Waclawczyk, T. & Koronowicz, T., Comparison of CI
CSAM and HRIC highresolution schemes for interface
capturing, J. Theor. Appl. Mech., Vol. 46, pp. 325345,
2008
