# Calculations The Flow Inside Different Types of Can Gas Turbine ...

Mechanics

Feb 22, 2014 (4 years and 4 months ago)

189 views

Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



Calculations The Flow Inside Different Types of
Can Gas Turbine Combustion Chamber

Arkan K. Al
-
Taie

University of Technology

Abstract

The paper presents a mathematical model for three
-
dimensional, swirling, re
-
circulating turbulent
flows inside

a diff
erent types of
c
an

gas turbine combustion chamber
. The present model is restricted to
single
-
phase, diffusion
-
controlled combustion, with radiation heat transfer. The mathematical model
comprises differential equations for: Continuity, momentum, stagnation

heat transfer, turbulence energy, its dissipation rate, and the mean square of concentration fluctuations.

The simultaneous solution of these equations by means of a finite
-
volume solution algorithm yields
the values of

the variables at all internal grid
nodes. the prediction procedure, composed of the
mathematical model and its solution algorithm, is applied to predict the fields of variables within a
representative can combustor; the results are compared with some work
ers. The predicted results give the
same trends as the measured ones, but the quantitative agreement is not always acceptable; this is
attributed to the combustion process not being truly diffusion
-
controlled for the experimental conditions
investigated by

some workers.

ةصلاخلا

ت
مممض

لا
مممحب

مممسارد
ممم ن

مممضاير

اممميرجل

مممثلاث

،
ِ
داعبلأا
، ر،مممض

امممفتلا

،

ممم
دامممعإ

ري دمممتلا

مممخاد

ممم مممف تخ خا مممنا
ممم رغ
مميب علا اارممتحلاا
.
مم
لا

ن

ا
دختممس ل

اممب،ي

ممم

امميرجلا مملاح
داممحأ

اارمممتحلاا خ ممن ر مم،لا
راممشتنا

َ

،

ممخلأا

اممم تنا راممبتعلاا رممظنب
ب رارمممحلا
عمممشلإا

ِ
خا
.
إ

ممم نلا

مممضايرلا

مممح مممضت

لادامممع لا

ممميتلاا

:
يرار تمممسلاا
،

ممملا
، خز

مممي ةلا ممميبلاثنلاا
،

لا
،زممميةرت

اممم تنا

لا

ِ
رارمممح

خاعشلإاب
،

ِ
يبار،ضلاا يةرحلا قا،لا
،

تشتلا دع
،

ِ
زيةارتلا ب ت دع ،س ت
.

كل ددح لا جحلا ير، ادختسا ت ، ر ة لا لاداع لا حل

داجيلإ

اارمتحلاا م رغل مي خادلا ،ام نلا دمنع مب ، لا اريغت لا يق
، يب علا
نرا ت ث

م ع اماي متلا حبلا ضعب جئاتنلا

مناة

جئامتنلا

يمل مبرا
ملإ

بمسب ا ما ابام،تلا دمح
إ

اارمتحلاا
سحب كل ايراشتنا اي ي ح ةي ل

هي ع ن ضت ا

بل ي علا حبلا
. يثحابلا ضع

List of Symbols

Latin Symbols

Symbol

Description

Unit

C
1
,C
2

Constant in Turbulence Model

-

C
D

Constant in Turbulence Model

-

C
R

Constant in Eddy
-
Break
-
Up Model

-

C
P

Specific Heat at Constant Pressure

J/Kg.K

C

Constant in Turbulence Mod
el

-

E
w

Black Body Emissive Power at the Wall
Temperature

W/m
2

f

Mixture Fraction

-

g

Time Mean Square of the Concentration
Fluctuation Defined by g=

2
f

-

G
k

Quantity in the Generation Term for k

-

h

Specific Enthalpy

J/Kg

h
o

S
tagnation Enthalpy

J/Kg

J

Radiation Flux in the Direction of Negative r

W/m
2

H
fu

Heat of Combustion (Net Calorific Value of Fuel)

J/Kgmole

J

Radiation Flux in the Direction of Negative r

W/m
2

K

Turbulence Kinetic Energy

m
2
/sec
2

K

Direction of Positive z

W/m
2

K
m

Thermal Conductivity for Mixture

W/m.K

L

Length

m

L

Radiation Flux in The Direction of Negative z

-



m
j

Mass Fraction of Mixture Species

-

m
fu

Mass Fraction of Fuel in The Inlet Fuel

-

m
ox

Mass Fraction of Oxygen in Th
e Inlet Oxidant
Stream

-

M

Radiation Flux in The Direction of Positive

W/m
2

N

Radiation flux in the direction of negative

W/m
2

P

Correct Pressure

N/m
2

P*

Guessed Pressure

N/m
2

P

Static Pressure

N/m
2

P
o

Stagnation Pressure

N/m
2

Pr

Prandtl Number

-

R

Gas Constant

J/Kg.K

R

Universal Gas Constant=8314.3

J/kgmole.k

R
x,

R
y
, R
z

A dependent Variable for Radiation Fluxes in
The Three Coordinate Directions

-

s

Mass Fraction of Stoichometric Oxidant/Fuel
Proportions

-

S
c

Scatteri

-

T

Temperature

K

Greek Symbols

Symbol

Description

Unit

The General Dependent Variable

-

Effective Exchange Coefficient

Viscosity

Kg/m.sec

t

Turbulent Viscosity

Kg/m.sec

eff

Effective Viscosity

Kg/m.sec

s
m
m
ox
fu

-

ox

s
m
ox

-

fu

fu
m

-

Stefan
-
Boltzman Constant

W/m2.K4

T,t

Turbulent Prandtl Number

-

T,l

Laminar Prandtl Number

-

k

Schmidt Number for k=0.9

-

Schmidt Number for

as in Equation (3
.21)

-

f

Schmidt Number for f=0.7

-

h

Effective Prandtl Number for Enthalpy=0.7

-

1
-
Introduction

The first multi
-
dimensional codes appeared in the 1970s, and aimed at the
solution of the full Navier
-
Stockes equations including finite rate chemistry.

Among
the earliest papers, some appeared in the first AGARD symposia on combustion and
fuels in gas turbines (AGARD, 1980). In subsequent meetings (AGARD 1987, 1993)
improvements were made on these methods and new ones appeared. A full review of
these tec
hniques will not be attempted here. Instead the general philosophy behind
these methods and their strengths and limitations will be highlighted.

