7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Eulerian Simulations of Polydisperse Flows using a Moments Model with a
Relaxation Approach for the Moment Transport Velocities
J.N.E. Carneiro, P. Dems V. Kaufmann, W. Polifk~e
Lehrstuhl flir Thermodynamik, TU Miinchen, Boltzmannstr. 15, 85748, Garching, Germany
carneiro~td.mw.tum.de, polifke~td.mw.tum.de
Keywords: Polydisperse Flows, EulerEuler, Moments Model, presumed Functions
Abstract
Polydisperse flows consist of a continuous phase and a population of particles with different sizes and different
velocities. In order to accurately describe the interactions between both phases, it is essential to take into account
effects of polydispersity (e.g. dependency of the drag force on the particle diameter). In the present work, a novel
approach is applied, which is derived from the integration of the MultiFluid equations over the spectrum of particle
sizes, leading to a formulation based on transport equations for the moments of the size distribution function. Closure
for the momentum interaction terms and the moment transport velocities is achieved by means of presumed Number
Density Functions (pNDF) and a concept using a force balance and the particles relaxation times. This idea is based
on the work of Ferry and Balachandar (2001), where the particle velocity is expressed as a first order expansion of the
fluid velocity scaled with the particle relaxation time. An extension for inertial particles was presented by Bollweg
et al. (2007) and its adaptation for the Moments Model context will be outlined here. The models are implemented
in the open source CFD library tool OpenFOAM and simulations are carried out using two test cases. The first
one represents a 1D configuration for which an analytical solution exists, where a distribution of particles relaxes
towards the continuous phase velocity. The second one is the configuration of Baessler et al. (2007), which consists
of a water spray generated by an ultrasonic atomizer. In this case, results for mono and polydisperse simulations
are compared to each other and to the experimental data. Velocity profiles for the polydisperse model show a very
reasonable agreement with the experimental profiles at different axial positions. Different distribution reconstruc
tion strategies using Gamma functions are also tested for both cases and its advantages and shortcomings are discussed.
Nomenclatu re
r relaxation time (s)
Subscripts
c continuous phase
cl dispersed phase
Superscripts
(k) order of moment average
t turbulent quantities
Roman symbols
D particle diameter (nt)
f size distribution function (#/m )
k turbulent kinetic energy (ni s )
M moments of the size distribution (nk3)
M momentum exchange term (N kgi)
P pressure (N nt)
u velocity(nisi)
Greek symbols
Dynamic viscosity (N s ni )
v kinematic viscosity (nr s)
p density (kg ni3)
o turbulent Prandtl number ()
r shear stress tensor (N n )
Introduction
Recently, alternatives to classical MultiFluid and La
grangian polydisperse models have been developed,
which make it possible to account for particle trans
port and population dynamics in practical systems (such
as bubble columns, fluidized beds and combustion en
gines) more efficiently. These methods rely upon formu
lif not stated otherwise, the term "particle" refers throughout the text
to either bubbles, drops or solid particles
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Mathematical Model
The moments of the size distribution function f(D) are
defined as:
lations involving moments (or some other information)
of the particle size distribution to determine the disperse
phase dynamics and its interaction with the continuous
phase flow. Several examples can be found in the litera
ture, such as in the work of Wacker and Seifert (2001),
Gharaibah et al. (2002), Beck and Watkins (2003) and
Fox et al. (2008). In the present work, a novel moments
approach is presented which consists of the integration
over the size spectrum of the particle equations written
in an Eulerian frame, resulting in transport equations for
the moments and its transport velocities. Extra model
ing effort is required for closure of the interphase ex
change terms depending on unknown moments, which
is achieved with help of a presumed Number Density
Function (pNDF) of the Gamma type. This assump
tion allows an easy, algebraic reconstruction of the dis
tribution from low order moments. The representation
of the effect of particles on the continuous phase (and
viceversa) still requires the treatment of sizedependent
particle velocities (also called "polycelerity" by Boll
weg et al. 2007) which appear in the integral exchange
terms. Here, a novel closure for the sizevelocity cor
relation based on the "Equilibrium Eulerian method" of
Ferry and Balachandar (2001, 2002) is presented, where
the particle velocity is expressed as a first order expan
sion of the fluid velocity scaled with the particle relax
ation time. It is expected to be accurate if the character
istic particle response times are sufficiently small. The
adaptation of this idea for the context of the Moments
Model allows the determination of moment transport ve
locities from reference velocities, which can be obtained
by the solution of proper momentum equations. There
fore, this method provides in combination with the pre
sumed distribution function approach a possibility of
analytical closure for the Stokes drag term and the mo
ments transport velocities, specially for flows where re
laxation times are sufficiently small to allow a first order
approximation for the exact velocity. It is pointed out
that in any method that aims at transporting the moments
of the size distribution function, care must be taken in
the reconstruction process, as the outcome of the solu
tion process might generate invalid moment sets; i.e.,
an unphysical combination of moments, which has no
corresponding distribution function (Wright 2007). This
will be discussed further here, and criteria for the con
sistent check of the validity of moment sets will be pre
sented. Two configurations are chosen for the validation
of the model: the first one represents a 1D configuration
to which an analytical solution exists; the other investi
gates the droplet dispersion in a real spray generated by
an ultrasonic atomizer.
D/ nf (D)
M(A)
A physical interpretation for each of the moments fol
lows''
~\"' : total number of particles (per unit cell
volume)
Af(1) : sum of particle diameters ( II )
S.\ ' li' : total surface area of particles ( II )
(Ad) : total volume of particles ( II ); or the local
volume fraction of the disperse phase
The moment transport velocities are defined for the
kth moment as:
u(D)Df (D
which are related to the integral over size and size
velocity correlation; f(D) and u(D), respectively.
Hence, u( )'s are the mean velocities at which the mo
ments are convected2. They will only be the same (and
equal to the dispersed phase velocity) in two limiting
cases: for a monodisperse distribution; or if all particles
have the same velocity. Otherwise, there is no reason,
e.g., for u(2) and u(3) to be the same, as they are surface
area and volume average velocities, respectively.
As shown by Carneiro et al. (2009), the moments
and its transport velocities appear automatically in the
derivation of the model as a result of the integration of
the Eulerian particle equations (or "MultiFluid Model"
equations) over the diameter spectrum, giving rise to the
moments formulation (shown next). The resulting set of
equations can be solved in addition to the conservation
equations for a turbulent gas phase. The main idea is to
minimize computational costs while being able to cap
ture the polydispersity of the flow through the correct
description of the evolution of the moments.
Moments Model. The moments transport equations (in
the absence of phase change, breakup or agglomera
tion) are given by:
+ V (Af u ))
0. (3)
2Implicit here is the assumption that particles with the same size have
locally the same velocity
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
The turbulent stress tensor is evaluated using the
transport velocity u( ) and the corresponding moment
Af C), yielding the approximate expression:
v f"~fD rdD
The momentum equation for the moment transport ve
locity u( ) reads:
fD M chlD V ( fD 7 cD (4)
with the momentum exchange term
M = VP +V 7 + MDray + TD. (5)
MDray is determined by:
2
Yu()I)
3
3J,
7 A [Mi) ylv'(Vu)+ Vu k
Due to the form of the drag coefficient used, the inte
gration of the particle equations causes the appearance
of unknown moments and the corresponding moment
transport velocities (Af C ), u~l 2) in the formulation.
SA simple way to overcome this issue is to assume a func
tional form for f(D), calculating unknown moments by
reconstructing the distribution from a few loworder mo
ments. Different methods of reconstruction are intro
duced here, using a Gamma type of presumed distribu
tion function. These reconstruction algorithms will be
outlined next, followed by the relaxation approach for
the moments transport velocities.
Presumed Gamma distribution. Let's assume that
f(D) is given by a Gamma distribution:
uc u
MDr ay
with p, the continuous phase density and the Stokes re
sponse time 7 being given by
PdD (7)
The turbulence dispersion term MTD is given by an
adaptation of the model according to Oliveira and Issa
(1998) (see Carneiro et al. 2009, for decl.nllol and will be
presented in integral form later on.
The disperse phase Reynolds stress term 7Lt is mod
eled by the Boussinesq approximation
r) = v,'(Vu+Yu2VuI)2kdI, (8)
3 3
where the disperse phase turbulent viscosity vLt and tur
bulent kinetic energy kd were introduced. A description
of the algebraic models used for these terms (in conjunc
tion with the k t turbulence model for the gas phase)
can also be found in Carneiro et al. (2009).
The integral of the momentum exchange term M
leads, for the pressure and shear stress contribution, to:
f (P)cl + f (VV~r xn( Te) lD=
At1 )VP + Af C)V 7e 19)
The integration of the drag and turbulent dispersion
terms leads to:
D
D,1e 
Co 9
f (D)
with r(q) = SOr zle'
be explicitly calculated by the general expression:
r(4nf( + k 'r
r()
Af f)
with the parameters p, q being given by
Mid 2)Mid) (Mid 13))
9 = A (j 2)Af(j) (Af(j 1))2
and l'n by
,if j EN
,if = 0
./ nr n.m
18#,Af M")(u 2
Pd
uc), (10)
Hence, this algorithm provides a simple away to re
construct the Gamma distribution function from any
(11) (loworder) consecutive moments. In the present work
two alternatives will be tested: \Il " \Il' ( 0) and
Af(1)M(3) (j 1).
~/TMTDIlD
18#c, Vye ""
 _cV M )
