1. technology and computing

' THE AMERICAN SOCIETY OF MiCHANICAL ENGINEERS
34 ath Si.; New York:, let 10017 '
. •
:.1
.•
..97-07309 t .
- ,1 ..
The Society shag not be reaponsibia foiliatamejrn er cipirtiens adViincecl in iipers Of dIsailslori at Meeting's S the &coyote:4 its Divisions Or
Sections or printed in its publicathorts. Discussion is printed only It the Paget is published in an ASME Jounial. Authorization to photocopy
material for Internal or personal use under circumstance not
within the fair use newts's:ins of the Copyright Act is (rented by 'ASME to
libraries and other users registered with the Copyright CleatInce Corner (CCC) :flailsaCtionsi Reporting Sr in/Ice
provided that the base fee of so.30, •
per page is paid directly to the DOC, 27 Congress Street Seem
' MA 01970. Re2eiteafcispecial permission Of but repioduCtion -should be addiessed
tonne
Copyright 01997 by ASME
MIghtsReserved
%' • e
•'"Printed in' LISA _
III 1 11 1 111119,11 II 11 1 111
CFD ANALYSIS OF UQUID SPRAY COMBUSTION IN A GAS TURBINE COMBUSTOR
Mark K. Lal
National Research Council Canada
Ottawa, Ontario, Canada
ABSTRACT
A numerical method is presented for predicting steady, threedimensional, turbulent, liquid spray combusting flows in a gas turbine
combustor. The Wain conservation equations for gas flow and the
Lagrangian conservation equations for discrete fuel liquid droplets were
solved. The trajectory computation of the fuel droplets provided the soma
terms for all the gas-phase equations. A standard k-e subrnodel was used
for turbulence. The combustion submodel used was a global local
aquarium morel, where chankal species (C 111,. 02 COP Hp, CO, 112
and NO approached local thermodynamic equilibrium with a rate
determined by a combination of local turbulent mixing and global chemical
kinetics times. The numerical methodology for gas-phase calculations
involved a staggered finite-volume formulation with a multi-block
curvilinear orthogonal coordinate computational grid, and the PISO
algorithm. This numerical code was applied to a model gas turbine
combustor similar to that of the Allison 570KP currently in use by the
Canadian Navy. The combustor was equipped with an advanced airblast
fuel nozzle. The calculations included the analysis of the internal passages
of the fuel nozzle. Through the numerical study at full-power and lowcruise operating conditions, a better understanding of the physical processes
of flow and temperature fields inside the primary zone was obtained.
Predicted hot spots corresponded to locations where deterioration of the
combustor liner has been observed in practice.
NOMENCLATURE
as...;
A's
Br
C„
C;,C.2
Cd
cd
cr
cPI
coefficients of the quartic equation, Eq. (49)
empirical constants, Eq. (59)
Spalding's mass transfer number, Eq. (25)
empirical constant in turbulent viscosity, Eq. (7)
empirical constants in the dissipation rate, Eq. (9)
aerodynamic drag coefficient, Eq. (15)
droplet liquid specific heat
specific heat at constant pressure
specific heat at constant pressure of species i
drag function, Eq. (14)
specific internal energy of liquid droplet
rate change of momentum due to phase interaction
generation term
specific Gibbs free energy
specific stagnation enthalpy of mixture, including sensible,
chemical and kinetic energy effects
specific enthalpy of mixture
specific enthalpy of liquid droplet
specific enthalpy of species!
specific enthalpy of vapour droplet
unit tensor
heat of formation of species i at standard conditions
hiD
K.
K8
equilib
ri um constant
co
turbulence kinetic energy
latent heat of vaporization
dissipation length scale, Eq. (18)
•
mass flow rate of liquid
droplet mass
md
•
total number of droplets
Nusselt number
Nu
•
molar concentration of species i
number of droplets per unit time along a droplet trajectory j
143
•
pressure
saturated vapour pressure of the fuel vapour
•
Prandd number, Eq. (31)
Pr
heat transfer rate at droplet surface
Oe
rate change of energy due to phase interaction
mean fomiation rate of species i due to chemical reaction
Ri
Reynolds number, Eq. (16)
Re
RN
univenal gas constant
droplet radius
source term induced by the computational grid
god
Schmidt number, Eq. (26)
Sc
Sherwood number, Eq. (24)
Sh
Presented at the International Gas Turbine & Aeroengine Congress & Exhibition
Orlando, Florida — June 2,5,1997
Downloaded From: http://proceedings.asmedigitalcollection.asme.org/ on 12/29/2014 Terms of Use: http://asme.org/terms
INTRODUCTION
Sauter mean diameter
specific entropy
temperature
liquid droplet temperature
reference temperature, Eq. (60)
Favre-averaged velocity vector of gas phase
liquid droplet velocity vector
gas turbulence velocity vector
volume of computational grid
molecular weight
local average molecular weight of all species exclusive of fuel
vapour
molecular weight of species i
CO -CO2 mole ratio, Eq. (50)
droplet position vector
mass fraction of species i
Due to recent advances in both computer hardware and the
development and application of mathematical modelling and numerical
methods, the last two decades have seen a major contribution to engineering
design and development of gas turbine combustors from the science of
Computational Fluid Dynamics (CFD) (So et al., 1986). CFD has proven
to be an important tool in reducing product development costs, optimizing
designs, and improving our understanding of physical processes in advanced
combustors, especially the combustor primary zone.
This paper outlines a methodology for simulating steady, threedimensional, turbulent, liquid spray combusting flows in a gas turbine
combustor. A CFD code based an this methodology is applied to a model
gas turbine combustor incorporating many geometrical features of the
special low smoke combustor of the Allison 570KF, as delivered to the
Canadian Navy in 1987. It should be noted that the combustor under study
has now been superseded in service by a more advanced design effusion
combustor. This combustor is equipped with an advanced airblast fuel
nozzle. As Fuller and Smith (1993) have shown that CFD results are
extremely sensitive to the fuel nozzk/swirler boundary specification, the
present CFD analysis includes the modelling of the internal filet nozzle
passages to obtain better prediction of conditions in the combustor primary
zone.
Preliminary three-dimensional (PD results for gaseous fuel
combustion in this combustor have been reported previously (Lai and
Cheney, 1995). Numerical simulations of liquid spray combustion in the
primary zone at full-power and low-cruise operating conditions are described
in the paper. It will be shown that predicted hot spots in the combustor
correspond to locations where deterioration of the combustor liner has been
observed in the some of the engines operated by the Canadian Navy
(Cheney and Lai, 1995).
Greek Symbols
ac ,aH ,a0 total number of atom per unit volume
V
gradient operator
Kronecker delta
at
timestep
turbulence energy dissipation rate
thermal conductivity
molecular viscosity
turbulent viscosity
P,
v,
stoichiometric coefficient of species i
mass density of mixture
liquid droplet density
Pd
rate
change of mass due to phase interaction
PS
pD
fuel vapour diffusivity in gas phase
equivalence ratio
variance of the Gaussian probability distribution
Prandtl number for H
Schmidt number for k
Jk
air
Schmidt number for chemical species
Schmidt number for e
at
summation
characteristic reaction time
rd
droplet relaxation time, Eq. (21)
eddy lifetime, Eq. (19)
interaction time, Eq. (22)
laminar conversion time, Eq. (39)
residence time, Eq. (20)
rt
turbulent conversion time, Eq. (40)
random number
MATHEMATICAL MODEL
Gas Phase Equations
-
(
For compactness, the conservation equations are written in vector
notation. The Favre-averaged conservation equations of mass, momentum,
chemical species, and energy are given by
9P0) = Ps
V.[puel - (P*1 1 ,)V 01 = Vq(11 +11 ,XVO) r
- (P + i(pk + (p +p i fiTO))i) +
v.(pOY, - Ent +lid/My)? =
P5
,
Pokr
SUDefSCriDte
•
(overbar) average of droplet surface and ambient values
transpose of tensor
equilibrium value
and
+ 111)430270 =
Os •
60 is the Kuonccker delta and species 1 is the species of which the spray
droplets are composed. All the phase interactions source terms
(( s, Fs , Q5 ) and the mean formation rate of specks i due to chemical
reaction ft; are defined later. The expression of the gradient operator V
depends on the particular coordinate system used (Lai, 1992). The values
of off and ay are assumed to be equal to unity.. Thermal radiation is
neglected.
To complete the equation set, the following relations are used
Subscripts
I, fuel
VqP0H
fuel vapour
liquid droplet
droplet surface
2
Downloaded From: http://proceedings.asmedigitalcollection.asme.org/ on 12/29/2014 Terms of Use: http://asme.org/terms
(5)
{ (24 « 4Re 213)1Re
and
Re < 1000
Cd
(15)
0.424
11 h • !CPO • k ,
(6)
where
and
p, = Co pe/e .
(7)
and
where
v.(pee
-
((p + pycia rk) = G - pc
+ p i )/ar re) = (C IG - C 2pe)etk
-
G = (1),1(V0) C70)TI - i(Pk
P,VO)f}:\70 .
Re = 2rplo • a
.
-
(16)
Each component of gas turbulence velocity d is assumed to follow a
1[73Tic and is calculated randomly from
Gaussian distribution with a =(2
The modelled equations for k and e and their empirical constants are given
by (Launder and Spalding, 1974)
(pok
Re 2 MOO
1
( 17)
(IC1) .
J = (Jo) siial(C)Eo0
(8)
is a random number selected from an uniform distribution in the range
di
1 ' Ern is the complementary error function. The Willies of the
complementary error function Erre are evaluated and stored in tabular form
at intervals 0.05, and are calculated by a linear interpolation at intermediate
values of
The characteristic eddy size is assumed to be the dissipation
length scale Id given by
(9)
1(1.
(10)
= Vie
klauld-Phase EauatIong
and the eddy lifetime re is calculated from
Following the practice of Gasman and loannides (1983) and Dukowicz
( 1980 ). the liquid phase is modelled in a stochastic manner as a spray of
computational droplets. Each computational droplet represents a number
of physical droplets of identical size, location, velocity, and temperature. For
each droplet the conservation equations are written in Lagrangian form and
in the same coordinate system as for the gas phase.
In an airblast fuel atomizer, fuel is supplied at low pressure (i.e. low
velocity) to a filming surface where the shearing action of a high velocity
swirling air acts to atomize the fuel. The liquid film is broken into ligaments
and then into many sizes of droplets. Droplet-droplet interaction has been
found to be important in the initial breakup region where coalescence and
collision as well as secondary atomization occurs. In the present study, the
drop breakup and atomization processes are not modelled and the liquid
spray is assumed to be thin permitting volume-displacement (Dukowicz,
1980) and other thick-spray (O'Rourke, 1981) effects to be neglected. The
liquid is assumed to enter the combustor as a fully atomized spray
comprised of spherical droplets.
The droplet mass is given by
= te/ 0.Th
2
(19)
The residence time r, is determined from the solution of a linearized
equation of motion of the droplet in a uniform flow (Gasman and loannides,
1983)
rr. = -r4 Ink - Verd i() +
1 ..
t d = (8/3)rpo /(pCd1U + a
where
did
- Od1 ) •
(20)
(21)
The interaction time in an eddy T i is given by
ti
=
>
te
min(ro ,;)
10
•
a - Odl td
(22)
10 • a - Odl td
To calculate droplet mass and heat transfers with the surrounding gas,
a uniform temperature model of a single droplet is used. The evaporation
rate is given by the Frossling correlation (Faeth, 1983)
MdEPer 3
The droplet trajectory and momentum equation are, in vector (cam
did =
dt
(12)
Ud
dme = - 2nrpDln(1 + R y )Sh
dt
(23)
Sh = (2 + 0.6Re 112 Sc "3 ) ,
(24)
=
(25)
where
and
dO
A (u a
.
di
_ od)
grid
(13)
fir
md
The source vector
contains curvature terms arising from the use of the
tad
coordinate system for the computational grid. The drag function D d ,
assuming that only the aerodynamic drag affects the droplets, is given by
(Faeth, 1983)
Dd
2
d1°
inr2C
g °di
(Ynsii
YaridA I
Yfuels) •
Sc = —
pint) .
and
(26)
The fuel vapour mass fraction is obtained from
Yneir =
(14)
+ Wdrll'v
-
(27)
where Wo is the local average molecular weight of all species exclusive of
fuel vapour and P, is the saturated vapour pressure of the fuel vapour at
the droplet temperature Td.
3
Downloaded From: http://proceedings.asmedigitalcollection.asme.org/ on 12/29/2014 Terms of Use: http://asme.org/terms
The rate of droplet temperature change is determined by the
conservation energy equation at the droplet surface
drn,,
d71,
m d Cil =
di
=
=
d
The heat transfer rate to the droplet is given by the Ranz-Marshall
correlation (1952)
= 2nri(T - Td )Nu
where
(29)
Nu = (2 + 0.6Re 112 Pr i8 )1n(1 + dill y
E
au W ,(n) - n;)
=
- ni) .
This ensures that the longer of the laminar conversion time r i and the
turbulent mixing time it controls the overall conversion from one species
to another (Spalding, 1971). r, is obtained by matching computed and
measured laminar burning speed for mixtures of n-octane and air over a
lunge of equivalence ratios rb, pressures, and temperatures (Kuo and Reitz,
1989)
Phase Interactions
The droplets and gas interact by exchanging mass, momentum, and
energy (Migdal and Agosta, 1967; Crowe et al., 1977). The actions of the
gas phase on the liquid phase are introduced through the drag function
terms in the liquid-phase momentum equations, and the heat transfer terms
in the liquid-phase energy and mass conservation equations. As in the PSICELL method (Crowe et al., 1977), the influence of the liquid phase on the
gas phase is incorporated by adding source terms into the relevant gas-phase
conservation equations.
At each computational cell traversed by the droplets, the rate of mass,
momentum, and energy transferred to/from the gas phase for the cell
volume V are, respectively
E {ill(m 4 ) b, - (fil d ),„„)i
Ps v = E
and
—E
(38)
tt •
(31)
Pr
Os =
(37)
For the present application,
(30)
t
and
(36)
is the local and instantaneous thermodynamic equilibrium value
of the molar concentration of species j.
= Oat local equilibrium, and
a y is the elements of the lacobian -RR I /8(WI n) • Assuming that
= 0 (i * j) and a, = -lit (i = j), then only one characteristic
reaction time t for all species needs to be defined and Eq. (36) becomes
(28)
+
+
.
.75x I 0 -8 TP ""expi( I -0.0810-1.151)15098/Ti .
(39)
t, is given by a modified form of the eddy-breakup model of Magnussen
and Hjertager (1976)
r, = (kit)/
min[4, 2(Yup+ Yco, *Yco *Yrr, )/(Ycji, - YC,tt p * Yo,* Yo)1 •
The values of n; are derived as follows.
The algebraic equations for the thermodynamic equilibrium
concentration of the 6 reactive species consist of the 3 atom balance
equations (C-H-0), the relation for the zero equilibrium fuel concentration,
and 2 non-linear equilibrium relations for the dissociation of CO 2 and H20.
These equations are
(32)
i(nd ud ),,, — (nd ud ).,k, ,
(rif(ind hd ). - (inti k
= I
(33)
(34)
+n
The summation is over all computational droplets which traverse the
computational cell. The subscripts in and out refer to values, respectively,
at entering and leaving the computational cell. ri d is the number of droplets
per unit time along a droplet trajectory j.
yn
2
nc.,) = xnejl,^co,
+ 2n, 0 + 2n,;, = y
2n(7
+ 2n,p +
nc.c, = ac ,
2
(41)
= an , (42)
2"Co , no
c = 2n0,+ 2nCO, + nCO = a0
•
(43)
Combustion Model
The global local equilibrium model of Abraham et al. (1985) is used.
Since this model allows H2 and CO and their equilibria to be included, the
proper heat release is approached for all equivalence ratios. Also, the model
uses only one characteristic conversion time for the achievement of such
equilibrium and thus is computationally efficient. A brief description of the
model is given below.
Following the practice of Reitz and Bracco (1982, 1983), the following
one-step global irreversible reaction is considered
n
= 0,
(44)
(45)
Kco, =
K H 20
. ,(ni;,onc.o) -I •
4.2 n CO
( 46)
The equilibrium constants are calculated from
CH +(02 + 3.76N 2)
Y
v CO, CO 2 + v H,0 H2 0CO
+ v CO + v H, H 2
+
376v 0 1 N 2 '
(35)
-RaTInKcp, gco, gco - "go,
The chemistry source terms in Eq. (3) are initially expanded up to first order
about local equilibrium
and
Tin K,,p = gco, *
After rearranging Ns. (41-46), we have
4
Downloaded From: http://proceedings.asmedigitalcollection.asme.org/ on 12/29/2014 Terms of Use: http://asme.org/terms
— :CO — 8 /4,0 •
(47)
(48)
a4X 4 • a3X 3 + a2X 2 + NA' + av = 0
T"/(T • Ao2 ) ,
p=
(49)
pD =
(59)
= AAI TLNT+ Al2) .
X a nco
.
where
,
(50)
The empirical constants O w A pr Am , A m, A u , A 22) are obtained from
Amsden et al. (1989). The ovetbar values of p. pD, and A denoting the
avenges of the droplet surface and ambient values are calculated at the
reference temperature Tie given by
= K„ ,„(ac - a0) ,
- (1 • K119 )a0 +
a3 = (1 •
2ac
a2
(51)
a 11 t2 aoKtip
ao
= rd
+ Kilo ) , ao = 2Kc.0
-2 2 .
=
nc:0a ac X(1 +
,
,
L =
no% = (Xit'co, ) -2 ,
(61)
Ad = ed • "v/pd .
XX„,0)1 .1 .(52)
=0.
(62)
Latent heat and vapour pressure tables are obtained from Vargaftik (1975).
Because P. varies acidly with T, its value is stored at intervals of 10K up
to the fuel critical temperature Tait . I. is stored at intervals of 100 K. The
values of ed are calculated at intervals of 100 K from Eqs. (61, 62). Their
intermediate values within 100 K intervals are obtained by a simple linear
interpolation. The value cd is approximated by the difference between
adjacent tabular values of e d , divided by 100 K.
Thermophyslcal Relationshipa
Br the gas phase, an ideal gas mixture is assumed to compute the state
relations. The specific Gibbs free energy go is used to obtain the
equilibrium constants at constant pressure for the combustion model
calculation, where
gi = hi - Tsi
- hd
where h. is the vapour enthalpy and the liquid enthalpy hd is defined by the
relation
no; = a„12(1 + XX,10)] -1 ,
/11 1. 2 = a ti n(1110 12(l
(60)
The liquid density pd is assumed to be constant. Since the latent heat
of vaporization I. is the energy required to convert a unit mass of liquid to
vapour at the equilibrium vapour pressure
The coefficients a, in Eq. (49) are calculated using the local values of gas
temperature and density. After the quartic equation fE.q. (49)1 is solved for X
by the non-iterative method of Salzer (1960), the equilibrium molar
concentrations ni' are computed from
neo
. = ac (1 •
(T - Td)/3
(53)
SOLUTION METHOD
The specific enthalpy hi is defined as
Grid Generation
hi = fcp.o dT • hi...
(54)
13
1
For each chemical species i, the specific enthalpy h, and specific entropy s,
are specified, respectively, as functions of temperature (Gordon and
McBride, 1971)
12
hi Wi /R. = a i r+ a2T 2/2 • a3T3/3 • a4r/4 • 0271/5+ ad (55)
and
s WI /R N
a l tar+ a 2 T + a 3 T2/2 • a 4 T3/3 • a 5 T4/4 + a
E Ya k .
;
1O
17
(56)
I;
76
5;
4
3
I
2
Fig. 1: Schematic Diagram of the Combustor
Showing Air Entry Features
Dome Alr —a-
(57)
Outer AlrOuter
To compute p. the equation of state is used
Inner
P = p7R.E(Yo lWi ) .
20
18
18
Vol Nozzle
Two sets of least square coefficients (a i ,...,a2 ) an tabulated for two
adjacent temperature intervals, 300-1000 K and 1000-5000 K, with the data
constrained to be equal at 1000 K. The values of go and hi are tabulated
at 100 K intervals. Their intermediate values within 100 K intervals are
obtained by a simple linear interpolation. After the mixture specific enthalpy
Is obtained from the solution of the energy equation, the corresponding
value of r is calculated via a linear search algorithm and the following
equation
A
14 15
! i
(58)
t
Inner Fuet —/
The molecular vErosity, the vapour diffusivity, and the thermal conductivity
are calculated, respectively, from the empirical correlations
Fig. 2: Schematic Diagram of the Fuel Nozzle and
Air Swirlers
5
Downloaded From: http://proceedings.asmedigitalcollection.asme.org/ on 12/29/2014 Terms of Use: http://asme.org/terms
outlet boundary of the computational domain was extended to the plane at
the outlet of the last set of cooling louvres, which was also the beginning of
the converging section of the combustor.
Figure 1 shows a schematic diagram of the combustor which was
analysed. Annotations on this figure indicate the positions of air inlet
features. The combustor uses many geometrical features of the low smoke
combustor of the Allison 570KF engine as delivered to the Canadian Navy
in 1987. Since the combustor is annular and periodic with 16 fuel nozzles,
a 22.5° sector has been modelled. The fuel nozzle is located centrally in the
sector and consists of a pressure atomized primary fuel nozzle and a main
airblast fuel atomizer with an inner and two outer air swirlers (Fig. 2). The
main fuel passages are not swirled. The inner air swirler is counter-rotating
relative to the outer swirlers and has a 10 ° swirl angle. The vanes for the
outer air miners are at 60 ° and are co-rotating. The combustor employs 4
types of cooling openings: splash cooling skirts (#6, #8, #10, #11, #13 and
#15), stacked rings (#9 and #12), liner wall cooling louvres (#1, #3, #5,
#16, #18 and #20), and rows of radial injection orifices entering the primary
(#7 and #14), intermediate (#4 and #17) and dilution (#2 and #19) zones.
One of the inner intermediate orifices, the outer intermediate orifices, and
the splash cooling skirts on the inner and outer liner walls are aligned with
the centre of the atomizer. The remaining cooling openings are staggered.
Numerical Procedure
Numerical solutions for the gas-phase equations were obtained using
the TURCOM-Bit computer code (Lai, 1992). This code used a staggered
finite-volume formulation with a curvilinear orthogonal coordinate
computational grid. The conservation equations given previously were
dtrcretized with the hybrid differencing scheme and then solved by the PISO
algorithm (Issa, 1986). Each set of discretized equations was solved
sequentially in an iterative fashion using a block iteration procedure. The
iteration procedure was organized to sweep alternatively across different
direction planes. in this sequence, the equations on each plane were solved
by the line-by-line method.
The conservation equations for each droplet were integrated in time
with a fourth-order Runge-Kutta method, using a timestep which was
dynamically adjusted based on the droplet velocity, grid cell size, and
interaction time. For numerical stability and accuracy, the timestep was also
restricted to be less than 10's. The equations were integrated starting from
the injection location at which the initial conditions were stochastically
prescribed. The integration proceeded until the droplet left the calculation
domain, evaporated to a negligible size, or reached the combustor walls.
During the integration, the gas-phase properties appearing in the droplet
equations were prescribed from the prevailing values at the nearest node of
the computational grid.
The overall solution was obtained by iterating between the calculations
of gas- and liquid-phase equations. First, a prescribed number of droplets
was introduced and their trajectories were calculated according to the above
procedure. For each computational cell traversed by the droplets, the phaseinteraction terms Eqs. (32-34) were computed and then inserted into the
gas-phase equations. Then, a prescribed number of iterations were
performed for the gas-phase equations. This cyclic process was repeated
until the gas-phase calculation converged.
(a)
(b)
Initial And 13g.mnsienanditIOne
In the absence of experimental data, the droplets were introduced at the
location near the exit plane of the outer fuel atomizer in the manner
described by Amsden et al. (1989). To form a spray from these droplets, the
mass flow rate of fuel injected M was apportioned among N droplets
according to the relation
Fig. 3: Computational Grid for the Combustor Sector:
(a) Plane Passing Through the Axes of Fuel
Nozzle and Combustor, (b) Perspective
= (rimd))
Figure 3 shows the computational grid used for the 22.5° sector of the
combustor. The curvilinear orthogonal coordinate computational grid was
a single-block structured grid consisting of several grids, each of which was
generated by an algebraic/numerical method. In the region outside the fuel
nozzle, the grid was generated numerically in diametral planes and was
rectilinear in the axial direction. The procedure used to generate the
curvilinear orthogonal grid involved solving the inverse Laplace equations
for the Cartesian coordinate positions of the grid points with multi-block
grid topology (Gasman and Johns, 1979). The multi-block orthogonal grids
were generated by subdividing a given geometry into a number of blocks.
Each block was gridded independently of the other blocks, but the
continuity of grid lines was maintained across the neighbouring block
boundaries. In the fuel nozzle, the grid was generated numerically in radialaxial planes and was symmetric about the centerline of the nozzle in the
circumferential direction. The inlet boundary started just downstream of the
swirl vanes and at the outlet planes of the fuel and air flow passage. The
for
j= I
N.
(63)
Each droplet was assigned different initial conditions of size, position, and
velocity, using a uniform random distribution. The initial temperature of the
droplets was assumed to be the delivery temperature at the atomizer. The
droplet radius was sampled randomly from a specified size distribution
where the most droplet mass occurred The initial position of the droplets
was varied randomly within a specified number of radial locations at the exit
plane of the fuel nozzle. The injection velocity Ud was chosen randomly
from velocities in the range
0
Od 0 .
(64)
The treatment of droplets impinging on the combustion chamber walls
is one of the most difficult aspects of modelling spray flames. The droplets
may evaporate, shatter, and/or reflect with reduced momentum after
impinging on combustor walls. In the absence of detailed experimental
verification, the droplets were assumed to undergo instantaneous
6
Downloaded From: http://proceedings.asmedigitalcollection.asme.org/ on 12/29/2014 Terms of Use: http://asme.org/terms
vaporization upon collision with the walls (El Banhawy and Whitelaw,
1980).
The combustor Walls were assumed to be no-slip, impermeable and
adiabatic. Boundary layer drag and heat transfer were calculated using the
wall function method (Launder and Spalding, 1974). Inlet conditions for
k were calculated using a turbulence intensity of I% of the resultant velocity
through an inlet opening in question, while the length scale used to obtain
a (Eq. 18) was chosen to he 3% of a representative dimension of the
opening. At the exit of the combustor, the normal velocity was calculated
from the integral mass balance and zero axial gradient was imposed for
other variables. Periodic boundary conditions were applied on the
circumferential end planes of the combustor sector.
On the outer liner wall, a hot spot is located around the region at the exit of
the splash cooling skirt (#15) and behind the liner wall cooling louvre (#16).
On the inner liner wall, a hot spot is seen behind the liner wall cooling
louvre (ID) and between the adjacent fuel nozzles. Referring to the dome
details, a limier temperature hot spot is located between the adjacent nozzles
and spreads to the lip of the outer swirler. These hot spots are mainly
created by the swirling and separation processes of the dome and outer air
swift flows, as will be discussed later. Good agreement between the
predicted positions of hot spots and the observed locations of the
deterioration of the real combustor is seen in Fig. 5.
RESULTS AND DISCUSSIONS
Table 1: Operating Conditions for CFD Calculations
Air Flow
Rate
(MP.)
AO
Temperature
CIO
Low.
Clone
0.274
Rill.
Power
0.997
Comiffion
Air
Pressure
(kill)
Fuel
Temperaiure
(K)
Rid FiClw
Rae
(kiin
410.0
7.5
350.0
OM
618.0
18.7
403.0
0.375
Numerical calcul lions were carried out at full-power and low-cruise
operating conditions (Table I). The thermochemical properties of Navy
distilbte (diesel) fuel were assumed to be equivalent to those of n-dodecane
(C 12l-1,). Using the specified villUti of fuel mass flow, compressor delivery
air mass flow, temperature, and pressure, the flow splits and corresponding
flow angles were calculated from the specified geometries of combustor,
wirier and atomizer. A uniform inflow velocity was calculated for each
inlet using the specified mass flow and a discharge coefficient of unity. The
total number of droplets N used was 1030. Numerical studies performed
using different values of SMD Skated that a value of I Opm gave the most
physically plausible results. Thus, this set of the results is presented here.
All the numerical results were obtained using 60x40x43 grid points,
which was the maximum number possible with the available computing
hardware. Further simulations using a liner grid and/or a higher-order
discretization scheme would he required to assess the effect of discretization
(i.e. numerical) errors on the results, but such errors are likely significantly
smaller than the mathematical modelling errors due to the complex physics
involved. The results are nonetheless believed lobe valid in a qualitative
sense and therefore useful for the engineering purpose of identifying
possible causes or the combustor liner deterioration referred to in the
Introduction.
The solution was assumed to have converged when the sum of the
absolute residuals in the gas-phase conservation equations, normalized by
the total inflow of the conserved quantity in question, had fallen below
10 -3 . To achieve convergence, the number of iterations required was at
least 3,000. Each cycle of the liquid-phase calculation was performed after
every 250 iterations of the gm-phase calculation. The CPU time on a
Power Indigo2 workstation was 50 s per iteration for gas-phase calculation
and about I hr per one cycle for liquid-phase calculation.
Figure 4 shows numerical result, for the gas temperature distribution
along the liner walls and dome details, The results for the two flow
conditions are very similar in distribution although different in magnitude.
Fig. 4: Predicted Gee Temperature Distribution Along
Dome and Liner Walls: (a) Low-Crulae,
(b) Full-Power
Figure 6 shows numerical results for the equivalence ratio and the
velocity vector distributions along the plane passing thmugh the axes of the
0072le and combustor. In the equivalence ratio plot, colour represents the
value of the equivalence ratio, defined as the ratio of the calculated fuel-air
ratio to the smichiometric value. In the velocity vector plot, each arrow
shows the magnitude and direction of the gas velocity while the colour
represents the calculated value of the temperature. Figure 7 shows predicted
particle traces through the primary zone and the dome region. The traces
represent streamlines of the flow and colour represents the calculated
temperature value. Figures Sand 9 show the fuel spray trajectories at two
operating conditions. In the plots, the colour represents the value of the
7
Downloaded From: http://proceedings.asmedigitalcollection.asme.org/ on 12/29/2014 Terms of Use: http://asme.org/terms
droplet size. The trajectories end when droplets have completely evaporated
and provide information about the distance traversed by the fuel spray.
(a)
11.1•SICO
t.
IS
ILO
-- 777-- frsta.
(k
"
7-' )•-•-' 91,"; I'
i
,. :e•nelleffiji s4
D
t
.,
..%, - '-.
• 1 A • . .
rrrra
(a) View From Combustor Ent of Dome and Liner Damage
-
-
tea
(b)
Torrenbas
_
DOS
SS
1210/0 ifla
1110130
IMO
Fig. 6: Predicted Distributions of Equivalence Ratio
and Velocity Vector for Full-Power Condition:
(a) Equivalence Ratio, (b) Velocity Vector
(I)) Borescope View of Horde end Dome Details
Fig. 6: Close-Up of Dome and Liner Walls of the Combustor
Hardware Showing the Damage Areas
Figure 6 shows that there are two distinct reaction zones in the
combustor primaty zone. A Central Recirculation Zone (CRZ) is created
along the nozzle centerline by the co- and counter-rotating airflows through
the duce passages inside the nozzle. The CRZ behaves as a toroidal vortex.
A Dome Interaction Zone (DIZ) is created by the sudden expansion of the
nozzle airflow discharging into the combustor. The DIZ is located in the
region between the dome wall, liner wall and CRZ. In the primary zone, the
toroidal swirls from the adjacent fuel nozzles interact both with the dome
details and with each other, producing a complex reacting nowt -mid. Since
the nozzle flow is strongly swirled, it generates large Slat and axial
pressure gradients, producing a CRZ which is large enough to separate the
Drz above and below the nozzle (see also Fig. 7). The CRZ extends to the
inner and cuter liner wall in the primary zone and terminates at the axial
plane of the intermediate holes (#4 and #17).
SAO
000A IMO 'Na. nal
8•000
Fig. 7: Predicted Distribution of Particle Trace at
Full-Power Condition: (a) Primary Flame Zone,
(b) Dome Region
8
Downloaded From: http://proceedings.asmedigitalcollection.asme.org/ on 12/29/2014 Terms of Use: http://asme.org/terms
5.0 pm
7.5 pm
10.0 pm
Droplet Radius Orderers)
Owlet Rens Colon)
0.0
10
40
10
10.0 pm
0.0
0.0
10.0
2.0
4.0
ILO
0.0
10.0
Flg. 9: Fuel Spray Trajectories at Low-Cruise Condition
for Different injection Droplet Radii
Fig. 8: Fuel Spray Trajectories at Full-Power Condition
for Different Injection Droplet Radii
Figure 7 shows that the primary jets (07 and #I4) are deflected
sideways by the CRZ and penetrate across the combustor. Inside the DIZ,
a complex secondary flow is observed. This flow causes the cooling
airflows from the dome wall (#9 to #I2) and the primary jets to mix with
the rest of the dome air in the DIZ and then to flow towards the dome
region of the adjacent nozzle. Due to the action of turbulent mixing with
the hot combustion products in the CRZ, this secondary flow creates a hot
spot along the dome wall and the rp of the outer swirler (see also Fig. 4).
These results indicate that the cooling airflows from the dome wall and the
primary jets ate inadequate to cool the dome wall and to maintain separation
of the very hot CRZ from the combustor liner walls.
Figures 6 and 8 show that the fuel droplets are convected downstream
along the shear layer fanned between the CRZ and DIZ. The main
combustion process occurs along the shear layer and the highest intensity
combustion region is located near the axis of symmetry. The silt of the
extremely fuel-rich region is small, which is desirable to minimize soot
formation. No droplet enters the CRZ and DIZ and most of the droplets
evaporate before impinging on the combustor liner walls. Since the droplets
are small, the droplet trejectories are seen to be greatly influenced by the gas
flow corning out of the fuel nozzle. Comparison between Fig. 8 and Fig. 9
shows that the liquid fuel traverses a much greater distance before
vaporizing when the combustor operates under the low-cruise condition.
Also, under the operating conditions of interest here, the residence time of
the liquid fuel droplets within the combustor is seen to be significant so a
single-phase calculation, which assumes that the fuel is fully vaporized,
would be of questionable value. Comparison between the gaseous fuel
combustion predictions (Lai and Cheney, 1995) and the liquid spray
combustion results presented here will be reported elsewhere.
CONCLUSIONS
A CFD method has been presented for predicting steady, threedimensional, turbulent, liquid spray combusting flows in a gas turbine
combustor. Numerical solutions at full-power and low-cruise operating
conditions were obtained for a model combustor similar to that used in the
Allison 570KF gas turbine. The calculations included the analysis of the
internal passages of the fuel nozzle. Predicted hot spots of the liner and
dome walls have been shown to correspond to areas where deterioration of
the combustor liner has been observed in practice. Flame patterns in the
primary zone have been qualitatively captured in the CFD analysis,
providing an understanding of the physical processes involved thereby
potentially enhancing the combustor design process. Further work is
requited to evaluate and improve the accuracy of the simple method used
to specify inlet velocities and mass flows, the submodels of turbulence and
chemistry, neglect of radiatizo, the spray injection model, and the numerical
differencing.
ACKNOWLEDGEMENTS
Financial support for the CFD program was provided by the CRAD,
Department of National Defence Canada (DND), Contract
No. 220794NRC09, under the Project Manager of Mr. P. Cheney and the
Scientific Authority of Mr. R.R. Hastings. The views presented are those
of the author and do not necessarily represent those of DND.
9
Downloaded From: http://proceedings.asmedigitalcollection.asme.org/ on 12/29/2014 Terms of Use: http://asme.org/terms
REFERENCES
Migdal, D., and Agosta, V.D., 1967, "A Source Flow Model for
Continuous Gas-Particle Flow", Journal of Applied Mechanics, Vol. 35,
pp. 860-865.
O'RourIce, PJ., 1981, "Collective Drop Effects in Vaporizing Liquid
Sprays", Ph.D. Thesis, Princeton University, and Los Alamos Ntaional
Laboratory Report LA-9069-T.
Patankar, S.V., and Spalding, DM., 1972, "A Calculation Procedure
for Heat, Mass and Momentum Transfer in Three-Dimensional Parabolic
Flows", International Journal of Heat and Mass Transfer, Vol. 15, pp.
1787-1806.
Ranz, WE., and Marshall, W.R., 1952, "Evaporation from Drops",
Chemical Engineering Progress, Vol. 36, pp. 141-146 and pp. 173-180.
Reitz, R.D., and Bracco, F.V., 1982, "Toward the Fommlation of a
Global Local Equilibrium Kinetics Model for Laminar Hydrocarbon
Flames", Notes on Numerical Fluid Mechanics, Vol. 6, Edited by N. Peters
and J. Wamatz,•Viewag-Verlag.
Reitz, RD., and Bracco, F.V., 1983, "Global Kinetics Models and
Lack of Thermodynamic Equilibrium", Combustion and Flame, Vol. 53,
pp. 141-143.
Salzer, H.E., 1960, "A Note on the Solution of Quartic Equations",
Mathematics of Computation, Vol. 14, pp. 279 281.
So, R.M.C., Whitelaw, J.H., and Mongia, H.C., (Editors), 1986,
Calculations of Turbulent Reactive Flows, AMD Vol. 81, ASME, New
York.
Spalding, ail, 1971, "Mixing and Chemical Reaction in Steady
Confined Turbulent Flames", Proceedings of the 13th International
Symposium on Combustion, The Combustion Institute, pp. 649-657.
Tolpadi, A.K., 1995, "Calculation of Two-Phase Flow in Gas Turbine
Combustors", Journal of Engineering for Gas Turbines and Power, Vol.
117, pp. 695 703.
Vargaftik, N.B., 1975, Tables on the Thermophysical Properties of
Liquids and Gases, John Wiley and Sons Inc.
Abraham, J., Bracco, F.V., and Reitz, RD., 1985, "Comparisons of
Computed and Measured Premixed Charge Engine Combustion",
Combustion and Flame", Vol. 60, pp. 309 322.
Amsden, A.A., O'Rourke, PJ., and Butler, T.D., 1989, "KIVA-II: A
Computer Program for Chemically Reactive Flows with Sprays", Los
Alamos National Laboratory Report LA-11560-MS.
Borman, G.L., and Johnson, J.H., 1962, "Unsteady Vaporization
Histories and Trajectories of Fuel Drops Injected into Swirling Air", SAE
Paper 598C.
Cheney, P., and Lai, M.K.Y., 1995, First Progress Report, the 15th
Meeting of the Defence Gas Turbine R&D Advisory Committee, Ottawa,
February.
Crowe, C.T., Sharma, M.P., and Stock, D.E., 1977, "The ParticleSource-In-Cell (PSI-CELL) Model for Gas-Droplet Flows", Journal of
Fluids Engineering, Vol. 99, pp. 325-332.
Dukcnvicz, 1.K., 1980, "A Particle-Fluid Numerical Model for Liquid
Sprays", Journal of Computational Physics, Vol. 35, pp. 229 253.
El Banhawy, Y., and Whitelaw, J.H., 1980, "Calculation of the Flow
Properties of a Confined Kerosence-Spray Flame", AIAA Journal, Vol. 18,
pp. 1503-1510.
Faeth, G.M., 1983, "Evaporation and Combustion of Sprays",
Progress in Energy and Combustion Science, Vol. 9, pp. 1 76.
E1., and Smith, CE., 1993, "Integrated CFD Modeling of Gas
Turbine Combustors", AIAA Paper 93-2196.
Gordon, S., and McBride, 13.J., 1971, "Computer Program for
Calculation of Complex Equilibrium Composition, Rocket Performance,
Incident and Reflected Shocks and Chapman-Jouquet Detonations", NASA
Report SP-273.
Gosman, A.D., and loannides, E., 1983, "Aspect of Computer
Simulation of Liquid-Fueled Combustors", AIAA Paper 81-0323, 1981 and
Journal of Energy, Vol. 7, pp. 482-490.
Gosman, A.D., and Johns, R.J.R., 1979, "A Simple Method for
Generating Curvilinear Orthogonal Grids for Numerical Fluid Mechanics
Calculations", Internal Report FS/79/23, Department of Mechanical
Engineering, Imperial College, University of London.
Issa, R.I., 1986, "Solution of the Implicitly discretised Fluid Flow
Equations by Operator-Splitting", Journal of Computational Physics,
Vol. 62, pp. 40-65.
Kuo,T.W., and Reitz, R.D., 1989, "Computation of Premixed-Charge
Combustion in Pancake and Pent-Roof Engines", SAE Paper 890670.
Lai, K.Y.M., 1992, "Calculation of Turbulent Reactive Flows in
General Orthogonal Coordinates", Technical Report IME-CFE-TR-001
(NRC No. 32128), National Research Council Canada, Ottawa, Ontario,
Canada.
Lai, M.K.Y., and Cheney, P., 1995, "CFD Analysis of the Allison
570KF Gas Turbine Combustor", I I th Symposium on Industrial
Application of Gas Turbines, Banff, Alberta.
Launder, B.E., and Spalding. D.B., 1974, 'Me Numerical
Computation of Turbulent Flows", Computer Methods in Applied
Mechanics and Engineering, Vol. 3, pp. 269 289.
Magnussen, B.F., and Hjertager, B.H., 1976, "On Mathematical
Modelling of Turbulent Combustion with Special Emphasis on Soot
Formation and Combustion", Proceedings of the 16th International
Symposium on Combustion, The Combustion Institute, pp. 710-729.
Martino, P.Di., Colantuoni, S., Cirillo, L., and Cinque, G., 1994,
"CF) Modelling of an Advanced 1600K Reverse-Flow Combustor",
ASME Paper 94-GT-468.
-
-
-
-
-
-
-
10
Downloaded From: http://proceedings.asmedigitalcollection.asme.org/ on 12/29/2014 Terms of Use: http://asme.org/terms