The turbulence model most frequently used in combustion application is the k
-

model (Kuo, 1986; Schetz, 1993).

It relates the turbulent fluxes to the mean velocity
gradients through a "turbulent" viscosity (Boussinesq assumption), which in turn is
related to the kinetic energy of turbulence k and a turbulent dissipation rate

. These
two variables are obtained fro
m the solution of transport equations derived from the
Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



original NS equations, and which are also partial differential equations. This increases
the complexity of the original numerical problem, but allows the "history" of the flow
to be accounted for the t
urbulent viscosity.

Several researchers pointed some of the inadequacies of these models to capture
the quantitative aspects of flows within combustors, among other reason for their
dependence on constants obtained from less complex flows (Maidhof and Jani
cka,
1993), and because the standard k
-

assumes isotropic turbulence which may not be
true in recirculating or swirling flows (Hwang, et al., 1993). The Boussinesq
assumption is frequently cast into doubt in combustion and other applications (Correa
and S
hyy, 1987; Schetz, 1993).

Turbulence and chemistry modeling are in reality coupled because reaction rates

are

highly nonlinear functions of the instantaneous values of the main variables, not
just their mean values. Therefore the mean reaction
-
rates (that

appear in the Reynolds
equations) are not the reaction rates evaluated with the mean variables. One
possibility is to evaluate the production rate by both an Arrhenius expression based on
mean values and an eddy break
-
up model to account for the effects o
f turbulence, and
use the smaller of the two (Khalil, Spalding and Whitelaw, 1975; Nikjooy and Peck,
1988). This approach is based more on intuition than on first principles, but they have
enjoyed some success and some argue this due to the assumption that

the combustion
rate is proportional to the inverse of the turbulent time
-
scale (Gran, et al., 1993).
Another approach is the use of Probability
-
Density Functions (PDFs) (Spalding, 1976;
Kahlil, 1982; Pope, 1991). They are being used in practical applicati
ons with more
frequency (Topaldi et al., 1996), usually in conjunction with the fast
-
chemistry model
(Shyy, et al., 1988; Alizadeh and Moss, 1993; Baron, et al., 1994; Di Martino et al.,
1994), but in all cases it requires greater computational effort.

It
is relatively common to assume that the flow is single
-
phase or locally
homogenous, which implies fast fuel
-
droplet evaporation or equivalently infinitely
-
small droplets (Serag
-
Eldin and Spalding, 1979). Some researches do not consider this
assumption to b
e appropriate in practical combustor applications, and therefore take
into account two
-
phase interaction. This means writing additional equations for liquid
phase and coupling them with the Reynolds equations for gas phase. Two approaches
are common, depen
ding on whether the discrete particles are followed along their
trajectories or the liquid phase treated as an interacting continuum: Lagrangian and
Eulerian (Swithenbank, et al., 1980; Sirignano, 1993; Hallmann, et al., 1993). Also
the different methods c
an be classified as deterministic or stochastic, depending on
whether the effects of turbulence are neglected or accounted for. Whatever approach
is used the complexity of the solution procedure is considerably increased. In addition
a fair amount of model
ization (and uncertainty) is introduced regarding droplet
-
size
distribution, droplet drag, heat
-

and mass
-
transfer coefficients, temperature within
droplet, etc.

As these models become more complex in an attempt to improve the accuracy of
the equations in
which they are based, greater demands are placed on the numerical
methods needed to solve them. The vast majority of applications are based on the
SIMPLE methodology introduced in the 1970s by Spalding and Patanker (Patanker,
1980; Kahlil, 1982). It is bas
ed on the sequential solution of the governing equations
together with a velocity
-
pressure correction. This approach offers a "simple" way to
add transport equations for variables as needed. This method was among the first to be
used in combustion applicat
ions, and still remains the most widely used (most papers
cited so far are of this type). Other researchers complain of poor convergence due to
the sequential approach and the explicit treatment of the chemical source
-
terms
(Chen,



1987
)
. Therefore they have

started using strongly
-
coupled and fully
-
im
p
licit schemes,
popular in non
-
combustion applications (Su, A., Chang,

J.
-
H., Ho, W.C., 2003). The
downside of these methods is given by their complexity and greater computer
-
storage
requirements.

In brief, the m
odeling effort necessary to solve the k
-

equations and the
demands it places on the numerical methods imposes some limitations on their
practical use:

The intensive computational effort required to solve the governing equations
puts a definite limit on th
e size of the physical domain that can be solved in a
given computational run. If the solution is required for a complete
configuration, it usually proceeds in two successive steps (McGurik and
Spencer, 1996). First, external flows into and around the prim
ary flow are
calculated. This may

be

done with CFD methods (conditions inside the
primary have to be assumed) (Karki, et al., 1992), or with one
-
dimensional
codes, sometimes they may be obtained from experimental data). The next
step is to calculate the in
ternal flow with previous mass
-
splits as inputs (no
attempts are apparently made to iterate between both calculations). In other
applications where more detailed information is required within the primary
-
flow (i.e., finer grids) a sequential multi
-
block a
pproach is used in that part of
the combustor (Crocker and Smith, 1993).

Almost all of the available solutions are for steady
-
state conditions. Even
-
state.
Transient simulations could be used
to predict weak extinction limits (Allan,
1996) or responses to variations in air and/or fuel flow
-
rates. Not surprisingly,
they are rarely attempted with conventional CFD methods.

CFD has become an indispensable tool in the design of modern gas
-
turbine
co
mbustors. Several researchers report fair to good quantitative results (Crocker and
Smith, 1993; Lawson, 1993; Baron, et al., 1994; Di Martino, et al., 1994). As
mentioned before they have been used successively in design applications to improve
existing d
esigns (Crocker and Smith, 1993; Lawson, 1994) or to predict the
performance of completely new combustors (Danis, et al., 1996). The main point here
is that these methods just have not reached the level of accuracy and ease of use that
would allow them to
be the sole analytical tool. Therefore in the design process there
is still a place for lower
-
order methods, based on basic principles and some empirical
information that can tackle problems not currently handled performance of a
combustor before applying
more sophisticated numerical tools or resorting to full
-
scale tests.

Thus, with the need
s

of combustion chamber designer
s