with od = 0.7.
O'nj
pjn " 1
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
It is also possible to reconstruct the Gamma distribu
tion using only two moments; maintaining either p or q
constant, for example. If the total number and volume
of particles (per unit volume) are chosen (that is, il" I''
and M(3)) to be transported, a very simple method can
be derived for the determination of the parameters p or
q. If q is kept constant throughout the domain, p is given
by:
p= () (16)
il""' q(q + 1)(q + 2)
On the other hand, for a constant p we have:
higher order might become zero much faster than the
lower order ones. The authors applied restrictions to the
moment transport velocities in order to prevent this sit
uation from occurring, while Watkins (2005) bounded
the Gamma and Beta parameters obtained in the recon
struction process. Bounding of the distribution param
eters is not necessary if the conditions outlined above
are met and this criterion was consistently checked in
all simulations carried out herein, with replacement of
the moments being performed in order to yield positive
p and q throughout the whole domain, simply by apply
ing the convexity condition given by Eq. (19) for M( )
or 11' '. Although this is less restrictive than positive
ness of all determinants presented above, it is a suffi
cient condition to generate physical distributions by the
Gamma approach used here. In fact, once the distribu
tion is reconstructed, the HankelHadamard inequalities
must automatically be valid. Nevertheless, it is illus
trative to show the evolution of the HankelHadamard
determinants for the test cases used here, which will be
done later on. The reconstruction algorithms using only
two moments are less problematic with no need for re
placement schemes here.
Relaxation approach. The "Equilibrium Eulerian
method" was introduced by Ferry and Balachandar
(2001, 2002). For very small response times, it is ap
propriate to express the particle velocity as a first order
expansion in its relaxation time in terms of the continu
ous phase velocity. By doing so, a computationally ef
ficient method was developed, which allows to readily
obtain the particle velocity fields for each size simply
by evaluating a correction of 0(3() on uc. This idea was
used by Bollweg et al. (2007) to determine minor par
ticle size velocities from a linear interpolation between
the continuous phase and reference particle velocities.
Here, we make use of this approximation in order to de
rive relations for the determination of moment average
velocities, thus allowing the development of a closure
for u(2
The first order expansion in 7 of the particle velocity
(around the continuous phase velocity) can be written,
without the gravity contribution, as:
A2 3A 3
3A
91~ 81 ()()1 12
Conditions for a valid moment set. The most straight
forward way to check the validity of a moment set is by
calculating the convexity condition:
M(") M( ) (ME1))2 > 0. (19)
For the Gamma distributions determined by the three
moment scheme, k 2 and 3 correspond to conditions
necessary to avoid negative values of p, which is clearly
unphysical. This criterion is a necessary but not suffi
cient condition for a valid moment set. Necessary and
sufficient conditions for existence of a distribution func
tion given a set of moments (M( ), k 0, 1, 2...21 + 1)
are nonnegative HankelHadamard determinants (see
Shohat and Tamarkin 1963):
Mt )
M(" )
ak,l
M(" 1)
M(" )
M(" )
M(" 1 1)
..M(" 1)
..M("+1+1)
..M(" 21)
Due1 + O(T )
Dt LYI/
u7=u7 o 7
u,
with k 0, 1; I > 0; Ak~l 0 for a monodisperse
distribution. The spatiotemporal evolution of moments
during a simulation often leads locally to invalid sets
of moments, especially so if the M(")'s are transported
with different velocities, u(") f u(l), k f 1. If this
is the case, ad hoc replacement of the moments or re
strictions to the distribution parameters are necessary in
order to lead to physical distributions. This was done,
e.g. by Bo and Watkins (2004), which reported seri
ous difficulties in the numerical simulations of a spray.
At the spray edges, specially at the front, moments of
In the above expression, note that u7=o = uc, since
at the limit D, 7 0 the particle follows the contin
uous phase perfectly. To close the model similarly to
Ferry and Balachandar (2001), we considered the ma
terial derivative Duc/Dt following a fluid particle and
the total (Lagrangian) derivative du/di approximately
equal, which was shown to be accurate for sufficiently
small 7. Furthermore, as proposed by Bollweg et al.
uc = 0.5 m/s
>X 20 cm
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
equations, independent of k. Therefore, only two mo
mentum equations are needed to determine these veloc
ities. It is clear that a closure problem also exists in
Eq. (25), since in order to determine a given moment
Mt ), M(" ) is needed. This will be handled simi
larly to the drag closure problem, with unknown mo
ments being calculated after reconstructing the distribu
tion function with available moments. The choice of
the additional transport equations for the loworder mo
ments will depend on the reconstruction algorithm cho
sen. The volume (or mass) of the distribution function is
a very if not the most important characteristic to be
ssion conserved. Furthermore, in many situations, specially
where the sizevelocity correlation is monotonic, the
transport velocities of the moments \Il ""\Il' lie be
i us tween uc and u(3), which makes the interpolation proce
dure between these velocities most appropriate. Hence,
TC"3) and u(3) are chosen here for the reference proper
(21) ties to be used in the interpolation procedure for u(") (in
this case, the reference particle size can be thought to be
the volumemean diameter). Closure for the moments
e ap transport velocities requires, therefore, solution of Eq.
mo (4), for k = 3. As a consequence, the drag term can be
e ex rewritten with help of Eq. (22) as:
Figure 1: Schematic representation of
given by Eq. (22).
the expre~
rum,
(22)
nter
onse
time
(23)
.(7)
(24)
(2007), we consider a linear approximation of du/d
ing reference properties, such that:
Due du u7=o uo
where the index 0 refers to reference values.
We identify here that this framework can also b~
plied to determine first order approximations of the
ments transport velocities. In order to do that, th~
pression for u7 is integrated over the size spect
yielding:
1 rc)
u ~f ul D'dD= u +I (uo 
~~Mt) o, To
Therefore, this closure corresponds to a linear i
polation between uc and uo (see Fig. 1) in the resp
times space, with the averaged particle relaxation
T(c") defined as:
,( ) = M ) T Dd.
T(c") can be in turn computed by introducing Eq
into (23), giving rise to:
~(l) MC )
Mt")
/ MoragdD = Pc C)I u)
uc), (26)
to1
with K1 Q. It is important to point out that the
above formulation of the drag term offers the advantage
that it explicitly depends on the velocity which is solved
for, u( )
Test cases
Numerical model. The implementation of the model
in the C++ class library OpenFOAM was done within
a steady state version of the twoPhaseEulerFoam
solver (see Rusche 2002) with an adapted SIMPLE al
gorithm. A brief description of the implementation of
the moments equations can be found in Carneiro et al.
(2008).
1D particle relaxation. This 1D test case consists of
a population of particles that have a higher inlet ve
locity than the continuous phase and will be deceler
ated towards the gas velocity along the axial position
(see Fig. 2). Particles have density p 998 kg/m3
It follows from the moments transport equations with
Eq. (22), that:
M + V (M")u,)
V [M ((u, u( ))]. (25)
D
Hence, all moments transport equations contain a con
vective term which depends on the continuous phase ve
locity and a correction term which depends on the order
of the moment (k) through the average response times
T(c ). The relative velocity (u, uo) is common to all
Figure 2: Schematic representation of the 1D test case.
r
L.
5E05 0 0001 0005 U
D (m)
+
T
~.
1*

:  *  
*
0 4
I__ _
I
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
1 8E+14
1 6E+14
14E+14
1 2E+14
S1E+14
S8E+13
6E+13
4E+13
2E+13
Figure 3: Distribution functions at three different axial positions: x: = 0 (left), x: = 2 cm and x: = 20 cm (right).
bottom), since u(3) decelerates from 1 (at the inlet) to 0.5
m/s, the equilibrium value must be twice as high as the
inlet value. However, it is apparent that for the Gamma
distribution with j 0, M(3) is not conserved, since
it is not used explicitly in the reconstruction process.
This can be a serious drawback for this reconstruction
method, because the nonconservation of the total mass
of the distribution is unacceptable in most of the practi
cal cases. Therefore, this reconstruction method will be
dropped for the following test case.
(water droplets) and the continuous phase (air e stan
dard conditions) velocity is uc,x, 0.5 m/s. The con
tinuous phase is considered not affected by the pres
ence of the particles. The inlet discrete size distribu
tion is constant (D 15 115pm; AD = 5pm), with
the following loworder moments: \Il''' = 7.05 10'
1/m M( ) 4.65 10s 1/m il' ' 36.15 1/m
and M( ) 3.15 10 The sizevelocity correla
tion is such that the inlet moment transport velocities
are: u) 0.79 m/s, u) 0.87 m/s, u) 0.94 m/s
and u/ = 1 m/s. It is important pointing out that this
test case also serves to show that, despite of the fact that
Gammabased reconstruction methods can poorly de
scribe TopHat or monotonic shaped functions (which,
eventhough outcome of the simulation conditions, will
very unlikely occur in a real spray, for example), the
physical properties of the distribution namely the low
order moments and its transport velocities can still be
well reproduced by the model.
Figure 3 shows the distribution functions obtained
with the different moment closure methods presented
here for three axial positions (x = 0, 2 and 20 cm) and the
corresponding analytical distributions. Generally, the
correct trend is captured by the model: as the particles
are decelerated by the continuous phase, they accumu
late upstream, and a small shift towards the bigger par
ticles of the spectrum occurs. While no substantial dif
ference can be seen between the Gamma reconstruction
algorithms for j = 1 and (p,q) constant, the j = 0 dis
tribution presented smaller peaks for all axial positions
(except at the inlet). Nevertheless, similar shapes were
obtained in all cases. Figure 4 (top) shows the variation
of the third moment transport velocity (ul 3)), which de
celerates to the continuous phase velocity with increas
ing axial position, reflecting the evolution the droplet
velocities. A very good agreement is obtained between
all reconstruction methods and the analytical solution,
showing the ability of the model to reproduce the overall
momentum exchange between phases. For M( ) (Fig. 4,

Gamma (mims)
Gamma (mom2)
Gamma (p const)
Gamma (q const)
Analytical
uc
0 05 ~ 1
x (m)
5 10 0 2
Figure 4: Axial evolution of u( ) and M( ) for the 1D
case.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
10005 01 015 02
x (m)
Figure 5: Axial evolution of the HankelHadamard de
termmnants.
The analysis of the first four HankelHadamard determi
nants (ao,l, alni, ao a and Ali,, normalized with their
inlet values) reveals 1ih.II Illi' remain positive throughout
the whole domain (see Fig. 5). This guarantees the real
izability of distribution functions for all sets of moments
obtained. It is also possible to observe that the axial vari
ation of the HankelHadamard determinants is similar to
that of the moments, increasing asymptotically to a con
stant value corresponding to the equilibrium distribution
function. It suggests that the situation of particle dec
laration, leading to its accumulation with the axial po
sition, is less problematic than the contrary situation of
particle segregation, which might induce invalid sets if
the determinants vanish.
Droplet dispersion by an ultrasonic atomizer. Figure
6 sketches the geometry of the test case. An axisymmet
ric domain is used to represent the experimental setup
of Baessler et al. (2007).
main air walls (adiabatic)
~I slmm !Ssmm, 5
clarrira r ax s
120mm 280 mm
Figure 6: Sketch of the geometry used in the test case.
The atomizer is located in the central part of a tube
with 71 mm diameter and 400 mm length. Fuel is fed
to an ultrasonic nozzle at 15 ml/min, which uses a
small mass flow of carrier air to improve droplet dis
persion. Main air flow is injected at a substantially
higher flow rate (285 In/min) than the atomizing air
(15 In/min), flowing through an annular passage 71 mm
Figure 7: Contours of the continuous phase velocity
magnitude.
long, mixing with secondary air flow and spray. All flu
ids were considered at 25 C and the outlet section was
prescribed with a pressure of 1 atm. Boundary condi
tions for the moments were calculated from the experi
mentally determined distribution measured at 8 mm be
hind the oscillating plate where droplets were generated:
Af(1) = 3.83 10s m2, 3 / = 1.43 10 ml and
Af3 = 6.36 10 Contours of the continuous phase
velocity magnitude obtained by the Moments Model us
ing the threemoment Gamma scheme are shown in Fig.
7. The two different air streams can be easily identified,
with higher velocities induced by the carrier jet, while a
strong recirculating region is created by the detachment
of the main flow in the annular region. The two streams
mix inside the tube towards the outlet of the domain,
For comparison, monodisperse simulations with a
TwoFluid Model (Rusche 2002) were also carried out
for different droplet diameters (D = 15, 30, 45 and 60
pm). The axial velocity profiles at two positions along
the zcoordinate (8 mm and 150 mm behind the oscil
lating plate) are shown in Fig. 8 with the experimental
profiles obtained from PDA measurements of Baessler
et al. (2007). Two distinct streams are observed near
atm the spray inlet, which tend to mix while traveling to
wards the outlet of the pipe, where a more uniform ax
181 velocity profile is observed. It is possible to notice
that a certain mean diameter is not necessarily appro
priate for different regions of the spray. For example,
while at z =128 mm the profile seems to be better rep
resented with D 30 pm (although the peak is slightly
better captured by D =15 pm), further downstream this
is not the case anymore, with the bigger diameters show
ing better results near the axis of symmetry. This illus
trates the need for a polydisperse approach which is able
to represent the population of droplets and the exchange
processes with the gas phase in a consistent manner. The
same axial velocity profiles are investigated for the Mo
ments Model using the presumed Gamma function ap
uc (Axial)
0.500 1.88 3.25 4.62 6.00
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010

0 005
5
4
20
20 000 00 0 15 0 0 0 25 003 0 E35
r(m)
for the mono and polydisperse cases: z = 128 mm (top) and
in the region near the symmetry axis.
M^(0) Norm
0.00100 0.00473 0.0224 0.106 0.500
r,
Figure 9: Visualization of the spray through the con
tours of the normalized zeroth moment (by the inlet
value).
The parameters p and q of the reconstructed Gamma dis
tributions for the three methods are shown in Fig. 10.
The variation along the centreline is mostly smooth, with
the value of q for the threemoment scheme reaching the
lower limit very near the inlet, while p reaches a peak
at the same region. For all schemes investigated, the pa
rameters (but specially q) are apparently very sensitive
to the local flow conditions, and further reconstruction
schemes will be investigated in future work.
0 01 0 015 0 02 0 025 0 03 0 035
r(m)
Figure 8: Axial velocity profiles in the radial direction 1
270 mm (bottom).
preaches which conserve the volume of the distribution
in the reconstruction process. Since the experimental
profiles represent a number average value of the droplet
velocities, it is natural that the transport velocity of \Il "' '
(u(n)) agrees better with the data. Near the spray inlet,
smaller droplets accelerate faster due to the higher con
tinuous phase inlet velocity, having a bigger contribu
tion to u(n), while u(3) carries mostly the contribution of
bigger droplets. The width of the peak is underpredicted
by the model, with smaller velocities obtained between
the highest value and the recirculation zone. This might
be a direct consequence of the interpolation procedure
of the relaxation approach, which in this case forces the
intermediate moment transport velocities to lie between
the reference velocity (u(3)) and uc. Future work will
investigate extensions of this type of closure. Further
downstream, where equilibrium with the gas phase can
be assumed, the moments transport velocities are nearly
the same and agree very well with the experimental pro
file.
The typical structure of the spray can be identified in
Fig. 9. The radial dispersion of the droplets (mainly
due to turbulent dispersion and drag) can be observed
inside the domain, with a consequent decrease on the
normalized number of droplets per unit volume 1 .\1 "')
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
HHDET1
0.00 005 .5 005 100
SGamma(mim3)
J 0005 001 0015 002 0025 003 0(35
r (m)
~~z
z(m)
Figure 10l: Parameters p and q of the Gamma distribu
tions for the three reconstruction methods.
Contours of the determinant a,1, which corresponds to
the convexity condition for k 2, are shown in Fig. 11.
Throughout the whole domain, positiveness of the de
terminants was observed. Note that also for the main air
stream in the annular region a pseudodistribution with
very few droplets was also defined at the inlet; this dis
tribution obviously has positive determinants and this
information is also carried through the domain. Fur
thermore, an interesting feature is observed here: aol
reaches its minimum values (close to zero) in the recir
culation zone, which is "out of the reach" of droplets
(see also Fig. 12). Eventhough it was not the case here,
it highlights that whenever a transition between regions
with and without particles occurs, a potential situation
for the generation of invalid sets exists.
The reconstructed distributions at the inlet and two
positions along the centreline are shown in Fig. 13 for
the three methods. All methods can describe the exper
imental shape very well, but the peak is better captured
by the threemoment scheme. In all cases, a decrease of
the distribution function can be seen (since il ""'' dimin
ishes), while the shape of the distributions remain posi
Figure 11: Contours of the HankelHadamard determi
nant ani.
Figure 12: Profile of the HankelHadamard determinant
a,1 at z = 128 mm.
tively skewed along the axial position for each particu
lar case. For the cases where p or q were kept constant,
similar distribution functions were obtained. However,
the position of the peaks predicted by the threemoment
scheme were shifted towards smaller diameters, with a
smaller presence of bigger droplets observed.
Conclusions
A Method of Moments based on presumed functions for
the particle size distribution and a novel relaxation ap
proach for the sizevelocity correlation was presented
and tested for different configurations. First, a one
dimensional test case consisting of a population of par
ticles being decelerated in a constant gas phase veloc
ity field was investigated. Different particle relaxation
times due to Stokes drag generate heterogeneous distri
butions with the axial position. Results for the model
using different reconstruction algorithms involving pre
sumed Gamma distribution functions were analyzed and
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
z = 220 mm
D (m)
z = 340 mm
D (m)
Figure 13: Distribution functions at three different axial positions: z = 0 (left), Jr = 220 mm and Jr = 340 mm
(right).
compared to an analytical solution. A very good agree
ment was obtained for the third moment profiles and its
transport velocity, with similar results for all approaches,
except the one for j = 0, which presented a volume
conservation issue, since for this scheme M(3 is not ex
plicitly used in the reconstruction algorithm. The model
was also applied in a configuration comprising a spray
generated by an ultrasonic atomizer, where droplet dis
persion occurs due to turbulence and drag. Results show
that the basic structure of the spray is well captured
by the model and a good quantitative agreement with
the experimental velocity profiles was obtained by the
threemoment Gamma scheme. However, further types
of closure for the sizevelocity correlation still need to
be investigated, specially to include effects of bigger re
sponse times in the model. Eventhough the sizevelocity
closure can be used with any type of momentmethod, it
is essential to investigate the influence of the determina
tion of unknown moments. The Gamma parameters for
the schemes used here were found to be very sensitive
to the local conditions and other distribution forms and
reconstruction algorithms (such as Beta or QMOM) will
be tested in the future.
Acknowledgements
The authors want to acknowledge the financial sup
port given by CNPq (Brazilian National Council for
Research and Technological Development) and GRS
(Gesellschaft flir Reaktorsicherheit).
References
S. Baessler, K. Moesl, and T. Sattelmayer. NOx emis
sions of a premixed partially vaporized kerosene spray
flame. Journal of Engineering for Gas Turbines and
Power Vol. 129, pages 695702, 2007.
J. C. Beck and A. P. Watkins. On the development of a
spray model based on dropsize moments. Proc. R. Soc.
Lond. A, 459:13651394, 2003.
Y. Bo and A. P. Watkins. Mathematical development
and numerical analysis of further transport equations for
the droplet size moment theory. 19th Annual Meeting of
ILASSEurope, Nottingham, UK. September 2004, 2004.
P. Bollweg, A. Kaufmann, and W. Polifke. Deriva
tion and application of a polycelerid method for poly
dispersed twophase flows. 6th international Confer
ence on Multiphase Flow, Leipzig, Germany, 2007,
2007.
J. N. E. Carneiro, V. Kaufmann, and W. Polifke. Imple
mentation of a moments model in OpenFOAM for poly
dispersed multiphase flows. In Open Source CFD ht
ternational Conference 2008, Berlin, Germany, 45 De
cernber 2008, 2008.
J. N. E. Carneiro, V. Kaufmann, and W. Polifke. Nu
merical simulation of droplet dispersion and evapora
tion with a momentsbased CFD model. In Proceed
ings of COBEM 2009 20th international Congress of
Mechanical Engineering November, 1520, 2009 Gra
inado, RS, Brazil, 2009.
J. Ferry and S. Balachandar. A fast Eulerian method for
dispersed twophase flow. htt. J. of Multiphase Flow, 27:
11991226, 2001.
J. Ferry and S. Balachandar. Equilibrium expansion for
the Eulerian velocity of small particles. Powder Tech
i,. I. . \, 125:131139, 2002.
R. O. Fox, F. Laurent, and M. Massot. Numerical sim
ulation of spray coalescence in an Eulerian framework:
Direct quadrature method of moments and multifluid
method. J. of Computational Physics, 227:30583088,
2008.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
E. Gharaibah, M. Brandt, and W. Polifke. Numerical
model of dispersed two phase flow in aerated stirred ves
sels based on presumed shape number density functions.
In Proc. of the GennanJapanese Workshop on Midti
Phase Flow, Karlsruhe, Gennan3', 2002.
J. P. Oliveira and R. I. Issa. Numerical prediction of par
ticle dispersion in a mixing layer using an Eulerian two
phase flow model. In Proceedings of FEDSM98 ASME
Fluids Engineering Division Sainier Meeting, Washing
ton, DC., 1998.
H. Rusche. Computational Fluid D3'nainics of Dispersed
TivoPhase Flows at High Phase Fractions. PhD thesis,
Imperial College of Science Technology and Medicine,
London, UK, 2002.
J. A. Shohat and J. D. Tamarkin. The problem of ino
inents. Providence, RI: American Mathematical Society,
1963.
U. Wacker and A. Seifert. Evolution of rain water pro
files resulting from pure sedimentation: spectral vs. pa
rameterized description. Atmospheric Research, 58:19
39, 2001.
A. P. Watkins. The application of Gamma and Beta num
ber size distributions to the modelling of sprays. 20th
Annual Conference of IIASSEurope, Orleans, France.
September 2005, 2005.
D. L. Wright. Numerical advection of moments of the
particle size distribution in Eulerian models. Journal of
Aerosol Science, 38:352369, 2007.
