7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Droplet Evaporation and Influence on Turbulence in a Swirl Combustor
K. Luo*, H. Pitscht and 0. Desjardinst
State Key Laboratory of Clean Energy Utilization, Zhejiang University, Hangzhou, Zhejiang 310027, P. R. China
t Center for Turbulence Research, Stanford University, Stanford, CA 94305, USA
Department of Mechanical Engineering, University of Colorado at Boulder, CO 80309, USA
zjulk@zju.edu.cn, h.pitsch@stanford.edu and Olivier.Desjardins@colorado.edu
Keywords: Direct numerical simulation, droplet evaporation, spray combustion, swirl combustor
Abstract
Direct numerical simulations of nheptane spray droplet evaporation and combustion as well as their interactions with
turbulence in a 3D model swirl combustor are performed to provide insights into realistic spray evaporation and
combustion. The variabledensity, lowMach number NavierStokes equations are solved using a fully conservative
finite difference scheme in cylindrical coordinates, with second order space and time discretization. Droplet
evaporation is described by an equilibrium model with corrections. Vapor combustion is modeled using an adaptive
onestep irreversible reaction. Interphase twoway coupling of mass, momentum, and energy are applied based on the
particleincell approach. It is found that to correctly predict droplet evaporation, the underlying grid size has to be at
least ten times the droplet diameter. The presence of droplets and combustion decreases the velocity spectrum in the
higher wavenumber region and enhances the central recirculation zones, but increases the turbulent velocity spectrum
in the lower wavenumber region and reduces the outer recirculation zones. As expected, the effect of preferential
concentration is stronger when the droplet Stokes number is closer to unity.
Introduction
Turbulent multiphase combustion is encountered in a
number of engineering applications such as internal
combustion engines and gasturbine aircraft engines.
The ability to perform accurate numerical predictions of
these systems is of paramount importance to improve
their design and efficiency. However, the underlying
physics of spray combustion are extremely complex.
The liquid phase undergoes primary and secondary at
omization. The resulting droplets are subject to evapo
ration, condensation, further breakup, or collision and
coalescence with other droplets. The resulting vapor
in the gasphase undergoes turbulent mixing supplying
the flame with unburned evaporated fuel. These con
current processes of liquid phase dynamics and evapo
ration, turbulence, as well as combustion interact and
strongly affect each other, which makes experimental
measurements, highfidelity simulations, and modeling
very challenging.
To understand these multiphysics phenomena and
their strong coupling interactions, direct numerical sim
ulations (DNS) are very helpful. In recent years, DNS
of spray combustion in simple configurations have been
performed by several groups (Domingo et al. 2005;
Reveillon and Vervisch 2005; Nakamura et al. 2005).
However, most studies are limited to 2D or simple con
figurations. In the present work, nheptane spray droplet
evaporation and combustion as well as their interactions
with turbulence in a 3D model swirl combustor are in
vestigated by means of direct numerical simulation to
provide insight into realistic spray evaporation and com
bustion. It should be noted that while the turbulent
physics are treated directly by solving the Navier Stokes
equations, the flow inside the liquid droplets and the flow
in the boundary layers around the droplets is unresolved.
These assumptions will be discussed in more detail be
low. Also a onestep chemical mechanism is assumed
for the combustion of nheptane. This is an acceptable
assumption, since for the questions addressed here, only
the heat release and its interaction with mixing and the
turbulence is important, and the details of the chemical
structure of the flame is of minor importance.
Mathematical models
In the present study, the NavierStokes equations for
the gas phase are solved in an Eulerian framework,
whereas the governing equations for the dispersed phase
are formulated in the Lagrangian sense. The droplets are
treated as point sources of mass, momentum, and energy
to the gas phase. This coupling approach is valid for a
mixture with droplet diameter smaller than the smallest
scales of the gas phase.
Governing Equations for the Gas Phase. The gas
phase is described using the variabledensity, lowMach
number NavierStokes equations. The influence of spray
droplets appears in terms of a mass source term Sm in
the continuity equation and a momentum source term S,
in the momentum equations. These equations read
9p Opui
at dxi
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Governing Equations for the Liquid Phase. Dis
persed droplets are known to follow the Basset
BoussinesqOseen equations. However, due to high den
sity ratio (PL/P > 1) between liquid nheptane droplets
and air, the Basset history term, the added mass effect,
and other unsteady drag effects can be neglected, since
all these terms are much smaller than the Stokes drag
force (Crowe et al. 1998). Thus, it is assumed here that
the droplet momentum equations can be described by the
dominant drag force and gravity only. Collisions and co
alescence are also neglected, as the volumetric loading
of droplets is small. Then, the equations for the droplet
displacement (Xd) and velocity (Ud) can be written as
(Maxey and Riley 1983)
dxd,i _
(1) and
dud,i fi
S (i Ud,i) + gi
Opui .' .. ..
at Oxj
Op i dj
where
aij = /( Oaz + O ]
2 au e
32i OJk
3' ) k
(3)
Three additional scalar transport equations correspond
ing to the mass fractions of fuel (F), oxidizer (0) and
products (P) are solved based on a onestep combustion
model that will be described later. They include a chem
ical source term ci, i E {F, O, P} as well as the same
mass source term as that in the continuity equation from
spray evaporation, 8.
In Eq. (9), gi is the gravitational force, and Td is the par
ticle time constant, defined as
pLd2
Td (10)
18P
where d is the droplet diameter, and PL is the liquid den
sity. fl is an empirical correction to the Stokes drag law
for larger droplet Reynolds numbers, which is here taken
to be (Kurose et al. 2003)
1 + 0.0545Red,slip + 0 1 I;,' ,lip (1 0.03Red,slip)
f1 b
1 + aReb,b
where
a = 0.09 + 0.077exp (0.4Red,slip)
a (pDF )+F+Sm (4)
Oxj Oxj
j pDo aO)z +o (5)
a (pD YPj + (6)
9xj O xj
In addition, a transport equation for the gas temperature
is solved
OpT CpTuj
P tCp 
8 8yj
O + WT + ST (7)
axyj axjJ
where ST is the energy exchange rate with the droplets,
and DT is the rate of heat release by combustion.
b = 0.4 + 0.77exp (0.04Red,slip) (13)
The local droplet Reynolds number based on the slip ve
locity between the droplet and the gas, Red,slip, is de
fined by
p u(xd) Ud d
Red,slip = (14)
and Red,b = is the local droplet Reynolds number
based on the blowing velocity ub due to evaporation.
Droplet Evaporation Model. For liquid droplet evapo
ration, there exist a variety of models (Miller et al. 1998;
Sirignano 1999; Sazhin 2006). In terms of complexity,
these models can be divided into six groups, i.e. constant
droplet temperature model, infinite thermal conductiv
ity model, finite thermal conductivity model, effective
conductivity model, vortex model and model based on
dpY, dpYYUj
at Oxj
+
apYo0 +pYoUj
at Oxj
nd
OpYp dpYpUj
9t Oxj
full solution of the NS equation. Recently, Sazhin
et al. (2007) have suggested an approach to numerical
modeling of droplet heating and evaporation with radi
ation based on analytical solutions of the heat conduc
tion equation of a single droplet. In the present study,
considering both the complexity and the accuracy, we
use the infinite thermal conductivity model with some
corrections, and assume that there is no temperature
gradient inside droplets. Heating and evaporation of a
singlecomponent droplet can be expressed in terms of
the masstransfer number defined by
Tysurf YF
BM = u Y (15)
1 _ypurf
where YF is the ambient mass fraction of fuel at the
droplet position, and ypurf is the mass fraction of fuel
at the surface of the droplet, obtained by assuming lo
cal equilibrium between the droplet and the ambient gas.
From the ClausiusClapeyron saturation law, the mole
fraction of fuel at the droplet surface can be expressed
as
Xsur Psat Patm [LvWp ( 1 1
tF P P R eb T )]
(16)
In this expression, P is the local pressure at the droplet
position, Patm is the atmospheric pressure, Tb is the liq
uid boiling temperature, R is the universal gas constant,
WF is the molecular weight of the fuel, and Lv is the la
tent heat of vaporization. The mass fraction of fuel at the
droplet surface yurf can be obtained from Xurt using
Xsurf
ysurf ___ F (17)
f Xysurtf +(1 surt)W/W (17)
Fp +(1 FA )W1WF
where W is the molecular weight of the mixture. Based
on the above assumptions, the following equations for
droplet temperature Td and mass md are obtained
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
This correction was shown to be very important (Miller
et al. 1998). To account for the effect on evaporation
by convection, the Nusselt and Sherwood numbers are
usually corrected. The popular correction A
Nu 2 + 0.552Red,slip Pr (22)
Sh 2 + 0.552Red,slip2 Sc (23)
and correction B
0.552Red,slip 2 Pr3
Nu = 2 + (24)
(1+ 1.232 2'
Red,slipPr3
0.552Red,s ip Sc
Sh = 2 + (25)
(1+ 1.232 4)2
Red,slipSSc
are tested in the present simulations.
Another important issue in droplet evaporation mod
eling is the choice of the physical properties used in the
above equations. It has been shown that evaporation rate
predictions are sensitive to the choice of property values
of gas and vapor (Law and Law 1976). The general ap
proach is to define a reference temperature and a vapor
mass fraction that are used to evaluate both the gas and
vapor properties
TR Td + a (T9 Td),YR Yur + (Y Yr)
(26)
where a is a constant coefficient. Given the vapor and
gas properties evaluated at the reference temperature, the
physical properties for mixture can be calculated using
the reference mass fraction weighted averaging proce
dure. For example, the heat capacity of mixture sur
rounding the droplet can be expressed as
Cp,m = YCpy,v + (1 YR)cp, .
F2 rndLv
(T Td) + dLv
Td mdCL
Sh md
Sh In( 1 + B),
3Sc rd
where Pr and Sc are turbulent Prandtl number and
Schmidt number, CL is the heat capacity of liquid, Cp,m
is the heat capacity of ambient mixture, T, is the gas
temperature at the droplet position, and f2 is a correc
tion for effect of droplet evaporation on droplet heating,
defined as
c 1'
3Prr1 dTd
2md
However, the above property calculation is generally
needed at every time step and can bring Nigiliik.ii com
putational expense in spray combustion cases, where a
large number of droplets is usually involved. To avoid
this, we follow Miller et al. (1998) and evaluate physi
cal properties based on the estimated wet bulb tempera
ture, which has been demonstrated to be efficiently and
sufficiently accurate. In this approach, we only need to
evaluate physical properties at the beginning of the sim
ulations. The liquid droplet density PL, heat capacity
CL and latent heat of vaporization Lv are assumed to be
constant and computed based on the estimated wet bulb
temperature, too.
InterPhase Coupling. Full mass, momentum and en
ergy coupling between gas and liquid droplets are in
cluded. The source terms Sm, Si, and ST due to spray
evaporation are written by assuming that each droplet n
dTd Nu c,,m
dt 3Pr CL
dmd
dt
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
in a cell of volume V acts as a point source to the gas
phase, leading to
V dt
Sm (,"Zd (28)
n
1 d
_, (,,) (29)
n
and
1 d en$dn) d
T ( T(" LT) htF ( ) (30)
n
where h is the enthalpy of formation of the fuel vapor,
related to Lv by
Lv (T)= h (L cp)T (31)
if we define the reference enthalpies for the liquid and
all species except the fuel to be null at T=0.
Vapor Combustion Model. For gaseous vapor com
bustion, detailed chemistry is generally preferable. Nev
ertheless, it is very computationally expensive. In the
present study, we focus on the global characteristics
of spray combustion. Thus the combustion process is
chosen to be represented by a onestep global reaction
model of nheptane. The choice is motivated by the
low spatial resolution requirements of such a model in
comparison to multiplestep models. However, in order
to retain as much as possible of the physical processes
found in realistic spray combustion, the approach pro
posed by FernandezTarrazo et al. (2006) will be fol
lowed, in which the coefficients of the Arrhenius law are
fitted in order to accurately reproduce the burning veloc
ity for the problem of interest in both the lean and rich
environments. This methodology has been shown also
to be able to reproduce the rate of strain at extinction.
The chemical source terms Cp, wo, wp, and CT in
Eqs. 4, 5, 6, and 7 are obtained from the following one
step irreversible reaction of nheptane:
vFF + voO vpP + Qt
Here Qt is not a constant any more, but is specified as
function of the equivalence ratio (FerandezTarrazo et
al. 2006). The preexponential factor A is set to be 9.7 x
108 m3/(mols), while the activation temperature TA is
another parameter specified as function of equivalence
ratio to accurately reproduce the burning velocity and
the strain rate at extinction.
Numerical Algorithms
The lowMach number NavierStokes equations are dis
cretized using a secondorder finite difference scheme in
space, and a secondorder semiimplicit scheme in time
in a staggered spacetime grid (Desjardins et al. 2008).
This scheme is an extension of the one from Morinishi
et al. (2004), and has excellent conservation properties.
Mass and momentum can be exactly conserved even on
nonequidistant grids. Kinetic energy is conserved ex
actly for constant density flows, and to the order of the
scheme for variable density cases. Governing equations
of droplets are advanced first, followed by the scalar and
momentum equations. The velocities are then corrected
by solving a Poisson equation to satisfy continuity. The
equations for species mass fraction and for the temper
ature are rendered especially stiff because of the chem
ical source terms. As a result, a fully implicit treatment
of the chemical source terms is required and has been
implemented for an accurate integration of those terms.
The droplets are described using a Lagrangian solver
that uses a secondorder, fully coupled RungeKutta
time integration of the droplet equations. The informa
tion from the gas phase is interpolated at the droplet po
sitions using a trilinear interpolation. The transfer of
the source terms back to the gas phase from the spray is
based on the commonly used PIC (particleincell) ap
proach (Crowe et al. 1977). Adaptive time stepping is
required to allow for an accurate integration of the equa
tions when the droplets become very small without sig
nificantly increasing the computational cost.
Single Droplet Evaporation
where the stoichiometric coefficients vi are taken
from the global nheptane oxidation reaction C7H16 +
1102 7CO2 + .TT.0O, and Qt is the heat release.
Then the chemical source terms can be written as:
o F
Ujp
CT
VFWF
voWo
VpWP
Qt
where
S(pYA )CpYo)x (
WF, Wo
TA
T
To study evaporationcombustion interactions, it is nec
essary to make sure that the droplet evaporation and va
por combustion models can correctly describe the corre
sponding physical processes. Although some validations
and comparative analysis of droplet evaporation mod
els have been presented in previous studies (Miller et
al. 1998; Sazhin et al. 2006), most of them were con
ducted in zerodimensional cases and the variation of
gas temperature at the droplet position was neglected.
However, in practical 3D unsteady simulations, the am
bient gas temperature and vapor mass fraction are not
constant, but vary both in time and space. In these cases,
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
the gas temperature T, in Eq. (18), which is used to
calculate interphase heat and mass transfers, becomes
a key issue for transient heating and evaporation of a
droplet. The local gas temperature obtained by inter
polation from surrounding grid points appears to be a
good choice. However, since the models for calculat
ing heat transfer between droplets and gas are typically
formulated in terms of farfield ambient gas tempera
ture or mainstream gas temperature, the employed lo
cal gas temperature must approximately satisfy this con
dition. Evidently, the local gas temperature is strongly
dependent on grid size and influenced by the cooling
effect from droplet evaporation, which imposes another
complication for droplet evaporation predictions. In the
present study, the heating and evaporation model of a
single droplet is validated against available experimen
tal data in 3D simulations initially based on constant
ambient gas temperature. Then the local gas tempera
ture obtained by interpolation is used and the gridsize
dependence is discussed.
4
t/d2, s/mm2
Figure 1: Time variations of normalized squared droplet
diameter of nheptane using the constant ambient tem
perature in the evaporation model ( Tg=741 K,
predicted; 0 Tg=741 K, exp.;  Tg=648 K, pre
dicted; Tg=648 K, exp.;  Tg=555 K, predicted;
E Tg=555 K, exp.; Tg=471 K, predicted; U
Tg=471 K, exp.).
Figure 1 shows the comparisons of the predicted time
variations of normalized squared droplet diameter of n
heptane with those obtained in the experiments (Nomura
et al. 1996). The predicted results agree well with the
experimental data, especially for the higher temperature
cases. Note that we use the constant ambient tempera
0.8
0.6 
0.4
0.2
0 1 2 3
t/d2, s/mm2
Figure 2: Time variation of normalized squared droplet
diameter of nheptane at P=latm and ambient tem
perature of 741 K for different mesh resolutions us
ing the constant ambient temperature in the evapora
tion model(0, exp.; A/do=l;  A/do=5;
 A/dn=10:  A/dn=20:  A/dn=100).
t/d2, s/mm2
Figure 3: Time variations of normalized squared droplet
diameter of nheptane at P=latm and ambient tem
perature of 741K for different mesh resolutions us
ing the local temperature in the evaporation odel(0,
exp.; A/do=l;  A/do=5;  A/do=10;
 A/do=20; A/do=100).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.5 1 1.5 2
t/d2, s/mm2
2.5 3 3.5
I
350
300 ,, ,,
0 1 2 3 4 5
t, S
Figure 4: Time variation of the droplet diameter squared
(top) and the droplet temperature (bottom) for decane
( predicted; U exp.). The experimental data are
from Wong and Lin (1992), and the conditions are:
Tg=1000K, Td=315K, do=2mm and Redo=17.
ture and neglect the variation of the local gas tempera
ture. In this case, the effect of the grid size on droplet
evaporation may be negligible, as shown in Fig. 2. How
ever, in practical simulations, the ambient temperature is
not constant, and a local gas temperature has to be used,
which obviously depends on the ratio of the cell size to
the droplet diameter. Figure 3 shows that when the grid
size is on the order of or larger than ten times the droplet
diameter, using the local gas temperature can yield re
sults consistent with the experimental data. Otherwise,
the predicted droplet evaporation will significantly devi
ate from the experiment data. This would indicate that
Figure 5: Effect of different corrections to the Nus
selt and Sherwood numbers on nheptane evaporation at
P=latm and ambient temperature of 741 K.
V0.4
0.2
1.5 2
t/d s/mm2
Figure 6: Effect of initial droplet position in a cell on
nheptane evaporation at P= latm and ambient tempera
ture of 741 K.
the cooling effect resulting from evaporation can only
influence the gas temperature within a range of about
ten times the droplet diameter, which is in good agree
ment with the recent analytical solutions of the heat con
duction equation to a sphere obtained by Sazhin et al.
(2007). It also suggests that the grid size has to be at
least ten times the droplet diameter if the local gas tem
perature is used for 3D evaporation prediction in prac
tical numerical simulations.
 B..
0.8
06
06
V0.4
0.2F
o Exp.
 No correction
 Correction A
r\
rr
CorcinB
n
I
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
The validation of the evaporation model is also ex
tended to other fuels, such as decane and hexane. Fig
ure 4 shows the time variations of droplet diameter
squared and the droplet temperature for decane. The ex
perimental data reported by Wong and Lin (1992) are
also plotted. The predicted results based on the applied
model are generally in good agreement with available
experimental data, which further validates the evapora
tion model used in the present study.
The Nusselt and Sherwood numbers are crucial pa
rameters for heat and mass transfer. Figure 5 presents
the effect of corrections A and B on droplet evaporation.
Both corrections can correctly predict nheptane evapo
ration, but the results from correction B are more con
sistent with the experimental data (Nomura et al. 1996).
In actual twoway coupling calculations, the gas
phase quantities at the droplet position are usually ob
tained by interpotation from surrounding grid points.
For certain interpolation scheme, different initial droplet
positions may influence the convergence of the solu
tions. In the present study, a trilinear interpoation
scheme is used. To check the effect of initial droplet
position in a grid cell on the droplet evaporation predic
tion based on the trilinear interpoation scheme, three
cases with different initial droplet positions in a grid cell
are tested and compared with the experimental data (No
mura et al. 1996). Here (0,0) means that the droplet is
put at the grid node, (0.5,0.5) denotes that the droplet
is put at the center of the cell, and (0.25,0,25) indicates
that the droplet is put at the location between the grid
node and the cell center. As demonstrated in Fig. 6, the
trilinear interpotation scheme is not sensitive to the ini
tial droplet position. The numerical solutions converge
to the experimental data for all cases.
DNS of a Model Combustor
Flow Configuration and Grid System. Most previous
DNS studies on spray combustion are limited to 2D or
very simple 3D configurations. In the present study, we
simulate a more realistic 3D swirl combustor. The flow
geometry of the model combustor is shown in Fig. 7.
The central swirl air of temperature 500 K is injected
through a pipe of inner diameter Din = 3.75 mm with a
mean axial velocity Uinj = 4.5 m/s and a mean swirl ve
locity Winj = 4.5 m/s. The secondary swirl air of temper
ature 500 K is injected through an annular pipe of inner
diameter Din = 5 mm and outer diameter Dout = 10 mm
with the same mean axial and swirl velocities as those
of the central one. This corresponds to a geometric swirl
number S=1.0. The combustion chamber is 40 mm wide
and 60 mm long, and the outside injection pipe is 10 mm
long. The flow Reynolds number based on the mean ax
ial velocity and outer diameter of the pipe is about 3000.
Figure 7: Geometry of the 3D model combustor.
The spray is assumed to have been fully atomized and
the resulting nheptane droplets are issued from the tip
of the wall regions between the central and the annular
pipes with temperature of 300 K. When the droplets are
issued, they are assumed to be in dynamic equilibrium
with and have the same velocities as the carrier air. This
leads to a spray cone angle of 900 and a low Weber num
ber limit. Thus, secondary breakup of droplets is not
considered in the present simulations. For droplet size
distribution, a commonly used lognormal distribution
with mean diameter 10 pm, maximum diameter 20 /m,
and minimum diameter 1 pm is used in the simulations.
In DNS of spray combustion, there are strict require
ments on grid resolution. On one hand, the grid size
has to be small enough to resolve both the Kolmogorov
length scale and the reaction zone thickness of the flame.
On the other hand, the grid size has to be around ten
times larger than the droplet size to get correct droplet
evaporation dynamics if the pointsource assumption of
droplets is used, as demonstrated above. To determine
the grid resolution, we have performed griddependence
studies. Figure 8 shows the influence of grid resolution
on the 1D axial velocity spectrum. The spectrum is
calculated based on a position in the shear layer region
with 10 mm distance to the nozzle and 10 mm distance
to the centerline for each case. It is observed that the
velocity spectrum becomes independent of grid resolu
tion when 384x 192x256 points are used. Considering
the resolution requirements for the flame, the grid res
olution with 768 points along the axial direction, 384
along the radial direction and 256 along the swirl direc
tion is chosen in the present study which consists of ap
proximately 75 million grid points. In order to capture
the important structures as much as possible, the mesh
is nonuniform. In the regions of interest, such as the
recirculation, strong shear, and near wall regions, finer
spacing is used, while in other regions a coarser spac
ing is used. This clustering is believed to resolve both
the turbulent and chemical scales in the interesting re
102
103
104
105
106
W 107
108 
109
1010
1011 _
102
103 104
Figure 8: Dependence of 1D axial velocity
spectrum on grid resolutions. k5/3;
 192x96x32 (Twophase);  384x192x128
(Twophase);  384x192x256 (Two
phase);  768x384x256 (Singlephase);
768 x 384 x 256 (Twophase).
gions. The accuracy of the mesh has been confirmed by
the DNS results.
In addition, it is also interesting to point out from
Fig. 8 that the presence of droplet evaporation and
combustion decreases the velocity spectrum in the
higher wavenumber region, but increases it in the lower
wavenumber region, compared with the single phase
case. This turbulence modulation due to reacting
droplets is consistent with the previous study of Yang
(2004).
Boundary Conditions and Computational Cases.
The flow is periodic in the azimuthal direction, and no
slip boundary conditions are used for all walls. The
downstream convective outflow condition is obtained by
solving a convection equation, allowing for a smooth
exit of all structures without perturbing the rest of the
flow significantly. To generate turbulent inflow con
ditions for the swirl air jets, separate pipe flow DNS
computations are conducted in the present study and the
data from these pipe simulations are used as the inflow
boundary conditions for the spray combustion DNS. At
first, hot air at a temperature of 1500 K is injected to ac
celerate the evaporation and trigger combustion. Then,
the temperature is reduced to 500 K. Data for post
processing are stored when the flames become approxi
mately stable.
The spray combustion DNS cases presented here sup
posedly represent the physics found in aircraft engines.
Different designs for gas turbine propulsion engine com
bustors have been proposed. Two designs that have been
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Parameter Turbine Model
Spray cone angle 90 90
Density ratio 0(103) 0(103)
Swirl number 1.0 1.0
Damkhler number Da 0(50) 0(50)
Karlovitz number Ka 0(0.1) 0(0.1)
Reynolds number Re 0(106) 0(103)
Stokes number St 2.5 10 0.1 0.4
Weber number We 0(50) 0(0.1)
Table 1: Comparison of major parameters between re
alistic gas turbine combustor and the model combustor
in the present DNS.
specifically designed to lower NOx emissions are the
socalled Lean Direct Injection (LDI) combustor and
the Richbum/Quickquench/Leanburn (RQL) combus
tor. To simulate these two cases, the droplet mass flow
rate is adjusted to achieve global equivalence ratios of
0.7 (lean) and 2.1 (rich), respectively. Here, the central
and secondary air jets rotate in the same direction (co
swirl). Recent experimental measurements (Hadefa and
Lenzeb 2008) show that counterswirl (flows rotates in
opposite directions) imposes different effects on spray
combustion. To investigate these effects, lean and rich
cases with counterswirl air jets are also simulated. To
further examine the importance of the incoming turbu
lence to spray combustion, a lean case with a laminar
tophat velocity profile inflow boundary condition is in
cluded. Additionally, a single phase cold flow case is
simulated for reference. To see how the model com
bustor in the present DNS matches realistic gas turbine
combustor conditions, Table 1 gives a comparison of
the main parameters used in the DNS and those typi
cally found in realistic engines. The main differences are
the flow Reynolds number, droplet Stokes number, and
droplet Weber number, all of which are chosen within
the capability of current DNS of spray reacting flows.
Results and Discussions
Instantaneous Vortex Structures. Swirlstabilized
combustion technology has been widely utilized in en
ergy conversion systems. The underlying mechanisms
and benefits are well documented in the literature and
depend mainly on the formation of largescale struc
tures in swirl flows, including a central toroidal recircu
lation zone which recirculates heat and reactive chem
ical species to the root of the flame and allows flame
stabilization and an outer recirculation zone (ORZ). Fig
ure 9 shows the threedimensional vortex structures for
different cases. In the reacting cases, there are vor
tex tubes originated from the KelvinHelmholtz insta
bilities due to strong shear in both swirl and axial di
rections between the CRZ and the ORZ. These vortex
tubes don't show the socalled helical dipoles (Okulov
and Fukumoto 2004), but form antiparallel vortex pairs
(Goto and Kida 2003). These selforganized structures
identify the inherent tendency of vortex tubes to align
themselves antiparallelly in turbulent flows, which is
consistent with previous numerical analysis (Goto and
Kida 2003). The selforganized vortex structures are
also observed in counterswirl flow, but the cone angle is
smaller. When looking at the singlephase case, obvious
differences appear. In the singlephase case, the vortex
structures are finer and more numerous. Since the CRZ
is not well developed, the cone angle of structures is
smaller. This also reflects the influence of droplet evap
oration and combustion on turbulence.
7'
4? '
Figure 9: Instantaneous vortex structures characterized
by Q criterion (Dubief and Delcayre 2000) for different
cases (Top left: lean coswirl case A; Top right: lean
counterswirl case B; Middle left: rich coswirl case C;
Middle right: rich counterswirl case D; Bottom left:
lean counterswirl with bulk inflow condition case E;
Bottom right: singlephase case F. Droplets are colored
by diameter).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Instantaneous Flame Structures. Previous DNS
shows that spray combustion is composed of premixed
and diffusion flames. To visualize the flame struc
tures in the present swirl spray configuration, the nor
malized flame index (Domingo et al. 2005) is used.
As seen in Fig. 10, for examples, the spray flame is
quite complicated. There are not only isolated diffu
sion and premixed flames, but also composite struc
tures, for instance, pocket diffusion flames enclosed by
pocket premixed flames, pocket premixed flames en
closed by diffusion flame sheets, and premixed flame
bands connecting diffusion flames. In addition, there
Figure 10: Instantaneous spray flame structures in
the lean coswirl case A ( : stoichiometric mix
ture fraction isoline; : diffusion flame isoline;
 : premixed flame isoline).
are local nonburning pockets in burning flames. All
these complex structures are related to the swirl fluid dy
namics and turbulent mixing, and bring igniliik.iiii chal
lenges for spray combustion modeling. In particular, the
coupled structures always consist of one thick layer and
another very thin one. It is impossible to resolve all these
thin layers in practical simulations. Although equiva
lence ratio, co or counterswirl and bulk inflow condi
tion have an influence on fluid dynamics, droplet dis
tribution and flame structures, the spray combustion is
typically characterized by lean premixed, rich premixed,
and diffusion flames for each case.
TimeAveraged GasPhase Velocity. Timeaveraged
statistics are obtained over a period of several charac
teristic time scales to provide information for the mean
gasphase characteristics. Averaging is also performed
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.8,
over the homogeneous azithumal direction to obtain a
twodimensional averaged field. Figure 11 presents
0.61
> 0.4
0 I I I I
0 100 200 300
Evaporation rate
0.8 i 1.5
> 0.4
x
2
1.5
1
0.5
0.
1
x
Figure 11: Contours of timeaveraged axial velocity for
different cases (Top left: lean coswirl case A; Top right:
lean counterswirl case B; Middle left: rich coswirl case
C; Middle right: rich counterswirl case D; Bottom left:
lean counterswirl with bulk inflow condition case E;
Bottom right: singlephase case F. The solid lines rep
resent locations of zero axial velocity).
the contours of timeaveraged axial velocity for differ
ent cases. In the singlephase case F, the CRZ is not
developed well. It has an irregular shape and smaller
. 4.
0 I l I I
0 0.5 1 1.5
d
Figure 12: Radial profiles of averaged evaporation rate
(top,x=0.25) and droplet diameter (bottom left: x=0.25;
bottom right: x=1.00)for different cases ( : case
A;  : case B;  : case C;  : case D;
 : case E).
size. On the contrary, the ORZ develops well with larger
size compared with reacting cases. This indicates that
the spray evaporation and combustion are able to en
hance the central recirculation zones, and reduce the
outer ones. Among the reacting cases, coswirl leads to
stronger CRZ, but smaller ORZ compared with counter
swirl. Although counterswirl helps to enhance turbu
lent mixing, it depresses the magnitude of the CRZ and
potentially causes instability. In addition, counterswirl
slightly increases the maximum values of axial veloc
ity. When comparing turbulent inflow case B with the
laminar tophat inflow case E, it is interesting to see that
both the CRZ and ORZ are similar inside the combus
tor, though the axial velocity values have obvious differ
ences.
TimeAveraged Liquid Phase. Some timeaveraged
x 4
X
1 0
2 3 4 5
x
r
r'
r
r
statistics for the liquid droplets are also obtained. Here
only the most interesting features are briefly presented.
Figure 12 shows the radial profiles of averaged evapo
ration rate and droplet diameter. Although the droplets
are distributed in a wider space, the evaporation occurs
in a narrower region and high equivalence ratio leads
to higher evaporation rates. It is also observed that the
profiles of droplet diameter at x = 0.25 are close to the
initial lognormal distribution, but they change signifi
cantly at the downstream location of x = 1.0, especially
for the rich cases. The larger droplet diameter is found
to be located in the spray boundaries. This means that
the droplets with larger diameter move towards the spray
boundary with increasing distance from the nozzle exit.
Previous numerous studies on muliphase flows (Squires
and Eaton 1991; Ling et al 1998) have demonstrated that
when the particle Stokes number is close to unity, pref
erential concentration will become important and lead to
high levels of particle dispersion. In the present study,
the droplet Stokes number is limited to below 0.4 due to
resolution requirement. This implies that the larger the
droplet diameter is, the stronger the influence of prefer
ential concentration, and the higher level of the disper
sion, which is consistent with the above observation.
Conclusions
nheptane droplet evaporation and combustion as well
as their coupling interactions with turbulence in a 3D
model swirl combustor are investigated by means of di
rect numerical simulation. The gasphase is formulated
in an Eulerian framework, and dispersed droplets are
tracked in a Lagrangian sense. An equilibrium evapo
ration model is used to describe the droplet evaporation,
and an adaptive onestep irreversible reaction is used for
gas combustion. Full interphase twoway coupling is
applied. The results show that the underlying grid size
has to be at least ten times the droplet diameter to get
correct droplet evaporation dynamics. Compared with
singlephase case, droplet evaporation and combustion
are able to decrease the velocity spectrum in the higher
wavenumber region and enhance the central recircula
tion zones, while increase the velocity spectrum in the
lower wavenumber region and reduce the outer recircu
lation zones.
Acknowledgements
The authors are grateful to NASA's support of this
project. K. Luo appreciates inspiring discussions
and help from Professor Parviz Moin, Professor Luc
Vervisch, Dr. Madhusudan Gurpura Pai, Dr. Edward
Knudsen, Dr. H. ElAsrag, Dr. Mehdi Raessi, and other
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
colleagues at Stanford University. Computational re
sources at Texas Advanced Computing Center (TACC)
are acknowledged. Funding from a Foundation for the
Author of National Excellent Doctoral Dissertation of
PR China (2007B4) and NSFC (50976098) are also ap
preciated.
References
Domingo P, Vervisch L. and Reveillon J., DNS analysis
of partially premixed combustion in spray and gaseous
turbulent flamebases stabilized in hot air, Combust.
Flame, vol. 140, pp. 172195, 2005
Reveillon J. and Vervisch L., Analysis of weakly tur
bulent dilutedspray flames and combustion regimes, J.
Fluid Mechanics, vol. 537, pp. 317347, 2005
Nakamura M., Akamatsu F., Kurose R. and Katsuki, M.,
Combustion mechanism of liquid fuel spray in a gaseous
flame, Phys. Fluids, vol. 17, pp. 172195, 2005
Crowe C. T., Martin S. and Tsuji Y., Multiphase Flows
with Droplets and Particles, New York: CRC Pres, 1998
Maxey M. R. and Riley J. J., Equation of motion for a
small rigid sphere in a nonuniform flow, Phys. Fluids,
vol. 26, pp. 883889, 1983
Kurose R., Makino H., Komori S., Nakamura M., Aka
matsu F. and Katsuki M., Effects of outflow from the
surface of a sphere on drag, shear lift, and scalar diffu
sion, Phys. Fluids, vol. 15, pp. 23382351, 2003
Miller R. S., Harstad K. and Bellan J., Evaluation of
equilibrium and nonequilibrium evaporation models for
manydroplet gasliquid flow simulations, Int. J. Multi
phase Flow, vol. 24, pp. 10251055, 1998
Sirignano W. A., Fluid Dynamics and Transport of
Droplets and Sprays, United Kingdom: Cambridge Uni
versity Press, 1999
Sazhin S. S., Advanced models of fuel droplet heating
and evaporation, Progress in Energy and Combustion
Science, vol. 32, pp. 162214, 2006
Sazhin S. S., Krutitskii P. A., Martynov S. B., Mason D.,
Heikal M. R. and Sazhina E. M., Transient heating of a
semitransparent spherical body, Int. J. Thermal Science,
vol. 46, pp. 444457, 2007
Law C. K. and Law H. K., Quasisteady diffusion flame
theory with variable specific heats and transport coef
ficients, Combust. Sci. Technol., vol. 12, pp. 207216,
1976
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
FemandezTarrazo E., Sanchez A. L., Linan A. and Ling W., Chung J. N., Troutt T. R. and Crowe C. T.,
Williams F. A., A simple onestep chemistry model for Direct numerical simulation of a threedimensional tem
partially premixed hydrocarbon combustion, Combust. poral mixing layer with particle dispersion, J. Fluid
Sci. Technol., vol.147, pp. 3238, 2006 Mecha., vol. 358, pp. 6185, 1998
Desjardins O., Blanquart G., Balarac G. and Pitsch H.,
High order conservative finite difference scheme for
variable density low Mach number turbulent flows, J.
Comp. Phys., vol. 227, pp. 71257159, 2008
Morinishi Y., Vasilyev O. V. and Takeshi O., Fully con
servative finite difference scheme in cylindrical coor
dinates for incompressible flow simulations, J. Comp.
Phys., vol. 197, pp. 686710, 2004
Crowe C. T., Sharma M. P. and Stock D.E., The particle
source in cell (PSICell) model for gasdroplet flows,
Int. J. Fluids Eng, vol. 6, pp. 325332, 1977
Sazhin S. S., Kristyadi T., Abdelghaffar W. A. and
Heikal M. R., Models for fuel droplet heating and evap
oration: Comparative analysis, Fuel, vol. 85, pp. 1613
1630, 2006
Nomura H., Ujiie Y., Rath H. J., Sato J. and Kono M.,
Experimental study on highpressure droplet evapora
tion using microgravity conditions, Proc. Combust. In
stit., vol. 26, pp. 12671273, 1996
Wong S. C. and Lin A. R., Internal temperature distribu
tions of droplets vaporizing in hightemperature convec
tive flows, J. Fluid Mech., vol. 237, pp. 671687,1992
Yang V., Liquid Rocket Thrust Chambers: Aspects of
Modeling, Analysis, and Design, New York: AIAA,
2004
Hadefa R. and Lenzeb, B., Effects of co and counter
swirl on the droplet characteristics in a spray flame,
Chemical Engineering and Processing: Process Inten
sification, vol. 47, pp. 22092217, 2008
Okulov V. L. and Fukumoto Y., Helical Dipole, Doklady
Phys, vol. 49, pp. 662667, 2004
Goto S. and Kida, S, Enhanced stretching of material
lines by antiparallel vortex pairs in turbulence, Fluid
Dyn. Res., vol. 33, pp. 403431, 2003
Dubief Y and Delcayre F., On coherentvortex identi
fication in turbulence, J. Turbulence, vol. 1, pp. 122,
2000
Squires K. D. and Eaton J. K., Preferential concentration
of particles by turbulence, Phys. Fluids, vol. 3, pp. 1169
1178,1991