in mind, the work
described in this paper has been devoted to the presentation application, and
validation of a mathematical model fo
r can combustors. The model simulates the
actual physical processes by means of differential equations for the dependent
variables. The simultaneous solution of these equations by means of a finite
-
difference procedure yield the values of the dependent var
iables at all internal grid
nodes. From the computed values of the dependent variables, and with the help of
algebraic
, relations, the values of the auxiliary variables are derived.

The model

presented here

is applicable to three dimensional,
single phas
e,
turbulent,
diffusion
-
controlled combustion,
swirling, recirculating flows
, with

Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



2
-

Governing Equation for Aero
dynamic and Thermal Predictions:

The finite difference procedure employed in this work is derived from the
pioneeri
ng work of Professor Spalding and Co
-
workers at Imperial College. It uses
the cylindrical polar system of co
-
ordinates to predict the complex three dimensional
swirling and reacting flow inside a combustor. The numerical combustion scheme for
solving the g
overning non
-
linear equations with arbitrary boundary conditions
coupled with mathematical models of turbulence and combustion has the pressure and
velocities as the main flow variables. The model is restricted to diffusion
-
controlled

effects are incorporated by means of a six flux model
(Swithenbank,
et al.,

1980).

A)
Equations of Continuity and Momentum
:

0
u
x
w
r
1
rv
r
r
1

…(1)

X
-

x
u
2
x
x
w
r
u
r
1
x
v
r
u
r
r
r
1
x
p
x
uu
wu
r
1
r
rvu
r
1

…(2)

R
-

x
v
r
u
x
r
v
w
r
1
r
2
v
r
1
r
w
r
r
r
1
r
v
r
2
r
r
1
r
p
r
w
x
uv
wv
r
1
r
rvv
r
1
2

…(3)

-
momentum

3. Density and Viscosity:
-

In general, the density is related to pressure, temperature and the composition of
the gas mixture through an equation of state:

=

⡐Ⱐ听⁭
j
)

…(5)

W桥牥 吠楳i 瑨攠瑥浰t牡瑵牥 潦o 瑨攠晬畩搠浩x瑵牥 a湤n 瑨攠m
j
's are the mass
fractions of the component species of the mixture. Variations of T and m
j
's either form
part of the problem specification or, in general case are obtained by the solution of
additional differential equations. In the latter case, it may be
necessary to introduce a
reaction model which allows the appropriate to be setup.

For the evaluation of viscosity, laminar and turbulent flows have to be
considered separately. For laminar flows, the viscosity is assumed to be a function of
temperature and

mixture composition:

=

⡔Ⱐm
j
)

…(6)

u
r
1
x
w
x
r
v
2
w
r
2
r
1
v
r
1
r
w
r
r
r
v
r
1
r
r
w
r
r
r
r
1
p
r
1
r
vw
x
uw
ww
r
1
r
rvw
r
1

…(4)



For turbulent flows, however, the problem is more complicated. The turbulent
contribution

t

to the effective viscosity is a function of local quantities such as
velocity gradients, etc. The evaluation of

t

and the t
urbulence model are discussed in
the next section.

4. The Turbulence Model and Effective Viscosity:
-

4.1. k
-

潤e氺
-

Several models of turbulence have been put forward by different authors. These
models differ in complexity and range of applicability; the
y also involve the solution
of different numbers of differential equations. The turbulence model incorporated in
this work is the high Reynolds number k
-

two equation model. This model requires
the solution of two differential equations, for the two turbu
lence properties: the kinetic
energy of turbulence k, and its dissipation rate

. The model is moderate in
complexity. It has been extensively used by many investigations and has proved to be
adequate over a wide range of flow situation. Here the governing

differential
equations are presented below in details (Launder,
and

Spalding, 1974).

Differential equations for turbulence
-
energy k and dissipation rate

used in
combustion are respectively as follows:
-



k
k
k
k
G
x
k
x
r
k
r
r
1
r
k
r
r
r
1
x
uk
wk
r
1
r
rvk
r
1

…(7)

k
2
2
C
k
G
k
1
C
x
x
r
r
1
r
r
r
r
1
x
u
w
r
1
r
rv
r
1


…(
8
)

4.2. The Effective Viscosities:
-

The turbulent viscosity

t

is related to k and

via:
-

t
=C
D

k
2
/

…(9)

䅮搬⁴桥⁥晦ec瑩癥 癩獣潳楴y 楳⁧楶敮i批:
-

eff
=

t
+

l

…(10)

W桥牥

l

is the laminar or molecular viscosity of the fluid. Often

t

is very

large
compared with

l

and,

eff

can be taken equal to

t

without introducing errors
(Spalding, D.B., 1976).

Recommended values for the constants appearing in the above equations are
from (Spalding, 1976; Versteeg

and

Malalasekera, 1995) as follows.

C
1
=1.
43

C
2
=1.92

C
D
=0.09

k
=0.9

…(11)

for the dissipation rate of turbulence is calculated from:
-

D
1
2
2
C
C
C

…(12)

W桥牥

is the von
-
karman constant. A value of

0.4 is assumed for

in this work
(Serag
-
Eldin,
and
, Spalding, 1979).

4.3. Length Scale:
-

A length scale of turbulence, l, can be defined as (Jones,
and
, Launder, 1972,
Versteege

and

Malalasekera, 1995):
-

Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



2
3
k
l

…(13)

Since the length s
cale is a more interpreted quantity than, for example, the
dissipation rate, its use with turbulence energy, allows the nature of turbulent flows to

5. Reaction Models:
-

Conservation equation for a chemical species J is defined

as (Khalil,
et al.,
1975,
Serag
-
Eldin,
and
Spalding, 1979, Versteeg,
and

Malalasekera, 1995):
-

x
j
m
j
x
r
j
m
j
r
1
r
j
m
j
r
r
r
1
j
R
j
um
x
j
wm
r
1
j
rvm
r
r
1

…(14)

W桥牥 m
j
is the mass fraction of the chemical species J, R
j

the mass rate of
creation of species J by chemical reaction per unit
volume and

j
is the exchange
coefficient.

In a multi
-
component system simplification can be introduced through the use of
the concept of a simple chemically reacting system. This assumes that fuel and
oxidant react chemically in a unique proportion, furth
ermore, the effective
diffusivities of all the chemical species are taken to be equal and the reaction is a
single step with no intermediate compounds or two
-
step.

