EPA-650/2-73-045
December 1973
Environmental Protection Technology Series
-------
EPA-650/2-73-045
A STUDY OF COMBUSTOR FLOW
COMPUTATIONS AND COMPARISON
WITH EXPERIMENT
by
R. F. Anasoulis and H . McDonald
United Aircraft Research Laboratories
400 Main Street
East Hartford, Connecticut 06108
Contract No. 68-02-0267
Program Element No. 1AB014
ROAP No. 21ADG-10
EPA Project Officer: D. W. Pershlng
Control Systems Laboratory
National Environmental Research Center
Research Triangle Park, North Carolina 27711
Prepared for
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
WASHINGTON, D.C. 20460
December 1973
-------
This report has been reviewed by the Environmental Protection Agency and
approved for publication. Approval does not signify that the contents
necessarily reflect the views and policies of the Agency, nor does
mention of trade names or commercial products constitute endorsement
or recommendation for use.
ii
-------
ABSTRACT
A computational procedure for calculating the coupled flow and chemistry
within combustion devices is presented. The procedure solves the time-
averaged Wavier-Stokes equations with coupled chemistry, including the effects
of turbulence and radiative heat transfer, using a novel field relaxation
method developed at the United Aircraft Research Laboratories. Although a
relatively simple turbulence model is employed in the procedure, modification
to this model is easily implemented within the framework of the computational
method. Computations of the flow and chemistry within a representative
furnace have been made using the procedure and are presented and compared
with experimental data.
This report was submitted in fulfillment of Contract 68-02-026? by
United Aircraft Research Laboratories, under the sponsorship of the
Environmental Protection Agency. Work was completed as of December 1973.
111
-------
CONTENTS
ABSTRACT iii
LIST OF FIGURES vi
LIST OF TABLES viii
ACKNOWLEDGMENTS ix
SECTION
I INTRODUCTION 1
II THEORETICAL ANALYSIS 1+
III COMPUTATIONAL ANALYSIS 33
IV RESULTS AND DISCUSSION 56
V CONCLUSIONS ?8
VI REFERENCES ..... 80
VII LIST OF SYMBOLS 85
-------
FIGURES
No.
1 Coordinate System for Radiant Flux 22
2 Burner Flow Model 28
3 Relaxation Procedures 35
U Vorticity Residual History 36
5 Vorticity Residual History 38
6 Computed Velocity Distributions 39
7 Extrapolation Procedure ^5
8 A Comparison Between Predicted and Measured Axial Decay of
The Velocity Maximum in a Turbulent Swirling Jet 60
9 A Comparison Between Predicted and Measured Axial Decay of
The Swirl Velocity Maximum in a Turbulent Swirling Jet 6l
10 Comparison Between Theory and Measurement of the Mean Velocity
Profile in a Turbulent Swirling Jet 62
11 A Comparison Between Predicted and Measured Axial Decay of
The Velocity Maximum in a Turbulent Swirling Jet 63
12 A Comparison Between Predicted and Measured Axial Decay of
The Swirl Velocity Maximum in a Turbulent Swirling Jet 64
13 The Comparison Between Theory and Measurement of the Mean
Velocity Profile in a Turbulent Swirling Jet 65
lU A Comparison Between Predicted and Measured Axial Decay of
The Velocity Maximum in a Turbulent Swirling Jet 66
15 A Comparison Between Predicted and Measured Axial Decay of
The Swirl Velocity Maximum in a Turbulent Swirling Jet 67
vi
-------
No.
FIGURES (Concluded)
l6 Comparison Between Theory and Measurement of the Mean
Velocity Profile in a Turbulent Swirling Jet 68
17 Radial Distribution of Velocity Vector 69
18 Radial Distribution of the Axial Velocity 70
19 Radial Distribution of Axial Velocity 71
20 Theoretical Streamlines for IGT Swirl Burner 72
vii
-------
TABLES
No. Page
I Coefficients of General Elliptic Equation 8
II Chemical Rate Constants 1?
Ill Auxiliary Reactions for Hydrocarbon Combustion Equilibrium
Model 20
viii
-------
ACKNOWLEDGMENTS
Dr. W. R. Briley contributed greatly to the successful development of
the FREP code described in the present report. In addition, Mr. W. A. Carey
constructed the code input/output procedures, including the overlay and
peripheral storage routines while Mr. R. C. Buggeln gave considerable
assistance with the running of the test cases. The work was supported by the
Environmental Protection Agency under Contract No. 68-02-026? with
Mr. David Pershing as the Contract Officer.
IX
-------
SECTION I
INTRODUCTION
Concern over air pollution and increased demands for higher efficiency
have placed considerable burden on present day designers of furnaces. Methods
for modeling furnaces have generally been less than completely satisfactory
largely due to lack of understanding of the fundamental flow processes which,
through heat, mass, and momentum exchange, directly influence pollutant
emission generation and combustion efficiency. For example, it has been
shown (Refs. 1 and 2) that stability and combustion intensity of flames as
well as residence time in furnaces are strongly influenced by swirling flow.
Residence times and flame stability and intensity are, in turn, related to
furnace performance and efficiency as well as to pollutant emission formation
(Refs. 3 and U). Thus, because of the strong influence of flow processes on
combustion characteristics, it is particularly important that the fluid
mechanics be properly taken into account in development of furnace modeling
techniques.
Procedures employed for modeling combustion processes have generally
been highly simplified, particularly in regard to flow modeling where
stirred reactor concepts and one-dimensional assumptions are employed (Refs. 5
through 10). Examples are the studies conducted by Fletcher and Heywood
(Ref. 5) and Hammond and Mellor (Refs. 6 and 7), who employed the stirred
reactor concept to model the flow in the primary zone of gas turbine combus-
tion chambers in order to assess the effect of residence time on combustion
behavior and to predict pollutant emissions. Although radiant heat transport
was neglected in both these studies, a sophisticated chemistry model, based
on quasi-global finite-rate hydrocarbon mechanism, was employed by Hammond
and Mellor to account for the fuel decomposition. Other studies included the
analysis of Roberts, et al., (Ref. 8), who, in order to predict nitric oxide
production in gas turbine combustors, divided the combustion chamber into
three regions: one corresponding to the central recirculation portion of the
upstream zone; a second representing the flow region surrounding the
recirculation zone which was interpreted to be a one-dimensional reacting
zone; and the third downstream zone modeled as a one-dimensional region. An
interesting result of the Roberts analysis was the small difference noted in
predicted nitric oxide levels in comparisons between coupled and uncoupled
hydrocarbon thermochemistry systems which were studied; thereby, suggesting
that a decoupling of the hydrocarbon chemistry from the nitric oxide chemistry
might be a reasonable approximation. No account was taken of radiant heat
transport in the Roberts analysis. An extension of the Roberts, et al.,
analysis, directed toward low power gas turbine operation, was developed by
-------
Mosier, et al., (Ref. 9)> where, by employing a more sophisticated finite-
rate hydrocarbon mechanism, trends were predicted which were in agreement
with experimental data. In another study, Edelman and Economos (Ref. 10),
in an attempt at developing a general analytical procedure for predicting
combustion behavior, applied a modular technique whereby various critical
combustion processes were treated on an individual basis or coupled as a
function of operating conditions. The Edelman approach may be critized in
its manner of accounting for recirculation (a stirred reaction is used),
disregard of radiative transport effects, and its inability to provide a
unified descriptien of a combustion chamber under a given set of operating
conditions.
As a general criticism, it seems that all of the foregoing methods are
lacking primarily in their ability to properly account for mixing and reverse
flow phenomena occurring in the recirculation zone of a furnace. However, a
computational method which does allow for calculation in recirculation
regions has been proposed recently by Gosman, et al., (Ref. 11), who solved
the governing equations by a relatively inefficient explicit point by point
relaxation procedure. Employing the Gosman method to demonstrate the
feasibility of making computations in the recirculation zones of combustion
chambers, sample computations were conducted at UARL in a representative
burner (Ref. 12). Highly simplified chemistry and turbulence models were
employed in these initial predictions but no account was taken of radiative
transport. The results obtained with this procedure demonstrated qualitative
agreement with experimental observations, and were very encouraging and, as
a consequence, led to development by UARL of an improved numerical code for
solving combusting flows containing recirculation zones. The UARL procedure
is a novel implicit computational scheme in which residuals are relaxed
simultaneously throughout the entire flow field rather than one at a time
which is characteristic of the explicit point methods. Because of this
feature the UARL code has proven to be considerably more efficient than the
point methods, especially as the number of grid nodes are increased.
A principal objective of the study under this contract was to further
develop and refine the UARL computational procedure in order to determine and
demonstrate the potential for predicting air pollutant emissions from
furnaces. For this purpose, computations of the flow field within a
representative furnace were undertaken as part of this contract and compari-
sons made against both hot- and cold-flow measurements. The computational
procedure, which is referred to as the FREP code (acronym for Field Relaxation
Elliptic Procedure), solves the time-averaged Navier-Stokes equations with
coupled chemistry, turbulence, and radiant heat transfer models. The
turbulence model consists of an eddy viscosity model derived from an extension
of Erandtl's mixing length hypothesis and a wall-flux model which is employed
-------
to resolve the flow field in the immediate vicinity of the wall where large
gradients are known to exist. (The turbulence model was partially developed
under a jointly sponsored FAA/WPAFB-APL contract, No. F33615-72-0-20^2.) The
radiant heat transfer model incorporated into the procedure is a four-flux
collimated model requiring solution of two differential equations to obtain
the flux of radiant energy in the radial and axial directions. The radiant
heat flux obtained by solution of these equations is included as an
additional term in the energy equation. The chemistry model is based on
equilibrium hydrocarbon chemistry coupled with finite-rate nitric oxide
chemistry. The intent in development of the flow models, described above,
was to minimize undue complexity and sophistication and to provide a
reasonably good framework within which refinements could be easily
implemented at a future time, if warranted by comparison with experimental
data.
-------
SECTION II
THEORETICAL ANALYSIS
II-A. APPROACH
The flew regime considered in the present study is a quasi-steady
gaseous phase turbulent flow system with combustion. Turbulent exchange
coefficients, defined by analogy with Newton's, Fourier's, and Pick's laws
for laminar flows, are employed to define momentum, enthalpy, and chemical
species fluxes, respectively. A knowledge of these turbulent fluxes is
required to close the governing system of equations. Values for the exchange
coefficients are determined by specification of a turbulence model together
with values of Prandtl and Schmidt numbers taken from knowledge of turbulent
flows of gases and gas mixtures. A chemistry model and a radiation model are
employed to evaluate chemical species rate expressions and radiation heat
transfer effects on furnace characteristics.
In order to solve the combustion problem, it is necessary, in addition
to specifying turbulent transport models, to employ a computational method
for solving the complex system of equations. The computational procedure
must be capable of treating the flow resulting from the interaction of
mixing and chemical reactions in the turbulent diffusion flame and from
sudden changes in flow properties which may be prompted as a result of the
combustion process, any or all of which could lead to eventual divergent
behavior. Taking into consideration these potential difficulties, a powerful
novel numerical procedure, developed at UARL, is applied under this contract
to the solution of the equations governing the flow and chemistry within a
representative furnace. Time-mean-average velocities, pressure, temperatures,
and concentrations are computed and comparisons are made against measurements
to allow for an evaluation of the procedure and possible improvements in
hypotheses and in furnace analytical modeling techniques.
II-B. GOVERNING EQUATIONS
Under consideration is the flow of a turbulent radiating chemically
reacting multicomponent mixture with heat and mass transport. The elliptic
time-averaged partial differential equations describing this combustion system
are based on the conservation laws of mass, momentum, energy, and chemical
species (Ref. 13). For simplicity, these equations are expressed in vector
notation below. Negligible coupling between diffusion and thermal gradients,
no body forces, and negligible bulk viscosity all are assumed. Fick's law is
-------
presumed valid throughout, implying equal binary diffusion coefficients for
each pair of species in the mixture. The resulting set of equations can be
written
Continuity :
V pv =0 (!)
Conservation of Species :
V- pv m, - V- ( I T: V mj) - r, - Irr = 0 (2)
i i
Conservation of Momentum:
[Vv + Vv -- 1- (V- v ) g] + V P =0 (3)
Conservation of Energy:
j) - V- ( Th CpVl) - V- (^ff V ( v2/2))+ QR -- o
There is also the thermodynamic relationships
"= (5)
H= (2 cp. m-,) dT + 2h| m. + -JL_- (6)
-------
which provide additional equations necessary to close the system. These
equations are written for variable fluid properties and embody mass diffusion,
chemical reaction, radiation heat transfer, and viscous dissipation. In the
energy equation, the kinetic heating terms,which assume importance only when
the swirl component of velocity predominates, have been neglected. In order
to solve the above equations in addition to boundary conditions, it is
necessary to specify expressions for turbulent exchange coefficients, I\, Fn,
and Meff, rate of generation of chemical species, ri (in this case only
nitric oxide), and radiant energy transport, Qy. In regard to the turbulent
exchange coefficients, because in the present analysis Prandtl and Schmidt
numbers are employed in the computations based on knowledge of turbulent flows
of gases and gas mixtures, only a turbulent momentum exchange coefficient
(i.e., Meff) need be specified. This is a consequence of the fact that the
enthalpy and chemical species turbulent exchange coefficients are related to
the turbulent effective Prandtl and Schmidt numbers through the effective
viscosity via the relations
P''« " r (7)
sc.f( • ^ (8)
To define the effective viscosity, fieff> a turbulence model is employed.
Likewise, a chemistry model and a radiation model are employed to define the
rate of production of chemical species and radiant energy transports,
respectively. These models are discussed in detail in subsequent sections.
At this time it is convenient to express the momentum equations in alternate
form in terms of vorticity and stream function. This is done in order to
make the subsequent mathematical treatment easier to follow and because it is
more convenient in the numerical computational procedure since it eliminates
pressure as a variable.
For the axisymmetric case, vorticity and stream function can be defined
in cylindrical coordinates according to the following relations:
-------
— - -^- (9)
ax dr I9j
and
(10)
Employing Eqs. (l) and (2), all of the governing differential equations,
Eqs. (l) through (U), can be converted into the following general form for
axisymmetric flow (see Ref. 11 for details)
ac
[ v -
(11)
= 0
where the coefficients a
-------
TABLE I
COEFFICIENTS OF GENERAL ELLIPTIC EQUATION
a^
6x
JL
6r
- b b. r
= 0
FUNCTION
mi
F
rvt
M.
r
'
>+
1-0
0
0
i r\
\ • u
I.O
r2
0
*,
r'eff
0
. -i.o
IVK.)
w
r2
1
b*2
r'eff
-i.o
«*•*!»
0
.^
r2
1
c*
I.O
I.O
I.O
i o
*
Meff
I.O
dc£
-n
a 2
Ktf
0 4 S
=— IG-o-T )--7r-vG-F)
r 2r
i a r ( f i ' ^ av 2/2
— * 1 r^6TT ' I c>— ^ \/
rdxL *• rr dx
+ 2 (-^ -\
' dx J J ~^~ ar [^e^f
0
a /-w.2, r a /u2+v2 \ dp
d*(p ' } 'lax \ 2 / ar
d / u2+v2\ dp i
ar \ 2 / ax J
- r2 So,
CJ
"~r~
8
-------
So»= -f- [ V (ir • v) • V(i • VM ) + V (ir v )
6 (12)
is excluded from the computation, an omission which may be expected to result
in only small errors. These terms can be included readily if circumstances
are encountered where they might be significant. In the numerical solution
procedure, all equations are solved in the form of Eq. (ll).
Outside of boundary conditions, to close the system of equations, there
remains to specify a turbulence model, a chemistry model, and a radiation
model to define effective viscosity, chemical species production, and
radiation heat flux, respectively. These models are discussed separately
below.
II-B.I. Turbulence Model
The flow in furnaces is known to be predominantly turbulent. In order
to account for this turbulent behavior in the solution of the governing time-
mean equations of motion, a turbulence model is introduced to define the
effective viscosity. An historical review and account of turbulence models is
available in the literature (see, for instance, Ref. lU). Prandtl (Ref. 15)
was perhaps the first to introduce a turbulence model when he postulated that
time-averaged shear stresses and time-averaged velocities are proportional as
in laminar flow and, furthermore, that the length scale which enters into
proportionality, which he called the mixing length, is proportional to the
turbulent shear region thickness. Prandtl's mixing length model has been
employed successfully by a number of investigators (for example, see Refs. 16
through 20) in a variety of problems primarily involving turbulent motion
along walls and in free turbulent flows. A disadvantage of the mixing length
model is that it is an equilibrium model (i.e., turbulence is presumed to be
produced and dissipated locally) and requires an ad hoc mixing length
distribution. Use of more advanced turbulence models requiring solution of
additional differential equations would serve to eliminate these limitations;
however, the development of these sophisticated models is in the formative
stages and successful applications of these techniques to practical flow
problems has not been clearly demonstrated.
-------
In view of the early stage of development of the more sophisticated
turbulence models, a mixing length model has been employed in the present
analysis. The mixing length model was partially developed for use with the
FREP code under a jointly sponsored FAA/WPAFB-APL contract, No. F336l5~72-C-
20^2. The formulation of the turbulence model employed in the analysis is
based on an evaluation of the mixing length from the work of Williamson
(Ref. 19) and Lilley (Ref. 20) with a correction added to the Williamson
model to account for the swirling flow. Williamson's model, which is derived
from the flow in ducts and channels, is employed in the main body of flow
away from the injection ports which are located near the axis of a represen-
tative burner (IGT) under consideration. Lilley's model, which includes a
correction for swirling flow, is directly applicable to free jet flows and
is employed in the vicinity of the injection ports. Transition between the
two models is accomplished by monitoring the cross over point of the mixing
lengths computed for each model at each axial station taking care to ensure
that a smooth transition occurs. If there is no cross over, the Lilley model
is used. In mathematical terms the mixing length expression, as employed in
this study, for effective viscosity in the rx and r9 directions for the case
of axisymmetric swirling flow takes the form (Ref. 21)
du
+ n I A/ V 11-
LAM
Mrx
where the mixing length, jt, is given by either Williamson's model
--=,.75 xe('-"" <«>
10
-------
or Lilley's model
um
where
X = 0.08 ( 1+ X. Sx) (17)
X = 0 60 (18)
S, = Gfl/(Gx- r^ ) (19)
- -. o.oi
,
and
= 1.0 + 5.0 S '/3 (2°)
The factor (1+\SSX) represents the stretching of the length scale due to
swirl, the mixing length parameter, \, approaching an acceptable nonswirling
value of 0.08 for the case of no swirl. Constants in Eqs. (15) through (20)
have been determined from comparisons of predictions with experimental data.
The turbulence model employed in this study represents then an extension of
Prandtl's isotropic mixing length model in that nonisotropy of the viscosity
is accounted for through definition of the variable r8 viscosity, Eq. (1*0,
and through linking of the rx and r8 shear through Eq. (l3)«
11
-------
In implementation of the effective viscosity, defined by Eqs. (13)
through (20), a difficulty arises in regions close to the wall if grid
spacing cannot be refined sufficiently to resolve the large flow gradients
occurring there. In the event grid resolution in the near wall region is
poor, errors in the local flow properties will occur and be propagated into
the main body of flow. To prevent this occurrence, a universal velocity
profile known to be accurate in the vicinity of the walls is imposed on the
solution at the near wall nodes. This well-known universal velocity distri-
bution, which is referred to as the 'law of the wall1, has the form (Ref. 22)
u
•f
+ = ' |n (i + y + ) + [(i . I _ c a) y+ - c. 1 e " ay ^ C2 (2l)
Ui L ^| * *• J
where
+ U
U - ^= (22)
(23)
The constants a, C^, and C% are empirically determined. With the proper
velocity distribution and, consequently, also the velocity gradient computed
at the near wall nodes from Eq. (21), the effective viscosity can then be
computed from
^ '
which is the boundary layer equivalent of Eq. (13) which is appropriate for
application near the wall. Implicit in the application of Eqs. (21) through
(2U) is the assumption that normal to the wall the shear stress in the
immediate vicinity of the wall is constant. Thus, specifying the velocity
and velocity gradient at the near wall nodes in this fashion is equivalent
to assigning a 'slip* velocity at the wall itself and imposes an artificial
12
-------
velocity distribution between the wall and its adjacent nodes. Flow
resolution in the narrow region very close to the wall, however, is of little
importance in the present task and, hence, the use of the law of the wall is
a small price to pay for the high accuracy attainable in the major portion of
the flow field. However, in the event accurate heat and mass fluxes at the
wall are desirable at a future time, refinements to the wall flux procedure
can be easily implemented within the framework of the computational
procedure.
II-B.2. Chemistry Model
Closure of the governing system of equations, Eqs. (l) through (6),
requires specification of the rate expression, r^, in the chemical species
equations, Eq. (2). However, because only nitric oxide chemistry is
kinetically treated in this investigation, only an expression for rate of
production of nitric oxide need be formulated. To develop the rate of
expression for nitric oxide requires knowledge of concentrations of species
entering into the nitric oxide reaction scheme which are dependent upon the
hydrocarbon fuel (methane, C!fy, in this case) combustion process. Under this
study, combustion of the hydrocarbon fuel is treated from chemical equilib-
rium considerations and, as a consequence, concentrations for products of
combustion to be employed in the nitric oxide rate expression are computed in
a straightforward manner from a chemical equilibrium analysis. Justification
for the hydrocarbon equilibrium model is based on the observation that, under
pressure and temperature conditions generally existing in burners, fuel
oxidation reactions can be expected to go to completion rapidly compared to
nitric oxide formation rates. Thus, the chemistry model employed in the
analysis consists of two parts: a hydrocarbon oxidation process which is
considered to be mixing controlled and a nitric oxide formation process which
is kinetically controlled. Also, within the analysis, because the nitric
oxide is considered a trace species (i.e., its' concentration is low enough
to have little influence on mixture properties), the nitric oxide and hydro-
carbon chemistry is uncoupled in the solution procedure. Hence, the nitric
oxide species equation is solved separately from the other equations and only
at the end of the iteration process after final converged values for all flow
variables have been determined. The two components of the chemistry model,
representing the equilibrium hydrocarbon chemistry analysis and the finite-
rate nitric oxide analysis, are discussed separately in detail below.
13
-------
II-B.21. Nitric Oxide Chemical Analysis
Solution of the species equations to account for convection, diffusion,
and production of nitric oxide requires definition of an expression
describing rate of creation of nitric oxide (see Eq. (2)). In order to
develop this expression, knowledge of the reaction mechanism by which nitric
oxide is formed is required. A generally accepted reaction scheme for nitric
oxide formation and decomposition in post flame gases is that proposed by
Lavoie, et al., (Ref. 23). It consists of the following six reactions:
N + NO •» N2 + 0 (25)
N +• Og s± NO + 0 (26)
N + OH * NO + H (27)
H + N20 =^N2 + OH (28)
0 + N20 ^ N2 + 02 (29)
0 + N20 -NO + NO (30)
The first two reactions, Eqs. (25) and (26), form the Zeldovitch
mechanism {Ref. 2k) and is considered to be the principal nitric oxide
formation reaction mechanism. These two reactions together with the third
reaction, Eq. (27), which assumes minor importance under fuel rich conditions,
form the extended Zeldovitch mechanism which is employed in the present study.
At low temperatures when nitric oxide concentrations are much greater than
equilibrium values, the fourth, fifth, and sixth reactions, Eqs. (28) through
(30), involving NgO as intermediary, may become important; however, because
overall reaction rates are so low, the net effect of these reactions is
probably negligible and, therefore, these reactions have been neglected in
the present analysis. The results of Bowman and Seery (Ref. 25), in
investigation of nitric oxide formation kinetics in combustion pressures, lend
support to this approach and it has been followed by other investigators
(Refs. 5 and 2?).
-------
With the reaction mechanisms for nitric oxide formation defined by the
extended Zeldovitch mechanism, Eqs. (25) through (27), it becomes possible
to develop an expression for the rate term, r^> entering in the chemical
species equation, Eq. (2). The mathematical development of the rate term
has as its basis the 'Law of Mass Action' which states that the rate of
chemical reaction is proportional to the active masses of the reacting
materials. Thus, referring to Eqs. (25) through (27), rate expressions for
nitric oxide and nitrogen can be written
dt = k.b[N2][o] + kzf[N] [o2] 4 k3f [N][OH]
- klf [N] [NO] - k2b [NO] [o] - k3b [NO] [H]
M
= klb [NJ [o] 4 k2b [NO][O] + k3b [NO] [H]
(32)
-klf [N] [NO] - kzf [N][OZ]- k3f [N][OH]
where [ j denotes moles per unit volume and k^ and k^ are the rate constants
(volume/mole-sec) in the backward and forward directions, respectively.
Because relaxation times for Eq. (32) are several orders of magnitude shorter
than for Eq. (31), it is a good approximation to assume steady concentration
for [Nj. Therefore, setting dN/dt =0, Eq. (32) can be solved for [NJ as
, = h,b[N2][o] 4 k2b[No][0]4 k3b[NQ][H] ^^
J klf H 4 k2f [o2] 4 k3f [OH]
Equation (31) can now be written
r~i . r.._i M . r i r i
(3»0
M [H]
[N]{k2f [oj+ k3f [oHJ-klf [NO]}
15
-------
where [N] is obtained from Eq. (33). Equation (3*0 is the desired rate
expression for [No] in terms of concentrations of N2, 0, H, ©2, and OH which
are obtained from an equilibrium hydrocarbon air combustion analysis,
discussed in detail below. The dimensionally correct expression for the
rate term, r^, to be employed in solution of Eq. (2) is
(35)
The rate constants governing the nitric oxide reaction scheme, Eqs. (25)
through (27), and utilized in the rate expression, Eq. (3*0, are employed in
the solution procedure in modified Arrhenius form, i.e.,
= AjTNj exp (-Bj/RT)
(36)
where data for the activation energies, Bj, the frequency factors, A,-, and
exponent, Nj, are taken from Refs. 27 and 28 and tabulated in TABLE II.
II-B.2ii. Equilibrium Hydrocarbon Air Combustion Analysis
In order to determine the state of the burnt gases and provide the
necessary information concerning concentrations of species entering into the
nitric oxide rate expression, Eq. (35), an equilibrium hydrocarbon-air
combustion model is employed in this investigation. Chemical equilibrium
is that condition under which the forward and reverse speeds of a chemical
reaction are equal and, hence, concentrations of reactants and products may
be determined. In the present analysis, the equilibrium hydrocarbon model
is implemented by first determining fuel-air ratios by solution of the
governing species diffusion equations without rate terms (i.e., assuming the
combustion process to be mixing controlled), followed by implementation of an
equilibrium hydrocarbon-air computational routine to establish coexisting
concentrations of products of combustion at the existing temperature and
pressure conditions. The equilibrium computational analysis employed for the
hydrocarbon-air combustion model is based on the work of Brinkley (Refs. 29
and 30). For a brief account of the approach, consider a chemically correct
reaction of a representative parafinic hydrocarbon with air
16
-------
TABLEn
CHEMICAL RATE CONSTANTS
k = ATNexp(-B/RT)*
REACTION
NO + N ^± O-f N2
N +O2 £=^ NO + O
N + OH ^f^NO + H
FORWARD RATE
A
3.10 x 1013
6.43 x 109
4.22 x 1013
B
334.
6250.
0
N
0
1.0
0
REVERSE RATE
A
1.36x1014
1.55x109
1.64 x1014
B
75400.
38640.
48600.
N
0
1.0
0
* UNITS ARE CM, CAL, °K, g-MOLE, SEC AND ATM
-------
Cn H2n+2 + b02+ 3.76 x bN2 "~C C02
(37)
+ d H20 + b x 3-76 N2
To solve Eq. (3?) for the concentration of the products of reaction, the
coefficients b, c, and d must be determined. For this purpose the required
number of additional equations for an atom balance can be written between C,
0, or H atoms resulting in an algebraic set of three equations which can be
solved for the three unknowns. The high temperature levels commonly
encountered in combustion devices justify the assumption that the hydro-
carbon-air reaction goes to completion as implied by Eq. (37); however, in
the usual case, the mixture of hydrocarbon fuel and air is not stoichiometric
as indicated. For the general case of fuel rich or fuel lean mixture the
products of combustion, in addition to C02, 1^0, and N2, may also include CO,
H2> 02» OHJ H» and ° atoms and radicals and other less populous constituents.
To solve for the concentration of these additional species, it is recognized
that additional equations are required to define additional unknown
coefficients. Therefore, auxiliary reactions for which equilibrium constants
are known are employed to provide the required additional equations. These
additional equations are written in terms of the unknown species concentra-
tions and the known equilibrium constants. As a simple illustration,
consider the case of two additional unknowns, as for example, in the
following reaction
+ 3-76 x bN2 *"° CO + C C02
dH20+ 6H2 + 3.76 X bN2
In addition to the reaction equation three equations can be written on the
basis of a C, 0, and H atom balance; however, there are five unknowns, a, b,
c, d, and e. The necessary fifth equation is obtained from an auxiliary
reaction for which the equilibrium constant is known. For example, the
water-gas reaction can be written
18
-------
C02 + H2 v H2 0 + CO (39)
from which the known equilibrium constant, k_, would be written as
.. . ad
0+0)
resulting in five equations from which the five unknowns coefficients can be
determined. The choice of additional reactions is somewhat arbitrary and can
vary depending on local fuel-air rates and temperature conditions, in order
to optimize the iterative solution procedure. In the present case, a test
is made of the oxygen concentration and depending on its magnitude relative
to the carbon and hydrogen concentration, three possible sets of additional
reactions are employed. These three cases are listed in TABLE III. Depending
on which of the three categories the oxygen concentration falls into, the
presumed primary products of combustion are indicated along with five
representative auxiliary equations which may be employed in the reaction
scheme for which equilibrium constants are known. The equilibrium constants
for the auxiliary equations are obtained from the data tabulated in
reference 31. Thus, with the equations obtained from the atom balance and
additional equations provided by the auxiliary equilibrium reactions, the
problem is completely specified and the unknown concentrations (some of which
are employed in the nitric oxide expression, Eq. (3^))» are solved for using
a Newton-Raphson numerical iterative procedure.
II-B.3 Radiation Model
The purpose of the radiation model employed under this investigation is
to define the radiant energy contribution to combustion heat transfer rates.
It is employed specifically in the analysis to define the radiant heat flux
term, Qp, entering in the energy equation, Eq.(U), through which coupling
between radiation, convection, and conduction modes of heat transfer occurs.
The radiative heat transfer process is inherently different from conductive
and convective heat transfer processes in that, in the two latter cases, the
energy is transferred by molecular collision and transport whereas, in the
radiative process, energy transfer is in the form of electromagnetic waves,
not requiring contact of the molecules. Radiation becomes of importance
19
-------
TABLED!
AUXILIARY REACTIONS FOR HYDROCARBON COMBUSTION EQUILIBRIUM MODEL
ro
o
OXYGEN CONCENTRATION
PRINCIPAL PRODUCTS
AUXILIARY REACTIONS
CASET
H2 , H20 , CO
C0+H20 * C02 + H2
2H20 5t 02 + 2H2
H20 = O + H2
1/2 H2 sfc H
H2 0 - PH + 1/2 H2
CASEH
H2 0 , CO , C02
CO+ H2 0 ^ C02+ H2
2C02 ^ 2CO + 02
1/2 CO+l/2 H2 0 -
1/2 C02 + H
1/2 C02 + 1/2 H20 =*
1/2 CO + OH
C02 - CO+ 0
CASETQ
H20, C02, 02
H20 ^ H2 + 1/2 02
1/2 H2 0 - H+ 1/4 02
H20 4- 1/2 02 ^ 2 OH
1/2 02 ^ 0
C02 st CO + 1/2 02
-------
relative to conduction and convection when there is appreciable conversion
between electromagnetic and molecular kinetic energy within the system, which
may occur, for example, if the medium is absorbing and emitting, as is
generally the case in combustion chambers.
The radiant heat transfer model employed under this analysis takes the
form of a discrete-flux model. This model has been employed successfully by
a number of investigators (Refs. 32 through 39) and is described in the
literature (Refs. 35 and 36). A basic hypothesis of discrete flux models is
that the net radiant flux passing through a unit area in each coordinate
direction, is presumed to consist of contributions from collimated forward
and backward fluxes. Thus, in the present two-dimensional axisymmetric flow
case, four fluxes are considered. The configuration of interest is
illustrated in Fig. 1 in terms of one coordinate direction only, representing
the two fluxes, K and L, in the forward and backward axial directions,
respectively. If the monochromatic intensity of radiant energy traveling in
the direction Q is represented by K^, then the total flux of radiant energy
passing forward through the unit plane which is perpendicular to the x-axis
is given by the integral
/X = CD * 0= 7T/2
2?T / K. cosfl sin0d0dX
J X
Similarly, the total flux passing backward in the negative x-direction
through the unit area is
GO
= f
;
B ~
27r f ~ 2 L cos0sin0d0dX (U2)
; *
and the net flux passing through the area is
NET AXIAL RADIATION FLUX = K-L
which has units of energy per unit area per unit time. Similar relationships
can be deduced for I and J fluxes passing forward and backward in the radial
direction, respectively. Thus, the net flux in the radial direction is
written
21
-------
FIG. 1
COORDINATE SYSTEM FOR RADIANT FLUX
dA
22
-------
NET RADIAL RADIATION FLUX = I-J
Employing Eq. (U3) and Eq. (UU), the radiant heat flux term, Qp, appearing in
the energy equation, Eq. (U), and requiring definition, can be written for the
axisymmetric flow case as
In order to derive relationships for the I, J, K, and L fluxes required
to define
reference is made to energy losses and gains incurred by a
collimated flux of radiation as it travels an incremental distance. Energy
losses are incurred by absorption and scattering and energy gains incur by
scattering from the backward flux and by thermal emission from the incremental
element. For an isotropic, gray (absorption coefficient independent wave-
length) medium this transport process can be described by the following
equations for each of the two axisymmetric coordinate directions:
— (rl) «r -(
dr (
— + KaaT
r
— (I+J+K+DJ (U6)
4 1
— (rJ) = r j(Ka + KS )j + KacrT
Ir ( r
4 - — (1+ J + K + L )
dK
dx
Ko-T + —
23
-------
— s (Ko + K*)L-Kfl<»-T4 - — (I+J+K
dx 4
where a is the Stefan-Boltzmann constant, T is absolute temperature, Ka is an
absorption coefficient, and Ks is a scattering coefficient. With the equa-
tions written in this form, it is seen that one-quarter of the total scattered
radiation is presumed to go into each of the four backward and forward
directions. Also, it is to be noted that the equations are consistent with
the solution I = J = K = L= crT1* so that in an equilibrium situation, aF*
represents the flux across any plane.
Before solution to the set of four radiation transport equations,
Eqs. (U6) through (^9), can be attempted, information regarding the mean
absorption coefficients, 1^, and scattering coefficients, Ks, is required.
However, because scattering can be considered to be of secondary importance
in the particle-free gaseous medium under study here, definition of
scattering coefficients becomes unnecessary and, hence, omitted in the
analysis. For an evaluation of the absorption coefficients, reliance is
placed on the considerable amount of data available on emissivity of gases
and mixture of gases (Refs. 37 and 38) from which values for absorption
coefficients can be extracted. For example, conversion between emissivity
data and absorption coefficients can be shown (Ref. 35) to be governed by the
following relationship
where LQ is a characteristic gas-thickness, in this case, taken to be a
representative grid spacing. The indicated extrapolation in Eq. (50) is due
to the fact that the concept of gas emissivity refers to an isothermal gas
and a radiating gas approaches isothermal conditions only in the limit.
Employing Eq. (50), absorption coefficients for air, water vapor, and carbon
dioxide are reported by Cess (Ref. 35).
Once the absorption coefficients are determined to implement the
radiation model, it remains to solve the radiation transport equations,
Eqs. (U6) through (1*9), in the computational procedure. As an alternative,
however, a considerable simplification can be realized by recasting these
equations in alternate form, whereby it will only be necessary to solve two
-------
additional equations instead of the four. For this purpose the following
definitions are made (the scattering coefficient, KS, is retained in the
subsequent mathematical development, for generality):
F = (I + j)/2 (51)
Q = I - J (52)
G = (K + L)/2 (53)
R = K - L (5*0
The manipulation of Eqs. (U6) and (kj) with Eqs. (5l) and (52) leads to
Q. I? *L (55)
(K.+K.+!) «r
and
d(rQ)
dr
2K0r(o-T4-F) + K8r(G-F) (56)
and the combination of Eqs. (55) and (56) can be written
dr Ka+K.+! dr
1 r
Similarly, Eqs. (Ub) and (H9), togetner witn Eqs. (53) and (5*0, can be
manipulated to yield
(57)
-------
R= "2 -£- (58)
KQ+K8 dx
4B. 2 Kn (aT4-G\ + K IF - 0) (59)
dx ° V
and the combination of Eqs. (58) and (59) is written
-^VKO (G-crT4) + i (G-F) (60)
v 2
dx KO+K, dx
Equations (57) and (60) represent the two differential equations of
interest. Solution of these two equations indirectly provides all the
information required to define the heat flux term, Q,R, in the energy equation,
For example, with distribution of F and G provided by solution of Eqs. (57)
and (60) in the iterative solution procedure, Q and R can be determined from
Eqs. (55) and (58), respectively. Then, to compute the desired I, J, K, and
L fluxes necessary for definition of the radiant heat flux term, QR, in
Eq. (^5), Eqs. (5l) and (52), and Eqs. (53) and (5^) can be rearranged so
that the following explicit expressions for the derived fluxes, in terms of
known values of F, Q, G, and R, may be obtained:
I = F + |d (57)
J = F - & (58)
K = G + |R (59)
L = G - |R (60)
26
-------
Equations (57) through (60) thus provide the necessary input information to
Eq. (^-5) in terms of forward and backward fluxes in each coordinate direction,
to enable definition of the radiant flux term in the energy equation.
II-C. BOUNDARY CONDITION
Solution of the governing partial differential equations in the form of
Eq. (ll) is contingent upon specification of boundary conditions for each of
the.dependent variables under consideration. Equation (ll) represents a
system of eight partial differential equations for the eight unknowns:
vorticity, stream function, stagnation enthalpy, fuel mass fraction, nitric
oxide mass fraction, swirl velocity, and the two radiation components
representing radiation fluxes in the axial and radial directions (temperature,
pressure, density, viscosity, air mass fraction, specific heat, and other
mixture properties are deduced from the eight dependent variables by means of
auxiliary relationships, as mentioned previously). Because of the elliptic
form of the governing conservation equations, boundary conditions for each of
the eight variables are required on all the boundaries of the solution
domain. These boundaries are indicated by dashed lines in Fig. 2 where a
schematic of the furnace under study is illustrated. The boundary conditions
for each of the eight dependent variables are given below in terms of normal
and tangential coordinates (n,s) to each of the boundary surfaces. With
exception of the right (downstream) boundary and the centerline, which are
treated separately, other boundary conditions are grouped according to
whether the boundary surface is a solid wall or an injection port. (Note:
The boundary conditions for the F and G radiation functions require special
attention and are given detailed consideration at the end of this section.)
Injection Ports:
e--i/r( *VD ^L_)
ds dn
v// = //ovn ds
H= Zhj mj + KE
m* =1 at fuel port (otherwise zero)
mNO=0
v. = c-tua at air port (otherwise zero)
27
-------
BURNER FLOW MODEL
RECIRCULATION ZONE
ro
oo
7
PRIMARY
AIR
INJECTION
FUEL.
-------
Solid Walls :
Centerline:
Right Boundary:
f I
f r •
* r
an
^ = CONST
dH
an
= o
f
an
am
NO
an
= 0
= 0
= 0
= 0
dn
\> ~ CONST
aH
dn
am f
an
v°
il
an2
dn2
a2H
an 2
= 0
= 0
= o
= o
= o
29
-------
d m,
dn
d rv.
* = 0
an2
The boundary conditions expressed above are formulated to simulate as
near as possible the condition existing in the furnace. At an air port a
tangential velocity component is imparted to the inlet flow of "Ct" times
the axial velocity. The quantity "C^" varies according to the degree of swirl
required to simulate the actual furnace operating conditions. The velocities
at the inlet ports are presumed to be known. Along the solid walls, the
boundary conditions reflect the velocity no-slip conditions, the imperme-
ability to mass transfer, and the local heat balance of the walls. In the
exit plane, the boundary conditions are formulated to present as little
constraint as possible to the flow and the boundary conditions along the
centerline are deduced from symmetry conditions.
The radiation boundary conditions in terms of the functions F and G
remain to be defined and are derived below. Only the boundary conditions for
function F need be formulated as the boundary conditions for the function G
can be deduced by analogy. Thus, considering the radial coordinate direction
only, the boundary conditions in terms of the radial outward flux I and inward
flux J can be written in general form as :
b,i+b2J=o (61)
where bQ, b]_, and b2 are coefficients to be evaluated from known conditions
at the boundary. These conditions fall into one of the following three
categories.
(l) OPAQUE WALL OF PRESCRIBED TEMPERATURE - In this case the radiation
flux at the surface is due to reflection and emission only (transmissivity =
0) and the boundary condition follows as
30
-------
(62)
from which bQ, b]_, and bg can be easily deduced to be
wV
(63)
V-i
The emissivity of the wall ew is known from experiment or obtained from tne
literature as a function of wall temperature and composition.
(2) OUTGOING RADIATION ESCAPES - In this case there is no emission or
reflection and the radiation flux is just equal to the incoming radiation,
J< n. Hence,
and it follows that
b, =o (65)
(3) SYMMETRY AXIS - At a symmetry axis inward and outward fluxes are
equal, consequently,
I'J (66)
31
-------
and, therefore,
b0.o
(67)
The boundary conditions, as developed above, are expressed in terms of I and
J fluxes, however, boundary conditions in terms of F are required. In order
to convert from boundary conditions in terms of I and J to those for F, Eqs.
(51), (52), and (6l) are manipulated to give
(68)
Employing Eq. (55) for Q, Eq. (68) becomes
The analogous expression for G in the axial direction is
:0 (70)
Substitution of the appropriate values of bQ, bls and bg in Eqs. (69) and
(70) is sufficient to define the boundary condition for the functions F and G
on the radial and axial boundary surfaces, respectively.
32
-------
SECTION III
COMPUTATIONAL ANALYSIS
Exact analytical solutions to the time-averaged Navier-Stokes equations
are rare due to their high order and coupled nonlinearity. As a result,
routine solutions to these equations, to a large extent, can only be obtained
by numerical methods. Numerical techniques have undergone significant
advances over the years and computational methods for obtaining numerical
solution to the Navier-Stokes equations are now available (Refs. 39 through
kk). These numerical methods, however, have not proven to be entirely
satisfactory for routine use. Early attempts at solving the Navier-Stokes
equations (Refs. 39 through 1*2), for example, suffered mostly from inability
to maintain computational efficiency and stability over a sufficiently wide
range of flow conditions. In a recent attempt to develop a practical method
of routinely solving the Navier-Stokes equations the method of Gosman, et al.,
(Ref. 1*3) employed a point relaxation method. However, a point relaxation
method is relatively inefficient and, in addition, because of the use of
one-sided differencing, the method may be expected to yield less accurate
solutions than second-order central difference schemes for the same number of
grid points. The recently developed pseudo-time-dependent computational
procedure reported by McDonald, et al., (Ref. UU) solves the vorticity and
stream function equations separately in the inlet region of a duct. The
necessity for evaluating and applying several weighting factors to promote
convergence within the procedure leaves some doubt as to the applicability of
the method to general problem solving. Furthermore, the McDonald procedure,
which applies the alternating-direction-implicit (ADI) algorithm of Peaceman
and Rachford (Ref. U5) to the solution of the vorticity-stream function
equations, employs values of the ADI iteration parameters which are not
permitted to vary in the iteration cycle; a factor which may be expected to
detract from the overall efficiency of the method. To overcome shortcomings
of previous methods, a novel numerical procedure for solving the time-averaged
Navier-Stokes equations has been developed under UARL Corporate-sponsorship.
In the present investigation the UARL computational procedure has been further
developed to solve the equations governing the flow in a furnace and details
of this procedure are presented below.
Ill-A. THE COMPUTATIONAL PROCEDURE
Direct numerical solution of finite-difference approximations to the
time-averaged Navier-Stokes equations in the present study is accomplished by
implementation of a relaxation procedure. In applying relaxation procedures,
33
-------
residuals computed at each point of the flow field are reduced (relaxed) to
a sufficiently small value by successively updating the approximation to the
solution current at that stage within the iterative solution cycle.
(Residuals can be defined as deviations of the current solution from the true
solution of the difference equations at a particular stage of the iteration
cycle.) By analogy with time-dependent methods relaxation procedures may be
classified as either implicit or explicit methods, depending on how the
differencing is done. Explicit methods are computationally simpler than
implicit methods because solution at existing iteration levels can be obtained
from known information at previous iteration levels. With implicit methods
the information necessary to advance the solution is contained implicitly in
the equations and cannot entirely be obtained from known information at the
previous iteration level. Consequently, solution by an implicit scheme is
obtained by the computationally more troublesome simultaneous solution of the
set of finite-difference equations usually either by matrix inversion or
Gaussian elimination. Although computationally more troublesome, implicit
relaxation schemes can be considerably more efficient (as shown below) than
explicit procedures, particularly when practical constraints dictate a large
number of grid nodes.
Relaxation procedures are frequently grouped as either point, line, or
field methods depending on where in the field residuals are relaxed
implicitly. A comparison of the three methods is indicated graphically in
Fig. 3 where a representative five-point molecule employed in the present
analysis is depicted. As can be seen, a solution is obtained at an individual
grid node, along a line of grid nodes or at all nodes throughout the field for
the point, line, and field methods, respectively. Hence, it can be concluded
from the previous discussion that in advancing from the point to the field
methods the computational procedure, although becoming computationally more
troublesome, can be expected to become progressively more efficient.
Recognizing the importance of an efficient computational procedure for
application to furnace problems, wherein large density of grid points might be
required and a large number of governing flow equations might be present, the
computational procedure developed at UARL for solving the time-averaged
Navier-Stokes equations and employed in the present study is a second-order
implicit field relaxation procedure. It employs a novel scheme for derivation
of the finite-difference equations which are solved by the Peaceman-Rachford
(Ref. l&) alternating-direction-implicit (ADI) algorithm. As evidence of the
efficiency of the field (implicit) relaxation procedure in comparison with
point (explicit) relaxation procedures, sample computations of the flow in a
simple axisymmetric chamber are presented in Fig. k. For an 11 x 21 grid
mesh approximately a threefold improvement in computer efficiency is
indicated by the UARL field relaxation procedure in comparison with the point
-------
RELAXATION PROCEDURES
DEFINITIONS: + REPRESENTATIVE VARIABLE
II RESIDUAL
n ITERATION IND*.X
FIG. 3
KJ
1,1+1
•i
1.1"!
KJ
*h—O
IUH
POINT METHOD
(EXPLICIT)
+ (OTHER N-LEVEL TERMS)
LINE METHOD
(IMPLICIT AT GRID POINTS ALONG A LINE)
n+l
• ft(P + (OTHER N-LEVEL TERMS)
FIELD METHOD
(IMPLICIT AT GRID POINTS OVER ENTIRE FIELD)
• R|P > (OTHER N-LEVEL TERMS)
-------
FIG. 4
0.50
0.0002 h-
0.0001 U
0.00005U
0.00002
>L_
VORTICITY RESIDUAL HISTORY
AXISYMMITRIC CHAMBER
11x21 GRID
R»-10
'iiiiiiiiHiiiiiiiniiiiniUHi
OSMAN PRIPCOOi
O - UPWIND DIFFERENCES
O ^CENTRAL DIFFERENCES
20
40
60 80
TIME (SEC)
100
120
140
-------
relaxation procedure of Gosman, et al. In Fig. 5 the results of the same
case are presented for a 31 x Ul grid mesh with resultant improvements in
efficiency by the field methods of one order of magnitude or more.
Representative velocity profiles for these sample cases are presented in
Fig. 6 at two axial locations establishing the identity of the two velocity
predictions.
A detailed description of the computational procedure employed in this
evaluation, particularly in regard to development of finite-difference
equations and the alternating-direction-implicit (ADI) solution procedure
employed to solve the difference equation, is presented in the subsequent
sections.
III-B. DEVELOPMENT OF FINITE-DIFFERENCE EQUATIONS
Any of nonlinear partial differential equations governing the flow in
burners may be represented by the general partial differential equation
developed previously, i.e.,
(ID
where the coefficients have been tabulated in TABLE I for each of the
dependent variables, , under consideration. It should be noted that the
stream function equation serves as an auxiliary equation relating vorticity
and velocity and is not solved alone sequentially in the general manner to be
described. For solution the stream function is solved simultaneously with
the vorticity equation. For the general numerical solution to Eq. (11) a
finite-difference approximation is employed based on a three-point central
differencing formula (see Fig. 3 for description of five-point molecule
employed in this study).
37
-------
FIG. 6
UJ
C/5
Q
55
LU
(E
o
QC
1.0
0.50
0.20
0.10
0.05
0.02
0.01
0.006
0.00
0.001
0.0006
0.0002
0.0001
0.00006
VORTICITY RESIDUAL HISTORY
AXI8YMMITRIC CHAMBIR
31 *41 GRIP
R»-10
-OOIMAN PHIP COOI
'Hiiiiiininiiiiiiiinniinii,
UARLFRIPCODE
I I
I
I I I I
678
TIME (MINI
10 11 12 13 14
-------
FIG. 6
(•» t-0 • aos
a
ee
1.0
0.1
0.0
0,4
0.2
0
0.33
COMPUTED VELOCITY DISTRIBUTIONS
AXIIYMMITHIC CHAMMA
Mc-10
0.20
0.40
an
u/u.
0.80
1.00
MOCIDUNI Of OCMMAN CT AL
IMP. 111
TTJO
U/U,
39
-------
The difference analogue to Eq. (ll) is written:
* i + i,j —" Ci -i,j *i-i,j ~ Ci,j + i *i,j+i
(U (71)
ci,H *i,j-,+ Cii*ii +
where the pseudo-linear coefficients are defined as
(72)
-------
{[Cy2
^r
(75)
\ i
; i r w •» "j '
Cvi (2) d b<^r I
•" a r •" ij ^7 '
V
and the generalized coefficient arrays, Cxi, Cx2, Cyl, and Cy2, which are
formulated to allow for nonuniform grid spacing, are defined according to the
following relationships :
Cxi (I) - Ax,/Ax 2 (Ax, + Ax 2)
Cxl(2) = (Ax2- Ax,)/Ax, Ax2
C xi (3) = Ax2/AX|(Ax,+ Ax2)
C x 2 (I) = 2/Ax2(Ax, + Ax2)
Cx 2 (2) - 2/ Ax, Ax2
C x 2 (3) = 2/Ax, (Ax, -f Ax 2 )
Ul
-------
and
Cyl (I) = Ar,/Ar2 ( Ar,+ Ar2)
Cyl (2) = (Ar2 - Ar,)/Ar, Ar2
Cyl(3) = Ar2/Ar, (Ar,-f Ar2)
Cy2 (I) = 2/Ar2(Ar,+ Ar2)
Cy2 (2) = 2/Ar, Ar2
Cy2(3) = 2/Ar, (Ar, + Ar2)
where
AX. = X:: - X: .
• 'I '-'I (79)
and
(80)
Tb numerically solve Eq. (7l)j an initial guess is made for $ and an
iterative field relaxation procedure is applied in which the governing set of
residual difference equations, written for each node in the field, are solved
implicitly by the ADI algorithm. In this manner, the residuals throughout
the field are relaxed simultaneously at each iteration, a factor contributing
greatly to the efficiency of the method. For development of the residual
difference equations, it is convenient to introduce the residual, RA, which
can be written at the n^k iteration as
-------
(81)
For solution to Eq. (71) it is seen that it is necessary to reduce the
residuals, defined by Eq. (8l), to some sufficiently small value at all points
in the field. For this purpose, appropriate adjustments are made to the
representative dependent variable, <£j_j, at all grid nodes simultaneously
aimed at reducing all the residuals in the field to zero. This procedure is
repeated at each iteration level until convergence is reached. The
iteration formula employed in the iteration process is based on the
following "chain rule" representation:
where the partial derivatives are computed from Eq. (8l) as
n n
;) = " Ci+i,j
n n
"= - c" ,
-'
1*3
-------
Equation (82) represents a set of finite -difference relaxation equations
which are solved to advance the solution to the (n+l) iteration level. The
term jj represents a collection of terms appearing because the coefficients
depend implicitly on <£ and because flow variables other than the dependent
variable may induce changes in the residual. Since the values of ^y are
not initially known, Eq. (82) is solved by a two-step procedure which is a
generalization of the well-known secant method for finding roots of trans-
cendental algebraic equations. During the first, or predictor step, inter-
mediate values of the dependent variable 0ij, denoted 0jt» are computed
from Eq. (82) by neglecting jj altogether. The information gained from the
predictor step is then extrapolated to determine $y which is then employed
to determine (n+l) level iteration values for the dependent variable in the
second, or connector, step. 3hus, writing d ^ = 0jj - 0D4 and employing
Eq. (71 ), the relaxation formula becomes for the predictor step:
n *
With Eq. (81+ ) written at each grid node, a system of linear equations is
obtained which is solved by an ADI algorithm, described below. From solution
of Eq. (8U), the intermediate value of the dependent variables, 0^, are
determined. Once these values are known, intermediate values for the
coefficients in Eq. (71) are computed and then employed for computing the
residual, R0i;j . Sufficient information is now available to compute values for
*ij fa advance) for displacement from the n-level to the (n+l) level by linear
extrapolation of known values of the residual at the n-level and *-levels .
For example, if one considers first a zero displacement in 0 at the n-level,
so that d0?j = 0, it follows that the residual computed from Eq. (71) would'
dust be Hfcj and, hence, by Eq. (8l), ** = R^. if a second displacement
is now considered from the n-level to the *-level, it is known from the
predictor step that the residual is R^,, and the value for *n. for this case
was considered to be zero. Hence, employing a secant extrapolation (see
Fig. 7) the desired values for *» to be employed in the corrector step are
given by
(85)
lilt
-------
FIG. 7
EXTRAPOLATION PROCEDURE
COUPLING AND LINEARIZATION TERM
(SEE EQ. (63) )
O
(/}
UJ
DC
-------
Hence, in the second step of the procedure, values for . . are computed
from the following expression (which is identical to Eq. (8U) with exception
of the *°.j term) derived from Eq. (71 )
"cm,i
n n-H n n-n . n n +1
- " c'-i *i-» " Ci ^.i*1
n " (86)
„ n+i
where Eq. (85) is employed to determine fc
n-H .n n+i n
n n-H n n+t . d " n
+ s-(-).r*
where Eq. (67) is employed to determine
-------
procedure. In the high Reynolds number cases, loss of efficiency was
determined to be due to the fact that coupling and linearization effects
enter only through the corrector step; consequently, at high Reynolds numbers,
when coupling and linearization effects are more significant, the predictor
step is not particularly accurate and the efficiency of the procedure was
found to suffer.
The code was subsequently developed to correct the computational
efficiency problems referred to above by solving the vorticity and stream
function equations simultaneously as a coupled pair. For example, the coupled
solution of the vorticity and stream function equations permits the vorticity
boundary conditions to be incorporated as part of the solution procedure,
thereby eliminating the extrapolation procedure altogether. In addition,
because of the coupling, it became possible to include extra stream function
implicit terms in the vorticity finite-difference residual equations (see
Eq. (8U)). The addition of these terms improved the imbalance between the
predictor and corrector steps of the solution procedure at high Reynolds
numbers because it added coupling terms in the predictor step. Thus, by
solving the vorticity and stream function as a coupled pair, refinements to
the numerical code are implemented which improve the efficiency of the
computational procedure and the range of flow variables to which it can be
applied.
III-C. ADI SOLUTION PROCEDURE
In order to iterate to the solution of the residual finite-difference
analogue to the governing elliptic equations representing the flow in a
furnace (Eq. (?l)) a system of linear algebraic difference equations must be
solved. The difference equations requiring solution are of the general form
(see Eq. (8U) or Eq. (86)):
n n+l n n+i n n + i
" Ci-.,J *i-i.j + C^ *ji " Ci+''i **'.i (87)
n n+' n
where i, j are grid nodes and n is the residual iteration level. Solution of
a pentadiagonal matrix is required to solve Eq. (8?) in this form, and in view
of the special structure of this pentadiagonal matrix the highly efficient ADI
algorithm of Peaceman and Rachford (Ref. MO can be used. The Peaceman-
Rachford algorithm solves the system of equations iteratively and at each
-------
iteration reduces the pentadiagonal matrix to a tridiagonal form which is
easily and efficiently solved by Gaussian elimination. To implement the ADI
procedure, it is convenient to separate the horizontal from the vertical
differencing terms in Eq. (87) in the following two ways
_ n , n-H ^n .n-H „ n , n+i ,n+i
- Ci-U *i-i,i + Cxij*ij - Ci+.,j *!+,
and
n n+i " n-H n n-fi n+i n-fi
Ci,j-i *i.j-i +Crij*ij -
where p may be regarded as a psuedo time step parameter referred to as the
iteration parameter" or the "acceleration parameter" (see section on
ITERATION PARAMETER). The Peaceman-Rachford ADI iteration method solves the
set of equations represented by Eq. (87) by solving Eqs. (88) and (89)
alternately by sets of rows or columns until convergence is reached
Accordingly, Eqs. (88) and (89) are written
n-H
M+l
(90)
M+2
(9D
-------
where M represents the ADI iteration index. Thus, Eq. (90) is employed to
advance the solution to the (M+l) iteration level followed by Eq. (91) which
is used for advancement from the (M+l) to the (M+2) iteration level. Values
at the (M-l) iteration level are presumed known either from the previous
iteration, or from an initial guess. Solution to Eqs. (90) or (91) is easily
and efficiently obtained by employing an algorithm for solving tridiagonal
matrices.
Convergence is determined in the ADI procedure by the following test
n+i.. M . n+t M-i
J-lft] Ji < c (92)
~
i <
J ~
I MAX
which dictates the maximum fractional change of 0 which is allowed in the
field must not exceed a prescribed value, el5 which is typically in the range,
0.001 to 0.005. An alternate criterion is sometimes employed when the value
of 0 at a particular node is very small relative to surrounding values. In
this case the following formula is employed
n+L M „ n-H, M-I
(93)
, <€
" '
MAX
- *2
MAX
where the change in the variable has been divided by the maximum value in the
field. Typically, e2 is set in the range 0.0001 to 0.00005.
m-D. ITERATION PARAMETERS
The rate of convergence of the ADI procedure is very sensitive to the
choice of acceleration parameters, p, and improper values will very likely
result in divergent behavior. Unfortunately, little is known of the correct
values for the acceleration parameter to be applied to the coupled nonlinear
equations requiring solution for the present problem and information must be
inferred from knowledge gained in simpler linear systems. The parameters
employed by Peaceman and Rachford for solving Laplace's equation are defined
according to
-------
-= b (a/b)
where a and b represent minimum and maximum eigenvalues of the tridiagonal
matrix, and m refers to the number of parameters to be employed. The rate of
convergence of ADI methods can be appreciably increased by the application of
several iteration parameters which are used successively in cyclic order.
Estimates of optimum choice of the number of parameters to be used can be
determined (Ref. 1^6) by finding the smallest integer m such that
2m
(8) < 4- (95)
where
8= S~T- I = 0.414 (96)
Ihis procedure for evaluating iteration parameters and determining the
optimum number to be employed was not found to be completely satisfactory in
the computational procedure, for cases where convection and source terms in
the governing differential equation were large relative to the diffusion
terms. Accordingly, a theoretical analysis was undertaken wherein a simple
time -dependent linear analogue of the governing differential equations was
investigated and compared against the analytical solution to the diffusion
equation, for which the parameters are known from Eq. (9U), for the purpose of
obtaining some insight into appropriate iteration parameters to be employed.
Consider the following linear differential equation
50
-------
For boundary and initial conditions specified as
(98)
solution becomes
r X- .X<\
^i *X, / _A _j . I ''O
^p ^ ^p ^ s v^D ***^p I crop I •^•"•••"•••'^^••••^ i
r 2 / (99)
exp [-a (x-x0) - a a (t-t0)J + f (t-t0)
The iteration parameters can be shown to be related to the time step in the
following manner (by differencing Eq. (97) and comparing against Eq. (88))
(x"x°)2 (100)
(t-t0)
Employing Eq. (100), Eq. (99) can be written
A<£ = A«#>0 ERF v^, exp [-a Ax
where
=<#»!- ^>0 (102)
Ax = x-x
0
-------
Three cases can now be considered in order to provide insight to appropriate
values for p by considering the governing equations to be either diffusion
dominated, convection dominated, or source dominated.
CASE 1. DIFFUSION DOMTHATED - If the diffusion terms in the governing
differential equation are much larger than the convection source terms, Eq.
(97) effectively reduces to the diffusion equation
(103)
a at
with solution
A*= A*0ERF y ^ (ioU)
where satisfactory values for p are given by Eq. (9*0 as determined by
Peaceman and Rachford:
CASE 2. CCHVECTION DOMINATED - For the case where the governing
equations are convection dominated the governing equation, Eq. (97) is
written
2a -. = no-
ax a at (105)
with solution
A00 exp [-a Ax (1+
(106)
-------
which is solved for p:
- a(QAx)2
CASE 3. SOURCE TERM - For the source dominated case, the governing
equation, Eq. (97), is written
with solution
— (109)
which when solved for p becomes
f Ax2
p -- (110)
By comparing Eq. (97) against the general elliptic partial differential
equation being solved under this investigation the following relationships can
be deduced
a = v (111)
53
-------
_L r L + _!_
2 L brc ax b^r ar
f=
By employing Eqs. (ill) through (113) it is possible to determine whether the
governing elliptic partial differential equation is diffusion, convection, or
source term dominated. Depending on the case, Eq. (9*0> (107), or (110),
respectively, is employed to evaluate the acceleration parameter. To solve
for p in these three equations, however, requires knowledge of A00 and
A0(AX is taken to be a representative grid spacing). For this purpose,
Eq. (l(A) can be employed (p and ot are known for the diffusion dominated case
from Eqs. (9^) and (ill), respectively). Thus,
(U»0
DIFFUSION DOMINATED
Equation (ll!+) provides sufficient information to determine p in Eqs. (9U) or
(107); however, to determine p in the source dominated case, Eq. (110), the
ratio A0/A0Q is not sufficient since values for both A0Q and A0 are
required. In this case, the assumption is made that across a grid spacing a
one percent change in the function is allowable in the given time step. Thus,
= 0.01 (115)
and A00 can be computed from Eq.
-------
III-E. ITERATIVE PROCEDURE AM) ORDER OF SOLUTION
Equations solved by the computational procedure under this investigation
include: vorticity, stream function, stagnation enthalpy, fuel mass fraction,
swirl velocity, two radiant heat transfer functions, and nitric oxide mass
fraction. Pressure, temperature, density, viscosity, other species concentra-
tions, and other flow variables are deduced from these dependent variables.
While solving these equations, in order to remain within a UNIVAC 1108
computer storage limitations, it was necessary to transfer data on and off
auxiliary storage facilities. Consequently, the above equations were solved
within the program in three sets, comprising groups of five or less. The
first group solved in the computational procedure consists of vorticity,
stream function, stagnation enthalpy, and fuel mass fraction, in that order.
The second and third groups consist of the following equations: swirl
velocity and two radiant heat transfer equations in the second group, and the
nitric oxide mass fraction equation in the third group.
Starting with an initial guess for all flow properties, a cycle of
iterative sequence of operations begins with computation of residuals for
each equation to be solved (with exception of stream function equation) . For
the predictor step, iteration is initiated by solving Eq. (8U) for all the
dependent variables by employing the ADI solution procedure. When convergence
is reached for all the equations, the secant extrapolation is implemented to
evaluate corrections for nonlinearity and coupling which were neglected in
the predictor step, but which are required for execution of the corrector
step. The corrector step is then implemented by solving Eq. (85) by the ADI
method to determine updated values for all the dependent variables. The
nonlinearity and coupling correction is introduced through the variable, 4>ij .
As stream function is considered an auxiliary equation, it is coupled with the
vorticity equation and solved simultaneously with the vorticity in the
corrector step as in the predictor step by employing Eq. (8U). Auxiliary
relationships are then invoked to update pressure, temperature, density,
viscosity, and other flow variables before start of the next iteration cycle.
Convergence is determined by the following test, as in the ADI procedure
< € (116)
n
i MAX MAX
where e is of the order of 0.001 or less.
-------
SECTION IV
RESULTS AND DISCUSSION
To validate the computational procedure comparisons were made with the
nonreacting turbulent swirling jet flow of Chigier and Chervinsky (Ref. 14-6).
In particular, it was recognized that for low and moderate degrees of inlet
swirl Lilley (Ref. 20) had obtained very satisfactory comparisons with
Chigier's data using a finite-difference boundary layer procedure together
with the turbulence model detailed in II .B.I. As was mentioned previously,
the appearance of a recirculation zone and the occurrence of large radial
velocities with high degrees of inlet swirl velocity, made calculations using
the boundary layer equations unacceptable in this flow regime. However, it
was felt that predictions made using the FREP code should agree with Lilley's
predictions (and the data) at low to moderate inlet swirl velocities. These
comparisons were performed under a Corporate-sponsored program but were
reported here in detail in view of their importance in validating the FREP
code. To demonstrate the applicability of the FREP code to predict the flow
field within a representative furnace, further comparisons were performed
under the present contract with the data of Larson and Shoffstall (Ref. U?).
Larson, et al., measured the mean velocity field in the turbulent swirling
jet of a furnace both with and without combustion.
IV-A. COMPARISON WITH THE DATA OF CHIGIER, ET AL. (REF. 1|6)
In Ref. U6 Chigier, et al., report an experimental and theoretical
investigation of the turbulent swirling jet. Subsequently, Lilley (Ref. 20)
compared Chigier's measurements to predictions of the jet made using a
turbulence model which he constructed together with a finite-difference proce-
dure to integrate the boundary layer equations. For low and moderate degrees
of swirl velocity the comparisons were very satisfactory and, consequently, it
was decided to perform a similar set of comparisons with Chigier's data using
the FREP code. Since the turbulence model was the same in both comparisons,
good agreement between FREP code predictions and Lilley's predictions and, of
course, the data, was expected. However, it became apparent in setting up the
cases that detailed initial profiles at the inlet plane were not available to
start the calculations and the subsequent behavior, in particular, the decay
of the potential core region, was very sensitive to the choice of the initial
profiles. In addition, it was not clear how Lilley had started his calcula-
tions, or what special procedure he used to calculate the potential core
decay. As a result of these considerations, an alternate comparative scheme
was adopted based on the similairty hypothesis.
-------
Basically, the similarity hypothesis argues that for turbulent flow far
from the orifice the jet properties would become invariant (similar) with
streamwise direction when expressed in a certain coordinate system. Chigier,
et al., (Ref. U6) investigated the similarity hypothesis and observed that
when expressed | , T| coordinates, where £ is x + a, a being the streamwise
distance from the effective origin of the flow and Tj being r/(x+a), then mean
profile similarity was obtained after approximately U to 10 jet diameters of
streamwise development. Accordingly, comparisons are presented in terms of
the similarity variables for three degrees of inlet swirl velocity correspon-
ding to a weak, intermediate, and strong coupling between the axial and
rotational velocity fields. In all cases a reference position x^ was adopted
downstream of which the flow could be considered to be similar. The subse-
quent decay of the axial velocity maximum, Um, was expressed as a fraction of
the axial velocity maximum at the reference location Umi, and similarly for
the swirl velocity maximum Wm. The velocity decay in the streamwise direction
is presented in terms of the ratio |i/| for the axial component and (£i/§)2
for the swirl component since Chigier indicated that on the basis of similarity
the decay would be linear when expressed in this manner. The comparisons
between the predictions of the FREP code and the data are presented in Figs.
8 through 16. A 28 x 20 grid was used in the computation and it can be seen
that a very acceptable level of agreement between prediction and measurement
has been obtained. For the high swirl case only at the last measuring station
was the flow development approximately similar. However, it did appear for
this high swirl case that the potential core region was negligible and that a
comparison in real variables might be warranted. Such a comparison is shown
in Fig. 17, where the axial development of the profile of the velocity vector
is presented. Obviously, further refinements to the inlet profiles could have
improved the comparison between prediction and measurement but that qualita-
tively the agreement is quite good. It was also evident from the results
that with this high degree of swirl velocity the radial velocity component
was quite large, nonnegligible, and insofar as the measurements are concerned,
unknown.
IV-B. COMPARISON WITH THE DATA OF D. H. LARSON AND D. R. SHOFFSTALL (REF. 1*7)
Under EPA Contract No. 68-02-0216 with the Institute of Gas Technology,
the flow field within a gas fired furnace has been mapped for both hot- and
cold-flow conditions. With cold-flow two degrees of inlet swirl velocity
were explored, designated minimum and high swirl. For the high swirl case a
large recirculation zone formed downstream of the inlet. The hot-flow data
were obtained by burning natural gas.
57
-------
IV-B.l Cold-flow comparisons
The comparisons between the predictions and the IGT measurements were,
like the Chigier comparisons, also influenced by the unknowns in the inlet
flow conditions. However, the details of the inlet flow are much more
likely to influence the low swirl calculations than the high swirl compari-
sons, since with large swirl velocities there is very rapid intense mixing
which tends to quickly eradicate the details of the inlet profiles. In this
regard it is perhaps worth noting that the calculation domain could have been
extended into the inlet pipe upstream of the influence of the expansion into
the'burner. At such an upstream location it would be presumably much easier
to specify a set of reasonable pipe flow conditions from which to initiate
the calculation. In view of the computer logic involved, for the present,
all the calculations have been initiated at the burner entrance. In Fig. l8
a comparison between the measured and predicted axial velocity profiles in the
turbulent jet with the minimum degree oif inlet swirl velocity is given. In
general, a very good degree of qualitative and quantitative comparison between
measurement and prediction is shown, with the theory indicating slightly too
slow a decay of the profile at the third and fourth measuring stations. In
Fig. 19 a similar comparison is shown for the highly swirling case and it can
be seen that although the inlet conditions have not been matched precisely
the profiles in the recirculation zone are predicted remarkably well. In
Fig. 20 the theoretical streamlines for the highly swirling case are
presented and from it it is evident that close to the inlet the flow is
changing very rapidly and in view of the streamline curvature the radial
velocity component is obviously quite large. It should perhaps be mentioned
at this point that the upper boundary condition applied for all of the
comparisons shown so far was that the variable second derivative with respect
to the radial direction was zero. This very weak boundary condition was
applied in favor of going out a very large distance to a solid wall or as in the
Chigier case to some distant location where the vorticity was zero.
IV-B.2 Hot-Flow Comparisons
In view of the difficulties associated with making measurements in hot
flows only a very limited amount of aerodynamic data was obtained by Larson
and Shoffstall (Ref. Uy) for the burner considered here. In making the
theoretical calculations it was assumed that the same inlet profile of fuel
and air would be obtained as in the comparable cold flow case and the magni-
tude of the inlet fuel and air velocities were simply scaled to give the
quoted fuel flow rate and overall stoichiometry. Insofar as the swirl
component was concerned some doubt was expressed as to the behavior of the
aluminum swirl generating blocks under hot flow conditions. As a consequence
of this two calculations were performed, one with the same relative swirl
distribution as in the cold flow case and the second with the swirl component
58
-------
increased by 1.333. With the same relative swirl distribution, only a very
small recirculation zone was obtained in the hot flow case; however, in the
hot case the flow in the region where the cold flow recirculation zone had
existed was of very low velocity, although not quite reversed. In the second
case where the swirl distribution was scaled by 1.333 a much larger
recirculation zone, similar in size to the cold flow recirculation zone, was
obtained. Comparisons with the data have been made only for the higher swirl
case.
Before proceeding to discuss the hot flow results, two points
concerning the hot flow calculations should be made. The first of these
was that, it was found beneficial to start the calculation from a cold flow
solution and under-relax the temperature field during the first few
iterations around the stream function - vorticity and the stagnation enthalpy
equations. The under-relaxation consisted of taking some fraction of the
computed local temperature together with one minus this fraction of the
previous local temperature, to determine the current local temperature. A
new density was then computed to be consistent with the under-relaxed local
temperature. The fraction of the new temperature increased smoothly as the
interation proceeded until after five outer iterations round the entire
system, the degree of under-relaxation was negligible.
As a second point it should be noted that in solving for the diffusion
of nitric oxide, the nonequilibriurn species production term in the governing
equations, ri? in Eq. (2), was nonnegligible and highly nonlinear. The
incorporation of this r.j_ term, or indeed any source term d^ containing the
dependent variable, into the present computational framework was straight-
forward. To account for such a term, it is necessary in the development
of the change in residual, Eq. (82), to allow for the contribution due to
the rate of change of the source term d^/b^, with the dependent variable <£.
This contribution is added to Cn, in Eq. (83), and the resulting linear
system solved using FREP. At this point a new rate of change of d^/b^ with
dependent variable <£ is evaluated throughout the field and, if necessary,
the linear system may be resolved. The process is repeated until the computed
dependent variable does not change to within a specified tolerance. This
process is obviously simply a variation of the Newton-Raphson method for
solving nonlinear equations. For the calculation of the nitric oxide
concentrations in the cases examined to date, the technique proved highly
efficient.
59
-------
FIG. 8
A COMPARISON BETWEEN PREDICTED AND MEASURED AXIAL DECAY OF
THE VELOCITY MAXIMUM IN A TURBULENT SWIRLING JET
MODERATELY SWIRLING JET, S = 0.134
g MEASUREMENTS OF CHIGIER (REF.46)
——" PREDICTION
1.0
0.8
0.6
<
DC
I M
Ul
0.2
0.2
0.4
0.6
INVERSE OF NORMALIZED AXIAL LENGTH.
1.0
60
-------
FIG. 9
A COMPARISON BETWEEN PREDICTED AND MEASURED AXIAL DECAY OF
THE SWIRL VELOCITY MAXIMUM IN A TURBULENT SWIRLING JET
MODERATELY SWIRLING JET, S = 0.134
• MEASUREMENTS OF CHIGIER (REF.46)
PREDICTION
INVERSE OF NORMALIZED AXIAL LENGTHip-L)
61
-------
FIG. 10
COMPARISON BETWEEN THEORY AND MEASUREMENT OF
THE MEAN VELOCITY PROFILE IN A TURBULENT SWIRLING JET
MODERATELY SWIRLING JET. S =0,134
ODA09 MEASUREMENTS OF CHIGIER (REF.46)
——— PREDICTION
1.0
D O A
O D
_L
_L
0.2 0.3 0.4
NORMALIZED RADIAL LOCATION, r/(x + a)
0.5
62
-------
FIG. 11
A COMPARISON BETWEEN PREDICTED AND MEASURED AXIAL DECAY
OF THE VELOCITY MAXIMUM IN A TURBULENT SWIRLING JET
INTERMEDIATELY SWIRLING JET, S = 0.416
A MEASUREMENTS OF CHIGIER (REF.46)
—_ PREDICTION
E
CC
>
o
o
LU
0.2 0.4 0.6
INVERSE OF NORMALIZED AXIAL LENGTH,
0.8
63
-------
FIG. 12
A COMPARISON BETWEEN PREDICTED AND MEASURED AXIAL DECAY
OF THE SWIRL VELOCITY MAXIMUM IN A TURBULENT SWIRLING JET
INTERMEDIATELY SWIRLING JET, S = 0.416
A MEASUREMENTS OF CHIGIER (REF.46 )
' ' PREDICTION
INVERSE OF NORMALIZED AXIAL LENGTH,/^Llf
\x + a
-------
FIG. 13
THE COMPARISON BETWEEN THEORY AND MEASUREMENT OF
THE MEAN VELOCITY PROFILE IN A TURBULENT SWIRLING JET
INTERMEDIATELY SWIRLING JET, S = 0.416
O D A 0 • MEASUREMENTS OF CHIGIER (REF.46)
—— PREDICTION
g"
o
o
LU
0.1 0.2 0.3
NORMALIZED RADIAL LOCATION,
0.5
r
(x + a)
-------
FIG. 14
A COMPARISON BETWEEN PREDICTED AND MEASURED AXIAL DECAY
OF THE VELOCITY MAXIMUM IN A TURBULENT SWIRLING JET
HIGHLY SWIRLING JET, S = 0.640
O MEASUREMENTS OF CHIGIER (REF. 46)
_____ PREDICTION
cf
cc
>~
O
O
111
INVERSE OF NORMALIZED AXIAL LENGTH.
(x; + a)z
(x + a)
66
-------
FIG. 15
A COMPARISON BETWEEN PREDICTED AND MEASURED AXIAL DECAY
OF THE SWIRL VELOCITY MAXIMUM IN A TURBULENT SWIRLING JET
HIGHLY SWIRLING JET, S = 0.640
O MEASUREMENTS OF CHIGIER (REF. 46)
PREDICTION
0.2 0.4 0.6
INVERSE OF NORMALIZED AXIAL LENGTH,
-------
FIG.16
COMPARISON BETWEEN THEORY AND MEASUREMENT
OF THE MEAN VELOCITY PROFILE IN A TURBULENT SWIRLING JET
HIGHLY SWIRLING JET, S = 0.640
• MEASUREMENTS OF CHIGIER (REF. 46)
PREDICTION
1.0
0.8
0.6
O
$
oc
O
3 0.4
ut
0.2
0.1
0.2
0.3
0.4
NORMALIZED RADIAL LOCATION,
(x + a)
0.5
68
-------
RADIAL DISTRIBUTION OF VELOCITY VECTOR
HIGHLY SWIRLING JET, S = 0.600
Q MEASUREMENTS OF CHIGIER (REF. 46)
—— PREDICTION
x/d = 4.1
x/d = 6.2
x/d = 8.3
x/d = 10.0
50
50
0 50
VELOCITY (FPS)
50
50
-------
RADIAL DISTRIBUTION OF THE AXIAL VELOCITY
MINIMUM SWIRL
O MEASUREMENTS OF LARSON ET AL (REF 47)
PREDICTION
1.0
0.8
0.6
0.4
0.2
8 o
a.
_i
^ -0.2
a
DC
-0.4 -
-0.6 -
-0.8 —
-1.0
7.6 CM STATION
17.8 CM STATION
La L
30.5 CM STATION
63.5 CM STATION
50
100
0 50
AXIAL VELOCITY (FPS)
50
50
O
oo
-------
RADIAL DISTRIBUTION OF AXIAL VELOCITY
HIGH SWIRL
O MEASUREMENTS OF LARSON ET AL (REF. 47)
PREDICTION
1.0
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1.0
LL
O
2
_i
<
Q
2.5 CM STATION
6>
O 7.6 CM STATION
_L
17.8 CM STATION
30.5 CM STATION
O
O
O
O
_Q
50
100
0 50
AXIAL VELOCITY (FT/SEC)
50
50
P
CO
-------
THEORETICAL STREAMLINES FOR IIT SWIRL BURNER
SWIRL NUMBER -0.8
COLD FLOW
T]
P
ro
o
-------
Turning now to the hot flow results the predicted distributions of axial
velocity are compared with the data in Fig. 21. It is immediately evident
that the comparison with data is poor and that obviously the assumed inlet
conditions were considerably different from those actually pertaining in
the experiments. Time did not permit further (arbitrary) changes in the inlet
conditions in order to better match the flow profiles at the 12.7 cms station.
Evaluation of the hot flow predictive capabilities of the present procedure
must, therefore, await the availability of data where the aerodynamic inlet
conditions are known in detail. In Fig. 22 the predicted swirl velocities
are compared with the measured data and, as in the case of the axial velocity
profiles, the comparisons are poor. In particular, it is apparent that
the degree of inlet swirl must have been much higher and much more localized
than assumed in the calculations on the basis of the cold-flow measurements.
In Fig. 23 the predicted radial profiles of temperature are compared with
the data and in this case the predicted profiles are uniformly several
hundred degrees Faranheight higher than those measured. Again, some of
this discrepancy can be attributed to the unknown inlet conditions but
almost certainly some of the lower measured temperatures must be the result
of the energy extraction taking place at the furnace wall, and this energy
extraction, for simplicity, was neglected in the present calculations. The
indications are, therefore, that radiative transport was not entirely
negligible, as was assumed, in performing the present calculations.
Finally, in Fig. 2U, the profiles of nitric oxide are compared with the
data and here, at least at the 12.7 cms station, the agreement is fair. As
it happens, the agreement at the 30.5 cms location is not nearly as good
and to some extent the Zeldovitch mechanism's known predisposition to under-
predict nitric oxide production is balanced against the over predicted
temperature levels. However, the maximum local production of nitric oxide
did not exceed 10 percent of the computed equilibrium nitric oxide
concentration.
In summary, little can be said of the results of the comparison with
data other than it should and will be redone for some case where the inlet
conditions are very carefully monitored. Finally, there is a suggestion that
radiative transport effects cannot be neglected in attempting to predict
pollutant levels.
73
-------
RADIAL DISTRIBUTION OF AXIAL VELOCITY
HIGH SWIRL
O MEASUREMENTS OF LARSON ET AL
-------
RADIAL DISTRIBUTION OF TANGENTIAL VELOCITY
HIGH SWIRL
O MEASUREMENTS OF LARSON ET AL (REF. 47)
— PREDICTION
-0
z
o
ro
1.0
0.8
0.6
0.4
0.2
8 o
Q_
-0.2
D
-0.4
-0.6
-0.8
-1.0
3.3 CM STATION
50
o o
CD
12.7 CM STATION
J_
100 0 50
TANGENTIAL VELOCITY (FT/SEC)
30.5 CM STATION
50
P
ro
NJ
-------
RADIAL DISTRIBUTION OF TEMPERATURE
HIGH SWIRL
O MEASUREMENT OF LARSON ET AL (REF. 47)
— PREDICTION
z
8
to
0)
I
w
1.0
0.8
0.6
t
"J n
2 °
I
<
DC
-0.4
-0.6
-0.8
-1.0
J L
P°
o
12.7 CM STATION
I I
30.5 CM STATION
J L
10
15 20 25
30 35 10 15 20 25
TEMPERATURE, 102 °F
30
35 10
15 20 25
30 35
40
P
NJ
CJ
-------
RADIAL DISTRIBUTION OF NITROGEN OXIDE
HIGH SWIRL
O MEASUREMENT OF LARSON ET AL (REF. 47)
PREDICTION
z
O
IO
to
01
I
1.0
0.8
0.6
0.4
P
u_
z 0.2
O
§ °
D-
_l
5 -0.2
Q
<
cr
-0.4
-0.6
-0.8
-1.0
3.3 CM STATION
_L
12.7 CM STATION
tr
o
o
o
o
O
O
30.5 CM STATION
50
100 0
50 100
NITROGEN OXIDE ppm
50
100 0
50
P
ro
-P..
-------
SECTION V
CONCLUSIONS
The aim of the present study was to evaluate the predictive capability
of a computational procedure for solving the time-averaged Navier-Stokes
equations for the flow within a combustion device. Turbulent swirling flow
with coupled chemistry and radiative heat transfer was allowed for and the
following conclusions were reached:
1. Routine calculations using a novel finite-difference procedure for
computing solutions of the time-averaged Navier-Stokes equations can be made
for representative furnaces. Such calculations can be performed for about
six hundred grid points to within a tolerance of one third of a percent in a
residual for approximately ten to thirty minutes of UNIVAC 1108 computer run
time, depending on the number of ancillary equations to be treated.
2. For nonreacting flows at low to moderate degrees of inlet swirl
velocity very satisfactory comparisons between measured and predicted flow
fields are obtained using a fairly simple turbulence model based on Prandtl's
mixing length hypothesis.
3. For nonreacting flows with a high degree of inlet swirl a
recirculation zone may be obtained downstream of the inlet. Even with the
simple turbulence model the comparisons between measured and predicted flow
fields in the region of the recirculation zone are in quite good agreement.
k. It has generally been accepted that in the recirculation zone itself
that the boundary layer type of approximations are not valid, thus requiring
the full Navier-Stokes equations to describe the mean flow behavior in this
region. However, it seems clear from both prediction and measurement that
the radial component of velocity is also large whenever the swirl velocity
becomes appreciable even if a recirculation zone is not obtained. Thus the
boundary layer approximations do not appear valid in flows with appreciable
amounts of swirl velocity and such flows require the full Navier-Stokes equa-
tions to accurately describe the mean flow behavior.
5. For the reacting flows compared to in the present study little
can be said about the predictions in view of the lack of information available
on the inlet flow conditions. It can be said, however, that hot swirling
flow calculations with equilibrium hydrocarbon combustion can be performed
for a very acceptable computational cost employing a reasonable mesh and a
fairly stringent convergence criterion by means of the UAEL FREP code.
78
-------
6. Highly nonlinear chemical kinetic reactions exemplified by the
extended Zeldovitch mechanism for the production of nitric oxide can be
incorporated into the present computational framework.
7. Although relatively untested in the present evaluation, it is felt
that further improvements in the turbulence modeling, radiative transport,
turbulence effect on chemical kinetics, and a better understanding of the
chemical hierarchy will be required before reliable trends can be obtained
by a prediction technique.
79
-------
SECTION VI
REFERENCES
1, Beer, J. M. and N. A. Chigier: Stability and Combustion Intensity of
Pulverized Coal Flames - Effect of Swirl and Impingement. Journal of
the Institute of Fuel, December 1969.
2. Beer, J. M. and W. Leucker: Turbulent Flames in Rotating Flow Systems.
' Paper No. Inst. F-NAFTC-7, North American Fuel Technology Conference,
Ottawa, Canada, 1970.
3. Beer, J. M. and J. B. Lee: The Effects of Residence Time Distribution
on the Performance and Efficiency of Combustors. The Combustion
Institute, 1965, pp. 11&7-1202.
k. Marteney, P. J.: Analytical Study of the Kinetics of Formation of
Nitrogen Oxide in Hydrocarbon - Air Combustion. Combustion Science and
Technology, Vol. 1, 1970, pp. 37-^5.
5. Fletcher, R. S. and J. B. Heywood: A Model for Nitric Oxide Emission
from Aircraft Gas Turbine Engines. AIAA Paper 81-123, 1971.
6. Hammond, D. C. (Jr.) and A. M. Mellor: Analytical Predictions of
Emissions From and Within an Allison J-33 Combustor. Combustion Science
and Technology, Vol. 6, 1973, pp. 279-286.
7. Hammond, D. C. (Jr.) and A. M. Mellor: Analytical Calculations for the
Performance and Pollutant Emissions of Gas Turbine Combustors.
Combustion Science and Technology, Vol. k, 1971, pp. 101-112.
T»
8. Roberts, R., L. D. Aceto, R. Keilback, D. P. Teixeira, and J. M. Bonne 11:
An Analytical Model for Nitric Oxide Formation in a Gas Turbine
Combustion Chamber. AIAA 9th Aerospace Sciences Meeting, Paper No.
71-715, New York, New York, 1971.
9- Mosier, S. A., R. Roberts, and R. E. Henderson: Development and
Verification of an Analytical Model for Predicting Emissions from Gas
Turbine Engine Combustors During Low Power Operation. ^Ist Meeting
Propulsion and Energetics Panel of AGARD, 1973.
8~0
-------
10. Edelman, R. and C. Economos: A Mathematical Model for Jet Engine
Combustor Pollutant Emissions. AIM Paper No. 71-71U, 1971.
11. Gosman, A. D., W. M. Pun, A. K. Runchal, D. B. Spalding, and
M. Wolfshtein: Heat and Mass Transfer in Recirculating Flows.
Academic Press, New York, New York, 1969.
12. Anasoulis, R. F.: Computations of the Flow in a Combustor. United
Aircraft Research Laboratories Report K110885-1, November 1971.
13 • Bird, R. B., W. E. Stewart, and E. N. Lightf oot: Transport Phenomena.
John Wiley & Sons, Inc., New York, New York, 1960.
lU. Launder, B. E. and D. B. Spalding: Turbulence Models and Their
Application to the Prediction of Internal Flows. Imperial College of
London, Technical Engineering Department Report No. TM/TN/A/18, 1971.
15. Prandtl, L.: Bericht Uber Untersuchingen Zur Ausgebildeten Turbulenz.
ZAMM, Vol. 5, 1925, p. 136.
16. Patankar, S. V. and D. B. Spalding: Heat and Mass Transfer in
Boundary Layers. Intertext Books, London, England, 1970.
17. Maise, G. and H. McDonald: Mixing Length and Kinematic Eddy Viscosity
in a Compressible Boundary Layer. AIAA Journal, Vol. 6, 1968, pp. 73~80.
18. McDonald, H. and F. J. Camarata: An Extended Mixing Length Approach
for Computing the Turbulent Boundary Layer Development. Proceedings
of the AFOSR-IFP-Stanford Conference on Boundary Layer Prediction, 1968.
19. Williamson, J. W.: An Extension of Prandtl's Mixing Length Theory.
Applied Mechanics and Fluids Engineering Conference, ASME, June 19&9-
20. Lilley, D. G.: Prediction of Inert Turbulent Swirl Flows. AIAA Paper
No. 72-699. AIAA 5th Fluid and Plasma Dynamics Conference, June 1972.
21. Beer, J. M. and N. A. Chigier: Combustion Aerodynamics. John Wiley &
Sons, Inc., New York, New York, 1972.
22. Walz, A.: Boundary Layers of Flow and Temperature. The M.I.T. Press,
Cambridge, Massachusetts, 1969.
81
-------
23. Lavoie, G. A., J. B. Heywood, and J. C. Keck: Experimental and
Theoretical Study of Nitric Oxide Formation in Internal Combustion
Engines. Combustion Science and Technology, Vol. 1, 1970, pp. 313-326.
2k. Zeldovitch, Ya. B., P. Ya Sadounikov, and D. A. Frank -Kamenet skll :
Oxidation of Nitrogen in Combustion. Academy of Sciences of USSR,
Institute of Chemical Physics, Moscow -Leningrad,
25. Bowman, C. T. and D. J. Seery: Investigation of NO Formation Kinetics
in Combustion Process: The Methane -Oxygen -Nitrogen Reaction, Emissions
from Continuous Combustion Systems. Plenum Publishing Company,
New York, New York, 1972.
2.6. Caretto, L. S., L. H. Muzio, R. T. Sawyer, and E. S. Starkman: The Role
of Kinetics in Engine Emission of Nitric Oxide. Combustion Sciences and
Technology, Vol. 3, 1971.
27. Baulch, D. L., D. D. Drysdale, D. G. Home, and A. C. Lloyd: Critical
Evaluation of Rate Data for Homogeneous Gas Phase Reactions of Interest
in High -Temperature Systems. Report No. U Department of Physical
Chemistry, Leeds University, United Kingdom, December 1969.
28. Campbell, I. M. and B. A. Thrush: Reactivity of Hydrogen to Atomic
Nitrogen and Atomic Oxygen. Trans. Faraday Soc. 6U, Part 5> 1968,
pp. 1265-127U.
29. Brinkley, S. R.: Computational Methods in Combustion Calculations.
Combustion Processes, Section C. High Speed Aerodynamics and Jet
Propulsion, Vol. 2, B. Lewis, R. N. Peace, and H. S, Taylor, Eds.,
Princeton University, Princeton, New Jersey, 1956.
30. Brinkley, S. R.: Calculation of the Thermodynamic Properties of Multi-
Component Systems and Evaluation of Propellant Performance Parameters.
Proceedings of the First Conference on Kinetics, Equilibrium and
Performance of High Temperature Systems, A. S. Bahn and E. E. Zukoski,
Eds. The Combustion Institute, 1960, pp. 7^-8l.
31. Stull, D. R. and H. Prophet: Janaf Thermochemical Tables. National
Bureau of Standards, U. S. Department of Commerce, NSRDS-NB537, 2nd
Edition, June 1971.
82
-------
32. Chen, J. C.: Simultaneous Radiative and Connective Heat Transfer in an
Absorbing, Emitting and Scattering Medium in Slug Flow Between Parallel
Plates. AlChE Journal, Vol. 10, No. 2, March
33. Viskanta, R. and R. J. Grosh: Heat Transfer by Simultaneous Conduction
and Radiation in an Absorbing Medium. Journal of Heat Transfer, ASME,
February 1962.
3k. Larkin, B. K. and S. W. Churchill: Heat Transfer by Radiation Through
Porous Insulations. AIChE Journal, Vol. 5, No. k, December 1959-
35. Cess, R. D.: The Interaction of Thermal Radiation with Conduction and
Convection Heat Transfer. Advances in Heat Transfer, Vol. 1, Academic
Press, New York, New York,
36. Hottec, H. C. and A. F. Sarofim: Radiative Transfer. McGraw-Hill Book
Company, New York, New York, 1967.
37. McAdams, W. H. : Heat Transmission. McGraw-Hill Book Company,
New York, New York, 195^.
38. Hadvig, Sven: Gas Emissivity and Absorptivity: A Thermodynamic Study.
Journal of the Institute of Fuel, April 1970, p. 129.
39. Allen, D. N. and R. V. Southwell: Relaxation Methods Applied to
Determine the Motion in Two Dispersions of a Viscous Fluid Past a Fixed
Cylinder. Quarterly Journal of Mechanics and Applied Mathematics, Vol. 8,
1955 •
UO. Thorn, A. and C. J. Apect : Field Computation in Engineering and Physics.
D. van Nostrand Co., London, England, 1961.
Ul. Friedmann, M., J. Gillis, and N. Liron: Laminar Flow in a Pipe at Low
and Moderate Reynolds Numbers. Applied Scientific Research, Vol. 19,
1968.
U2. Runchal, A. K. and M. Wolfshtein: Numerical Integration Procedure for
the Steady State Navier -Stokes Equations. Journal of Mechanical
Engineering Science, Vol. 11, No. 5, 1969-
U3. McDonald, J. W., V. E. Denny, and A. F. Mills: Numerical Solution of the
Navier-Stokes Equations in Inlet Regions. Journal of Applied Mechanics,
ASME, Vol. 39, No. U, December 1972.
83
-------
kk. Peaceman, D. W. and H. H. Rachford, Jr.: The Numerical Solution of
Parabolic and Elliptic Differential Equations. Journal of the Society
of Industrial and Applied Mathematics, Vol. 3» 1965.
H5. Ames, W. F.: Numerical Methods for Partial Differential Equations.
Barnes & Nobels, Inc., New York, New York, 1969.
U6. Chigier, N. and A. Chervinsky: Experimental and Theoretical Study of
Turbulent Swirling Flow Jets Issuing From a Bound Orifice. Technion,
Israel Institute of Technology, Department of Aeronautical Engineering
FAE Report No. U6, November 1965.
^7. Larson, D. H. and D. Shoffstall: Aerodynamic Control Over Emissions of
Nitrogen Oxides and Other Pollutants From Fossil Fuel Combustion.
Institute of Technology, .Project No. 8933, EPA Contract No.
68-02-0216. Final Report Vol. I & II, 1973-
U8. Heap, M. P. and T. M. Lowes: Development of Combustion System Design
Criteria for the Control of Nitrogen Oxide Emission from Heavy Oil and
Coal Furnaces. International Flame Research Foundation Ijmuiden,
December 1972, EPA Contract No. 68-02-0202.
-------
SECTION VII
LIST OF SYMBOLS
A Frequency factor (see Eq. (3?))
a,b Minimum and maximum eigenvalues
bO»bl>b2 Radiation boundary condition coefficients (see Eq. (6l))
a<£jttyj <=$}
-------
ka Absorption coefficient
ks Scattering coefficient
k Thermal conductivity of mixture; kinetic energy of turbulent
motion of mixture; rate constant
JL Mixing length
L Heat of vaporization
L0 Characteristic gas thickness (see Eq. (50))
m Mass fraction; mass per unit droplet radius class:
iteration parameter
M Molecular weight
n Normal direction
Pr Effective mixture Prandtl number
P Pressure
r Radius: production term in species equation; droplet radius
Q, Heat of combustion
QR Radiation flux
R Residual: gas constant
S Tangential direction
Sc Effective mixture Schmidt number
Effective gaseous Schmidt number associated with
turbulent kinetic energy ( U
Sx Swirl number (see Eq. (20)
BUJ Group of terms in vorticity equation (see Eq. (13))
86
-------
t Time
T Absolute temperature
u,v Velocities in x and r coordinate directions, respectively
v Time-averaged velocity
w Swirl velocity
x,r Coordinates axially and radially, respectively
y Distance in radial direction
Yog Mass fraction of oxidant
e Convergence criterion; emissivity
6.. Khonecker delta (equals unity if indices are same,
or zero if different) (see also, Eq. (77))
IY Effective exchange coefficient of species i ( pDj_)
Fh Effective exchange coefficient of heat ( k/Cp)
Ij Effective exchange coefficient of particle class ( pDj)
8 Circumferential direction
X Mixing length parameter (see Eq. (l8))
H Absolute viscosity of mixture
v Kinematic viscosity
§ Normalized vorticity ( cu/r)
p Density of mixture
CT Viscosity weighting parameters (see Eq. (2l));
Stephan-Boltzmann constant
T Shear stress
-------
(U
Subscripts
eff
g
h
i
3
k
1
m
n,s
P
t
w
Superscripts
Representative dependent variable:
overall coupling and linearization term
Stream function
Vorticity
Imparts turbulent characteristic
Gas or gaseous phase
Upper
Gas phase component in mixture; grid node
Liquid phase component in mixture; grid node
Turbulent kinetic energy
Lower
Maximum
Normal and tangential direction
Reference
Particle or liquid phase
Tangential or swirl
Wall
Transpose of tensor
Refer to u+ and y+ identified in Eqs. (23) and (2U)
-------
BIBLIOGRAPHIC DATA
SHEET
1. Report No.
EPA-650/2-73-045
S.^ecipient's Accession No.
4. Title and Subtitle
A Study of Combustor Flow Computations and Comparison
with Experiment
"5. Report Date
December 1973
6.
'. Author(s)
R. F. Anasoulis and H. McDonald
8. Performing Organization Kept.
No.
'. Performing Organization Name and Address
United Aircraft Research Laboratories
400 Main Street
East Hartford, Connecticut 06108
10. Project/Task/Work Unit No.
ROAP 21ADG-10
11. Contract/Grant No.
68-02-0267
12. Sponsoring Organization Name and Address
EPA, Office of Research and Development
NERC-RTP, Control Systems Laboratory
Research Triangle Park, North Carolina 27711
13. Type of Report & Period
Covered
Final
14.
15. Supplementary Notes
16. Abstracts
The report presents a computational procedure for calculating the coupled flow and
chemistry within combustion devices. The procedure solves the time-aver aging
Navier-Stokes equations with coupled chemistry, including the effects of turbulence
and radiative heat transfer, using a novel field relaxation method. Although the
procedure employs a relatively simple turbulence model, the model can be easily
modified within the framework of the computational method. The flow and chemistry
within a representative furnace have been computed, using the procedure; and
the computations are presented and compared with experimental data.
17. Key Words and Document Analysis. 17a. Descriptors
Mathematical Models
Combustion
Combustion Chambers
Aerodynamics
Navier-Stokes Equations
Chemical Reactivity
Turbulence
Nitrogen Oxides
Furnaces
17b. Identifiers/Open-Ended Terms
Coupled Flow
Computational Procedures
Radiative Heat Transfer
17c. COSATI Field/Group 20M, 13B
18. Availability Statement
Unlimited
19. Security Class (This
Report)
UNCLASSIFIED
20. Security Class (This
Page
UNCLASSIFIED
21. No. of Pages
22. Price
FORM NTIS-35 (REV. 3-72)
THIS FORM MAY RE REPRODUCED
89
USCOMM-DC 1 49R2-P72
------- |