These suppositions enable various concentrations to be determined by solving
equation for th
e variable f, defined as the mass fraction of fuel in any form.

The mixture fraction, f, is defined as:

ox
fu
ox
f

…(15)

inlet
ox
ox
s
m

…(16)

inlet
fu
fu
m

,

W桥牥:
-

s
-

inlet
fu
m
a湤n

inlet
ox
m

a牥 瑨攠浡獳m晲ac瑩潮猠潦o晵f氠a湤n潸yge渠楮i

-

e煵q瑩潮⁷楬o⁢ ⁵獥搠a牥:

x
f
x
r
f
r
1
r
f
r
r
r
1
wf
r
1
vf
r
r
r
1
uf
x
f
f
f

…(17)

W桥牥:
-

f

-

f
eff
f

…(18)

䅮搠

f

is the Schmidt number for mixture fraction is assumed to be constant
and equal to 0.7 (Serag
-
Eldin,
and

Spalding,

1979).



6. Combustion Models:
-

The
re are three combustion models described below (Khalil,
et al.,
1975, Serag
-
Eldin,
and
Spalding, 1979):
-

Model 1:

the first model postulates a physically controlled, one
-
step reaction, with
fuel and oxygen unable to coexist at the same location. The only s
pecies equation to
be solved is that for the mixture fraction f, this equation has no source.

Model 2:

the combustion model which will be employed here is a relatively simple
one, which assumes "instantaneous reaction" with allowance for concentration
fluc
tuations.

The instantaneous
-
reaction assumption implies that whenever fuel and oxidant exist
together at a point, chemical reaction will proceed instantaneously to completion in a
single step, and producing complete combustion product. As a consequence of
this
assumption, fuel and oxidant may not be present in the same place at the same time.
However, because the flow is turbulent, the concentrations of fuel and oxygen
fluctuate with time so that at the same point, fuel and oxygen may be present
separately,

at different times. Thus, the local fluctuation of the concentration will
affect the local time
-
mean value of the unburned fuel mass fractions, and hence the
rate of reaction. Consequently, an estimate of the magnitude of the fluctuation and its
effect on

the rate of reaction is required. The approach adopted here is that of
reference (Spalding,D.B.,197
6
), where a single transport equation is solved for the
time
-
mean of the square of the concentration fluctuation; and assumptions regarding
the shape of the

instantaneous concentration
-
time profile, a hypothetical distribution
of the instantaneous concentration with time is derived. In this work a "battlement
-
shaped" profile, (Serag
-
Eldin, M.A., Spalding, D.B., 1979) is assumed.

Model 3:

in contrast to models

1and 2, a finite reaction rate is introduced and the
chemical process is not physically controlled (as, for example, in the case of premixed
fuel and oxidant mixtures). The rate of the chemical reaction will be the influencing
factor, solution of a source
-
free equation alone is not sufficient to determine all the
species mass fractions. In addition, a conservation equation for one of the species
must be solved. This could be any one of the three: m
fu,,

m
ox
or m
pr
, the computation
of the source of fuel for

kinetically influence reactions depending on types of flow
Laminar or Turbulent.

a)

Laminar

For Laminar flow, a form of the bi
-
molecular Arrhenius expression for the
reaction rate (Versteeg, H.K., Malalasekera, W., 1995):
-

T
R
E
exp
m
m
T
z
R
ox
fu
2
l
,
fu
2
1

…(19)

W桥

-

z⁩猠 桥⁰牥
-
ex灯湥湴na氠晡c瑯爠t湤n䔠楳⁴桥⁡c瑩癡瑩潮⁥湥牧y.

Turbulent

The rate of reaction which depend on Eddy Break
-
Up model of Spalding
(Spalding, 1976, Magnussen,
and
Hjertager , 1976 )

k
g
C
R
2
1
R
EBU
,
fu

…(20)

W桥牥:
-

C
R

is a constant
and equal to 1.5, g represents the local mean square
concentration fluctuation (Magnussen,
and

Hjertager , 1976).

Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



7. The Stagnation
-
Enthalpy Equation:
-

The stagnation
-
enthalpy, h
o
, is defined as:

2
2
2
o
w
v
u
2
1
h
h

…(21)

fu
fu
p
H
m
T
C
h




W桥牥吠楳i瑨攠瑥浰t牡瑵牥Ⱐ
fu
is the heat of combustion, and C
p

is the specific heat
at constant pressure of the gas mixture. The specific heat C
p

is a function of
temperature and mixture composition. For incompressible flows, the kinetic energ
y of
the mean motion is usually small compared with the specific enthalpy, and can be
ignored. In these cases h
o

is identical to h

The stagnation
-
enthalpy is given by (Serag Eldin,
and
Spalding, 1979):
-

h
o
h
o
h
o
h
o
o
o
S
x
h
x
r
h
r
1
r
h
r
r
r
1
uh
x
wh
r
1
rvh
r
r
1

…(

)

I渠瑨攠a扯癥 e煵q瑩潮o

S
h
, represent the sum of all the source terms for h
o
.
Depending on the particular application, S
h

may include: source term due to radiation,
heat transfer to wall boundaries, and heat exchange with the second phase.

While

h

is
the exchange coefficient an
d equal to

eff
/

h
,

h

is the effective Prandtl number for
stagnation enthalpy is assumed to be constant and equal to 0.7 (Serag Eldin,
and
Spalding, 1979).

-

The effects of radiation in the mathematical model are accounted for, by ref
erence
to the six
-
flux model of radiation. The differential equations describing the variations
of the fluxes are (Swithenbank, J., et al., 1980):
-

N
M
L
K
J
I
6
S
aE
r
J
I
S
a
r
rI
dr
d
c
c

…(

)

N
M
L
K
J
I
6
S
aE
r
J
J
S
a
r
rJ
dr
d
c
c

…(

)

N
M
L
K
J
I
6
S
aE
K
S
a
K
dx
d
c
c

…(

)

N
M
L
K
J
I
6
S
aE
L
S
a
L
dx
d
c
c

…(

)

N
M
L
K
J
I
6
S
aE
M
S
a
M
rd
d
c
c

…(

)

N
M
L
K
J
I
6
S
aE
N
S
a
N
rd
d
c
c

…(

)

W桥牥:
-

††
I
-

††
J
-

††
K
-

††
L
-

††
M
-

.

N
-
radiation flux in the direction of positive

.



a
-

absorption coefficient.

S
c
-

Scattering coefficient.

E=

T
4
, black body emissive power at

the fluid temperature.

-

Stefan
-
Boltzmann constant=5.669E
-
08 W/m
2
.K
4
.

8.1. Composite
-
Flux Equations:
-

The composite
-
fluxes R
y
, R
x
and R
z

are defined as (Swithenbank, et al., 1980):
-

)
J
I
(
2
1
R
y

…(

)

)
L
K
(
2
1
Rx

…(
31
)

)
N
M
(
2
1
R
Z

…(
32
)

They are employed to eliminate I,J,K,L,M, and N from the previous equations to
yield three second
-
order ordinary differential equations

Z
Y
X
c
x
x
c
x
c
x
c
R
R
R
2
3
S
E
R
a
r
R
S
a
1
r
R
S
a
r
r
r
1
)
dx
dR
S
a
1
(
dr
d

…(

)

r
R
r
1
S
a
1
r
x
R
r
1
S
a
1
x
)
dr
dR
r
1
c
S
a
1
r
(
dr
d
r
1
Y
c
Y
c
Y

Z
X
Y
Y
R
R
R
2
3
c
S
E
R
a

…(

)

r
Z
R
c
S
a
r
r
r
x
Z
R
c
S
a
1
x
r
Z
R
c
S
a
1
r
1

Y
R
X
R
Z
R
2
3
c
S
E
x
R
a

…(

)

-

-
e湴桡汰n e煵q瑩潮o

X
, R
Y

and

R
Z
in the source term of the latter.

The contribution of the radiation flux to the stagnation
-
enthalpy source term is given
by

E
3
R
R
R
a
2
S
Z
Y
X
h

…(

)

9. Solution Procedure:
-

The previous governing partial differential equations for mass, momentum,
energy, and species equation are elliptic in nature and as mentioned before, can be
convenient
ly put in general form:

S
x
x
r
r
1
r
r
r
r
1
w
r
1
v
r
r
r
1
u
x

…(37)

,

, and S

. as explained above for each equation

The equations are first reduced to finite difference equations by integrating
over finite control volumes and then solved by a procedure de
scribed in details in
Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



previous section for both orthogonal and non
-
orthogonal, two
-

three dimensional
equation.

It will be sufficient for the purpose of this section to summarize the main feature of
the solution procedure.

The numerical scheme follows the
well known SIMPLE pattern and can be
summarized as follows:

Iterative one which starts from given initial conditions for all variables and converges
to the correct solution on the completion of a number of iterations.

Each iteration performs the following
steps:

1
-

the momentum equations are discretized and linearized, leading to an algebraic
equation system for each velocity component, in which the pressure, the other
velocity components, and all fluid properties are treated as known (values
from the previous

iteration are used).

2
-

The improved velocity field and the prevailing density are used to calculate
new mass fluxes through the control volume faces and to invoke the mass
-
conservation equation; the result is the pressure
-
correction equation. Upon
solving i
t for p

, the mass fluxes, control volume
-
center velocities, density, and
pressure are corrected using the expressions given above. Only an

-
fraction
of p

(0.1
-
0.8); smaller values is used for highly
non
-
orthogonal grids.

3
-

ditional transport equations are solved in the same manner to obtain better
estimates (e.g. for energy, turbulence quantities, species concentration, and

4
-

If necessary, fluid properties are updated (e.g. the density, viscosity, and
Pran
dtl number may be a function of temperature or turbulence quantities).

The above steps are repeated until the residual level in each equation becomes
sufficiently small. Usually a reduction by three to four orders of magnitude
corresponds to a convergence
error of the order of 0.1%.

10
-
Geometry of the Combustion Chamber in the Present Work:

The geometry of the combustion chamber is displayed in Figure(1)
.

T
he
combustion chamber walls are cylindrical, having a 0.21 m internal diameter and a 2.0
m length. T
he inlet wall supports to coaxial
burner, which
comprises

an annular vane
-
swirler of 0.0781 m outer diameter and 0.0425 m inner diameter, surrounding a
0.0195 m diameter fuel tube. The swirler includes ten vanes forming 45
o

angles with
the burner axis is e
mployed to swirl the primary air before entering the chamber .

Six radial pipes, each of 0.0254 m diameter, carry the dilution air into the
chamber at a distance of 0.24 m downstream of the inlet wall. The pipes are spaced 60
deg apart. Since the dilution

air is introduced into chamber at discrete positions on the
circumference of the cylindrical walls, it results in an asymmetric flow pattern and
hence dictates the three dimensional treatment of the problem.

The flow of the
dilution air through each of th
e six ports is indepen
den
tly measured and controlled, so
as to insure equal distribution of the air among the ports.

The cooling of all chamber
walls is achieved by means of separately controlled air (other than the dilution or
primary air).

The experimen
tal conditions adopted are tabulated in Table (1)

(Serag
-
Eldin and
Spalding, 1979). Three conditions were investigated designated: Flame A, Flame B,
Flame C. Flame B exhibits the same primary air and fuel flow rates as Flame A, but
with a larger dilution a
ir flow rate, whereas Flame C exhibits approximately the same
ratio of primary air to fuel flow rates as Flames A and B, but with gr
e
ater primary air



and fuel rates, and smaller dilution air flow rates. The fuel used was
natural gas
.

While the
experimental

conditions for the second model

(McGurik

and
Palma,1993)
are tabulated in Table (2)

and shown in Figure (2)
.

Table (1)
:

Experimental Conditions (Serag
-
Eldin and Spalding, 1979)

Item

Flame A

Flame B

Flame C

Primary air flow
rate (kg/s)

0.039

0.039

0.0
562

Dilution air flow
rate (kg/s)

0.041

0.061

0.0352

Fuel flow rate
(kg/s)

0.00155

0.00155

0.00219

Overall fuel to air
ratio (kg/kg)

0.0194

0.0155

0.0240

Fuel to primary air
ratio (kg/kg)

0.0397

0.0397

0.039

Mean wall

temperature

553

556

583

Table
(2):Experimental Conditions(McGurik

and
Palma,1993)

Combustor with swirler

Swirler airflow
rate(kg/sec)

0.543(17%)

Primary air flow
rate(kg/sec)

1.04(33%)

Dilution air flow
rate(kg/sec)

1.6(50%)

Fuel flow rate(kg/sec)

0.083

Combustor without swirler

Swirler air flow
rate(kg/ec)

-

Primary air flow
rate(kg/sec)

1.04(40%)

Dilution air flow
rate(kg/sec)

1.60(60%)

Fuel flow rate(kg/sec)

0.083

Length of combustion chamber=0.21 m

Diameter of combustion chamber=0.074 m

Primary hole diameter=10 mm(6 holes
)

Dilution hole diameter=20 mm(6 holes)

Swirler number=0.45

Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



11
-
Results and Discussion:

The finite difference grid employed here exhibits (30*22*12) grid nodes in the
x,r, and

direction
s, respectively. The distribution of the grid nodes is revealed in
Figure (
3
).

The integration domain extends only 0.73 m in the downstream direction, since
this length was considered sufficient to give fairly well developed profiles at the
downstream end
of the integration domain. Since the flow pattern in the three
dimensions repeats itself in each 60 deg sector of the combustor cross section, the
integration domain covers only one 60 deg sector.

From Figure (3) it is seen that the grid nodes locations ex
tended beyond the B
1

and B
2

-
boundaries of the integration domain; this is because the

-
boundaries
Figure(1):Grid
G
e
ne
ration for
T
hree

D
imensional
C
ombustion
C
hamber

(Type (1))
.

Figure(2):Grid
G
eneration for
T
hree
D
imensional
C
ombustion
C
hamber

(Type(2))
.



conditions are cyclic ones, i.e. the given boundary condition is that all the

-
planes
separated by multiples of 60 deg are corresponding planes (exhibit i
dentical profiles).
Thus, in order to implement these boundary conditions implicitly in the solution
procedure, the

-
boundary values are stored at

-
planes outside the integration
domain, which are so oriented with respect to B
1

and B
2

as to "correspond"
to internal
solution planes.

Figure(3):
G
rid Generation for
T
hree
D
imensional
C
ombustion
C
hamber

(Type(1))
.

Since the computations yield predicted three
-
dimensional distributions for all
the dependent and auxiliary variables,
it is practically impossible to display all the
predicted results in here. Thus only some selected profiles are displayed and discussed
for convenience.

Figures (4
-
5
) display predicted results for flame A. Figure (4) shows the velocity
vectors in the

-
pla
ne passing through one of the dilution air ports. In one of the two
-
diagrams, the vectors represent the velocity component in the

-
plane, in both
magnitude and direction, whereas in the other diagram the vectors are all of equal
length and indicate the ve
locity component's direction only. The two diagrams are
required in order to indicate both the velocity magnitude and direction, in regions
where the velocity magnitudes are small. Figure(5) displays the predicted temperature
contour lines in the same

-
pl
ane as the
velocity
vector diagram.

Figures (
6
-
7
) show the predicted results from flame B that correspond to those
of Figures (
4
-
5
) for flame A, whereas Figures (
8
-
9
) show the ones for flame C.

Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



Comparison of the three sets of predictions reveals interestin
g information
regarding the effect of decreasing the proportion of dilution air on the profile shapes.
For example, when proportion of dilution air decreases from its value for flame B to
its value for flame A, and then even further to its value for flame
C.

i)

The dilution air jet penetrates radially less into the chamber.

ii)

The proportion of dilution air recirculating upstream decreases until finally
all the dilution air flow forwards and the central recirculation zone
diminish.

iii)

Unburned fuel exists for a long
er distance downstream.

iv)

The maximum temperature in the field rises from 1784 K to 1972 K, and
its position is shifted further downstream.

v)

The circumferential non
-
uniformities in the temperature profiles at section
downstream of the dilution air ports.

vi)

The
radial temperature distribution downstream of the dilution air ports
and in the dilution air

-
plane changes from one where the minimum
temperature is close to the center line and the maximum temperature is
near to the wall to one where the maximum tempera
ture is near to the
center line and the minimum temperature is close to the wall.

Figure (13) display the corresponding comparison between the predicted and that
measured temperatures by (Serag
-
Eldin, Spalding, D.B., 1979) for the same model
using above bu
t the fuel used is 100%

methane while for model used in the above
research of fuel composition is

(93.63%Methane,3.25%Ethane,1.78%Nitrogen,0.69%
Propane, and others 0.65%).

The cause of disagreement of the results for flame A and C with the measured
values

is the neglect of the effect of chemical kinetics on the rate of reaction. Indeed,
referring to Figure (
10
) for example, it is noticed that in the region near the inlet wall
and close to the center line, the fuel is predicted to burn under condition of ri
ch fuel
mixtures and high rates of mixing, for which case the rate of reaction may well not be
diffusion
-
controlled. Thus, if this model is to be applied to similar conditions then the
effect of chemical kinetics has to be incorporated.



Figure(4a):Velocity Vectors Showing Direction only,

for Flame A.

Figure(4b):Velocity Vectors Showing Magnitude and Direction ,for
Flame A.

Figure(5):Temperature Contour Lines in Dilution air

-
Plane for Flame A.

Figure(
6
a
):Velocity Vectors Showing Direction only, for Flame B.

Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



Figure(
6
b):Velocity Vectors Showing Magni
tude and Direction ,for Flame B.

Figure(
7
):Temperature Contour Lines in Dilution air

-
Plane for Flame B.

Figure(
8
a):Velocity Vectors Showing Direction only, for Flame C.

Figure(
8
b):Velocity Vectors Showing Magnitud
e and Direction, for Flame C.



Figure(
9
):Temperature Contour Lines in Dilution air

-
Plane for Flame C.

Figu
re(
10
):Temperature Contour for Flame A in Downstream Planes.

Figure(
11
):Temperature Contour for Flame B in Downstream
Planes.

Figure(12):Temperature Contour for Flame C in Downstream
Planes.

Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



The second model for combustion chamber was studi
ed by (Heitor, et al.
,
1986
;

Bicen et al. 1990
;

McGuirk,
and
Palma, 1993), but take into account heat transfer to
the wall with pressure equal to the pressure at the end of the diffuser as shown in
Figure (
2
).

The calculations were restricted to a 60 deg se
ctor of the flow inside the
combustor for the same reason mentioned above. A numerical mesh of 60*32*22
. T
he grid spacing in the axial and radial
directions was changed smoothly to minimize the deterioration of t
he formal accuracy
of the finite difference scheme due to variable grid spacing and in such a way that a
higher concentration of nodes occurred near holes and the swirler. In the azimuthal
direction, the grid was regular except around the primary hole axis

and on the edges
of the domain.

The calculations were performed for two different flow conditions (Table (
2
)).
The discussion of results is spilt into calculations of combustor flow with and without
swirler.

Combustor Flow with Swirl:

The gas temperature

contours for the two longitudinal and the two transverse
planes are illustrated in Figure (
1
4
) and Figure (
1
5
), respectively. The fuel rich
mixture prevailing in the primary zone of the combustor results in relatively low gas
temperature in this zone. The

subsequent admission of air downstream brings the
mixture closer to the stoichiometric conditions resulting in hotter regions, in particular
near the wall. The information indicates the need to include extra holes (smaller)
between the primary and dilutio
n holes to handle the heat transmitted to the wall in
various zones. The transverse section at the primary jet plane shows higher
Figure(13):Comparison between the Predicted Temperatu
re and M
easured by
(Serag
-
Eldin and
Spalding, 1979)



temperatures to the sides of the sector than those in the core. At the intermediate jet
plane, regions of hot gases exist near

both walls.

Figure(14):Temperature Contour Lines View in the Longitudinal Planes
Containing the Axes of the Holes.

Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



Figure (15): Temperature Contour Lines View in the Cross Sectional

Planes Containing the Axes of the Primary and Dilution Holes.

Figures (
1
6
-
1
7
) show the velocity vectors for can combustor at two pla
nes: two
longitudinal at cross
-
sectional planes; that contains, the primary jet and dilution holes
respectively. The main features of the flow field, such as flow reversal in the primary
zone, the jet penetration, and the flow acceleration toward the exit
of the combustor,
are all illustrated in Figure (
1
6
-
1
7
).

The velocity vector plots in Figures (
1
6
) and (
1
7
) give the global picture of the
mean flow pattern inside the combustor under the present conditions (Table (
2
)). The
salient features are: i) Figure
(
1
6
), the toroidal vortex in the combustor head (caused
by the joint effect of the primary jet impingement and shearing, at the upstream side
of the jet, between the jet and the cross flow) and ii) in Figure (
1
7
b
), the two contra
rotating vortices in the c
ross section containing the dilution holes axes caused by span
wise deflection. These vortices are typical of jets in cross flow geometrics and
probably do not appear at the planes of the primary jets Figure (
1
7
a
) because of the
strong jet impingement.



Figure (17):Velocity Vector View in the Cross Sectional Planes that Containing

Primary and Dilution Holes.

Figure(16):Velocity Vector View in the Longitudinal Planes Containing
the
Primary and Dilution Holes.

Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



Figure (
1
8
) shows the resulting temperature distribution. Cold air occupies most
of the combustor, including the

annular paths and the primary region up to the dome.
As the fuel
-
air ratio immediately following is on the rich side, the maximum
temperature is reached only after air from dilution brings the mixture closer to
stoichiometry (or alternatively burns the ex
cess fuel). As more air comes into the
primary region, the temperature goes down. The resulting mixture finally leaves at a
temperature roughly intermediate between inlet and maximum.

Combustor Flow without Swirl:

The calculations were made with n
o flow through the swirler (both the azimuthal and
axial velocities were set to zero).

Figure (
19
) show the velocity vector plots for the flow conditions mentioned
above there is a single vortex in the head a recirculating region is larger than in case
of
swirling flow. In the intermediate region there is another toroidal vortex which is
eliminated if there is flow through the swirler. This is simply because the swirler flow
increases the axial momentum near the walls, downstream of the primary jets.

The co
mbustion flow with swirl also shows the correct shape of the turbulence
kinetic energy distribution across the centerline of the combustor. The maximum
values exist around the primary jet, especially near the impingement point.
Calculation without swirl pr
edicted higher values of turbulence kinetics energy and
stronger back flow than that the calculated values with swirl flow. This explains for
higher levels of swirl reduce the turbulence intensity in the center of the combustor as
shown in Figure (
20
).

Z(m)

Figure(18):Temperature Distribution along the Centerlines of the
Combustion Chamber.



Figure(20): Turbulence Kinetics Energy Distribution along the Centerlines

of the Combustion Chamber.

Figure(19):Velocity Vector (without swirl flow) View in the Two
Longitudinal Planes that Containing the Primary and Dilution Holes.

Journal of Babylon University/Pure and Applied Sciences/ No.(5)/ Vol.(18): 2010



12
-
References:

-
AGARD(1980) Combustor modeling, conference proceeding 275 propulsion and
Engineering panel 54th
(B) specialists' Meeting, DFVLR, Cologne, Germany, 3
-
5 October 1979.

-
AGARD (1987) Combustion And Fuel in gas turbine combustors, Conference
Proceeding 422. Propulsion and Energetic Panel 70th Symposium, Chania, Crete
Greece, 19
-
23 October 1988.

-
AGARD (19
93) Fuels and Combustion Technology for Advanced Aircraft Engine,
Conference Proceeding 536. Propulsion and Energetic Panel 81th Symposium,
Fluggi, italy, 10
-
14 May 1993.

-
Alizadeh, S. and Moss, J.B., "flowfield Prediction of NO
x

and Smoke Production in
Ai
rcraft Engine", see AGARD(1993).

-
Allan, J. (1996) "CFD Indications of Premix Flame Stability in a Gas Turbine
Combustor", ASME Paper 96
-
GT
-
182.

-
Baron, F., Kanniche, M., and Mechitoua, N. (1994) " Influence of Inlet Conditions
on Flow and Combustion Chara
cteristics of a Model Annular Combustor ",
ASME paper 94
-
GT
-
281.

-
Bicen, A.F., Tse, D.G., Whitelaw, J.H., (1990) "Combustion Characteristics of a
Model Can
-
Type Combustor", Combustion and Flame, 80, 111
-
125.

-
Chen, T.R.A. and Murman, E.M. (1985) "A finite
-
Volume Method for the
Calculation of Compressible Chemically Reacting Flows" AIAA Paper 85
-
0331.

-
Correa, S.M., and Shy, W., (1987) " computational Models and Methods for
Gaseous Turbulent Combustion" Progress in Energy and Combustion Science 13,
249
-
292
.

-
Crocker, D.S., and Smith, C.E. (1993) "Numerical Investigation of Enhanced
Dilution Zone Mixing in a Reverse Flow Gas Turbine Combustor", ASME Paper
93
-
GT
-
129.

-
Danis, A.M., Burrus, D.L. and Mongia, H.C. (1996) "Anchored CCD for Gas turbine
Combustors D
esign and Data Correlation "ASME Paper 96
-
GT
-
143.

-
Di Martino, P., Colantuoni, S., Cirillo, L., and Cinque, G. (1984) " CFD Modeling of
-

Flow Combustor", ASME Paper 94
-
GT
-
468.

-
Gran, I.R., Melaaen, M.C., and Magnussen, F. (1993)
" numerical Simulation Fluid
flow and Combustion in gas turbine combustors " ASME Paper 93
-
GT
-
166.

-
Hallmann, M., Scheurlen, M. and Wittig, S. (1993) "Computation of Turbulent
Evaporating Sprays Eulerins vs. Lagrangian Approach", ASME Paper 93
-
GT
-
333.

-
Hwa
ng, C.C., Zhu, G., Massoudi, M. and Ekmann, J.M. (1993) " A comparison of
the Linear and Nonlinear K
-

Turbulence Models in Combustors", Journal of
Fluids Engineering , 115, 93
-
102.

-
Jones, W.P., Launder, B.E., (1972), "The Prediction of Laminarization wit
h a Two
-
Equation Model of Turbulence" Journal of Heat and Mass Transfer, 15, 301
-
314.

-
Karki, K.C., Oechsle, V.L., and Mongia, H.C. (1992) "A computational Procedure
for Diffuser
-

Combustor Flow Interaction Analysis " Journal of Engineering for
Power , 114
, 1
-
7.

-
Kahlil, E.E. (1982) "Modeling of Furnaces and Combustors", Abacus Press.

-
Kahlil, E.E., Spalding, D.B., Whitelaw, J.H., (1975) "The Calculation of Local Flow
Properties in Two
-
Dimensional Furnaces", Journal of Heat and Mass Transfer,
18, 775
-
791.

-
Kuo, K.K. (1996) "Principles of Combustion", John Wiley and Sons.



-
Launder, B.E., Spalding, D.B., (1974) "The Numerical Computation of Turbulent
Flows" Computer Methods in Applied Mechanics and Engineering, 269
-
289.

-
Lawson, R.J. (1993) "Computational Mode
ling of an Aircraft Engine Combustor to
Achieve Target Exit Temperature Profiles", ASME Paper 93
-
GT
-
164.

-
Magnussen, B. F. , Hjertager, B. H., (1976) " On mathematical models of turbulent
combustion with special emphasis on soot formation and combustion",
In 16th
Symposium (Int'l.) on Combustion. The Combustion Institute.

-
McGuirk, J.J. and Spencer, A. (199
6
) "Computational Methods for Modeling Port
Flows in Gas Turbine Combustors ", ASME Paper 95
-
GT
-
414.

-
McGuirk, J.J., Palma, J.M.L.M., (1993) "The Flow In
side A Model Gas Turbine
Combustor: Calculations", Journal of Engineering for Power, 115, 594
-
602, July
.

-
Nikjooy, M., So, R.M.C., Peck, R.E., (1988) "Modeling of Jet and Swirl
-

Stabilized
Reacting Flows in Axi
-
symetric Combustor" Combustion Science and
Te
chnology, 58, 135
-
153.

-
Patankar, S.V., Spalding, D.B., (1972) "A Calculation Procedure for Heat, Mass and
Momentum Transfer in Three
-
Dimensional Parabolic Flows"Journal of Heat and
Mass Transfer, 15, 1787
-
1806.

-
Pope, S.B. (1991) "Combustion Modeling usin
g Probability Density Functions", in
Oran, E.S., and Boris, J.P., eds. (1991) "numerical Approaches to Combustor
Modeling ", AIAA Progress in Astronautics and Aeronautics, Vol. 135, 349
-
364.

-
Schetz, J.A. (1993) Boundary Layer Analysis, Prentice Hill.

-
Ser
ag
-
Eldin, M.A., Spalding, D.B., (1979) "Computations of Three
-
Dimensional
Gas
-
Turbine Combustion Chamber Flows", Journal of Engineering for Power,
101, July, 326
-
337.

-
Shy, W., Correa, S.M., and Braaten, M.E., (1988) "Computation of Flow in A Gas
Turbine C
ombustor" Combustion Science and Technology, 58, 97
-
117
.

-
Sirignano, W.A. (1993) "Computational Spray Combustor ", in Chung, T.J., ed.
(1993) Numerical Modeling in Combustion , Taylor and Francis, 453
-
542.

-
Spalding, D.B., (1976) "Mathematical Models of Tu
rbulent Flames:Review",
Combustion Science and Technology,13, 3
-
25
.

-
Swithenbank, J., Turan, A. and Felton, P.G. (1980) " Three Dimensional Two=Phase
Mathematical Modeling of Gas Turbine Combustors", in Lefebvre , A.E.(Ed.)
Gas Turbine Combustor Design Pro
blems", Hemisphere Publishing Corporation,
1980.

-
Su, Ay, Chang, Ja
-
Hwan, and Ho, Wu
-
Chi, (2003) "Numerical Simulations of Coal
-
Fired Swirl Rectangular Combustor", Industrial Technology Research Institute.

-
Topaldi, A.K., Hu, I.Z., Correa, S.M., and Burrus
, D.L. (1996) "Coupled Lagrangian
Monte
-
Carlo PDF
-
CFD Computation of Gas Turbine Combustor Flowfields with
Finite
-
Rate Chemistry" ASME Paper 96
-
GT
-
205.

-
Versteeg, H.K., Malalasekera, W., (1995) "An Introduct
ion to Computational Fluid
Dynamics The Finite Volume Method", Longman Group Ltd.