!

!

y

NASA/TM--1998-208823

JF-o ?

f

d_W c_p

Improved Modeling of Finite-Rate

Turbulent Combustion Processes

in Research Combustors

Thomas J. VanOverbeke

Lewis Research Center, Cleveland, Ohio

i

November 1998

https://ntrs.nasa.gov/search.jsp?R=19990009344 2018-05-11T16:59:19+00:00Z
The NASA STI Program Office... in Profile

Since its founding, NASA has been dedicated to

the advancement of aeronautics and space

science. The NASA Scientific and Technical

Information (STI) Program Office plays a key part

in helping NASA maintain this important role.

The NASA STI Program Office is operated by

Langley Research Center, the I_ad Center for

NASA's scientific and technical information. The

NASA STI Program Office provides access to the

NASA STI Database, the largest collection of

aeronautical and space science STI in the world.

The Program Office is also NASA's institutional

mechanism for disseminating the results of its

research and development activities. These results

are published by NASA in the NASA STI Report

Series, which includes the following report types:

TECHNICAL PUBLICATION. Reports of

completed research or a major significant

phase of research that present the results of

NASA programs and include extensive data

or theoretical analysis. Includes compilations

of significant scientific and technical data and

information deemed to be of continuing

reference value. NASA's counterpart of peerreviewed formal professional papers but

has less stringent limitations on manuscript

length and extent of graphic presentations.

TECHNICAL MEMORANDUM. Scientific

and technical findings that are preliminary or

of specialized interest, e.g., quick release

reports, working papers, and bibliographies

that contain minimal annotation. Does not

contain extensive analysis.

CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

CONFERENCE PUBLICATION. Collected

papers from scientific and technical

conferences, symposia, seminars, or other

meetings sponsored or cosponsored by

NASA.

SPECIAL PUBLICATION. Scientific,

technical, or historical information from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific

and technical material pertinent to NASA's

mission.

Specialized services that complement the STI

Program Office's diverse offerings include

creating custom thesauri, building customized

data bases, organizing and publishing research

results.., even providing videos.

For more information about the NASA STI

Program Office, see the following:

• Access the NASA STI Program Home Page

at http:l/www.sti.nasa.gov

• E-mail your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA Access

Help Desk at (301) 621-0134

• Telephone the NASA Access Help Desk at

(301) 621-0390

Write to:

NASA Access Help Desk

NASA Center for AeroSpace Information

7121 Standard Drive

Hanover, MD 21076

.t

ilJi
NASA/TM--1998-208823

Improved Modeling of Finite-Rate

Turbulent Combustion Processes

in Research Combustors

Thomas J. VanOverbeke

Lewis Research Center, Cleveland, Ohio

National Aeronautics and

Space Administration

Lewis Research Center

November 1998
j

NASA Center for Aerospace Information

7121 Standard Drive

Hanover, MD 21076

Price Code: A13

Available from

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22100

Price Code: A13

!rl 1 i_
ii

Abstract

The objective of this thesis is to further develop and

test a stochastic model of turbulent combustion in

recirculating flows. There is a requirement to increase

the accuracy of multi-dimensional combustion predictions.

As turbulence affects reaction rates, this interaction must

be more accurately evaluated. In this work a more

physically correct way of handling the interaction of

turbulence on combustion is further developed and tested.

As turbulence involves randomness, stochastic modeling is

used. Averaged values such as temperature and species

concentration are found by integrating the probability

density function (pdf) over the range of the scalar. The

model in this work does not assume the pdf type, but solves

for the evolution of the pdf using the Monte Carlo solution

technique. The model is further developed by including a

more robust reaction solver, by using accurate

thermodynamics and by more accurate transport of elements.

The stochastic method is used with Semi-Implicit Method for

Pressure-Linked Equations. The SIMPLE method is used to

solve for velocity, pressure, turbulent kinetic energy and

dissipation. The pdf solver solves for temperature and

species concentration. Thus, the method is partially

familiar to combustor engineers. The method is compared to

benchmark experimental data and baseline calculations. The

baseline method was tested on isothermal flows, evaporating

iii

sprays and combusting sprays. Pdf and baseline predictions

were performed for three diffusion flames and one premixed

flame. The pdf method predicted lower combustion rates than

the baseline method in agreement with the data, except for

the premixed flame. The baseline and stochastic predictions

bounded the experimental data for the premixed flame. The

use of a continuous mixing model or relax to mean mixing

model had little effect on the prediction of average

temperature. Two grids were used in a hydrogen diffusion

flame simulation. Grid density didn't affect the

predictions except for peak temperature and tangential

velocity. The hybrid pdf method did take longer and

required more memory, but has a theoretical basis to extend

to many reaction steps which cannot be said of current

turbulent combustion models.
ii! I !i
TABLE OF CONTENTS

Chapter 1 - INTRODUCTION

i.i - Motivation for the Thesis

1.2 - Progress in Simulation of Combustion

1.3 - Turbulent Combustion

1.4 - Stochastic Turbulent Combustion Models

1.5 - Contribution of the Present Work

Chapter 2 - BASELINE SOLUTION METHOD FOR CONVECTIONDIFFUSION PROBLEMS

2.1 - Governing Equations

2.2 - Discretization of the General Equation

2.3 - Solution of Pressure

2.4 - Grid System and Axisymmetric Equations

2.5 - Boundary Conditions

2.5.1 - Inlet Conditions

2.5.2 - Axis of Symmetry

2.5.3 - Outlet

2.5.4 - Solid Walls

2.6 - Spray Combustion Modeling

2.7 - Spray Solution

Chapter 3 - NUMERICAL SIMULATION OF COMBUSTION

3.1 - Infinite rate Combustion Models

3.2 - Finite-rate Combustion Models

3.3 - Turbulent Combustion Models

3.3.1 Eddy-breakup Combustion Model

3.3.2 Probability-Density-Function Models

iv

1

1

3

7

9

i0

13

13

17

21

24

24

25

25

26

26

27

31

36

39

43

48

49

51

V

Chapter 4 - HYBRID COMPOSITION PDF SOLVER 64

4.1 - Monte-Carlo Solution Method 68

4.1.1 - Convection/Diffusion Phase 72

4.1.2 - Molecular Mixing Phase 75

4.1.3 - Combustion Phase 81

4.2 - Further Development of Hybrid Pdf Method 83

Chapter 5 - COMBUSTOR MODEL VALIDATION 87

5.1 - Combustor Model Validation Data 88

5.2 - Spray Model Validation 94

5.2.1 - Isothermal Open Swirling Flow 94

5.2.2 - Evaporating Spray Model Validation 107

5.2.3 -Spray Combustion Calculations 136

5.3 - Monte-Carlo Pdf Model Calculations 178

5.3.1 - Non-Swirling Diffusion Combustor 179

5.3.2 - Swirling Diffusion Combustor Case 200

5.3.3 - Premixed Combustor Case 226

5.3.4 - Swirling hydrogen diffusion flame 238

"Chapter 6 - CONCLUDING REMARKS 258

References 266

Appendix 1 Gas Phase Source Terms 276

Appendix 2 Spray Source terms 277
iEl ! _,
1

Chapter 1

INTRODUCTION

i.i Motivation for the thesis

Improving the numerical simulation of combustor flows

is a never-ending endeavor. There is a continuing demand to

increase gas turbine performance and decrease emissions.

These demands do not complement each other. The demand for

increased performance is typically met by increasing cycle

pressure and final temperature. This is allowed by

continual advances in metallurgy. However, increasing

pressure and temperature usually increases emissions. Thus,

combustor designs are being evaluated to test for increased

performance with lower emissions. One method to help in

designing gas turbine combustors is Computational Fluid

Dynamics and combustor models. These models are benchmarked

or validated using existing engine or test data. CFD

combustor model predictions are accurate or calibrated

within similar engine designs. Thus, new concepts cannot

be totally evaluated by current CFD models. The number of

prototypes in an engine development program may be reduced

by using CFD combustor model predictions. Initial

computational designs must be tested on engine stands. CFD

can be also be used to correct design deficiencies found in

prototypes. Commercial gas turbines are very difficult and

costly to instrument and test. The use of CFD models

provides data at all locations in a combustor field. Thus,
the nature of CFD complements engine test programs.

Newer engine designs are going to higher peak

temperature and pressure to improve both fuel economy and

thrust to weight ratio. Severe reductions in emissions are

being required by governments. Some European countries are

considering instituting fees based upon the total amount of

pollutants emitted by aircraft flying through their

airspace. Older aircraft producing more pollution would be

more heavily taxed. The evaluation of proposed engine

designs is causing a demand for increased accuracy in CFD

combustor models as well as a wider range of validity in

those models. Engine designers are tolerating more complex

combustor models in a desire for more accurate CFD

predictions. Computer resources to allow more accurate

modeling is becoming more affordable with time. Combustion

simulations were once entirely done on mainframe computers

or supercomputers. Running many combustion simulations on

supercomputers can cost thousands of dollars. Workstations

have been developed to the point that taking into account

the amount of time jobs are in queues, the total turn-around

time is comparable for high-end workstations and

supercomputers. The cost of doing calculations on

workstation is significantly less than running on

supercomputers. To further reduce turn-around time of

workstations, computers may be clustered together to give

performance similar to dedicated super-computers. Many
workstations see relatively little use at off hours. With

suitable networking, these under-utilized workstations may

be used to perform complex 3-D predictions at small

additional cost.

1.2 - Progress in the Simulation of Combustion

Much of the model improvement in the last decade has

been in the area of gridding. Curved boundaries and round

dilution jets are modeled without stair-stepping the

geometry. Unfortunately, generating grids for complex

combustors can take up most of the time to perform the

calculation. One way around this is to overlap simpler

grids. This technique essentially allows imbedding grid for

a flow feature into another simpler grid which poorly

defines the flow feature. A good example is a cylindrical

grid for a rod imbedded into a cartesian grid for a

rectangular duct. Another technique being developed is the

use of unstructured grid solvers. In unstructured grids,

the grid is generated by breaking the flow geometry into

many smaller general elements. If more refinement is

needed, the elements are further divided. The data storage

for fluid elements is quite different than for structured

grid solvers. The unstructured grid method incurs

additional work to keep track of neighboring grid cells.

This technique is only in its infancy. It is unknown which

solution method will prove best for combustor design.

Over the past twenty years, most elliptic flow
calculations have been done using variants of the SemiImplicit-Method for Pressure-Linked-Equations developed by

Patankar I. The usual implementation is to assume low Mach

number flow. In this flow regime, pressure is largely a

relative variable, that is, pressure is fairly constant.

Changes in density are due to changes in molecular weight

and temperature. This solution technique is referred to as

pressure based method. Industry regards the solution method

as robust, but, as rather dissipative. If flow is at an

angle to the grid, calculated flow quantities may be

dissipated. A great deal of effort was expended in the

eighties to develop improved differencing methods to improve

accuracy of pressure based schemes. Skew differencing was

developed to take into account the angle of the flow to the

grid. Quadratic upwind interpolation and second-order

upwind differencing were also developed. Theoretically

these schemes are more accurate, but, the schemes are not as

robust as first-order upwind differencing schemes.

Convergence is much more difficult for improved schemes.

Some of these schemes blend in varying amounts of upwind

differencing to improve robustness.

Density based solvers are thought by some individuals

to be less dissipative than the pressure based solvers.

However, these solvers must include artificial dissipation

terms for robustness. Density based solvers are efficient

in the solution of higher Mach number or compressible flows.

!Elli
5

The convergence of the solver greatly deteriorates with

decreasing Mach number as the flow becomes incompressible.

Solution techniques using density based solvers for low

speed flows are being developed, such as NASA Lewis Research

Center's ALLSPD2 code. The ALLSPD code uses

preconditioning, gauge pressure and pseudo-time stepping to

maintain good convergence rates at low Mach number flows.

This solution technique is in the development phase.

Another technique to improve convergence rates is to

use multi-grid solvers. Joshi and Vanka 3 performed multigrid calculations for gas-turbine type geometries. The

method was also demonstrated for 3-D geometries for studying

hot gas ingestion for short take off and vertical landing

aircraft 4. Multi-grid was also used in a numerical ramjet

study by Vanka and Krazinski 5. A four-step combustion

mechanism with combination eddy-breakup and arrhenius

reaction rate submodel was used in the study. A single grid

system was used in the solution of species concentration due

to the strong influence of the source terms in the chemical

kinetic equations.

The overall chemistry or reaction kinetics of fuels can

be very complicated. Also, combustion in practical gas

turbine combustors is highly turbulent. Improved turbulence

models for variable density flows are the least developed of

all turbulence models, and typically have been calibrated

for isothermal incompressible flows. Turbulence is usually
6

ignored by researchers who study fundamental combustion

processes the most, the experts in chemical reaction

kinetics. These experts predict reaction rates accurately

over the range of combustor temperature, pressure and

species concentration. Unfortunately, this is accomplished

by employing dozens of chemical reactions with widely

varying time characteristics. These kinetic reaction

systems must be solved by specialized solvers using small

time steps, predictor-corrector steps, etc. Many reactions

involve species the engineer isn't concerned with. Not

enough effort is being spent on producing reduced chemical

reaction mechanisms that will give acceptable engineering

predictions. Laminar combustion experts have yet to develop

simpler kinetic reaction schemes which they feel are

accurate enough to predict exhaust emissions for different

engine concepts. Typical combustor models employ five to

ten species and as few chemical reactions as possible. The

chemical kineticists typically insist on twenty or more

reactions for hydrocarbons and would like to see as many as

fifty reactions. The time steps needed to accurately solve

species concentrations for these kinetic mechanisms can be

extremely small, much shorter than the pseudo time step

taken in typical combustor calculations. A good time step

for the LSENS 6 general kinetics and sensitivity solver is

of the order of ten micro-seconds. The code employs

variable time steps to improve solution speed. The solution

IF!li_
7

efficiency has been substantially improved from the original

solver. This code is one of the most highly developed and

efficient chemical kinetics solver existing today. However,

it is still too computationally expensive to use for solving

complex multi-dimensional combustor flows. Techniques used

to accurately solve laminar one dimensional combustion are

too complex and costly at the present time to use for the

solution of multi-dimensional turbulent recirculating flows.

Typical design models use a very small number of reaction

steps to predict performance. Engine emissions are usually

predicted by post-processing, adjusting constants or

employing empirical corrections.

1.3 - Turbulent Combustion

The prediction of combustion is greatly affected by

turbulence. The length of current gas turbine combustors is

fairly short. Fuel and oxidant must be thoroughly mixed

within a short distance at fairly high velocities.

Combustor velocities are greater than laminar flame

velocities. Combustor exit velocities in the newer gas

turbine engines are getting into the high subsonic realm.

The overall mixing of fuel and oxidants in a combustor is

greatly augmented by turbulence.

Turbulence causes changes in the instantaneous values

of velocity, species and temperature. Reaction rates are

nonlinear functions of concentration and temperature. Thus,

the averaged reaction rate isn't the same as the reaction
8

rate of averaged species concentration and temperature.

Previous modeling has employed many ad-hoc corrections to

account for the interaction of turbulence on combustion.

One such correction is to include a factor to convert

laminar reaction rates into turbulent reaction rates.

Turbulence produces eddies or modules with varying

concentrations of fuel and oxidant. These eddies go through

a breakup or decay process before the eddies completely mix.

As turbulence continually generates new eddies in the flows,

this is an ongoing process. Thus, turbulence causes

incomplete mixing which causes an overlap of fuel and

oxidant species concentration particularly near

stoichiometric conditions where there is the greatest

chemical reactivity. The overlap in average reactant specie

concentration makes turbulent combustion resemble finite

rate combustion. Current turbulent combustion methods do a

good job of predicting combustor exit temperature patterns

and major species concentrations. Minor species and

intermediates, especially carbon monoxide, are not well

predicted. Unfortunately, minor species can cause

significant pollution levels. The level of engine emissions

is being monitored and future emissions have been regulated

due to serious environmental considerations. The

development of an American supersonic transport was canceled

in the seventies due to concerns of immense nitrous oxide

production at high altitudes. Lakes and forests are being

7!|_ll
adversely affected by current fossil fuel combustion. The

combustion process must be better predicted to aid in the

design of lower emissions combustors.

1.4 - Stochastic Turbulent Combustion Models

A more physically correct model of the turbulencechemistry interaction is to treat combustion as a stochastic

process. This is the idea behind current assumed shape

probability density function models. These models assume

certain shapes for the probability of a conserved scalar.

The probability density function may be inferred by solving

modeled equations for the means and variances of variables.

The major assumption in presumed shape pdf modeling is that

reaction rates are fast in these models. This can be

partially relaxed by including another independent variable,

typically for reactedness. The additional variable is

usually assumed to be uncorrelated to the conserved scalar.

Unfortunately, the additional variable is usually

correlated. Including more species and chemical reactions

adds complexity and increases the number of assumptions.

Even allowing for incomplete or partial combustion in

assumed shape pdf models causes predictions of high reaction

completeness as the simplified combustion mechanisms

employed typically used are only valid at high temperatures.

Another stochastic method involves solving for the

evolution of probability density function. This modeling

eliminates the chemical reaction rate closure problem
I0

inherent in traditional combustion modeling. Stochastic

treatments of species transport are said to exactly include

the chemical reaction rate term. This has been known for

quite some time, but, the capability to perform stochastic

calculations of practical combustor flowfields is being

developed at the experimental or research level. Actual

engine combustor design predictions have to be done in a

short period of time to be effective in the combustor design

process. Also, industrial predictions have to be cost

effective. Stochastic modeling demands large expensive

computer resources. As computer power continues to increase

and the cost of doing combustion simulations becomes more

acceptable, the use of stochastic combustion models will

increase.

1.5 - Contribution of the Present Work

Assumed shape pdf models have been used for over a

decade and development of these models has largely

• stagnated. The next step in pdf modeling is the composition

pdf, which requires roughly an order of magnitude increase

in memory and computer time. Pdf modeling allows a

substantially more robust and accurate method of calculating

minor species for turbulent combustion, albeit at

substantially increased cost.

Newly developed hybrid pdf models have predicted

excellent agreement with experimental data for parabolic or

one-way combustor flow. The method is called hybrid as

![Ii i
ii

stochastic methods are used to solve for average species

composition, density and temperature, while traditional CFD

techniques are used to solve for other flow quantities.

Full stochastic simulations of all flow quantities is beyond

current capabilities. The hybrid method retains the

traditional CFD techniques which engineers are used to,

while it treats turbulent chemistry interactions in a much

more physically correct manner.

The purpose of my work is to incorporate and further

develop a composition and enthalpy pdf model into a

computational fluid dynamics model. The hybrid pdf model

uses individual species particles to solve for the evolution

of composition and temperature pdfs. The hybrid pdf model

consists of consecutively solving for velocity, pressure,

turbulence length and time scales, species transport and

chemical reaction. A Monte Carlo solution technique is used

to predict species transport, mixing and species reactions.

Many improvements to the original pdf model were made in the

work. The pdf model was improved to more properly calculate

species transport. The original model used an incorrectly

coded upwind differencing scheme to calculate species

transport. A more accurate scheme was implemented. The

original model had major thermodynamic simplifications.

This was corrected to better predict actual combustor data.

The original model had an extremely simple combustion

solver. This solver worked with the original reaction rate
12

constants, but, didn't for the reaction models in this work.

An improved species solver was implemented. Pdf results are

properly compared against good CFD results. Some pdf

predictions are performed against CFD predictions done years

ago, and claiming the improved results are due solely to pdf

modeling. This ignores many improvements which have been

made to CFD models over the years. Pdf predictions have

also been compared to initial or other poorly refined

predictions. The hybrid pdf model will be tested against a

proven baseline turbulent combustion model and experimental

data for various research combustor flowfields. The

baseline model uses the pressure based method to solve for

velocity, pressure, mean species, species fluctuation and

turbulence scales in elliptic flow fields. The k-E

turbulence model is used. A fuller description of the model

development in this work is covered in chapter 4.

UI -Ii-
13

Chapter 2

BASELINE SOLUTION METHOD FOR CONVECTION-DIFFUSION PROBLEMS

The Semi-Implicit Method for Pressure Linked Equations

or SIMPLE method of Patankar is used in conjunction with the

hybrid composition pdf method in this thesis. The SIMPLE

algorithm method is used for solving momentum, overall mass

conservation, turbulent kinetic energy and dissipation rate

of turbulent kinetic energy. The pdf method used here

solves turbulent combustion and species transport. This

chapter covers a rather short discussion on the SIMPLE

method, boundary conditions used and source terms for

combustion species.

2.1 - Governing Equations

The gas phase equations can be put in the following

general form for cartesian geometries:

0t

(2.1)

where p is the density, U is the axial velocity, V is the

radial velocity, F_ is the appropriate diffusion coefficient

and S_ is the source term for the variable 9. For example,

the diffusion coefficient for the laminar momentum equations

would be _i, the laminar viscosity. Turbulence introduces

additional complexity into the transport equation. The way

this is traditionally done for incompressible flow is to

perform a Reynolds decomposition of velocity into a mean and

fluctuating components:
where,

U = "_' + U I

Co+de

D = l'm I [ Ud_

:]_0"-*oot o

14

(2.2)

(2.3)

Substituting in the Reynolds decomposition of axial and

radial velocity into the momentum equation and time

averaging results in terms such as:

pu---I_, pv---i_,and pu--I_7, e tc

which are known as Reynolds stresses. The first two terms

are normal stresses. These can be i_corporated into the

pressure term. If density varies, the dilatation or

divergence of velocity isn't zero. The dilatation can also

be included in the pressure term:

2 (aUk_ + 2pk (2.4)

where _c is the turbulent coefficient of viscosity. The

Reynolds stresses are unknown. The Reynolds stresses must

involve modeling at some point due to problems with closure.

Lower level modeling involves simple correlations for the

stresses. Higher level or second moment modeling involves

solving modeled equations for the transport of Reynolds

stresses. These equations are formed by a sequence of

multiplication of the momentum equations, Reynolds

decomposition and time averaging. These equations have
15

unknown triple correlation terms, which must then be

modeled. Lower level turbulence modeling has proved

suitable for prediction of turbulent combustion. The nonnormal Reynolds stresses are modeled using a generalization

of Boussinesq's eddy viscosity concept or the gradient model

assumption:

dx_ )

Following the Prandtl-Kolmogorov definition of eddy

viscosity, the viscosity coefficient is composed of length

and velocity scales. The preferred method of supplying

these scales is the two equation k-e turbulence model. The

k-_ model was developed by Jones and Launder v. The term k

is known as the turbulent kinetic energy:

(2.6)

The term £ is known as the turbulent kinetic dissipation

rate. If the characteristic length is defined as k312/e and

velocity as k In, the turbulent diffusion coefficient can be

defined as:

k2 (2.7)

Equations can be formed for k and £. However, these

equations have been modified to better predict turbulence.

The modeled equations are:
and

a--{ 0 (pk) + a[p Uk-Fk-- _ -_

+ pVk-l" k : p,=G - p8

a-_ a (pc)+ a_[pu=-r"_]a=l+

--_y[pVe-l".-_] = C_C.Ge/k-C2pe2/k

where G is the turbulence generation term:

16

(2.8)

(2.9)

(2.10)

The k-_ model allows for the variability of length and

velocity scales with history effects. That is, the

viscosity isn't just a function of local gradients or

scales.

Some of the non-normal turbulent stresses can be

incorporated into diffusion terms in the general equation by

defining a effective diffusion coefficient which is the sum

of the laminar and turbulent diffusion coefficients.

(2 .ii)

The left over turbulent terms are put into source terms for

the momentum equations. For example, the axial momentum

source term for Cartesian flows is:

a_ au_ alp av/ (2.12)

![! ! _:
17

2.2 - Discretization of the general equation

Except for simplified problems, the general equations

do not have an analytic solutions. The equations must be

solved numerically. The

equation is discretized by

using the control volume

approach. The general

transport equation is

integrated over the dashed

control volume for the node

P shown in figure 2.1. The

values at the faces of the

control volume are presumed

to preside over the whole

face. This gives:

W

N

W

I

n j

I

P e'1

1

1

1

1

S '

!

S

E

figure 2.1 Control Volume

(pM, p-p,5) _xA.v

+F._.-D m (_z-_p) -

at (2-13)

where the source term has been linearized and the convection

coefficients at the control volume faces are given by:

Fe = (pU) e_y, (2.14a)

Fo = (p_ ,,_y, (2.14b)

F n = (pV) aaX, (2.14C)
18

F, = (p$0,x_x, (2.14d)

where _x and Ay are the width and height of the control

volume. The diffusion coefficients are given by:

(2.15a)

(2.15b)

(2.15c)

(2.15d)

where the term (_x) e is the distance between the points P

and E, etc.

The equation is further developed by manipulation of

the continuity equation. Integrating the continuity

equation over the same volume gives:

(pp_pO) ax&y + Fe _ Fw + F= - F s = 0 (2.16)

At

If the integrated continuity equation is multiplied by #p

and subtracted from the integrated general equation, the

following equation is obtained:

If!_li
_t

+Fo (_.-__)-DB(__-__)-Fw (_w-_p)+D_ (__-_w)

+F= (_=-_p) -D_ (_-_p) -F n (_8-_P) +Ds (_p-_s)

= (Sc+S_p) _x_y

19

(2.17)

This equation involves intermediate values at the cell

faces. The goal is to discretize the above equation into an

equation involving _p and its neighboring nodes:

a _ = a_ + a_w + a_/_H + a_s ÷ b (2.18)

This may be done by forming interpolations involving the

adjoining nodes for the intermediate values at the cell

faces. However, if any of the above coefficients is allowed

to become negative, the answers may become non-physical. To

ensure physically realistic results the hybrid differencing

scheme is used in this work. The hybrid differencing scheme

is a combination of central differencing and upwind

differencing. Hybrid differencing uses more accurate

central differencing when the absolute value of the ratio of

the convective flux to the diffusive flux (or Peclet number)

is less than two. When the Peclet number is greater than

two, the upwind differencing scheme coefficients are used,

which is physically realistic at high convection rates.

The use of hybrid differencing is illustrated with the

following expression for fluxes on the east side of the

control volume:
2O

Fe(¢e-_p)-DE(_E-¢p) = a_(_p-¢E) (2.19)

If the absolute Value of the Peclet number is less than two,

central differencing is used for the term #., thus:

which gives:

(2.20)

Fe for IPel<2 (2.21)

a E = D e - -_-,

If the Peclet number is greater than two, the convective

flux from node P overwhelms the diffusive flux. The

diffusive flux is ignored and _. is taken as #p, and:

F e(_P-_P) = aE(¢p-¢E) (2.22)

which gives:

(2.23) a_ = 0, for Pe > 2

If Pe is less than -2, then the convective flux from node E

is much larger than the diffusive flux, which is ignored.

Thus:

Fe (_E-_P) = a_ (CP-¢E) (2.24 )

which gives:

a E = -F o, for Pe < -2.

The general value for a_ has been worked out to be:

(2.25)

(2.26)

a_ = [-Fe,D e - -_-, Fo01

l!! i
where the symbols I _ stands for the largest of the

quantities within it. The formula for aw is:

A w = IFw, Dw+- Fwl _, 0

(2.27)

The formulas for a_ and a s are similar.

rest of the coefficients are:

The formula for the

ap = a E + a w + a_ + a s + a_ - SpaXny,

a_- p_AxAy

at "

(2.28a)

(2.28b)

21

00

b = ScAXAy + ap_p (2.28c)

2.3 - Solution for Pressure

An offset or staggered grid is used to solve the

momentum equations. The nodes for axial velocity occur on

the east and west faces of the scalar control volume. Thus,

the control volume for axial momentum is offset from the

•scalar control volume. Thus the discretization equation for

axial velocity is:

aeU e = _a_U_ + b + (Pp-PE)Ae (2.29)

where the change in pressure with axial distance has been

pulled out of the source terms. A, is the area term on

which the pressure is acting. A shorthand expression is

used for neighbor coefficients and velocities. A very

similar discretization equation for radial velocity can be

formed.
22

In the pressure solution process, the solutions of

velocity and pressure is divided into two components. The

current velocity and pressure solutions are denoted _,

and P'. Due to an inexact solution there is a pressure

correction P' which causes corresponding velocity

corrections U', V'.

U = U* + U I, (2.30a)

V = V* + V/, (2.30b)

p = p- + p/

A velocity correction equation may be formed by

plugging in _ and _ into the discretized momentum

equations. Subtracting the equation for _ from the

equation for U gives a velocity correction equation:

The first term is unknown, but, disappears at convergence,

so is ignored. This gives:

(2.30c)

UIe - Ae _IP/-P/_P E, (2 .32)

a o

Plugging this in the equation for U gives the velocity

correction equation:

Ae I I

Go = U: +

as

(2.33)

The continuity equation is integrated over the control

IF1 ! i-
23

volume for node P. The resulting equation involves

velocities at the faces of the control volume. The

equations for velocity correction are then substituted into

this equation. A discretized for P' equation is formed by

rearranging the terms:

a_D_ = a_D E'+ a_D_ + A_D_ + asp _ + b, (2.34)

where b includes the mass imbalance due to _ and _ and the

change in density with time.

This procedure for the solution of convection-diffusion

problems is a very brief description of the Semi-Implicit

Method for Pressure Linked Equations. The SIMPLE algorithm

and its variants are the most used solution methods for low

Mach number combustor flows. A much more general and

thorough development of the SIMPLE method is given in

Patankar [i].

Computational models using the SIMPLE algorithm and

hybrid differencing have been widely used and have proven

robust for combustor flows. The accuracy of the hybrid

differencing scheme is believed to deteriorate when the flow

does not align with the grid. Correa and Shyy et al. 8 have

studied second order upwind differencing and the QUICK

differencing scheme of Leonard 9. Using simple simulations,

it was demonstrated that one scheme was not generally

superior for all cases. Higher order schemes are generally

thought to have problems with robustness. However, Correa

and Shyy stated that convergence rates of the more accurate
24

numerical schemes were comparable to the hybrid scheme in

the context of curvilinear coordinates. It was also found

that with an appropriate grid distribution, the difference

between finite difference operators could be quite small.

2.4 - Grid System and Axisymmetric Equations

The grid system used in this work is orthogonal. The

grid setup is Practice B cited by Patankar [i]. Practice B

first sets up the control volumes and then places the scalar

gridpoints in the center of the control volume. The method

used in this work uses plain averages rather than the more

exact geometric interpolations or harmonic averages for

physical quantities such as viscosity, etc. at cell faces.

The boundary control volumes are of infinitesimal thickness,

which eliminates half control volumes for scalar values.

Axisymmetric geometries are solved in this thesis. The

equations are slightly different than for Cartesian

geometries. The general equation for orthogonal

axisymmetric geometries is:

ot ÷ pv - =s (2.35)

The source terms and diffusion coefficients for the various

momentum and scalar equations are given in Appendix i.

2.5 - Boundary conditions

Inlet conditions are read in and set the main program.

Initial values are set in subroutine INIT. A restart

lilI i
25

capability is also in the main program. Each variable has

its own calculation subroutine and entry point in the

boundary condition subroutine. The Boundary conditions are

implemented in subroutine MODPRO. The program must be

modified for different geometries.

2.5.1 - Inlet conditions

All inlet conditions are specified. A large portion of

this work is to validate computational algorithms. This

requires that inlet conditions be as well specified as

possible, so that quality of numerical predictions can be

judged. Inlet velocities, temperatures, species mass

fractions, and turbulence kinetic energies are specified

based on actual measurements, if possible. If there are

multiple inlets, the inlet pressure is set to be the same

for both.

2.5.2 - Axis of Symmetry

Radial velocity is zero at the axis of symmetry. The

radial gradient of all variables, except radial velocity is

zero at the axis of symmetry. Mathematically:

(4,) = 0 (2.36) r= 0

The axis of symmetry is located on the j=l gridline. The

axisymmetric boundary condition produces the result:

_i,j=1 = _i,j=2, except for velocity V (2.37)

This boundary contribution is implemented each iteration.

During the iteration process the neighbor coefficient for
26

the axis of symmetry is zero as the area term is zero.

2.5.3 - Outlet

At the exit of the combustor the axial gradient of all

variables except for axial velocity is assumed to vanish.

Thus, the neighbor coefficient a_ is set to zero at the

exit. Thus, exit variables, except axial velocity, need not

be specified.

For axial velocity, an exit velocity correction due to

mass conservation is applied every iteration. Generally,

this greatly helps the rate of convergence.

2.5.4 - Solid walls

The law of the wall is used for velocities parallel to

solid walls. The laminar sublayer extends to y+ of 11.63.

The value of y+ is calculated by using the density, laminar

viscosity, turbulent kinetic energy and appropriate length

for the grid cell adjacent to the wall. Thus for axial

velocity:

pc_ v_ _y/2 (2.38)

Y_ = _i

If y+ is less than 11.63, the shear stress is laminar and is

given by:

zw = (2_i_y) (2.39)

Otherwise the cell is in the turbulent flow regime and shear

stress is calculated using the law of the wall:

where K is taken as 0.4187. Shear stresses are included in

!_! 4
27

contributions:

Sp = -i.0 x103o ,

S c = 1.0 x l0H°g

Walls are treated as adiabatic.

(2.43a)

(2.43b)

Thus, the neighbor

coefficients for wall influences for heat transfer and

diffusion of species is set to zero.

2.6 - Spray Combustion Modeling

Kp i l (2.4 0)

rw = In (9y +) J

the momentum equations as source terms.

For turbulent kinetic energy, the neighbor coefficient

corresponding to the wall is set to zero and the source term

is:

Sk = _cG - [2C_/¢pkH/2y+IA_, y+< 11.63 (2.41a)

S k = _G- [2C_/4pkJ/41n(9y *)(KAy)I, y+zll.63 (2.41b)

The turbulent dissipation rate next to the wall is

calculated by assuming that the rates of generation of

turbulent kinetic energy and dissipation are equal so that

the value of £ is:

e = Kay , next to the wall (2.42)

The value of 8 next to walls is fixed by the forcing the

source term coefficients to overwhelm all other neighbor
28

Spray modeling has been under development for about

three decades. Many review type papers on spray combustion

modeling have been written. These include, Faeth's paper in

1977 I°, and Law's paper in 198211, which cites improvements

to original droplet evaporation modeling. The

sophistication of spray models being used for combustor

calculations is increasing. The earliest spray models

treated spray vaporization as a single droplet in quiescent

flow. The difference in velocity between the droplet or

fluid phase and the gas phase velocity was ignored. In a

198312 paper, Faeth cites models which include

consideration of the difference between the droplet velocity

and gas phase velocity. In Sirignano's paper in 198313 ,

fluid flow within the droplet itself was incorporated in a

spray model.

One of the earliest studies of spray combustion was

performed by Onuma and Ogasawara 14. Some conclusions were

that for the spray burner used, the region where the

droplets exist is limited to a small area above the burner

nozzle, and that it was reasonable that droplets do not burn

individually, but, that the vapor cloud formed by the

vaporization of the droplets burns like a diffusion flame.

Many past calculations have successfully ignored spray

vaporization or the fluid phase. Fuel was treated as being

completely in the gaseous phase. However, as the demand for

more efficient and lower emissions combustors has increased

I[I1i
29

along with the needed computer resources, spray modeling is

being developed and implemented in the prediction of gas

turbine combustors.

According to Law, the basic droplet combustion model

was formulated by various people, including Spaulding 15.

This model has since been termed the d_-Law because the

model predicts a linear relationship between the square of

droplet diameter and time. Many improvements have been made

to the d;-Law. These include modeling of droplet heating,

variable specific heat in the gas phase, and the multicomponent fuels.

Actual sprays consists of a fairly continuous

distribution of droplet size and droplet number density.

Spray modeling employs discrete droplets in the prediction

of spray combustion. Onuma and Ogasawara [14] thought that

modeling a single characteristic droplet would be sufficient

to accurately calculate spray combustion. The

characteristic droplet is determined by some sort of

averaging process regarding initial droplet size, position

and velocity. The effect of multiple droplets or groups is

found by taking multiples of the characteristic droplet

effects. Most gas turbine spray modeling assumes that

droplets do not interact, that the spray is dilute. Thus,

the calculation of spray vaporization is treated as a

function of individual droplet characteristics and the

surrounding gas phase. Current models track numerous
3O

droplet groups. A range of different size, position and

velocities may be used to approximate a droplet

distribution. The most computationally expensive droplet

models even model turbulent droplet dispersion. This is

done by incorporating the effect that a random turbulent

eddy would have on the droplet. Enough droplets must be

computed to accurately predict overall spray.

Consideration is being given to modeling non-dilute or

dense sprays, which exist in fuel atomizers. Faeth's paper

in 198716 includes a section on the modeling of droplet to

droplet interactions or dense sprays. Tsai and Sterling Iv,

show that the interaction between droplets increases as the

distance between droplets decreases. There is difficulty in

measuring dense spray properties. Non-intrusive optical

methods are usually used. When multiple droplets are

detected in a measurement volume, the measurement is

typically invalidated.

The fuel break up process is non-uniform. In a section

on the design and development of gas turbine combustors,

A.H. Lefebre _8 writes that the fuel atomization is normally

accomplished by spreading the fuel into a thin sheet and

then increasing its surface area until it becomes unstable

and disintegrates into threads or ligaments. These threads

then break up into rows of droplets. The ligaments and

droplets are highly nonuniform. The droplets are of various

sizes and shapes and are not generally spherical. Droplet

it|I i_
31

vaporization modeling assumes perfectly spherical droplets.

It simply is too difficult to calculate non-spherical

droplets in combustors or measure them for that matter. The

break up of fuel stream into droplets is generally handled

by empirical functions, rather than any physical modeling.

2.7 - Spray Solution Method

A spray combustion model including most of the above

mentioned characteristics was modified. Simplified internal

droplet fluid flow and temperature distribution is handled

by the simplified vortex model of Tong and Sirignano 19.

The spray model is a simplification of the multicomponent

spray modeling done by Raju and Sirignano 2°. The model

assumes a dilute spray. Constant properties are assumed

within the droplets. The droplet trajectories are solved

using discretized Lagrangian equations. The turbulent

dispersion of droplets isn't considered, as it would add at

least an order of magnitude time to spray computations.

Sprays consist of many thousands of droplet particles

each having their own velocity, drag forces, position,

droplet diameter, and distributions of temperature and

vortex strength. Due to limited numerical capabilities, a

limited number of droplets are calculated. Each

characteristic droplet has its own axial position xk, and

axial velocity 9, which is a function of time. The

subscript k refers to the kth droplet group. The change in

droplet position with time is given by:
32

@xk - Uk , (2.44)

at

@Yk _

at v,, (2.45)

az k

___ = wk (2.46)

The droplet velocities are affected by drag effects caused

by differences in velocity with the surrounding gas phase:

.@Uk _ 3 C_Re. (__Uk) " (2.47)

at Z6 p_

at 16 p_

(2.48)

at _6 pF_

(2.49)

Where the subscript g denotes the gas phase, r_ denotes the

radius of the _h droplet and Rex is the droplet Reynolds

number taken with respect to the difference in velocity with

the surrounding gas phase:

ae_=2_-_P[c_-._) ' + cT_-v_)=+c_-_)=] _/= (2.50)

The coefficient of drag is a function of droplet Reynolds

number :

IF! I _-
33

ae_/31 24 I+

CD- Re k 6

(2.51)

Molecular viscosity is evaluated by Sutherlands's equation:

T3/2 (2.52)

_g(T) = 1.4637x10 -6 T+I20

A reference temperature is used for calculation of viscosity

in the vicinity of the droplet:

1 2 T

rzef - 3 Tg + -_ ks (2.53)

where the subscripts refers to the droplet surface.

The droplet regression rate is a composite function of

droplet Reynolds number, Re k. Three ranges are used to

approximate the function:

d(:2)_ 2 1[!Rek]1 2f(B ),

dt Pk[ _ ]

Rek> 20 (2.54a)

d(r 2) _ _l I/3]Re_°771n(l+B k) dt _[i÷ (1÷_ek) l<Rek<l (2.54b)

d(r2) - _--l[l+(l+Re k)I/3]In(l+B_),

dt Pk

Rek<l (2.54c)

where B k is the Spaulding transfer number. The function
f(B_) is obtained from the solution of Emmon's problem.

The internal droplet temperature is calculated from:

a% k_ [ _% a%]

a--F: _-'cp,,,,_._["_j_ +(_+c(t)')-6&j

where k I refers to thermal conductivity in the droplet or

liquid phase, and C(t) is given by:

(2.55)

34

c (t) - _v r ': -aE (2.56)

represents the streamline of a Hill's Vortex in the

droplet. The boundary conditions for the droplet

temperature are:

a_#_ _.r_l _aT#

aa 17[ k, jr.-_

_=0 (2.57a)

where ik is the latent heat of vaporization of the fuel.

The initial droplet temperature is set to the droplet

injection temperature. The Spaulding transfer number is:

B_. C;,(T,-Th,) = (Y,,,-Y_)

Ik,of_ (i-Y_)

(2.58)

Ik,,_ is the effective latent heat of vaporization modified

for heat loss to the droplet interior:

!T! i
l,,.t::

35

(2.59)

where m_ is the droplet vaporization rate.

fraction at the droplet surface is given by:

I Wa -I -I

The fuel mass

(2.60)

where _ is the molecular weight of the gas excluding the

fuel and Z is the mole fraction which is found by assuming

phase equilibrium at the droplet surface, and using the

Clausius-Clapeyron relation:

rill11] (2.61)

where Pn is the normal atmospheric pressure, P is the

prevailing pressure, I a is a constant and _ is the boiling

temperature for the droplet. More accurate expressions than

the Clausius-Clapeyron relation involving additional

parameters can be found in Reid et al. 2_

The spray provides a mass source term in the continuity

equation. The droplets also affect the momentum equations

through drag forces. Thus, appropriate source terms must be

included in the gas-phase equations. The source terms are

listed in Appendix 2.
36

Chapter 3

NUMERICAL SIMULATION OF COMBUSTION

The field of combustion demands knowledge in computational fluid dynamics, turbulence modeling, chemical

kinetics modeling of combustion, real gas effects, combustor

concepts and turbulent-chemistry effects. A variety of

combustion models exist. Combustion models cannot be

implemented as easily as changing damping coefficients for a

turbulence model. Many of these turbulent combustion models

are application specific, that is, the combustion model was

developed for a specific combustion regime. The combustion

regime is decided by the length scale of combustion, which

is dependent on the rates of mixing and the rates of

chemical reactions. Much combustion model development of

late has involved combustion for jet flames where the combustion scale is of the order of the smallest scales of

turbulence. The chemical reaction rate is very fast

compared to the mixing rate. Unfortunately, combustion is

distributed in practical combustors due to the use of slower

reacting, more complex fuels and high turbulence levels.

The expertise in computational combustion is divided

into two main groups. One group very closely models

thermodynamic and transport properties. Turbulence is

generally ignored by this group. If there is transport, it

is laminar, and species diffusive velocities, viscosities,

etc. are precisely modeled. Equilibrium, partial

I!!li
37

equilibrium, or complete system of chemical reactions may be

solved. Equilibrium assumes infinite time for chemical

reaction. Species concentrations do not change after

reaching equilibrium. Reactions are still proceeding, but,

forward reactions are perfectly balanced by reverse

reactions. The final chemical species for equilibrium

combustion calculations are found by minimizing the

Helmholtz or Gibb's free energy depending on whether it is a

constant-volume or constant-pressure process. Finite rate

combustion involves time stepping through a large number of

ordinary differential equations, which are usually stiff.

Small time steps must be used to retain numerical stability.

The combustion process is viewed as a collection of a number

of elementary reactions involving likely combinations of

reacting species. Each of these elementary reaction steps

has its own constants for activation energy, pre-exponential

constant, the species concentration dependence and possibly

some extra temperature dependence. Many of these elementary

reactions are reversible and many of the species are in

numerous competing chemical reactions. Some reaction rate

constants may be calibrated using shock tube experiments.

Other reactions are currently unknown, so combustion rates

are inferred. When the kinetic scheme duplicates shock-tube

data over a range of conditions the scheme is said to be

accurate or validated. Not too surprisingly, there are

multiple schemes which have been semi-validated for the
38

combustion of simple fuels. Improving the accuracy of the

model is done by adding more reactions, which increases the

applicability of the kinetic combustion model over wider

temperature, pressure and species concentration ranges.

There are multiple competing paths in which fuel and oxidant

are burned. The influence of these competing reaction paths

is highly affected by temperature. At different temperature

ranges, different reaction paths can take precedence.

Including multiple reaction paths extends the validity of

the kinetic model. The current state of chemical kinetic

modeling is that the combustion of hydrogen and some of the

lower hydrocarbons is quite accurate. Laminar flame speed

and heat release are well predicted for the very simple

fuels. The kinetic combustion modeling of heavier fuels

used in practical combustors, such as kerosene, is still in

a development stage. Kinetic modelers will not advocate

combustion schemes for heavier fuels until the combustion

models for the lighter fuels are sufficiently validated, as

the heavier fuels break down into lighter fuels and other

intermediate species during the combustion process.

The fine degree of kinetics combustion modeling does

not lend itself to calculation of multi-dimensional

turbulent combustion flow. They are simply too

computationally expensive. Thus, the implementation of

kinetic reaction submodels into multi-dimensional CFD models

has only involved hydrogen or methane with no turbulent-

!il I i
39

chemistry interaction. That is, mean quantities are used in

the reaction rates, there is no turbulence correction.

Currently, the solution of species transport for

recirculating flows using CFD is such that the calculation

of individual species without chemical reaction is

inaccurate. That is, the transport of one species can

easily be off by over ten percent, a level of accuracy CFD

calculations should meet to have predictive validity in the

design of low-emission combustors. Designers would like to

be able to predict within a few percent. Various combustion

models have been implemented to predict combustor flows in a

reasonable time.

3.1 - Infinite Rate Combustion Models

The simplest CFD combustion model is the fast or

infinite rate reaction model. This type of scheme does a

good job of predicting the exit temperature profiles of

combustors. In this scheme, once the fuel and oxidant are

mixed, the reaction rate is such that the reactants are

immediately converted to products. The reactants are

consumed in strict proportions. A simple example is the

stoichiometric burning of methane. The reaction is referred

to as stoichiometric as it assumes stoichiometric burning of

CH 4 + 2 _- CO 2 + 2 H20

the fuel and oxidant. One molecule of methane is consumed

with two oxygen molecules to produce one molecule of carbon
4O

dioxide and two of water.

react with 64 Kg of oxygen.

each Kg of methane.

related.

rate term is eliminated.

conserved scalar is:

In units of mass 16 Kg of fuel

Four Kg of oxygen react with

The reaction rate source terms are

Rfu = Rox/4 (3.1)

Using the general transport equation for species,

equations may be written for the mass fractions of oxygen

and fuel. Subtracting a constant specific portion of the

transport equation of oxidant to that of fuel eliminates the

chemical reaction source term, leaving a conserved scalar.

Using a conserved scalar eliminates calculating troublesome

reaction rates. In the previous example, if we form a

variable composed of the fuel mass fraction minus one

quarter of mass fraction of oxygen the troublesome reaction

The general equation for this

f = YF-(F/O)stYo (3.2)

The term (F/O)st is the fuel to oxidant molecular weight

ratio for stoichiometric burning process. Other conserved

variables could be formed. Another popular conserved

variable is the sum of all fuel fractions, both unreacted

and reacted. This conserved variable is often referred to

as the mixture fraction. The solution of the conserved

variable defined in (3.2) gives the fuel and oxidant mass

fractions. If the conserved variable is positive, the fuel

!!IIi
41

mass fraction is the same as the conserved variable, and the

mass fraction of oxygen is zero. If the conserved variable

is negative, the fuel mass fraction is zero and the mass

fraction of oxygen is 4 times the negative of the mixture

fraction. The sum of the mass fractions of all species is

one so the mass fraction of products is simply:

Ypr = 1.0 - YF- Yo (3.3)

An atom balance is done to figure the individual product

species fractions. If the oxidant is air, the above

equations must be slightly modified to include nitrogen. In

this infinite rate combustion model, the width of the

reaction zone is zero. Combustion occurs in a flame sheet.

There is no overlap of oxidant and fuel. Fuel exists on one

side of the combustion sheet, and oxygen on the other.

An enthalpy equation is solved from which temperature

can be found from the various mass fractions. Radiation is

typically ignored as research combustors are usually

operated at lower temperatures and pressures. The equation

of state is used to solve for density:

p =PM/RuT (3.4)

The equation uses the universal gas constant and the mixture

molecular weight (inverse of sum of number of moles of each

species divided by each species molecular weight). The

mixture molecular weight in terms of mass fractions is:
42

M = (re°x- + mfu + __mpr )-I (3.5)

M:u M;r

Various improvements can be included, such as allowing for

variable specific heat with temperature. As specific heat

increases with temperature, this change alone can improve

the prediction of temperature by a couple of hundred

degrees. The use of infinite reaction rates is known to

give good predictions of exit combustor temperature

profiles, and major species concentrations. Partially

reacted species or nitrous oxides are not predicted by

infinite rate combustion models. Rather than assume

complete chemical reaction, equilibrium combustion can be

modeled. At equilibrium the species concentrations are

fixed, but, the rate of forward chemical reaction is just

balanced by the rate of back or reverse chemical reaction.

Equilibrium modeling allows for the presence of minor

chemical species, allows reactant species to overlap and a

broadening of the reaction zone width. Carbon monoxide is

usually overpredicted with rich fuel mixtures. Equilibrium

combustion models significantly overpredict chemical

reaction and product species at low temperature. Formation

of some species such as nitrous oxide is a function of time.

The equilibrium model does a very poor job of predicting

nitrous oxides.

One model that allows for overlap of reactant species

at low temperature is the finite rate combustion model.

It!I [
43

3.2 - Finite Rate Combustion Models

Finite rate combustion means that fuel and oxidant can

co-existing at the same time. A conserved scalar is still

conserved in finite rate combustion modeling, but, the

concentration of fuel or oxidant cannot be solely obtained

from the conserved variable as before. An additional

equation must be solved for that includes a finite rate

reaction source term. The calculation of the finite rate

reaction term is usually based on arrhenius reaction rates.

Arrhenius was the first to include the Boltzmann factor into

chemical reaction rates. At room temperature oxygen and

hydrogen can co-exist with no reaction. Energy must be

supplied to begin the combustion process. Arrhenius

postulated that collisions between fuel and oxidant do not

automatically produce combustion. The collision between

molecules must be strong enough that an energy barrier is

overcome. The Boltzmann factor takes this into account.

The equation:

Ea

k = exp(--Ru T)

(3.6)

is called Arrhenius law. A typical Arrhenius reaction rate

may look like

dC_u -E a

d-----_ = -A (Cfu)a (Cox)bex p (_uuT ) (3.7 )

where C is the molar concentration. By employing the
44

perfect gas law, the expression may be algebraically

manipulated to use mass fractions which is more commonly

used in CFD. A is the frequency number to take into account

the collisions between molecules. Sometimes a dependency on

temperature to some power is also included. For fundamental

or elementary chemical reactions, the exponents a and b are

whole numbers from the molar coefficients in the chemical

reaction equation. Fundamental chemical equations are

reactions that actually occur. The combustion of hydrogen

with oxygen isglobally represented by:

H 2 + 1/2 _ -H20

The actual combustion process occurs as a series of

reactions of which only a few are:

OH + H_ -. H20 + H

H+ 02 " OH+ 0

O +H 2 -" OH+ H

OH + OH + M _ 0 + H20

OH + HO 2 _ 02 + H20

H+ 02 + M " M + HO 2

There are numerous elementary reactions to get the final

product, H20. In the process there are many reaction

radicals which are also in many elementary competing
45

reactions. To further complicate things, the intermediate

species, H202, is also important. R.J. Kumar n, and S. W.

Kim 2_ have successfully used a set of 24 elementary

reactions involving 9 species to predict hydrogen-oxygen

combustion. This set of elementary reactions was reduced

from a set of 48. At different temperature ranges different

reaction paths take dominance. There is wide disparity in

the reaction rates which causes great difficulty in their

solution. Some kinetic reaction schemes are validated by

duplicating calculated flame structure, laminar flame

speeds, and chemical composition using large numbers of

fundamental reactions. Computational models used to do this

are one-dimensional, include modeling of individual species

laminar diffusion velocities and use very small time steps.

In multi-dimensional recirculating flow global reaction

models are typically used. The global reaction models are

only calibrated over limited temperature, species

"concentration and pressure ranges. In a 1972 paper, Dryer

and Glassman 24 used a temperature range of 1030 to 1230 K

for oxidation of carbon monoxide and methane. The most

often cited global reaction schemes are from Westbrook and

Dryer 2s. The reaction schemes were designed to capture

rich and lean flammability limits, information on flame

temperature and burned gas composition and flame speed.

They do not accurately describe the chemical structure of

the flame. Different models can produce different results
46

with the same reaction numbers. Coffee 2_ recommends higher

frequency numbers as his results were slightly different

than Westbrook and Dryer's, unless he adjusted the reaction

numbers. There is further controversy in regards as to what

phenomena should have primary importance in modeling.

Coffee 27, in another report, states that heat release

should be used to calibrate reaction numbers, not the

criteria used by Westbrook and Dryer. The frequency number

given in this paper is orders of magnitude higher.

If one wants to calculate intermediate species such as

carbon monoxide, at least one additional equation must be

solved. Nikjooy et al. 2_ gives an example with the

variables and source terms. A four step scheme has been

used by Srivatsa 29. This scheme is from Hautman et al. 3°

and allows for additional species of Hydrogen and ethylene.

This scheme includes product concentration in the reaction

rates. The scheme was validated over 960 to 1540 K for

propane.

An early study of global reaction schemes was done by

Abdalla et al. n One conclusion was that "inaccuracies

could arise from the unrealistic simplification of chemical

kinetics and the variable effects of turbulence". Duterque

et al. n also came up with a global scheme for methane,

propane, benzene and isoctane validated for a well-stirred

reactor model. The authors stated that the model could be

used in the predictions of burners and combustors.

If!I!T
47

In a more recent paper, Westbrook and Dryer 33 discuss

detailed kinetic modeling for various hydrocarbon fuels

along with a short section on global modeling. They note

that using a single reaction scheme leads to overprediction

of products and temperature. In 1986, Paczko et al. 34 gave

a four-step reaction mechanism for methane and a five step

reaction mechanism for propane. The propane mechanism

includes two reactions for the initial break up of propane.

A three-step mechanism is recommended by Seshadri and

Peters 3s for stoichiometric methane-air flames. The

mechanism was tested for a counterflow diffusion flame. In

1987, Kiehne et al. 36 discuss an improved eight-step

mechanism for propane. The model is based on the original

work of Hautman et al. The model better predicts ignition

delay and combustion at higher pressures. In 1988, Jones

and Lindstedt 37 proposed a four-step global mechanism for

the combustion of alkane hydrocarbons up to butane. The

scheme was tested on methane and propane flames, along with

diffusion flame data for a methane-air flame. Good

agreement was predicted for flame speed, flame thickness,

and species profiles. Bilger and Starner 38 constructed a

four-step mechanism from a set of 58 elementary reactions.

The model uses steady-state approximation for the oxygen

atom rather than a partial-equilibrium approximation. The

mechanism was used to predict flame structure for the

counterflow diffusion flame. Good predictions were obtained
48

for minor species over a range of flame stretch.

Some of these schemes involve the implementation of

source terms for calculating product species. This

researcher found problems when calculating product species

for multi-dimensional recirculating flows. The frequency

numbers used in arrhenius rates are fairly large. At higher

temperatures, arrhenius reaction rates can predict rapid

reaction even though reactant concentrations are quite low.

The arrhenius reaction rate is only accurate for an

extremely short time step, not the typical time steps used

in implicit calculations. Thus, product concentrations can

be overpredicted at high temperatures. Most often,

overpredicted concentrations are simply overwritten during

the iteration process.

3.3 - Turbulent Combustion

All practical gas turbine combustors operate in the

turbulent regime. Thorough mixing and burning must occur

•within rather short axial distances in current gas turbine

combustors. Turbulence for chemical species means that

there is variation with time in the chemical species and

temperature. A consequence of this is that even with

infinite reaction rates, fuel and oxidant can exist at the

same location, but, not at the same time. Where this

occurs, the time averaged concentration of fuel and oxidant

will not be zero. Many different turbulent combustion

models have been worked on in the last two and a half

If!1 i
49

decades. Many of these models do well within the range of

their assumptions.

3.3.1 - Eddy Break Up Combustion Models

One of the earliest turbulent combustion models is the

eddy-breakup model proposed by Spaulding 39. The original

application for the model was for premixed combustion. It

was observed that changes in velocity and fuel to air ratio

had little affect on combustion in premixed combustors. The

use of so-called laminar arrhenius reaction rate combustion

models didn't predict this. The controlling mechanism

wasn't temperature related, but, was controlled by mixing

rates. In the eddy-breakup model combustion is viewed as a

mixture of burned and unburned eddies. The combustion

process is controlled by the mixing and exchange of energy

from burned eddies to neighboring unburned eddies. The

eddies are stretched and folded around each other until

local gradients are strong enough that mixing occurs. This

is very similar to the cascade of turbulent kinetic energy

from large eddies to smallest eddies where the energy is

finally dissipated. It was further proposed that combustion

was also proportional to reactant concentration. If oxygen

concentration was low, as if most of the oxygen were already

consumed, the reaction rate should be lower. The initial

combustion model did a good job of qualitatively predicting

the correct hydrodynamic phenomena for premixed combustors.

The eddy-breakup model was modified over the next few years.
5O

In another paper 4°, it was proposed that reaction rates

should be proportional to the concentration fluctuation not

the reactant concentration. An equation for concentration

fluctuation was proposed and successfully tested for a

turbulent jet. In another paper, Mason and Spaulding 4_

tested the combustion model. From a later review 42, it was

reported that the improved more physically realistic model

predictions were worse than previous modeling. In 1976,

Magnussen and Hjertager 43proposed an another eddy-breakup

model where the reaction term was proportional to the

minimum of related constants times appropriate fuel,

oxidant, or product concentration. The model performed well

results, but, the constants had to be increased for a

different case. Isothermal turbulence modeling was also

developing. The k-e turbulence model was developed and has

proven successful. Using the k-£ turbulence model, the form

of Spaulding's eddy-breakup model is:

(3.8)

The reaction

A measure of the rate of eddy breakdown is the reciprocal of

eddy lifetime. A measure of eddy-lifetime is the ratio of

turbulent energy, k, over the dissipation rate, E.

The eddy-breakup rate was developed for premixed

combustors but has been applied to other combustors. In

work at Garrett research 44 the eddy-breakup model was

combined with an Arrhenius type reaction rate.

li!! i
51

rate was taken as the minimum of the Arrhenius rate or the

eddy-breakup rate. Multiple step chemical reactions were

also implemented. The combination of eddy-breakup model and

arrhenius model is easy to implement. The model was

superior to earlier models.

3.3.2 - Probability Density Function Models

The eddy-breakup model has been long regarded as a

semi-empirical combustion model, even though it has proven

successful. In turbulence modeling, the current thought is

to include as much of the physical process as possible to

improve the models. Thus, there has been great effort to

model various turbulent terms in the averaged Naviar Stokes

equations for chemical species. One method is to

incorporate a correction factor for turbulence which can be

applied to the laminar Arrhenius reaction rate. Employing

Reynolds decomposition and the perfect gas law, a typical

turbulent reaction rate may look like:

-BP2 (_u+Ylfu) (_ox+Yi)exp( -Ea

_fu = 2(T+TI) 2 Ru(_+T/)) (3.9) Ru

Using _ = EJR u and some algebraic manipulation, the

exponential part of the above expression is:

exp( -Ta -Ta Ta T '/_ ) (3.10)

÷ TI)= exp(--_--)exp( T 1 + T'/T-
52

The second exponential may be expanded in an infinite

series. This infinite series is only convergent if _ is

not large compared to the averaged temperature and if the

variation in temperature is not large compared to the

averaged temperature. This is not probable in practical

combustors. An additional problem with this method is

modeling the correlation term Y'_uY'ox . One way around the

reaction rate closure problem is to treat turbulent reacting

flow as a stochastic process 4s. That is, the solution of

average scalar quantities involves the weighting of the

scalar with the probability of thescalar.

1 = /oiP(f) df (3.11)

Yi = Yi (f) P( f) df (3.12)

Here it is implicitly understood that the probability

density function is also a function of time and position.

The first stochastic methods used in multi-dimensional

turbulent combustor flow involved the use of two-parameter

probability density functions. The pdf is completely

specified by two parameters, the mean and variance of a

conserved scalar. These are referred to as the first and

second moments.

One of the basic assumptions of this model is that the

li|I_
53

reaction is very fast or near equilibrium. This allows the

use of a single variable, mixture fraction, to specify the

mixture. Physically, the fast reaction assumption means

there is little or no overlap of reactant species at any

instance in time. The incorporation of stochastic modeling

of turbulence allows the overlap of reactant species over

time. Much like the eddy-breakup model, reactants occur in

separate parcels. Over time these parcels move around so

that the time average of reactant concentrations overlap.

This occurs with mixtures near stoichiometric conditions

causing a broadening of the reaction zone. As actual gas

turbine combustion zones are distributed, this modeling

better portrays actual combustion.

Typical two parameter probability density functions are

rectangular or battlement pdf, Gaussian pdf, clipped

Gaussian pdf and Beta pdf. Initially, rectangular or

battlement probability density functions were used due to

their simplicity. Gaussian pdf's are more appropriate in

highly turbulent flows, but, are more complicated to

calculate. Gaussian pdf's are determined by their mode or

most probable value and standard deviation, not the mean and

variance which is found by solving their CFD transport

equations. There is no direct way of solving for the

Gaussian pdf from the scalar mean and variance. The pdf

must be specified so that other scalar averages, such as

temperature and density can be solved. The most efficient
54

way of solving for the Gaussian pdf is to construct a table

lookup of mode and standard deviation versus mean and

variance. Thus, when the mean and variance of conserved

scalar are found during the iterative process, an

interpolation of the table is done to find mode and standard

deviation which specifies the gaussian pdf.

One problem with Gaussian pdf's is the tails which

extend to infinity. These tails can correspond to

physically impossible states. The clipped gaussian was

developed to take care of this. Delta functions are added

to take the place of clipped tails. The solution of clipped

gaussian pdf's is slightly more difficult. A pdf which

resembles the gaussian pdf is the Beta pdf, which was

proposed by Richardson 4_. The Beta pdf has been

successfully used by a number of researchers 4v. Unlike the

Gaussian pdf, the Beta pdf can be directly calculated from

the mean and variance.

The conserved variable or mixture fraction is typically

non-dimensionalized for calculations. The mixture fractions

at the fuel and oxidant inlet streams represent the range of

mixture fraction in the combustor.

8= [YF- (F/O)s_Yo]M- [YF- (f/O) s_Yo]a (3.13)

[zF - (F/o)s Zo] - [YF - (F/o)stzo]

_{

The mean mixture fraction is found by solving the averaged

![I:Ii
equation for the mixture fraction.

fraction is often referred to as g:

55

The variance of mixture

_ 1

g = 8 ''_= f(8-8) 2Pr (8)48

0

(3.14)

where the overbar refers to Favre or density averaging

rather than the normal time averaging. A modeled transport

equation is solved for the local variance. The formula for

F(x) -

beta pdf is:

F (a+b) x a-1(l-x) b-1

F(a)F(b)

where 0 <x<l,

and a>0, b>0

(3.15)

the formulas to find a and b from mean and variance are:

a = _[ _(i-_) I], (3 16)

g

b - a(i - 8)

8 (3.17)

The major limitation of this method is that computed

mixtures are close to equilibrium. Equilibrium assumptions

tend to overpredict CO concentrations. Another choice in

combustion models is the laminar flamelet model. In

reviewing experimental work done by Tsuji and Yamoka 4e,

Bilger 49 noticed that measurements of temperature and

species concentration, including CO appeared to be functions

of mixture fraction for laminar diffusion flames despite

differing shear rates and positions. That is, there exists
56

a state relationship to the mixture fraction. Thus, a table

of species concentration and temperature versus mixture

fraction can be constructed from the experimental data.

This is very useful in computations as the table can be used

to find thermodynamic quantities during the iteration

process. This eliminates calculating reaction rates and

enthalpy. Instead of assuming infinite rate or chemical

equilibrium, actual experimental data is incorporated into

numerical predictions. Jeng 5° used this in his thesis to

perform buoyant free jet combustion calculations using the

GENMIX sl computer program. Jeng also modified an

equilibrium program to check the experimental concentration

data. The program was modified for rich combustion. Very

high carbon monoxide concentrations were predicted,

characteristic of high equivalence ratio equilibrium

combustion calculations.

This thermochemical alternative to chemical

equilibrium was further developed into the flamelet model by

Liew et al. s2 In this model, reaction rates are very high

compared to the mixing rate. An instantaneous snapshot of

the combustion process would show a flame with finite

thickness or a flame sheet. The interaction of turbulence

with this flame surface involves movement and wrinkling of

the flame surface. In a geometrically similar type

combustor, laminar diffusion flame data is measured. These

data are then used to provide unique thermochemical

li[li
57

relationships solely as a function of mixture fraction.

These relationships are then coupled with an assumed shape

pdf to account for turbulence interaction with combustion.

The actual flame is viewed as an ensemble of laminar

diffusion flames. This method has produced substantial

improvements in prediction in fuel-rich regions of an open

combustor, particularly in CO concentration.

The laminar diffusion flamelet idea was further

developed by Liew and Bray s3 for highly stretched

combustion. It was found that if flames were substantially

stretched, the state relationship between thermochemical

values and conserved scalar was affected. As the flame is

stretched or instantaneous scalar dissipation is

substantially increased, it was noticed that temperature

decreased. Also, the mixture fraction at which maximum

temperature occurs, increases with increasing stretch. The

region surrounding this maximum temperature is the region of

greatest chemical reaction. Beyond a certain level of flame

stretch, the imbalance between local heat loss, local active

species loss, and smaller heat production due to flame

stretch leads to flame quenching or blowoff. The

thermochemical or state relationships now involve two

variables, mixture fraction and flame stretch or average

instantaneous scalar dissipation. Conceptually, a library

of flames is measured or calculated as a function of two

variables. Turbulence effects are still handled by an
58

assumed shape pdf. However, the pdf is now a function of

two variables, it is a joint pdf. The shape of joint

probability density functions are much more difficult to

deduce. The simplest joint pdf is to assume that the two

variables are completely independent. That is, the joint

pdf is a product of pdf for mixture fraction and another pdf

for scalar dissipation. Liew and Bray used a Beta pdf for

mixture fraction. They assumed that the pdf for scalar

dissipation was log-normally distributed at high Reynolds

number. Equations are then formed for the Favre average

scalar dissipation and other quantities needed to specify

the pdf for scalar dissipation. This model allows the

prediction of oxygen penetration into fuel jets where there

is large strain. The chemistry is laminar and can be

calculated very accurately. Drake and Blint 54 found that

finite rate kinetics have an important effect on flame

structure. The use of partial equilibrium modeling for

carbon monoxide concentration in low temperature fuel-rich

flames was shown to be inadequate. The local effects of

flame stretching correspond to a turbulent ensemble to

average the flamelet characteristics.

Flamelet models are based on a one-dimensional

approximation of the flame. The curvatures of the various

flamelets are neglected. Coherent vortices which are

present in many mixing layers are neglected. Borghi 5s, in

his review of turbulent combustion modeling, states that

]i[li
59

flamelet modeling should be confined to cases where the

chemistry is very fast and the Reynolds number is not too

large. Peters 56 did a review of the laminar flamelet

approach which describes the method more fully. The method

is an attractive way of treating the prediction of minor

species in turbulent flow calculations at reasonable

computer costs.

Janika and Kollman 57 used a two-variable pdf scheme in

their studies of turbulent H2-air flames. The two pdf

variables are mixture fraction and reaction completeness.

This method assumes partial equilibrium for its H2/air

chemistry. The kinetic mechanism selected consists of fast

shuffle reactions which are considered to be in partial

equilibrium:

H + O2 - OH + O

(RI)

O + H2 - OH + H (R2)

H 2 + OH', H20 + H (R3)

20H , H20 + 0

and slow three-body recombination reactions.

H+ OH + M-. H20 + M

(R4)

(R5)

H+ 0 + M- OH+ M (R6)

H+H+M-H 2 +M (R7)

The authors stated that the partial equilibrium scheme was
not good below 1200 K. To reduce the number of variables,

combined variables are used. This was proposed by DixonLewis et al. 58 The combined variable for hydrogen is:

3_

c;_: c_ ÷ 1 _2Cox÷--Co ÷ cH (3.18)

2 Mow M o 2 M H

The reaction rate for the combined variable is independent

of the fast shuffle reactions (RI-R4).

the hydrogen combined variable are:

6O

The source terms for

(_2 = -MH2(2(_5 + 2_6 + 2_7) (3.19)

The combined variable varies between its unburned and

equilibrium values. Thus, a nondimensional reaction

progress variable is:

,11- C;2 - C_2,u (3.20)

_U

Transport equations are used to solve for mixture fraction

and reaction progress variable. The solution of enthalpy,

mixture fraction, and reaction progress variable, coupled

with partial equilibrium relations completely specifies the

instantaneous thermochemical state. Turbulence is handled

by a combination of assumed shape probability density

functions. The mixture fraction and combined variable are

assumed uncorrelated so that the joint pdf is the product of

the pdf of mixture fraction and another pdf for the reaction

progress variable. Transport equations for the mean and
61

variance of the mixture fraction and progress variable are

solved. The transport equation for mean reaction progress

variable must include a reaction term for the recombination

reactions. A beta pdf can be used for the mixture fraction.

Three dirac delta functions can be used for the pdf progress

variable pdf corresponding to unburned, equilibrium and mean

of the reaction progress variable:

P(n) = c18(_) + _8(_-_) + _8(n-l) (3.21)

The coefficients of the dirac delta terms are determined

from simple functions of the mean and variance of the

progress variable. From Janika and Kollman's $9 1982 paper:

cI = -_--_._, (3.22a)

c 2 = 1 - 8,2 ; (3.22b)

(

c3 - _:_ (3.22c)

where favre averaging is actually used.

The scheme was used to calculate the jet diffusion flame of

Kent and Bilger 6°. The reaction progress variable varied

between 0.88 and 1.0 across the calculation at x/D of 20.

The progress values increased further downstream. Janika

and Kollman state that the more elaborate modeling used in

the calculation didn't significantly improve prediction of
62

stable components and temperature. However, the modeling

did considerably help in the prediction of monatomic oxygen

and reasonable NO concentrations.

Correa et al. 61 include a couple of CO reactions in

the above scheme for calculation of CO/H2/-air turbulent

diffusion jet flame. The calculation method successfully

predicted experimental measurements of superequilibrium

concentrations of OH. A year later, Correa 62 published

another paper where a progress variable was added for CO

concentration. The partial equilibrium equations used in

the above models are only valid above at least 1200 K. The

correlation assumptions made about joint probability density

functions is questionable. Pope and Correa 63 found that

joint pdf of velocity, conserved scalar and reaction

progress variable were correlated. The kinetics models for

higher molecular weight fuels are much more complicated.

Including additional variables into the joint pdf for

breakup of the heavier fuels into lighter components is even

more questionable and untested. Borghi in his 1988 review

paper generates partial pdf shapes for some situations. At

best, these shapes only qualitatively resemble Gaussian like

profiles. The next step in turbulent combustion

calculations is to actually remove some of the assumptions.

Some recent papers were written about the current status

of Combustion modeling. This include Hukam Mongia's

paper 64 which covers some of the calculations done at the

_,_! i
63

Allison Gas Turbine Division of General Motors Corp. Correa

and Shyy [8] covers the combustion modeling effort done at

General Electric's Corporate Research and Development

Center. Much of this effort was to develop their CONCERT

code. There are also quite a few texts on combustion which

include: Kuo's 6S text on the principles of combustion,

Williams TM reference on combustion theory, Bartok and

Sarofim's _v book on fossil fuel combustion, Chigier's 68

reference on energy, combustion and environment, and Oran

and Boris '69 reference on CFD simulation of reactive flow.
64

Chapter 4

HYBRID COMPOSITIONPDF SOLVER

Many assumptions have been made in developing the

previous assumed shape pdf models. Implicit in those

assumptions was that chemistry is very fast. Even

correcting this by using partial equilibrium solutions for

jet flames is only valid at higher temperatures and thus, a

high degree of reactedness. Another assumption was the

shape of the pdf. The focus of this work will be to further

develop and test a pdf method which does not make these

assumptions.

For this work, the pdf is initially assumed to be a

function of time, position, temperature or enthalpy, and

multiple species concentration. This is referred to a

composition joint pdf. The pdf is called joint because the

pdf is a function of more than one variable. The hybrid

composition pdf uses traditional CFD techniques to obtain

solutions of velocity, pressure, turbulence kinetic energy

and dissipation. The composition pdf handles the solution

of species concentration and temperature. The hybrid

composition pdf solution is thus a combination of

traditional CFD methods and stochastic modeling. If the pdf

is assumed to be also a function of velocity, the solution

method is vastly different from traditional CFD techniques.

Thus, the composition pdf represents a good intermediate

step to complete pdf modeling, which is in its infancy. In

il ! i:
65

1989 Anand et al. v° used the joint velocity pdf method to

calculate the turbulent flow behind a backward facing step.

In this calculation only the pressure field was needed from

another solver. Preliminary work is being done to calculate

all quantities, including pressure in the pdf solver. In

1993, Anand, Pope and Mongia presented predictions for a

swirling flow. The method did not incorporate a radial

dependence for pressure, although an externally imposed

pressure gradient could be incorporated into the equations.

The mean swirl number for the flow was reported to be 0.09.

As it has been shown that turbulent combustion models

are rather simplistic and full stochastic methods are not

developed, the hybrid pdf represents an excellent testbed

for testing predictions of turbulence effects on combustion.

Various pdf transport equations can be written. One

form for the composition pdf is:

N

=

N N

(4.1)

The terms represent the time rate of change, mean

convection, chemical reaction, turbulent convection, and a

mixing term, respectively. The symbol over P denotes Favre

averaging.
= pP/F (4.2)

66

ei. j is the scalar dissipation:

_i.j = Da_iaa_j (4.3)

where D is the diffusion coefficient. The use of favre or

density weighing eliminates some terms in the averaged

equation. The notation <xly> denotes the mathematical

expectation of a random function x conditioned on y.

Terms on the left side of the pdf evolution equation do

not need modeling. That is to say, they don't have

additional unknown variables which demand assumptions or

relations to produce closure models. The time rate of

change and mean convection are seen in traditional CFD

modeling. The last term on the left side is the evolution

of the pdf in transform space due to chemical reaction.

This term does not employ modeling assumptions about various

turbulence quantities such as _' and T'. The pdf method

can use the laminar kinetic reaction rate modeling. As was

mentioned previously, using regular moment methods to solve

for chemical reaction produces correlations that have to be

calculated or modeled. These correlations correct the

laminar chemical reaction rate for turbulence. All of these

correlations involve modeling assumptions. The pdf method

does not include these assumptions. This is why pdf

modeling is promoted as exactly modeling chemical reaction.
67

The terms on the right side of the pdf evolution

equation must be modeled. The turbulent convection term is

modeled by gradient diffusion:

_<ual_k>p , Dc_ ' (4 4)

where D_ is the turbulent diffusion coefficient, which is

related to the eddy viscosity by the turbulent Schmidt

number. The next term that needs modeling is the molecular

mixing term. A proper molecular mixing model causes the

scalar variance to decrease with time while the mean remains

constant. Some molecular mixing models are: the

coalescence/dispersion model attributed to Curl _I, modified

Curl model of Janicka et al. v2, and the quasi-Gaussian

model of DopazoV3. The general form of the

coalescence/dispersion model is:

n n

i=i j=1

M( PD

(4.5)

where T is the probability that particles will undergo a

transition in concentration between one particle with

concentration of _' and another particle with concentration

_" to produce two particles with concentrations _ and _'+_"-

respectively. In Janicka et al.'s paper, a general

equation for transition probability is developed using the

concentration mixing equations developed by Curl for droplet

mixing. Using the simplest approximation, the transition
probability is:

68

II 1

0 otherwise, "

This approximation assigns equal transition probabilities

for each point in the range _' to _". The time scale I; for

this process is taken as the integral time scale, _=k/s.

Traditional finite difference techniques are not

feasible for solving this pdf evolution equation as the pdf

is a function of a large number of variables (five species

plus temperature or enthalpy). Instead, a Monte Carlo

method will be used to solve for the pdf evolution.

4.1 - The Monte Carlo Solution Method

The Monte Carlo method is easy to implement, even if

the number of variables is large. The evolution of the pdf

entails the movement of particles in physical space as well

as the scalar space (_ space). The movement of the

particles is governed by the pdf evolution equation. The

solution to the pdf evolution equation is a finite set of

particles or events for each grid cell. The actual pdf is a

continuous function. Increasing the number of particles

does not result in a continuous solution. Another problem

with Monte Carlo is the level of solution accuracy. The

ensemble average a solution particle property (species

concentration, temperature, enthalpy or density) for a grid

cell is an estimate of the true mean since the expected

IF1li-:
69

value of the ensemble average is the mean. It can be shown

that the solution accuracy is proportional to the square

root of the number of particles. Thus to get an order of

magnitude increase in accuracy the number of particles must

be increased by two orders of magnitude.

Mathematically, the particles for the Monte Carlo pdf

are represented by a set of delta functions. The

distribution function is the integral of the probability

density function. The Monte Carlo distribution function is

a sum of heavy-side functions. An example from Pope's 1985

paper TM is shown in figure 4.1. A Gaussian pdf and a set

of 20 normally-distributed random numbers are compared. The

mean is 0.5 and the standard deviation is 0.01 for both

probability density functions. It can be seen that the

Monte Carlo distribution function approximates the

continuous distribution function in a step-wise manner. The

delta functions are represented by finite spikes. The

actual height of the spikes is infinite and the integral of

each spike is 1/20.

Using fractional time-steps, the pdf solution is

divided into 3 phases. The time derivative is first

discretized as:

= /At

then the pdf evolution equation can be written as:

(4.7)
70

DISTRIBUTION FUNCTION

0.1 0.3 0.5 0.7 0.9

4

2

0

Probability Density Function

/ %

/ \

/ \

/ \

ILi1112111

0.1 0.3 0.5 0.7 0.9

Figure 4.1 Comparison of continuous and discrete pdf's and

distribution functions

ii1 IF
71

Using approximate factorization, the above equation is

recast as

(4.8)

-(l+At_ (I+AtM) (I+_tR) P = (4.9)

where C denotes the convection/diffusion operator, M is the

molecular mixing, and R is the chemical reactions. Using

the above expression the different processes may be

performed consecutively:

(i) convection/diffusion between grid cells:

= (4.10)

(2) molecular mixing within each grid cell:

P_ = (I+AtMgP*, (4.11)

(3) chemical reaction in each grid cell:

_,÷i = (I+AtR)_. (4.12)

In the Monte Carlo pdf calculation, each particle has

its own species mass fraction and temperature-enthalpy.

Each particles has the same mass. Particle position within

gridcells is not modeled. In the Monte Carlo simulation a

continuous pdf is replaced by N x m delta functions,
72

N

in_ 8 (_ __m) x8 (_,-_n)) "'_ (_m-_")) (4.13)

each product of the m delta functions representing one event

of an ensemble of N sample events. The pdf is a function of

position and time.

The Monte Carlo pdf formulation/solution is initialized

with a SIMPLE method solution. The pdf solution technique

could be started with an unconverged flowfield, but, the

time needed to get a converged solution would be much

higher. Each of the initial particles in a grid cell have

identical properties. Inlet particle properties are set to

the inlet boundary conditions. The first step in the

iteration technique for the pdf solution is calculating the

exchange of particles between grid cells by convection and

diffusion.

4.1.1 - Convection/Diffusion phase

Replacing the continuous pdf with the Monte Carlo pdf,

the convection/diffusion step of the evolution equation can

be written as:

 atP, + Ba,P* = (D@.P*)

re-writing this for the x-y coordinate system:

The convection/diffusion step does not have any source

terms.

(4.14)

(4.15)

The SIMPLE scheme solves the transport of species

!_!I J
73

semi-implicitly. That is, very large time steps can be

taken. Solution information travels very quickly in this

system. The Monte Carlo pdf scheme uses the evolution of

particles to calculate scalar properties such as

concentration, temperature and density by ensemble

averaging. The Monte Carlo scheme does not implicitly

calculate particles or particle properties. Particle

properties are calculated from existing neighbor particles

in time and space. Thus, it is incorrect to use implicit

differencing in Monte Carlo simulations. Explicit

differencing must be used, which limits the time step that

may be taken. It should be obvious that the time step must

small enough that any particle does not travel further than

one grid cell in any direction. In the SIMPLE scheme the

weighing of neighboring grid cell properties is based on

coefficients calculated from the convection and diffusion

using an implicit formulation. The coefficients used in the

Monte Carlo solver must be figured using explicit formulas

as the exchange of particles is explicit. For example,

assuming I-D diffusive heat transfer we have as a general

formula:

apTp=aB[fTs+ (l-f) T_] +aw[fTw+ (l-f_ T_] + [a_- (l-f) a B- (l-f) a_ T_

(4.16)

where
74

ap = fag + fa w + a °

(4.17)

and a_ and aw are calculated as before. The letter f is the

weight of implicit differencing. Generalizing this for 2-

dimensional convection and diffusion with f=O for explicit

differencing :

a_p=a_+a_a_ON+a_Os+ [a__a _aw_as_a. ] _ (4.18)

where

ap = a ° = pc_xAy (4.19)

At

thus, the time step must be small enough so that the

weighing for the previous time step is non-negative (which

introduces non-physical solutions). A maximum step is

calculated for each iteration by finding the cell with the

smallest allowable time step and using that value.

Dividing the above equation by ap and using the Monte

Carlo pdf with i,j notation, the convection/diffusion step

is:

p_,j o o o eP_,j (4 20)

In this step the new Monte Carlo pdf is made of a fixed

number of particles. This is a combination of particles

from neighboring cells and the particles from the same cell

(particles from the previous time step). The particles from

each cell are randomly chosen. Only whole numbers of

particles are exchanged. In the conversion from a floating
75

point number into an integer number some fraction is lost.

The coding is such that the remainder is added in at random

iterations so that the exchange process is more accurate.

If the remainder were not randomly added in, it is thought

that the iteration process may become cyclic, at least with

a small number of particles per cell. The time step is

reduced so that the possibility of extra convection or

diffusion exchange particles does not cause the number of

particles at the previous time to go non-negative (ie.

exchange more particles than you actually have). Of course,

using high numbers of particles per cell cuts down on the

truncation error. The Monte Carlo scheme does not keep

track of where the particles are within each cell. The

solution quantities to a SIMPLE method calculation (as most

CFD solutions) represent the means and sometimes the modeled

variance for individual cells. The variance represents the

degree of mixing within a grid cell.

4.1.2 - Molecular Mixing Models

In the second fractional time step, molecular mixing is

performed using the continuous mixing model of Hsu. The

modified Curl model has non-gaussian higher moments. Curl's

model assumes complete mixing between chosen mixing pairs.

In the modified method of Curl, the Monte Carlo simulation

randomly chooses N_ pairs of particles based on the

following formula:
- 8t N - 8t8 N

N_ = 0 ._--_ = 0 ._--_ ,

76

(4.21)

where C is 6.0. The pairs are allowed to mix as:

S= ( t+8 t) = ASm (t) + (I-A) Sn (t) (4.22a)

S. (t+8 t) - At. (t) + (l-A) Sm (t)

where A = 0.5_, with _ a random variable uniformly

distributed on the interval [0,I].

particles remain unchanged:

$=(t+St) = $_(t)

(4.22b)

The remaining N - 2N=

(4.22c)

A plot of the mixing process in composition space is shown

in figure 4.2. The value of A used in this plot is 1/6.

This model is discontinuous as the sample particles change

properties abruptly regardless of how small the time

interval is made. Reducing the time interval reduces the

number of particle pairs chosen for mixing, but, does not

affect the total property change for the chosen particles.

A. Hsu and J.Y. Chen vS modified this process to become

continuous. In this model all particles are chosen to

undergo mixing and the time interval is incorporated into

the mixing level term A. Thus, the jump in physical

properties mixing pair particles undergo decreases with the

time step, and thus, the process becomes continuous. There

are N/2 particle mixing pairs for each cell. A is now
77

O

¢-

O

.B

.m

o

0

%(0

. Cn(t+_t)

"0

_m(t+at)

o. ,m(t )

0

0

• 1 composition space

Figure 4.2 Plot of Monte Carlo Mixing process

defined as:

A = c'_st, --E-' (4.23)

where C" = 2.0. Hsu showed that this method resulted in

more appropriate higher order moments, specifically flatness

and superskewness in his calculations. In a 1990 paper,

Chen and Kollman v6, use the modified Curl method for a twostep hydrogen combustor calculation and conclude that the

model did predict cross correlation properties reasonably

well for turbulent shear flows, despite the problems with

non-Gaussian character in the higher moments.
78

Other mixing models have been suggested and employed.

Pope_ suggested that the sampling probability of particles

should be biased in regards to the elapsed time since the

last mixing event. This entails an additional scalar for

the elapsed time since the last mixing event for each

particle and the use of an appropriate weighing function for

age. Another mixing model is the relax-to-mean model. In

this model, the cell mean is calculated, then each

individual particle's variance from the mean is reduced

exponentially depending on turbulent mixing and time step.

The formula for this is:

_n( t+8 t) =_+ (_n(t) -_) .exp (-0.5C'_8 t) (4.24)

This mixing model was used in A. Norris and A. Hsu's v8

Monte Carlo pdf calculation of Masri et al v9 CH4/CO/H2/N 2-

air pilot stabilized jet flame. A somewhat similar model

which uses a linear relaxation is the IEM model developed by

Villermaux, 8° and by Yamazaki and Ichigawa 8_. The formula

for this is:

at + uj axj _z_

+ _i

(4.25)

The model name comes from exchanges by Interaction with the

Mean Value. The term _m is an exchange time. This model

treats combustion simultaneously with mixing. This model

has been used by Correa and Braaten 82 to test reduced

chemistry mechanisms against full chemistry mechanisms for

Ill | 7
79

methane-air combustion in a Partially Stirred Reactor (PaSR)

model. This model was chosen in preference to pair exchange

models because it is easier to employ on parallel computers.

The simulations take many hours to perform on a single

processor. Doing calculations in parallel, that is,

dividing the numerical work between multiple machines will

result in greatly improved turn-around time. The PaSR

represents a step into pdf modeling which can use full

laminar chemistry for a single grid cell or reactor.

The PaSR pdf calculation admits a fixed number of inlet

particles per mixing step. At the same time, the same

number of particles is taken from the PaSR. The PaSR has

two other simpler reactor models as mixing limits. With no

mixing the reactor is simple plug flow. The other limit is

infinite mixing which corresponds to a fully stirred

reactor. Actual multi-dimensional flowfield calculations

can be thought of as a collections of Partially stirred

reactors. The exchange of particles between grid cells

amounts to exchange of particles between neighboring

reactors. The mixing frequency inside each reactor is

governed by local turbulence mixing. The results of Correa

and Braaten's PaSR calculations varied with the frequency of

mixing. Specifying low mixing frequencies resulted in

blowoff for a reactor with separate fuel and oxidant inlets.

Even with steady-state combustion, concentrations,

particularly NO and CO, varied with mixing frequency even
8O

using full-up chemistry. As reducing CO and NO emissions is

of prime concern in engine design, this is somewhat

troubling. Taking individual gas samples amounts to

measuring an average. This also disturbs the flow, which

affects accuracy. Laser Raman Spectroscopy can take nearly

instantaneous concentration measurements without disturbing

the flow, but, the measurements must undergo substantial

corrections. Thus, it may be some time before data are

calibrated well enough to fully validate the pdf method.

In a later paper, Sanjay Correa 83 compared the IEM

mixing model with Curl's model and the modified Curl's model

for premixed CO/H 2 combustion in a PaSR. At a frequency of

316 Hz, the results were similar. From figures 4-6 of his

paper, the difference in temperature between IEM and the

pair exchange models was on the order of 50 degrees. The

difference in predicted CO and OH concentrations was close

to a factor of two. Oxygen predictions were not affected by

the mixing model at this frequency. At a frequency of 1000

Hz, the pair exchange models and the IEM model produced

nearly identical results. As Correa and Braaten predicted

in-flame mixing frequencies of 1000 Hz and larger, Correa

stated that the choice of mixing models was not critical in

the distributed reaction regime of lean premixed combustion,

as long as the turbulent mixing frequencies were above 1000

Hz. This thesis will employ benchmark data from research

combustors. These combustors operate at lower temperatures

IllI_
81

and less intense combustion than practical combustors.

Thus, the IEM model is probably not appropriate to use for

this work. Research combustor data is used in this work as

experimental data from practical combustors is sparse.

Combustion was also effected in a CoalescenceDispersion Model used for infinite reaction cited by J.Y.

Chen and W. Kollman B4. If chemistry is very fast compared

to mixing, as in the flamelet combustion modeling,

combustion occurs within a very limited mixture range close

to stoichiometric. If combustion isn't treated

simultaneously with molecular mixing, combustion can be

vastly under predicted. In scalar space, if only the final

mixing concentration is examined, particles can pass through

a combustible mixture during the mixing phase without

reacting. This problem does not occur with distributed

combustion which is considered for this thesis.

4.1.3 - Combustion phase

Combustion is proceeds after the exchange and mixing of

particles. The pdf method does not have to employ

corrections to reaction rates due to turbulence. The

unmodified laminar reaction rate is used in the pdf

evolution term. Each particle has its own species

concentration and enthalpy. Each particle then undergoes

chemical reaction without having to worry about modeled

turbulence correlations, turbulent flame speed or other adhoc turbulent combustion rate. The laminar reaction rate is
82

the actual reaction rate. However, as mentioned previously,

there is much controversy about true reaction rates. Even

laminar combustion involves tens of reactions with diverse

time constants. Including "exact" kinetics into a pdf

scheme would result in 100,000 "exact" solutions per

iteration for a rough grid of 50 by 80 with 250 particles

per grid cell. Semi-implicit turbulent combustion

calculations take thousands of iterations to converge a few

orders of magnitude. Explicit pdf calculations require at

least ten thousand iterations for elliptic flows. The

current use of pdf methods must inherently include

simplification combustion modeling to give realistic

calculation times. Some of these approximations have been

already mentioned. If the number of scalars is kept low it

is possible to construct thermodynamic state tables. The

table is then interpolated to find thermodynamic data during

the iteration process. For instance, J. Y. Chen et al. 85

have constructed state tables for methane combustion using 4

and 5 scalar reaction mechanisms and two time increments.

Reaction submodels incorporated in the Lewis pdf model have

included infinite reaction rate, a two step global reaction

model for hydrogen by Chinitz, an Intrinsic Low-Dimensional

Manifold (ILDM) method of Maas and Pope 86. The ILDM method

uses dynamical systems approach to reduce the number of

degrees of freedom to simplify the chemical kinetics system

and construct a state look-up table. The method was
83

successful in the calculation of Norris and Hsu [78]. The

ILDM method for heavier, more complex fuels more difficult

and is currently being worked on. Also, the table look-up

method for higher molecular weight, more complex fuels

requires more reaction steps. This is beyond current

computer memory limits. A single finite-rate reaction

submodel is used in this work.

In 1990, Chen and Kollman 87 compared simplified

chemical reaction models for hydrogen-air combustion. At

conditions near equilibrium, the chemistry models all

produced similar results. Only one of the models was

reported to be able to predict flame blowout, but, the

blowout velocity was incorrect. They concluded that further

improvements in the turbulent mixing model and the numerical

scheme were required to obtain quantitatively correct

predictions. Also needed are data from simultaneous

measurements of radical species and mixture fraction.

4.2 - Further Development of Composite PDF Method

One contribution of this work is further development of

the composition pdf and testing the model for finite-rate

combustion in recirculating flows. The composition pdf has

been primarily used for parabolic or one-way flows with

modifications of equilibrium type combustion models.

Equilibrium or fast kinetics combustion models should be

used where the rate of mixing is much smaller than the

combustion rate. These models are valid for low molecular
84

weight fuels as hydrogen and methane at high temperatures.

Practical combustors use higher molecular weight fuels.

Burning is distributed in practical combustors. There is

significant overlap in reacting species, particularly at low

temperature.

The SIMPLE method of Patankar and its variants is the

most widely used solution method for practical combustor

flows. The composition pdf method is combined with the

SIMPLE method. Momentum, overall mass conservation,

turbulent kinetic energy and dissipation are solved using

the SIMPLE method. The composition pdf model used in this

work was initially developed by Hsu and co-workers at NASA

Lewis Research Center. The first release of this

composition pdf model was with the Lewis RPLUS 2-D combustor

code 8s. RPLUS is used to predict high speed combustor

flows. As Mach number is lowered the solution convergence

rate deteriorates. It is an extension of the Lower-Upper

Symmetric Solver used by Yoon and Jameson 89 to predict

species composition. The first dissemination of the

development model was at a workshop in 1993. The model

released at the workshop was set up to solve the flow

configuration of Burrows and Kurkov's supersonic hydrogen

combustor 9°. The model was set up to use a global five

species, two step reaction submodel originally proposed by

Rogers and Chinitz 91 for hydrogen-oxygen combustion.

Many changes were made to the original composite pdf
85

model. The implementation of the finite rate combustion

model didn't calculate steady state combustion using

suggested reaction numbers for methane. A different species

composition solver was implemented. The species

concentration solver uses the Newton-Raphson technique. In

the original coding, an incorrect implementation of upwind

differencing was used to calculate the number of particles

to exchange between cells. Pure upwinding is regarded as

highly diffusive in CFD. So as to fairly compare the

predictive capability of the combination pdf model against

the baseline combustion model, hybrid differencing was

implemented in the composition pdf for the calculation of

exchange coefficients. The thermodynamics of the base pdf

model were highly simplified. Constant Cp was assumed. Cp

changes markedly over the range of combustor temperatures.

Predictions of exit temperature were adjusted in development

model by adjusting Cp. The baseline combustion model used

in this work has variable Cp. Variable Cp was also

incorporated into the thermodynamics of the composition pdf

model. This necessitated the use of enthalpy rather than

temperature as one of the pdf variables. The development

combination pdf model used a combination of temperature

solver in Monte Carlo solver and enthalpy solver in RPLUS

code to solve for temperature. The model implemented here

properly solves for temperature completely in the pdf

solver.
86

Trial calculations are done using the SIMPLE method as

it is more economical than the hybrid pdf method. The best

solution is used to begin hybrid pdf solution and is also

used to compare predictions and measurements.

r! I i_
87

Chapter 5

COMBUSTORMODELVALIDATION

Experimental data used to evaluate modeling should be

selected based on the quality and quantity of experimental

data, especially initial or boundary condition data.

Experimental data should include measurement of multiple

chemical species, velocity, temperature and turbulence

quantities. If comparison is performed with limited data,

the quality of a CFD solution may not be adequately

evaluated. CFD solutions are usually repeated changing

initial conditions, gridding, turbulence modeling, etc.

until an optimum comparison is given. With limited data, an

optimum solution may be found for the wrong reason. A good

prediction of one species does not necessarily mean that all

other species are equally well predicted. Similarly, a good

prediction of one velocity profile does not mean that other

velocities profiles are equally as good. Good predictions

at one station in a combustor does not mean that all

stations are of the same quality.

Suitable and accurate data must be used to ascertain

the appropriateness and accuracy of numerical predictions.

Past measurement techniques have used hot wires and pressure

probes to determine velocities. Over the past decade,

optical measurements of velocity have become standard.

Measurements require some kind of access. Probes can be

handled by access from the rear of the combustor. Lasers
88

demand optical access through the combustor. This is

usually handled by installing heat resistant windows in the

combustor. The windows degrade faster with increasing

temperature, heat transfer and pressure. Thus, almost all

benchmark quality data for combustors have assembled using

research-type combustors. These combustors incorporate the

major features of practical commercial gas turbine

combustors, but, their designs are vastly simplified.

Running an actual combustor is very expensive and very

difficult to instrument. Most research type combustors are

axisymmetric, except for the possible addition of dilution

jets. The dilution jets are required to cool the combustion

liner and lower the overall equivalence ratio of the

combustor. The solution of combustors with dilution jets

requires 3-D models. Combustors flows undergo substantial

change in density. The flow through a combustor without the

flow expansion caused by combustion can be substantially

different than with combustion. Measurements of some cold

combustor flows have shown reverse flow at the exit of the

combustors. Most CFD models assume one-way flow at the

solution exit plane. One way exit flow may be enforced by

installing flow restrictions at the combustor exit.

5.1 - Combustor Model Validation Data

In 1985 an evaluation of data for parabolic turbulent

reacting flows was published. The editors were Strahle and

Lekoudis 92. The work was sponsored by the Air Force Office

!il 1i
89

of Scientific Research. The purpose was to conduct a

program similar to that run at Stanford in 19689_, except

on turbulent reacting flow as opposed to turbulent boundary

layers. A committee was set up to help in the evaluation of

the data. At the final committee meeting it was recommended

that a computational effort not be initiated at the time.

One reason was that most combustion theories are application

specific and cannot be readily used for flows of different

character or chemistry from those for which they were

developed. Some problems noted in the data were: lack of

completeness, incomplete specification of initial and

boundary conditions and excessive experimental uncertainty.

This work did not include data for complex reacting

turbulent recirculating flows.

One of the major problems with early combustion data is

the limited number of exPerimental quantities measured.

Many of the experiments were performed to characterize only

parametric variations of a combustor type, not to validate

CFD models. In a paper by Baker et al. 94, only velocities

and stresses are reported. In another paper by Lockwood et

al. 9s, only flame length and mixture fraction were

reported. McDannel, Peterson and Samuelson 96 gave species

and temperature measurements for a premixed flow. The flame

was stabilized by a recirculation zone formed by a high

velocity jet opposed to the direction of bulk flow.

Velocities were varied in a the study. E1 Banhawy,
90

Sivasegaram and Whitelaw 97 took measurements of

temperature, velocity and species concentration for a

premixed combustor in 1983. The measurements are presented

in contour form. In premixed combustors the fuel and

oxidant are combined upstream of the combustor inlet. In

this configuration, combustion is anchored by a hot

recirculation zone caused by a sudden area increase. This

geometry is known as a step combustor. Two different step

sizes were studied. Premixed combustors aren't used in

practical gas turbine engines, due to stability problems.

The concept has been under study for quite some time because

the concept reduces engine emissions. Smith, Giel and

Catalano 98measured hydrogen-air combustion in an

axisymmetric recirculating combustor flow. A lean mixture

was used. Due to its low molecular weight, hydrogen

diffuses much more rapidly than oxygen or water.

Preferential diffusivity wasn't included in the numerical

model, so this flow was not calculated. There is a wide

disparity in the inlet velocities, which causes the

recirculation zone. The hydrogen has a bulk velocity of

.914 m/s, while the air inlet velocity was 102 m/s.

Temperature and species concentration for the reacting case

was measured by laser raman spectroscopy. However, not

enough measurements were taken at each location for high

confidence levels. The temperature profiles are quite rough

for the first three measurement stations.

ii!Iil
91

Lightman et al. 99 performed measurements on a bluffbody combustor at Wright-Patterson Aero Propulsion

Laboratory. Measurements were for combusting and noncombusting flows. They reported that most of the features

of the cold flowfield also were present in the combusting

flow. However, in a later paper I°°, high speed pictures of

the combustor flow revealed the existence of large scale

structures. By examining laser Doppler anemometer

measurements and measurements of luminosity it was concluded

that the flow is comprised of non-luminous regions and flame

turbules. The combusting flow was unsteady. Schefer et

al. TM also saw coherent structures in their study of

another combustor flow stabilized by a bluff body. The flow

was characterized by large scale entrainment of air into the

recirculation zone. Pdf's of Velocity and velocity

correlation were measured. Bimodal pdf's were observed in

regions of high shear in the recirculating zone boundaries

and in the downstream stagnation region. This bimodality

was associated with the alternate passage of large-scale

turbulent structures consisting of unmixed fuel and air

through the measurement volume. In 1987, Sivasegaram and

Whitelaw I°2 published a paper studying oscillations in

axisymmetric dump combustors. Another study was published

in 1989 by $ivasegaram et al. I°3 Some unsteady combustion

predictions have been done, but, the practicality of the

work is somewhat controversial.
92

Samuelson and his students have produced a number of

papers regarding the design and experimental measurements of

axisymmetric combustors. The most complete 2-D database is

an University of California , Irvine, Combustion Laboratory

paper by R. Charles I°4. A parametric series of

measurements were made by modifying the combustor inlet

conditions including the fuel nozzle. Unfortunately, only

axial and tangential velocity measurements are reported.

The radial velocity is unknown, although it may be inferred

from continuity for non-combusting flows. In general, the

number measurements taken in the radial direction number

around ten. This number is too small to accurately validate

or benchmark combustor models. Axial velocities measured 1

cm from the inlet show considerable backflow. Azimuthal

velocity generally shows a peak near the outside wall at

this axial location. Peak azimuthal velocity at other axial

locations is much closer to the combustor centerline. Inlet

data from this configuration would have to be set up based

on flow geometry, bulk flow and experience of the combustion

modeler.

J.C. Pan 1°s did a parametric study of turbulent

confined premixed flames. The flame was stabilized by

conical bluff bodies. In the study, the size and angle of

the c0nical bluff bodies, the equivalence ratio, the inlet

velocity and turbulence were varied. The data given in the

reference were not in tabular form. Velocity and

Ill I!
93

temperature pdf's were given, velocity vector diagrams,

contours of velocity and reynolds stresses and centerline

profiles of velocity and temperature. Some cases used very

lean mixtures close to the blow-out limit. The data may be

good for validating flamelet combustion models.

Some three dimensional combustor flows have also been

measured. Heitor and Whitelaw I°6 measured velocity,

temperature and species in a combustor having dilution jet

holes. The combustor exit is constricted to one side. The

fuel was gaseous propane. The scalar field was dominated by

the presence of the dilution jets. Air to fuel ratio was

varied. Combustion efficiency was low for some of the runs.

In a 1990 paper, Bicen, Tse, and Whitelaw I°_ reported up to

98 % combustion efficiency for a similar can-type combustor.

The primary zone of the combustor was fuel rich.
94

5.2 - Spray Model Validation

The baseline combustion model was tested in various

stages to benchmark its predictive capability. The first

stage was predicting an isothermal open combustor flow. The

spray vaporization model was then tested. Finally,

calculations of spray combustion were performed.

5.2.1 Isothermal Swirling Flow

The isothermal swirling combustor flow of D. Bulzan I°8

was used to validate the baseline combustor model. This

combustor is part of an on-going effort to provide experimental data to validate spray combustion models. The experimental configuration is a Parker Hanifin research simplex

air-assist atomizer surrounded by a co-flowing air stream.

The configuration used here is unconfined, that is, it does

not have an outer wall. This was done to simplify taking

experimental measurements. A schematic drawing of the

combustor is shown in fig. 5.1. Most unconfined or free-jet

combustors are notorious for poor mixing, and require significant axial distances for complete mixing. The atomizer

and the co-flow air had 45 degree swirl in the same direction. This swirl causes the formation of a large recirculation zone which promotes rapid mixing and stabilizes the

flame for the reacting case. The recirculation bubble is

longer than the laser measurement system can traverse. The

model was able to handle this negative outflow condition by

suitable treatment of the exit boundary conditions as sug1_ :11I
95

gested by P. George

Huang. Normally, exit

boundary conditions

assume zero-gradient or

one-way flow where the

flow is not significantly changing with

axial distance. The

exit velocity is unknown, but, is needed

in the calculation of

Swirlers

pressure. The solution

is part of the iterarive process. The exit

velocity is set to the

last calculated upstream velocity plus a

Figure 5.1 Schematic drawing of the

combustor. Dimensions in mm.

velocity correction due to mass conservation. Thus pressure

for the exit computation cell can be calculated. In the

process of calculating the upstream velocity, the downstream

contribution to the exit computational cell is set to zero

without adding any significant error, since the downstream

contribution is either zero or very small. For the case

where the exit velocity is negative for a computational

cell, the exit neighbor coefficient is the upstream coefficient and the contribution is very significant and cannot be
96

ignored. This is actually an unspecified inlet condition,

which must be calculated. This is handled in the iteration

process by setting the boundary cell velocity equal to the

velocity of the downstream adjacent cell, which is calculated. Then the upstream or boundary neighbor coefficient is

added to the non-constant part of the source term. Then the

neighbor coefficient is set to zero (the boundary coefficient isn't included in the sum of neighbor coefficients).

To summarize, first neighbor coefficients are calculated in

the normal way. Then the (negative) exit velocity coefficient is added to the non-constant part of the source term.

Then the exit velocity is treated as if it were zero in the

difference equation. Pressure is still calculated the same

as before. The mass correction at the exit boundary is only

made where the exit velocity is positive.

The flowfield in the inlet region is quite complicated.

The combination of the blunt body atomizer with the highly

swirling co-annular flow produces two recirculation zones.

The first has already been mentioned. The second is a small

toroidal vortex adjacent to the blunt atomizer between the

outward swirling flow and the reverse flow along the centerline. The unconfined combustor entrains surrounding ambient

air. A velocity-vector diagram from a calculation is shown

in figure 5.2. The centerline is on the right side of the

figure. The flow is from bottom to top. The small vortex

near the bottom center-line is somewhat larger than what was

IF!I-
97

experimentally measured. Calculation of the flow-field is

difficult. One problem is reverse flow at the first measurement location. The current model can handle moderate

amounts of reverse flow at the combustor inlet, as it uses a
98

pressure based algorithm. Another problem is that the initial flow is rapidly expanding into an open environment.

Entraining flow is difficult to accurately predict using

pressure based algorithms. The strong swirling flow, which

produces the very large recirculation bubble, could not be

contained by boundary conditions intended for entraining

free-jet flow. Using wall-type boundary conditions did

contain the flow, but, the size of the calculated recirculation zone depended on placement of the non-physical wall.

Using a larger radial calculation domain produced a larger

recirculation bubble, with most of the positive axial flow

occurring close to the wall. The experimental measurements

show an axial velocity peak which gradually moves outward.

A solution to this was to specify velocities on the outer

radial boundary. Dr. Bulzan measured some velocities at

radial distances of 143 mm and 168 mm from the combustor

centerline. This data was obtained directly from Dr. Bulzan. The measurements at 143 mm were used for this calculation. At this radius, flow was entrained to an axial distance of about 130 _. Measured values of velocity at 5 mm

downstream of the atomizer were used as initial conditions

for the calculation. Predicted flow fields varied greatly

by using different initial values of turbulent kinetic

energy and dissipation. Normally, the k-£ turbulence model

is believed to produce only small change due to rather large

change in initial turbulence quantities. Typically, the use

!!l11
99

of poor initial turbulence quantities is overcome within a

few gridpoints. Thus, the k-_ model is perceived as a good

model to use as predictions are not unusually deDendent on

initial dissipation or length scale, which are usually

unknown. Unfortunately, the initial structure of the flowfield is quite complicated, is unconfined, and involves very

high shear. The atomizer causes large velocity gradients.

This causes a huge generation of kinetic energy and rapid

dissipation. Using the full measured turbulence energy and

a length scale of 30 % of the diameter of the atomizer and

co-flowing air passage smeared out the small secondary

vortex and blocked the reverse flow aiong the centerline.

Also, the predicted central recirculation zone was over a

factor of two wider than measured at the calculation exit.

Another calculation was done using peak turbulent

kinetic energy around an order of magnitude smaller than

measured and specifying a length scale of 30 % of the width

of the total calculation domain. This resulted in very good

predictions of the velocity field. Evidently, the low initial values allowed lower initial dissipation and which

later built up enough to predict the correct amount of

spreading.

Due to convergence problems, initial swirl was slowly

ramped up to the full value. The total number of iterations

was 6000. Calculated versus measured axial velocities are

shown in figure 5.3. The comparisons are done at axial
I00

distances of 20, 50, and i00 mm. The axial velocity predictions show quite good agreement with the measurements. The

velocity peaks are slightly underpredicted. Figure 5.4

shows comparisons of radial velocity. Radial velocities are

underpredicted.

Figure 5.5 shows comparisons of tangential velocity.

Similar to axial velocity profiles, the calculated tangential velocity profiles are very good. Peak velocities are

only slightly underpredicted. Plots of turbulent kinetic

energy are shown in figure 5.6. The comparison is poor for

the lower half of the profile at 20 mm, where the turbulent

kinetic energy is vastly underpredicted. This part of the

flow represents shear between the second vortex and the

backflow along the centerline from the large recirculation

zone. The local axial velocity peak is also underpredicted

in this shear region.

Predicted profiles at other locations are very good.

These calculations show that juggling initial turbulence

values has a significant effect on computational results for

this flow. As reducing initial turbulence values reduced

mixing and produced good computational results, a modification to the turbulence model was considered. The k-£ model

assumes isotropic turbulence and is not regarded as optimal

where swirl, body forces, or curvature exist. In the past,

some researchers have resorted to a curvature corrections to

improve predictions for strongly swirling flows I°9.

NI11
I01

0.125

0.100

E 0.075

v

¢n

a

<

tr 0,050

0.025

0.000

[1011experiment omm O

20 mm

100 mm

0

0

0

°

0

-5 0 5 5 5

AXIAL VELOCITY (m/s)

Figure 5.3 Calculated and experimental

axial velocity at 20, 50, and i00 mm.

0.125 "1

0.100

E 0.075

o9

tn

,<

tr 0.050

0.025

0.000 5

-5

I

calculated I 50 mm

experiment [ O

20 mm )

)0

0

0

0

0

0

0

0

0

I

0 5

I00 mm

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

D

9

5

RADIALVELOCITY(m/s)

Figure 5.4 Comparisons of radial velocity

at 20, 50, and I00 mm.
0.125

0.100

E 0.075

09

D

5

<

rr 0.050

0.025

10texperiment I,Ommo 0

o

D

0.000 i i

-5 5

20 mm

O

0 °

-- | ,

0 5

Figure 5.5

TANGENTIALVELO_TY(m/s)

Comparisons of

velocity at 20, 50 and i00 mm.

100 mm

O

!

)

)

(

5

tangential

102

0.125

0.100

E 0.075

09

a

<

n- 0.050

0.025

0.000

-5

IT calculated experiment

20 mm

• 0 0

O0 0

I

0 5

50 mm

O

O

O

O

%

(

C

C

OCD O

_ o

C

I

5

100 mm

O

O

I

(

C

%0 ___,

5

TURBULENTKINETICENERGY

Figure 5.6 Calculated and

turbulent kinetic energy.

experimental

!E!! '::
103

Several modifications were tested by Sloan et. al. n° for

thirteen cases. Unfortunately, the lack of initial conditions prevented drawing absolute conclusions about the

improvements. This work does not employ curvature or Richardson number corrections. Instead, a simpler method of

changing the turbulent viscosity coefficient based on shear

was tried. Recently, a modification for high shear was

proposed and tested by Shih et. al. nl This development

is based on invariant theory in continuum mechanics. The

modification is done so that negative normal stresses are

not calculated in any situation of rapid distortion. This

involves damping of the coefficient for turbulent viscosity,

which is normally treated as a constant.

are:

k2

vt= C_

The equations used

(5.1)

2/3 (5.2)

c_--_÷_

= Sk/z (5.3)

s =q2si_sij (5.4)

si_ = (u_,j÷u_._)/2 (5.5)

The predicted viscosity coefficient ranges from a high of

0._21 to a low of 0.0178 at the edge of the co-flowing air.
104

The results of this calculation are shown in figures 5.7-10.

The axial velocity profile at 20 mm is similar to the

results with the regular k-£ model, but, the 50 and i00 mm

profiles do not show as much radial spreading. The radial

velocity profiles are slightly reduced which accounts for

the lower spreading rate. This is also reflected in the

tangential velocity profiles. The peak turbulent kinetic

energy for the 20 mm axial location is improved, but, at all

other locations predictions are worse with much less radial

spreading. It appears that modification of the turbulence

can be made to improve predictions in one area, but, unfortunately reduce accuracy in other areas.

These calculations were done with hybrid differencing

which is known to be overly diffusive, especially for calculation of round free jets. Higher order differencing should

be used before judgements about the validity of various turbulence modeling improvements.

iilI !_
105

0.125

0.100

E 0.075

o9

D

tn

<

n" 0.050

0.025

0.000

experiment o

100 mm

O

O

o

20 mm

O

%

0

!

-5 0 5 5 5

AXIAL VELOCITY (m/s)

Figure 5.7 Axial velocity profiles at 20,

50, and i00 mm with modified turbulence

model

0.125

0.100

E 0.075

09

rn

<

rr 0.050

0.025

0.000

)

9

O

£

-5 0 5

50 mm

O

)

D

O

O

O

O

O

O

,g

0

0

O

3

I

5

I00 mm

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

O

9

9

5

RADIALVELOCITY(_s)

Figure 5.8 Radial velocity profiles at 20,

50, and I00 mm.
106

0.125

--]calculated I 50mm lOOmm 0 0

0 experiment 0 I0

O.100 20 mm 0

0.075

o'3

_D o

<

n- 0.050 C

o° _

o

o

o

0.025 0

0.000 -- ' _ _0 J

-5 0 5 5 5

TANGENTIAL VELOCITY (m/s)

Figure 5.9 Tangential velocity profiles at

20, 50, and i00 mm with modified turbulence model

0.125

0.100

E 0.075

D

<

cc 0.050

0.025

=l ou. experiment Iserum O

O

O

0.000 -----J

-5 5

100 mm

0

0

0

_0_m t° °o

o

o

0 o

'_ J O J

0 5 5

TURBULENT KINETIC ENERGY

Figure 5.10 Turbulent kinetic energy at

20, 50 and i00 mmwith modified turbulence

model

1ii 1:_
107

5.2.2 Evaporating Spray Model Validation

The methanol vaporization data of McDonell et. al. n2

was used to validate the spray model. McDonell's data were

specifically taken to validate spray models. The experimental apparatus consists of a research atomizer injecting

methanol spray downward into a 457 mm square duct with coflowing air stream at an approximate bulk velocity of 1.0

m/s. Measurements were made of axial and radial gas-phase

and droplet velocities, methanol vapor concentration, and

droplet fluxes. No reverse flow was measured in this noncombusting case. The droplet data were measured by Phase

Doppler Interferometry (PDI). The spray droplet data were

partitioned into i0 different size groups due to limitations

of the instrument. Data at an axial distance of 7.5 mm were

used as initial conditions for the calculation. Gas phase

velocities were measured to a radial distance of 88.0 mm.

Initial gas phase velocities beyond this distance were

ramped to the bulk velocity of 1.0 m/s. Spray droplet data

were specified at radial distances out to 13.0 mm in 1 mm

increments. The total number of droplet groups measured at

this axial location was 102 as the large droplets were not

found near the center line.

According to McDonell, the error in vapor concentration measurement was about i0 %. The error in droplet

measurements is higher, especially in the initially dense
108

spray. The sum of the measured methanol vapor and droplet

fluxes only accounted for 23 % of the known methanol flow at

the first axial measurement location. The fraction of the

measured methanol flux goes up to 85 % at the last measurement location. The PDI measures the number of fringes

crossed and does not accurately measure flows in very dense

sprays. Droplets are going through the measurement zone at

various angles, velocities, shapes, and sizes. The measurements are invalidated if there are multiple droplets detected in the measurement zone. As a rough approximation, the

corrected droplet count rate is usually taken as the sampling rate. The error in measuring small droplets in dense

sprays is larger as there exists a higher probability that

multiple droplets are in the measurement zone. The ratio of

the sample rate to the validated rate is around 5 for the

smallest droplets at the first axial measurement location.

This ratio drops to i.I for the largest droplets in the

outer less dense region of the spray.

Four sets of spray droplet input data were used in the

predictions. For the first calculation the droplet data

rate was incremented at each radial location so that the McDonell's flux rate at that radial location was achieved.

Then, an overall final correction was done for the flux

rates so that the total methanol flow rate was correct.

This droplet rate correction seemed to bias the numbers

towards the larger droplets. Larger droplets evaporate more

IrlI
109

slowly than smaller droplets as the surface to volume ratio

is smaller for the large droplets. The calculated methanol

vapor flux was lower than that measured by McDonell for the

last two stations, as shown in figure 5.11. In the second

prediction, the sampling rate was simply multiplied by an

overall value to give the total methanol flux. This slightly increased the vaporization rate. For the third prediction, the actual probe area values were obtained from McDonell. These probe volumes changed from location to location. This information allowed for a more accurate measurement of droplet flux. The calculated vaporization was similar to that obtained in previous predictions. In the fourth

prediction, the overall mass correction was weighed for each

droplet size according to the ratio of the sampling rate to

the validated measurement rate. This calculation had the

largest number of small droplets and predicted the highest

vaporization rate. It was hoped in using different initial

spray distributions that the velocity field predictions

would improve. As shown by figures 5.12-15 this did not

happen. The use of different droplet data did not drastically change the axial velocity profiles. Thus, the error

in predictions does not appear to be caused by the initial

spray distribution, or the spray model. Some of the error

may be due to the initial radial gas phase velocities.

McDonell's data for these velocities is sparse. Near the

atomizer, measurements were taken only very 4 mm. In a
ii0

0.050

o

"_ 0.040

C7)

rr

0

< 0.030

>

,_1

LU

D

LL

LL

0 0.020

0

_J

09 0.010

<

droplet dist. 1

" " " 2

__ . -- I! iv i! 3

" " " 4

....... low v-velocity

0 experiment

/

/

/

/

/

/

/

/

/0

/

/

/

/

/

0.000 l i = I i

0 20 40 60 80 100

AXIAL DISTANCE (mm)

Figure 5.11 Mass flow rate of vaporized Methanol

IF!11_
111

5

Or)

E

4

3

i

0

0

._l

LU

>

..n 2

<

X

<

0

i

• droplet dist. 1

0 ...... 3

_ ....... lowv-velooity

0 experiment

0 0

0

I I I

0 20 40 60

RADIUS (ram)

Figure 5.12 Axial Velocity at 15 mm

I l

8O 100
112

5

4

E

V

>" 3

0

_J

W

>

2

<

X

<

0

droplet dist. 1

.... N m g 2

-- -- • " " " 4

....... low v-velocity

0 experiment

1

0_

O0

0 i i i

0 20 40 60 80 100

RADIUS (ram)

Figure 5.13 Axial Velocity at 25 mm

I I

IIl 1i
113

5

4

_9

O

o

ILl

>

_j 2

<

m

X

<

droplet dist. 1

_ _ _. w w w 2

--'-- " " " 3

....... low v-velocity

O experiment

0 I 1 I

0 20 40 60

RADIUS (ram)

Figure 5.14 Axial velocity at 50 mm

I

8O

J

100
114

5

4

droplet dist. 1

.... w w I! 2

- -- - -- " " " 3

_ . w w M 4

....... low v-velocity

0 experiment

o

_ 3 )"'" o

9

I.

0 I I [ I

0 20 40 60 80

RADIUS (mm)

Figure 5.15 Axial velocity at i00 mm

I

IO0
115

distance of 12 mm, the reported radial velocity increases to

over 3.2 m/s and then decreases to -1.547 m/s. The turbulence is extremely high. Axial velocities vary as much as

5% between the axial-radial and the axial-azimuthal velocity

sweep measurements. In the predictions, the initial high

radial velocity components caused a rapid decay in radial

and axial velocity components. To check the radial velocity

components, another calculation was done with much smaller

radial velocity component near the atomizer. This produced

axial velocity profiles with higher peaks. Also, the peaks

occurred along the centerline, as is shown by most of the

measured velocity profiles. With the reduction in radial

velocity, the predicted velocity profile almost perfectly

matches the measured velocity profile at 50 mm, but, at i00

mm the predicted profile is flatter. The predictions appear

to be suffering from too much diffusion. Various curvature

effects can be applied to coefficients in the turbulence

model, but, with the rather small amount of curvature in

this flow, these corrections would do relatively little to

improve predictions.

Predicted methanol vapor concentration is shown in

figure 5.16 for an axial distance of 15 mm. The contours

are remarkably similar for all of the predictions at this

axial location. Experimental measurements were not given at

this location. A comparison of experimental and predicted

contours at 25, 50 and i00 mm is given in figures 5.17-19.
116

V

3.0

2.5

Z

O 2.0

!--

<

l--

z

LU 1.5

O

Z

O

O

co 1.0

09

<

0.5

0.0

f._L _ droplet dist. 1

.... m u m 2

_>.>N. ...... 3

"%,-,_- ...... 4

"%', \

",_,,, \ ....... low v-velocity

_

_ "._.,_

-

1 I I I !

0 5 10 15 20 25 30

RADIUS (mm)

Figure 5.16 Methanol concentrations at 15 mm

l!l! i_
117

4.0 (

3.5

n"

I--

z

LU 2.0

O

Z

O

0 1.5

or)

09

< 1.0

0.5

0.0

O

O

O

O

\O

- ",..._% \

6-_ \

o .'.%

o "-_

- 0 "'_

- Oo, \

_ %'_,

, , _---_

0 10 20 30

droplet dist. 1

.... " " " 2

B m N 3

-- -- - - • " 4

....... low v-velocity

O experiment

I I

40 50

RADIUS (mm)

Figure 5.17 Methanol concentration at 25 mm
118

0.0

\

\

\

\

\

0 \

0 \

o \

0

I I 1

droplet dist. 1

.... m w m 2

-- . -- • m • 3

__ . s m w 4

....... low v-velocity

O experiment

0 10 20 30

RADIUS (mm)

4O

Figure 5.18 Methanol concentration at 50 mm

50

iil IT
119

,.0I

o_.° 3.0 \\

g

- ,,,_

2.5 ,

uJ 2,0

O

Z

O

O 1.5 -

09

or)

"< 1.0 -

0.5

0.0

0

droplet dist. 1

I l m 2

I N W 3

" - " 4

low v-velocity

experiment

",I

0 2O 40 6O 80

RADIUS (mm)

Figure 5.19 Methanol concentration at i00 mm

I

100
120

All of the predicted profiles are flatter and wider than the

measured profiles, with the highest peak shown by droplet

distribution number 4 which has the largest number of small

droplets. The predicted methanol contours are diffusive,

suggesting a deficiency in the numerical differencing model.

Correcting this deficiency is beyond the scope of this work.

The outer portions of all of the predicted contours are

similar. McDonell states that some of the methanol concentration measurements are at a saturation state along the

center-line. Predicted methanol concentration initially

drops along the centerline. The maximum calculated concentration is never as high as the peak experimental measurements. This suggests that the spray model may be somewhat

deficient in regards to the saturation concentration. An

off-center peak is shown by all of the predicted contours at

i00 mm while the measurements show off-center peaks at 50

and 100 mm. In the predictions, an axial location is eventually predicted in the flowfield where all droplet groups

are continually spreading outward. Even if a droplet group

is started with a negative radial velocity, the calculated

droplet trajectory will eventually cross the centerline, be

reflected in the calculation and become positive. Thus, the

calculation should predict an off-center methanol concentration peak as the evaporating methanol droplets move further

from the center-line. Some predicted and measured droplet axial velocities at axial distances of 15, 25, 50, and

I[I!
121

75 mm are shown in figures 5.20-23, respectively, for the

third prediction (probe area taken into account, no individual weighing of error measurement). The data for the other

predictions should be almost identical, as the largest

difference should be the gas-phase velocities which was

shown to be almost identical for all the changes in spray

input. The predicted velocities are shown as curves while

the measurements are shown as symbols. The predictions were

done with initial droplet sizes of 7.4, 15.5, 25.5, 35.5,

45.5, 55.5, 68.0, 83.0, 98.0, and 120.0 _m in diameter to

correspond to the ten size ranges used in the measurements.

The profiles at 15 mm are very good and remarkably smooth.

This tends to verify the initial conditions used for the

calculation. The predictions for the large droplets at the

next axial measurement location are also very good. The

prediction for the 31 to 40 _m droplet groups show some

oscillation. Some of the droplet groups are crossing each

other. One group with lower axial velocity is crossing the

path of a droplet group with higher axial velocity, resulting in a discontinuous profile when plotting the velocities

of the droplet groups. The further downstream, the more

opportunity for the predicted droplet groups to cross paths.

However, all of the droplet groups do approach the gas phase

velocity as a function of their size and speed. This tends

to smooth the velocity profiles, canceling the first effect.

Also, this effect should not have much of an effect on the
122

E

v

>-

i

O

O

W

>

_J

<i

X

<

20

16

12

8

4

0

droplet size (lim)

gas phase pred.

.... 11.0-20.0 "

31.0-40.0 "

----51.0-60.0 "

O

J

0 20 40 60 80 100

RADIUS (mm)

Figure 5.20 Droplet and gas-phase axial velocity at 15 mm

ill| i:
123

E

>-

F-u

O

O

LU

>

.J

<

X

<

15

12

9

6

3

droplet size (_m)

gas phase pred.

.... 11,0-20.0 "

b._ 31.0-40.0 "

b-- 51.0-60.0 "

....... 76.0-90.0 "

0 gas phase data

[] 11.0-20.0 "

A 31.0-40.0 "

_'. _ sl.o-6o.o"

- _,/'Z_/k\_A <_><_", V 76.0-90.0 "

,/

_ It,x\

o, I!1

_0_ 0

0 I I I

0 20 40 60 80 1O0

RADIUS (mm)

Figure 5.21 Droplet and gas-phase axial velocity at 25 mm
124

E

v

>-

O

O

_J

w

>

<

X

<

10

8

6

4

2

0

V ,,

_7 _7 _',,,

_7",

v

- _\\ ,,

<><>\\

- °Ok> []

o .. 0\%

l I I l

0 20 40 60 80

RADIUS (mm)

O

Figure 5.22 Droplet and gas-phase velocity at 50 mm

I

100

lil ]i:
125

E

v

>-

0

0

W

>

6

5

4

3

<

2

<

O0 \

\

\

\

\

\

q

0

0

I I I

0

0 20 40 60 80 100

RADIUS (mm)

Figure 5.23 Droplet and gas-phase axial velocity at 75 mm
126

concentration prediction, which is more important than the

prediction of each droplet group's velocity. At an axial

distance of 50 mm the predicted velocities for the two

smallest groups shown are nearly the same as the gas-phase

velocities while the measured profile shows the same behavior at less than 28 mm radius but shows a second outer peak

with nearly the same velocities as measured for the 31-40

and 51-60 SMD droplet group measurements. In checking the

predictions more closely, it was found that no droplets in

this size range were found beyond a radius of 18.5 mm at

this axial location. The initial measurements at 7.5 mm do

not show a secondary peak. There is no gas-phase fluid

dynamic structure or effect that would cause this. The

only possibility is droplets from another droplet group. In

the predictions, only the smallest droplets significantly

evaporate, while the larger droplets only lose a small

fraction of their size. The smallest droplets eventually are

treated as totally evaporated when they reach a certain

size. This is done because the fixed droplet time step

could produce non-physical results. The use of smaller

droplet time-steps would allow tracking the droplets longer

but this would increase the number of computations and

memory needed dramatically. The predictions do not show any

droplets losing enough mass to be included in a lower size

droplet group within the first 100 mm. The predictions show

that some of the original 15.5 _m droplets are reduced to

I!I1!
127

10.5 _m in size when they pass the last measurement location

of i00 mm. Larger diameter droplet groups are reduced in

diameter even less. The largest droplets only lose about 2

_un in diameter in i00 mm. The predictions are started using

i0 different droplet sizes. Each inlet droplet represents

an average of a very large number of droplets spanning at

least 10 Bm in diameter. Conceptually, a group of droplets

near the smaller end of a droplet bin size range would drop

into the next size measurement bin at a downstream location.

The secondary peaks measured later in the flowfield are most

likely part of a group of droplets having larger diameter

and velocity near the injector. This effect could be duplicated in a calculation by increasing the number of initial

sizes. The smaller size would evaporate enough to drop into

the next size range during the simulation, providing data

for the smaller-diameter group at that location. However,

the prediction of these velocities is not a requirement for

accurate calculation of methanol concentration. The predicted droplet is still there, producing vapor, and other

effects. It's just that the spray model works with a discrete size and the measurements deal with a range of sizes.

The production of methanol vapor is largely unaffected if

additional droplet groups are added with the appropriate

averaging for the diameters. More droplet groups would help

along the calculation axis of symmetry where calculated

droplet trajectories eventually all diverge. The spray
128

model interpolates between droplet sources for individual

grid cell source terms. With no or few source terms near

the axis of symmetry, the predicted methanol concentration

can be depressed.

There is another difference between predictions and

measurements. The predictions show large droplets up to and

beyond a radial distance of 60 mm at an axial distance of 50

mm. The measurements were stopped at an radial distance of

48 mm. The measurements at this location show a fairly high

droplet rate. From these numbers, it may be deduced that

there are many droplets beyond 48 mm. This same defect is

also seen at an axial distance of 75 mm, where the calculation shows droplets up to 80 mm radius and the measurements

are cut off at 60 mm. This defect shows up also in bar

graphs of the droplet flux.

In the measurements, the number of droplets crossing

the probe area are counted. A radial profile of these

numbers is somewhat Gaussian. In the predictions, the

droplet flux must be converted into a number of droplets

entering the calculation domain in a certain time-step. The

number of droplets crossing an area must be converted into a

set of droplet groups with specified numbers, sizes and

velocities. The Lagrangian predicted radial location of

droplet groups changes as the groups move downstream. The

predicted radial location does not vary uniformly as does

the experimental measurements which are taken at specific

FI1-i-
129

locations. Also, in the predictions, the number of droplets

in a group does not change, unless that group evaporate s (is

taken out of the calculation). Thus, the measurements of

density and the predictions of droplet location cannot be

directly compared. In this case, the predicted data will be

processed to resemble a normalized droplet density to compare to the measured droplet fluxes. This was done by

dividing by a pseudo-area for the droplet group and then

non-dimensionalizing this by the total number of droplets in

that size range. A minimum area was specified so that a

unrealistically large number would be calculated at small

radii. Also, if more than one group was calculated at

nearly the same radius, the groups were sometimes combined

in the plots. The measured data are also normalized. That

is, the sum of all normalized droplet densities adds up to

1.0. Accordingly, the data obtained are a somewhat graphical representation of the calculated position of the droplet

groups and the fraction in each group. The droplet distributions at 15, 25, and 75 mm are shown in figures 5.24-31.

Both the measured and predicted droplet densities are somewhat Gaussian and show spreading. The numerical predictions

generally keep the same shape through the calculation until

droplet trajectories cross each other. At this point, the

predicted bar graphs become somewhat ambiguous. Droplet

group fractions should probably be combined in some of

plots. In the predictions, the effect of an increase in
0.35

0.30

>-

F0_

Z

UJ 0.25

a

FLU

.J

n 0.20

O

n"

E3

a 0.15

LU

_N

.J

<

0.10

nO

Z

0.05

I

t

I

1

1

L

[--- -_

I

I--I

0.00

0 10

axial distance

-- x=15mm

.... x=25mm

--'-- x=75 mm

,--'-I

-, I

- i

L.--,

I

i

I

'--'-t

_ L . --, "l . __ .

20 30 40 50

RADIUS (mm)

Figure 5.24 Measured droplet density data

for Ii - 20 micro-meter droplets

0.35

130

0.30

z

uJ 0.25

a

II.U

-.I

n 0.20

O

ha

a 0.15

LU

_N

-J

<

0.10

nO

z

0.05

axial distance

x= 15 mm

.... x= 25mm

x=75 turn

]1, I '

ii 1;

i:i

I;i

ii,

III

0.00 _ ' I I r I I

0 10 20 30 40 50

RADIUS (ram)

Figure 5.25 Calculated Droplet density for

10.5 - 20.5 micro-meter droplets

I[I :I];
131

0.5

0.4

(/)

z

III

a

w 0.3 .J

Q.

O

rr

a

tn

w 0.2

U

J

<

rr

Oz 0.1

0.0

m

i i

--7i I

!

- -- "'J I

I

axial distance

-- x= 15ram

--'-- x=75 mm

L___.

I

I

I

I

L-- I

0 10 20 30 4O

RADIUS (mm)

Figure 5.26 Measured droplet density for

31 - 40 micro-meter droplets

0.5

0.4

z

uJ

a

i-.-

uJ 0.3 _,J

o,.

0

no

o

uJ 0.2

N

,_1

<[

nO 0.1 z

0.0

axial distance

-- X = 15.'_

.... X = 25 rrtm

_'-- X=75 mm

-t

I

I

I.

t i JI

5;i---JS

0 10 20 30 40

RADIUS (mm)

I I

Figure 5.27 Predicted droplet density

for 30.5 - 40.5 micro-meter droplets
0.5

axial distance

>-

I- 0.4

O9

Z

UJ

a

FLU 0.3 ,--I

O.

O

ha

LU 0.2

N

-J

<

nO 0.1 Z

L--

I--iI

L, IlL

0.0 r _.'-- -,-_P -;z I

0 10 20 30 40

RADIUS (mm)

-- x=lSmm

.... x=25mm

--'-- x=75 mm

I

I

I

I

I

I

[°_

I -'-P50 60

Figure 5.28 Measured droplet density for

51 - 60 _m droplets

0.5

0.4

09

z

w

a

0.3

-J

13.

O

a

a

LU 0.2

N

_.1

CE

0 0.1 Z

axial distance

x=15mm

.... x--25 mm

_'-- x=75 mm

r,

'1

I,

'1

--- li "i i -I ii!"

i ii!'

i i i,_i

-._ I I

!

o I

I

i

_.

-

0 10

o.o -_ _-'--'-' I t I I I

20 30 40 50 60

RADIUS (ram)

Figure 5.29 Predicted droplet density for

50.5 - 60.5 Bm droplets

132
133

0.5

>-

P-- 0.4

O9

z

uJ

a

FLU 0.3 ..J

13_

O

n"

a

t'n

w 0.2

N_

._1

<

nO 0.1 z

axial distance

x=15mm

.... x=25mm

x=75mm

•-'11 i '--'

I -! _"

, I

I 1

i

I it... I

-i I

7 I-r" i i

,3 -_4 I l ,. l

7O

[I I

0.0

0 10 20 30 40 50 60 80

RADIUS (mrn)

Figure 5.30 Measured droplet density for

76 - 90 _m droplets

0.5

0.4

Z

W

r_

lt,U 0.3 ...J

ft..

O

rr

D

a

w 0.2

N

,<

nO o.1

z

o.o

axial distance

x= 15ram

.... x=_Smm

--'-- x=75 mm

i_ I

I

I -

I

I

I

t

I

I

I

I

i

I

_ I

.... i

I

I

_._ -i I

,i

I-I_ - I

20 40

RADIUS (mm)

i

I

i

I

I

l._

I

c. 71 I

60 80

Figure 5.31 Predicted droplet density for

75.5 - 90.5 _m droplets
134

droplet density is seen as a decrease in distance between

droplet groups, not as an increase in the number of droplets. An increase in droplet density is intuitively contradictory for this flowfield, and numerical predictions of

such effects most likely represent slightly incorrect initial velocities. However, the somewhat inaccurate calculation of numerous droplet densities does not have a large

effect on the prediction of methanol vapor concentration.

It would be computationally impossible to model each droplet. The numerical predictions do not show as much spreading as do the measured droplet densities, except for the

largest diameter droplets. The spray prediction model does

not take into account interaction between turbulent eddies

and droplet groups. Droplet turbulent dispersion models

take this effect into account. In this model, droplet

interactions with the turbulent eddies are modeled by varying the instantaneous gas phase velocity surrounding the

droplet by using statistical fractions of the turbulent

fluctuation velocity calculated from local turbulence quantities. This would tend to increase the droplet spreading.

However, to do this the number of droplet groups must be increased by at an order of magnitude in order to maintain

meaningful results. This would increase calculation time

also by an order of magnitude. This is impractical for this

work as testing turbulent combustion models will also substantially increase calculation time.

IF!I'_
135

The largest diameter droplet density graphs show similar spreading for both experimental and numerical predictions. The larger droplets are affected less by gas-phase

drag, having a lower surface to volume ratio. Thus, the

larger droplets tend to spread further. The full extent of

large droplet spreading was not found since the measurements

were terminated too early.

The current droplet spray model produces good results

for methanol concentration. Droplet velocities can be surprising well predicted. There is an error in the amount of

droplet spreading. However, there isn't a large enough

variance from experimental measurements to justify a turbulent dispersion spray model, with its substantially increased computational costs.
136

5.2.3 - Spray Combustion Calculations

Both spray modeling and spray measurement capabilities

have markedly improved in the last 15 years. In a 1980

paper, E1 Banhawy and Whitelaw II_ assumed two spray distributions. Also, the droplet vaporization model used in the

calculations did not take into account the effect of combustion. The capability now exists to measure droplet spray

distributions for size, number flux, and velocity. Spray

combustion data is being assembled from simplified researchtype combustors. Gaseous data is also needed to validate or

benchmark spray combustion models, as the gas phase must be

accurately predicted in conjunction with the liquid phase.

Thus, measurements of gaseous velocity, species concentration and temperature must also be taken. There is some

progress at measuring species and temperature using nonintrusive optical methods, but, this is a developing and

specialized discipline. Measurements of temperature and

•species can have substantial error. Spray measurements at

the inlet of the combustor have the greatest impact on

numerical predictions. These measurements are also the most

difficult. The flow is optically dense close to the spray

nozzle. Thus, the more accurate spray measurements are at

substantial distances from the spray nozzle. The accuracy

of inlet conditions for spray combustion leaves room for improvement. The gaseous species concentrations may have to

be inferred from the difference of liquid fuel delivered to

1F1II
137

the combustor and measurements of fuel spray flux.

A partial database for spray model validation is given

in Edwards and Rudoff n4. An optically dense zone of roughly 25 mm diameter and 25 mm long is shown in one of the

figures. Spray data are given at an axial station of 50 mm.

Velocity vector diagrams for the gas phase and spray show

that the reverse flow along the combustor centerline reaches

no closer than 25 mm from the spray nozzle. Droplet data at

the first two measurement locations are marked as biased. A

central recirculation zone is caused by swirling inlet flow.

The flow is not contained. Maximum velocities within the

burning flow are of the order of 15 m/s. The rest of the

flow appears quiescent. This type of open, recirculating

flow is particularly difficult to accurately predict. Phase

Doppler Particle Analyzer or PDPA is used in the measurements. In a 1988 paper _!5, the PDPA measurement method was

compared to another method. It was found to agree within

15%. In a 1993 paper 11_, a correction scheme is demonstrated. It was concluded that this correction scheme was

difficult to employ in swirling 3-D flows. A new correction

scheme is being worked on.

In another 1993 paper, Presser et. al. _I_ present some

non-burning and burning data for swirl numbers of 0.0, and

0.53. Laser sheet pictures are shown for swirl numbers of

0.0, 0.39, 0.53, and 0.76. The gas phase flow velocities of

the open combustor are low, giving a moderate sized central
138

recirculation zone. The maximum axial velocity is about 3

m/s. Maximum axial spray velocity is over 15 m/s. Droplet

axial velocities seem to be all positive for a swirl number

of 0.53. Increasing swirl increases the radial displacement

of the spray. Temperature and species concentration measurements are not given.

In another 1993 paper, Ghaffarpour and Chehroudi n8

present temperature and droplet data for a hollow cone

kerosene spray. Droplet characteristics were reported for

mean velocity and mean droplet diameter. Negative mean

axial droplet velocities are shown along the combustor

centerline at distances of 40, 55, and 70 mm for the nonburning case. The first measurement location was 15 mm from

the injector. Mean axial droplet velocities at this station

were all positive, although the velocities were lower for

the burning case. The burning case shows considerably more

radial spreading of the fuel spray. The only gas phase

velocities reported were from the non-burning case, which

indicated a central recirculation region. Measurements of

burning and non-burning droplet flux are nearly identical,

indicating negligible droplet vaporization and/or burning.

Conversely, temperatures within the hollow cone area were

fairly high, about 500 degrees C at 15 mm and 600 degrees at

25 mm. The authors infer that this is due to a downward

movement of the central recirculation zone for the burning

case, but this flow needs to be properly seeded to measure

IL|1 !
139

velocities in this region. Before this case can be used for

spray model development, the gas phase velocities must be

accurately measured and droplet statistics should be specified for various droplet size groups, not just for mean

diameters. Using mean quantities limits the number of droplet groups at the 15 mm location to about sixteen.

The combustion case of Bulzan [108] was used to benchmark the combustion spray model. The flow configuration is

identical to the isothermal case studied earlier, except

that the injector was supplied with liquid heptane fuel.

Temperature, gaseous and liquid velocities were measured.

Data were meticulously taken and repeated to minimize experimental error. Great pain was taken to eliminate 3-Dimensional effects, such as locating the swirler for the coflowing air-stream 140 mm upstream of the injector face.

The combustion experiment does not have combustor walls.

This makes it simpler to get optical measurements, but makes

it difficult for the combustion modeler, as the flow entrains ambient air. Runs were made with entraining boundary

conditions, but this modeler found it best to specify boundary conditions at measured radial locations. The average

radial velocities are rather low but show large variation.

This modeler was not able to duplicate such high turbulence

with low velocity. Initial calculations were performed

using 70 axial gridpoints and 55 gridpoints in the radial

direction. Measurements taken 2.5 mm from the face of the
140

spray injector were used as initial conditions. These data

are much closer to the combustor spray nozzle than reported

by others. The droplet groups are differentiated by droplet

size ranges and location. Close to the injector, droplet

measurements were taken at 1 mm intervals out to a radius of

about 13 mm, where the data rate became negligible. The

droplet data were separated into i0 size ranges or bins.

These bins are then characterized by average diameters. The

average group diameter range is from 6.893 to 122.65 micrometers. These data were reduced to 79 droplet groups as

some velocity data were unreliable due to the scarcity of

data, or the amount of fuel spray in the group was judged as

insignificant. Obviously, the computer time needed to

calculate the spray and its influence on the gas phase will

roughly scale with the number of groups. Data were taken as

a number flux at specified radius. For calculation purposes, the flux data were multiplied by an appropriate area to

get the total number of drops in a group. The total spray

mass was calculated and this was significantly different

than the measured fuel flow rate to the injector, which was

originally taken as 0.38 g/sec. The difference was input as

a combination of fuel and combustion products. Originally,

it was planned to take probe samples for discerning gaseous

species concentrations. Bulzan now feels that sample probes

would distort the flowfield producing inaccurate results.

Measured values of temperature are low in the dense spray
141

due to droplet impingement on the probe. Gas sample probes

would also be greatly affected by the spray. The only data

with which to infer gaseous mass fractions are the total

rate of fuel flow and temperature measurements. The temperatures in this area were inferred from neighboring temperatures not exhibiting droplet cooling. To come up with

appropriate inlet mass fractions, a sub-program was used to

calculate temperature as a function of increasing fuel mass

fraction. A fixed final ratio of unburned to burned fuel

was used while keeping the total enthalpy constant. The

numbers were iterated/modified upon to get the total overall

fuel rate. As initial rough grid calculations produced

higher than measured temperatures, the heating value of the

fuel was dropped by i0 % for better agreement. Initial

calculations were done assuming no correction to droplet

number flux. In McDonnel's thesis, it was concluded that

most of the droplets were not measured in the area of dense

spray close to the injector. If multiple drops are detected

within the measurement zone, the data are discarded. The

measurement of mass fractions were determined to be more

accurate than that of droplet flux, so the droplet flux was

corrected. The droplet flux closest to the injector was off

by a factor of almost 4. Bulzan didn't report species

concentrations which would have enabled correcting inlet

droplet number flux. Also, the data are further complicated

by combustion, widely varying temperature, and very rapid
142

evaporation. Several 70 by 55 grid calculations were made

with different spray flux corrections. It was found that

the best agreement was found by assuming that about 35% of

the fuel and fuel in combustion products were gaseous.

Temperature was significantly overpredicted. Since the

rougher grid calculations were made, it was found that the

total fuel flow rate was actually 20% lower. Succeeding

calculations used a finer grid and a corrected overall fuel

flow rate. The finer grid essentially doubled the number of

radial gridpoints in the fuel injector zone. A velocity

vector diagram of a portion of the flowfield is shown in

figure 5.32. The flowfield structure is simpler than for

the isothermal case. In the combusting case, the recirculation bubble is closed within the confines of experimental

measurement limits. Axial velocities are all positive 200

mm from the combustor injector. In the isothermal case,

reverse flow existed out beyond the last measurement station, and modeling proved to be rather difficult. In addition, the small toroidal recirculation bubble seen near the

fuel injector and co-flowing air stream in the isothermal

case isn't measured or predicted for the combusting flow

case. Calculation of entrainment is still a problem for the

combusting case.

Measured centerline axial velocity next to the fuel

injector is negative. This reverse flow is part of a large

central recirculation bubble caused by high inlet swirl.

Ht ! i_
143

, , , , , i r I

..... , i

, , , , u , J

...... i

....... r

....... r

....... l

r

r

i

i

f

I r

I f

r r l

i I I

i r I

I I

; I

i I

I I

i I

I I I

J

l

I

II/

I/

l I I I I I r illrt_,,,i,

I I l ; r rl,_,,,,,+,

I

r s,,, ........

i r , ...........

/ , , ................

// I ,, ...... , , I,_I,1_LI_II

/ I"" .... ,_ "_ _ _lllllllllll

lllillililll

t ' ..... I IUIIlIIIIII

I ...... I IUIIIIIIIII

I .... iillllllllll

\\

Figure 5.32 Velocity Vector Diagram

II ,-

ii

6\\\\\\\\\\\\k \ ,', _\\\\X\_\\k\\\z, :, I

llililllllli

llilllllilll

llllilllllll

lllliill[lll

illllillll]]

lllilillTHl

lillllllllll
144

Theoretically, the reverse flow should end very close to the

fuel injector. The closest Dr. Bulzan was able to perform

measurements was 1.0 mm from the injector. At this location, the axial velocity was non-negative, but the spray

measurement data were in question. It may be possible to

perform calculations starting from the injector face, but

the fuel stream breakup would then have to be modeled, which

cannot be done at the present time. In this work it was

desired to see how good the spray model was, without coupling it to additional modeling to develop or back calculate

suitable droplet distributions. It was decided to start

with the best possible measured droplet distribution and see

how well the combustion spray model predicted successive

axial spray distributions.

The numerical simulations started 2.5 mm from the

injector face. Plots of the experimental and numerically

predicted axial velocities are shown in figure 5.33-38

corresponding to axial distances of 5.0, i0.0, 20.0, 50.0,

I00.0, and 200 mm from the fuel injector. In taking the

experimental measurements a complete traverse was taken

across the combustor. The measurements on one side of the

combustor are plotted with negative radii. Numerical predictions were 2-D axisymmetric, thus, the numerical results

aresimply reflected for comparison purposes. The axial

velocity predictions are quite good beyond the region of

fuel injector flow. Surprisingly, there is little

FII!
145

4o

3O

_" 20

g

8

o

-10

-2O

0

0 0

0 0 O0

0 0

0

0

O, L 0 _

70 x 55

O 95 x 74

O

O exper

-30 = (_ I I

-0.050 -0.025 0.000 0.025 0.050

RADIUS (m)

Figure 5.33 Axial Velocity at 5.0 mm.

E

0

0

,._1

U.I

>

.,_1

x

<

3O

20

10

-10

0

oo

o

0

0

0

o,,_'_0 0 0/_::_

I

O 70 x 55

O ..... 95x55

OG_O O exper

I I I

-20

-0.070 -0.035 0.000 0.035 0.070

RADIUS (m)

Figure 5.34 Axial Velocity at i0.0 mm
146

Ev

>-

0

o

.J

III

>

.J

_<

X

<

10

-5

-10

0o 0

0 o

0 0

0 0

0

0

o 0

%P

-15 , i

-0.08 -0.04 0.00

RADIUS (m)

Figure 5.35 Axial Velocity at 20.0 mm.

70 x 55

95 x 74

© exper

I I

0.04 0.08

4

v

g

_, o

-4

O

0 0

&.__y ---7ox55

"_--J O ..... 95x74

O O

O_30 o exper

,--6 ! i i 1 [

-0.15 -0.10 -0.05 0.00 0.05 0.10

RADIUS (m)

Figure 5.36 Axial Velocity at 50.0 mm.

I

0.15

iil i ]
147

E

v

>-

0

0

.J

W

>

._J

<

x

<

-1

0

o7 _o ol \<_

o/ \ o o t \_o oo/ \ o / \o

\oOj <%___

k,o o,,/, - _

-3 I I I

-0.15 -0.10 -0.05 0.00 0.05 0.10

RADIUS (m)

Figure 5.37 Axial Velocity at i00.0 mm.

70 x 55

95 x 74

0 exper

1 I I

0.15

4

_2

1

0

-0.15 -0.10 -0.05 0.00 0.05 0.10

RADIUS (m)

Figure 5.38 Axial Velocity at 200.0 mm.

70 x 55

..... 95x74

0 exper

O0 O0

0 0

O0/ O

I I I I I 1

0.15
148

difference between numerical predictions despite the modeling changes.

The inlet reverse flow velocity specified next to the

centerline was about -23.0 m/s. The predicted centerline

velocity at an axial distance of 5 mm is about i0 m/s, while

the measured centerline velocity is about -30.0 m/s. Experimental centerline axial velocity increases while the simulations predict a decrease in centerline velocity. This

poor prediction pattern is also seen for peak axial velocity. The numerical simulations predict peak axial velocity of

approximately i0 m/s at an axial distance of 5 mm, while

measured peak axialvelocity is almost identical to the

inlet peak velocity of 35 m/s. The numerical simulations

are poorly predicting a critical zone of the combustor.

This area of the flow contains the highest concentration of

fuel spray. The gaseous velocities should be correctly

predicted in order that spray predictions are also correctly

,predicted. Many additional calculations were done modifying

turbulent inlet conditions and modifying the turbulence

model in an attempt to better predict reverse flow along the

combustor centerline and peak axial velocity. These simulations did not appreciably improve the prediction of axial

velocity close to the injector, but mostly shifted the axial

velocity peak in the radial direction.

Experimental and predicted peak axial flow velocities

at an axial distance of i0 mm are smaller than at the 5.0 mm

l[!l_
149

station. Predicted peak axial velocities are roughly onethird of the measured velocities. The simulations predict

almost complete merging of the air-assist fuel injector flow

with the surrounding co-flowing air-stream. The simulations

predict a larger radial displacement of the fuel-air injector flow than measured. Predicted reverse flow centerline

velocities are less than half of measured velocities.

The experimental data taken at an axial distance of 20

mm shows distinct peaks due to the fuel-air injector flow.

The simulations are nearly identical except in the central

recirculation zone. The simulations predict a single velocity peak. Qualitatively, predicted axial velocities are

fairly good, except in the central reverse flow region.

Predicted reverse flow velocities are about half of the

measured velocities.

At 50 mm, the axial velocity predictions seem to improve. Single velocity peaks are shown for predictions and

measurements. The simulations slightly underpredict peak

velocity. The quality of the axial velocity predictions in

the central recirculation zone are much improved from previous axial stations. The finest grid shows the lowest centerline velocity.

Peak velocity is well predicted at an axial distance of

I00 mm. There is little difference between the two numerical simulations. The width of the predicted axial velocity

profiles are noticeably thinner than the measured profiles.
150

This implies the simulations may be incorrectly predicting

the total mass flow in the axial direction. This is probably due to the complexity of calculating an open recirculating flow which is entraining ambient fluid.

At 200 mm, the simulations predict nearly identical

profiles. Centerline axial velocities are all positive.

The recirculation zone closed somewhere between i00 and 200

mm downstream of the fuel spray nozzle. Axial velocity is

well predicted from a radius of 5 cm to 15 cm. The simulations predict parabolic-type axial velocity profiles that

are symmetric about the peaks. The experimental profile is

noticeably non-symmetric around the peak. Large axial

velocity gradients are measured about the combustor centerline. The simulations predict a much lower gradient.

Qualitatively, the large velocity deficit at the combustor

centerline is well predicted. The velocity deficit along

the centerline will disappear much more quickly in the

experimental profile than in the predicted profiles. The

local minimum experimental axial velocity isn't located

exactly at the combustor centerline. The experimental axial

velocity profile is slightly shifted to the right. This is

probably due to a slight asymmetry in the flow.

The radial velocities are shown in figures 5.39-44.

The peak radial velocity is underpredicted by both simulations. However, the radial velocity peak is better predicted than the axial velocity peak. The fuel-air flow

li!I i
151

5O

4O

E 3O

0

0

...i 20

u.l

>

,_1

<

10

<

rr

0

0

A

f

/I I

©

70 x 55

O 95 x 74

cO O exper

O

JI

llk\N

I I

I I

-10

-0.050 -0.025 0.000 0.025 0.050

RADIUS (m)

Figure 5.39 Radial Velocity at 5.0 mm.

E

v

O

O

UJ

>

<

D

<

m

30

2O

10

0

O

%

o

O

Co

O

0 0 0

70 x 55

..... 95x74

O exper

O

I I I

-10

-0.070 -0.035 0.000 0.035

RADIUS (m)

Figure 5.40 Radial Velocity at i0.0 mm.

I

0.070
152

E

O

O

LU

>

.J

<

a

<

n"

20

15

10

0

70 x 55

95 x 74

O exper O

O

O

O

0 0 0

0

0

0

0

I I I I

-0.04 0.00 0.04 0.08

RADIUS (m)

Figure 5.41 Radial Velocity at 20.0 mm.

3

2

E

v

O

g

LU

>

" 0

<

O

<

n"

-1

%

oo _o 0

O O O

oOA

3/ \%

0 0

0

0 0

0 0

d_ o_

70 x 55

.... 95 x 74

0 exper

I I ]

-2 . I t l

-0.15 -0.10 -0.05 0.00 0.05 0.10 0.15

RADIUS (m)

Figure 5.42 Radial Velocity at 50.0 mm.

111'IT
153

E

v

>-

0

O

O

.J

W

>

-J -I

<

D

<

-2

O

o °

oo % o o °o

0

©

0

%

I I I I

-3

-0.15 -0.10 -0.05 0.00 0.05

RADIUS (m)

70 x 55

95 x 74

0 exper

I I

0.10 0.15

Figure 5.43 Radial Velocity at i00.0 mm.

1.0

0.5

vE 0.o

0

_ -0.5

I.U

_J

<

a -1.0

<

!1:

-1.5

0 0

0 O0 0 O0

0 0

0

0

0

O

70 x 55

..... 95x74

O exper j

I I I l I

-0.10 -0.05 0.00 0.05 0.10 0.15

RADIUS (m)

Figure 5.44 Radial Velocity at 200.0 mm
154

expansion due to combustion seems to occur more in the

radial direction. The simulations predict a larger radial

expansion of the fuel-air injector flow than measured. The

predicted width of the fuel-air injector flow is nearly

twice as wide as measured. The predicted radial over-expansion destroys the large radial velocity peak at successive

axial locations. At i0.0 mm, the predicted radial velocity

peaks are about twice as far from the centerline as the

experimental peaks. The predicted peak radial velocity is

less than one-fifth of the experimental peaks. Radial

velocities beyond a radius of 4 cm are well predicted.

At 20.0 mm, both the numerical simulations predict a

single radial velocity peak. The flow from the fuel-air

injector has completely merged with the co-flowing stream in

the simulations at this axial location. The experimental

data shows separate velocity peaks. The experimental velocity peaks are about four times higher than predicted peaks.

Radial velocity in the co-flowing stream is well predicted.

The experimental data shows a single radial velocity

peak at the 50.0 mm station. Flow originating from the

fuel-air injector and the co-flowing air-stream have combined. Qualitatively the simulations correctly predict the

radial velocity at this station. The simulations underpredict the velocity peak by about 50 percent. The radial

velocity in the central recirculation zone is off by at

least a factor of four. The negative velocities near the

F[ 1i;
155

combustor centerline indicate the recirculation zone is

starting to close.

The radial velocity profiles for the entraining flow

and outer co-flow are well predicted for the 5.0, i0.0, and

20.0 mm stations. However, the portion of the outer profile

that is in agreement, grows progressively smaller with axial

distance. At 50.0 mm, the agreement in the outer flow is

the smallest. The width of the predicted profile is underpredicted. There is the thought that the initial radial

velocities may be over-specified, causing an over-rapid

broadening of the injector flow. This researcher did a

calculation specifying reduced radial velocities for the

injector flow. This modification remarkably improved the

calculation of the peak velocity flow from the injector.

However, there is little evidence to support such a change

in radial velocities, other than it works. Another possibility of error is the flow direction to the computational

grid. The injector flow is roughly 45 degrees to the grid

lines, which is known to cause accuracy problems. It would

be interesting to implement an improved or higher order

numerical differencing scheme in the model.

At i00 mm, the numerical radial velocity predictions

are all negative. That is, there is inflow across the whole

combustor, no expansion. The experimental data shows some

expansion at larger radii. At the 200 location, the predicted radial velocities are very close to zero. The mea-
156

sured velocities at this station are generally larger.

There is significant non-symmetry for the i00.0 and 200.0

experimental data near the centerline. The velocity should

go to zero near the centerline but doesn't. Part of the

problem is the sign of the radial velocity. It is believed

that the flow is slightly off to one side, which causes

problems in interpreting the sign of the radial velocity

component. The experimental axial velocity profile exhibited sharp gradients near the centerline for the i00.0 and

200.0 stations, thus, there should be significant radial

velocity in this area.

The plots of tangential velocity are shown in figure

5.45-50. As in previous figures, the part of the flow

agreeing least with the experimental data is the injector

flow. The numerical predictions at 5.0, i0.0, and 20 mm

show a large drop in tangential velocity for the injector

flow. The inner profile is worst predicted at the i0.0 mm

station, where it is a quarter of what it should be. The

predicted inner profile actually gains some strength at the

20.0 mm station. The experimental data for the injector flow

exhibits quite a bit of variability and asymmetry. Between

the 5.0 and i0.0 station, the experimental velocity peak

changes sides! The difference in injector flow velocities

from one side to the other is of the order of 25%.

The outer velocity profiles are well predicted at the

5.0, i0.0, and 20.0 locations. Then the width of the

Irl1!
157

15

.. 12

E

v

O 9

O

_J

UJ

>

<

6

z

LU

Z

<

3

©

oo

0

0

0

cO o°,';

0

-0.070 -0.035

70 x 55

95 x 74

O exper

0.G30

RADIUS (m)

,', O

I % _ 0

,I

0.035

J

0.070

Figure 5.45 Tangential Velocity at 5.0 mm

10

_. 8

E

v

m

O

O

LU

>

J

_<

IZ

LU

Z

<

I--

6

70 x 55

..... 95x74

O exper

O

O

O

O

O

O

oO O O3 o_oO° _ o oc_ /_o

-0.035 0.(_0 0.035

RADIUS (m)

Figure 5.46 Tangential Velocity at i0.0 mm

I

0.070
158

6

i

O

O

w

> 4

<

i

lz

LU

C

Z

< 2

k70 x 55

95 x 74

O exper

O

0 0 0

0

oO o o

0 0 000

%o h,

//oooooo /t°

/

; \°J

-0.070 -0.035 0.000 0.035 0.070

RADIUS (m)

Figure 5.47 Tangential Velocity at 20.0 mm

4

g

g 3

Z

_<1

70 x 55

..... 95x74

O exper

0 i

-0.15 -0.10 -0.05 0.15

O

O

O°Cb

o o:P

O

)O

)O

O

\

0.00 0.05 0.10

RADI_S (m)

Figure 5.48 Tangential Velocity at 50.0 mm

[_ : r- Ji!I i
159

5 70 x 55

O

O 95 x 74

0 0% 0 exper

E ,

0 a

0

_J

w

>

_.1

_ 2

z

u.l

I--- 1

0

V

-_.15 -0.10 -0.05 0.00 0.05 0.10

RADIUS (m)

Figure 5.49 Tangential Velocity at i00 mm

I

0.15

.-. 4

(n

E

v

O 3

9

LU

>

...J

<

_- 2

z

ILl

(9

Z

<

I--

1

0

70 x 55

.... 95 X 74

O exper

O

O O

O

o

oooO¢ I Y I\ooo _J

-0.10 -0.10 --0.05 0.00 0.05 0.10

RADIUS (m)

Figure 5.50 Tangential Velocity at 200 mm

J

0.15
160

profile is off, but predictions significantly improve! At

50.0, i00.0, and 200.0 mm the prediction of tangential

velocity is good. The velocity peaks are only slightly

under-predicted. The width of the profile are underpredicted, caused by a sharper drop in tangential velocity than is

measured. The inner profile is underpredicted at the i00.0

and 200.0 locations. This was also seen for the axial

velocity profile. The experimental data indicates stronger

radial inflow near the centerline. Viscous effects are

acting on the on the strong velocity gradient near the

combustor centerline. The measured turbulence is quite high

at this location. At 350.0 mm, the experimental axial

velocity profile only shows a slight deficit along the

centerline.

Plots of temperature profiles are shown in figure 5.51-

56. The agreement on the left side of figure 5.51 is good.

The left side agreement is noticeably better for the 5.0

"through 20.0 mm stations. It appears that initial combustion is better on the left side for these stations. This is

the side that the flow seems to shift or expand away from.

The predicted profiles have rather sharp ridges at the outer

edges of the high temperature region that are probably nonphysical. The depressed temperatures due to fuel droplet

impingement on the thermocouple probe are evident. The

model predicts lower temperature in the fuel spray zone.

Both experimental and predicted temperature profiles show

i!l! i_
161

21 [

1800

!

I

I I

0

70 x 55

95 x 74

exper

%" 1500

w

_ 1200

LLI

O.

LIJ

I-- 900

0

0

O O 0

O

O

O O

O

O

O

O

O

6OO

O

O

3OO

-0.070 -0.035 0.000 0.035

RADIUS (m)

Figure 5.51 Temperature at 5.0 mm

0.070

2100

1800

1500

LIJ

i--

< 1200

¢t"

W

13-

LU

I-- 9OO

60O

/I

/ I

r I I

!

I

!

0

0 0

0

(b

0

0

0

3O0

-0.070 -0.035 0.000 0.035

RADIUS (m)

Figure 5.52 Temperature at I0.0 mm

0.070
162

2100

1800

1500

w

re

D

I.-

< 1200

rr

LU

13-

W

I.- 9OO

6OO

0

0

0

300 l

---0.08 -0.04 0.00 0.04

RADIUS (m)

Figure 5.53 Temperature at 20.0 mm

0.08

2100

1800

v 1500

w

n"

E)

I-

< 1200

n"

ILl

n

W

I.- 90O

6OO

I I t

I |

1 I

% I

\ I

0

0

300 i

-0.10 -0.05 0.00 0.05

RADIUS (m)

Figure 5.54 Temperature at 50.0 mm

i

0.10

lrl ! i
163

2100

1800

1500

UJ

I-

.< 1200

¢t"

UJ

0.

ILl

I-- 90O

6OO

0

D

?

0

300

•...-0.150 0.000

RADIUS (m)

o

-0.075 0.075

Figure 5.55 Temperature at 100.0 mm

70 x 55

95 x 74

exper

I

0.150

2100

1800

1500

LU

_C

D

I--

< 1200

W

OW

_- 900

60O

30O

-0.15 -0.10 -0.05

0 0

0 0 o

_noO°,°j , ,

0.00 0.05 0.10

RADIUS (m)

Figure 5.56 Temperature at 200.0 mm

I

0.15
164

very large gradients at the edges where combustion is not

occurring. The finer-grid calculation predicts higher

temperatures than the coarse grid. At 50.0 mm, the experimental temperature profile is fairly flat across the combustion zone. The coarser grid exhibits two small humps, while

the fine grid has very prominent peaks. The peaks are much

smaller at the I00 mm station. The width of the temperature

profile is slightly underpredicted at the i00.0 mm station.

The gradient at the edge of the high temperature profile is

overpredicted. This is probably due to the underprediction

in radial velocity. At 200 mm, the experimental and predicted profiles show a rather bell shaped curve. The experimental profile shows a little more spreading, again probably due to the radial velocity.

Fuel Droplet Comparisons

The droplet axial velocity comparisons are shown in

figures 5.57-60. Droplets were not measured past the 50 mm

station. The experimental data are shown for negative radii

in the figure, and the fine grid numerical predictions are

shown for positive radii. The solid line is the gas phase

velocity shown for comparison purposes, as the droplets will

eventually relax to the gas phase velocity due to drag. The

experimental gas phase velocity plots were not smoothed.

The different symbols denote the different size droplet

groups. Only 6 of the 10 droplet sizes are plotted. At 5.0

mm, the numerical calculation of spray droplet axial

IllI ',:
165

40

o

30

0

[] A

0.000

RADIUS (m)

gas

0 6.9

[] 15.4

23.8

¢ 40.8

v 66.2

iU 97.2_0

0.010 0.020

Figure 5.57 Droplet Axial Velocity at 5 mm

35

25

15

g 5

_1

<

-5

<

-15

0.00

RADIUS (m)

I

-25

-0.04 -0.02

<>

Z_

L_ V

Z_D

[]

gas

O 6.9

[] 15.4

Z_ 23.8

¢ 4O.8

V 66.2

0.02 0.04

L_

Figure 5.58 Droplet Axial Velocity at i0 mm
166

E

20

10

i

O

o, 0

LU

>

_J

<

X

<

-10

-20

-0.050

00_

0 0[3

,-:,0 O O

<> O_

O0 _vv_ 0

z_

gas

0 8.9

D 15.4

Z_ 23.8

<> 40.8

_7 66.2

[2 97.2

I I I I

-0.025 0.000 0.025 0.050

RADIUS (m)

Figure 5.59 Droplet Axial Velocity at 20 mm

o9

O

o,

LU

>

_<

15

I0

gas

O 6.9

[] 15.4

Z_ 23.8

o_$ <> 40.8

O [] _7 66.2

x

v 0

-5 j l J J

--0.10 -0.05 0.00 0.05 0.10

RADIUS (m)

Figure 5.60 Droplet Axial Velocity at 50 mm

i_i| i
167

velocity is slightly underpredicted. This is caused by the

under prediction of the gas phase axial velocity near the

vicinity of the fuel injector. The larger difference between the predicted gas phase and droplet velocities results

in higher drag forces, and a quicker reduction in droplet

velocities. The smaller droplets generally track the gas

phase velocity better. The shape of the droplet velocity

profile is good, except, the droplet groups show greater

radial spreading than experimentally measured. The experimental data at 5.0 mm indicates many particles near the

combustor centerline. Many of these particles show negative

axial velocities! The predictions do not show the same

patterns. There are few predicted droplet groups near the

centerline, and only one with negative axial velocity.

Gaseous fuel is provided by evaporating droplets. The

combustion of this fuel results in expansion and possible

increase in reverse flow velocity. The experimental data

shows an increase in negative axial velocity at this location, while the predictions show a decrease. The spray

model is not picking up on possible phenomena to explain

this effect. In setting up the droplet group input data,

droplet groups with negative axial velocities were ignored.

It would be useless to input this data. The spray flowing

towards the injector significantly contributes to gaseous

fuel mass fractions, but without doing some species measurements, it is unknown.
168

At I0.0 mm, The droplet axial velocities have decreased further. The predicted gas phase velocity does not

even show a peak corresponding to the air-assist stream from

the nozzle. The difference in predicted velocity for the

two phases is much higher than the experimentally measured

velocity difference. This will lead to further loss in

droplet axial velocity. The predicted droplet velocities

exhibit a much sharper peak than measured and show greater

radial spreading. The experimental droplet data shows

grouping of the data according to droplet size. Droplet

groups having the same average diameter seem to plot along

smooth curves. This isn't as obvious in the numerical

predictions, especially, later in the flowfield, due to

scarcity of particles.

The number of predicted droplet groups decreases as

the droplets are burned. There are only 5 droplet groups

plotted at the 20.0 location, and two at the 50.0 mm station. This is far too few show droplet size grouping. It

is also questionable if this small number could cause any

kind of local gas phase velocity peak. The experimental

data shows good correlation between the injector gas velocity peak and the placement of the droplets. The simulations

don't adequately predict the air-assist nozzle flow. To

help validate the spray model, it would be better to have

significantly more droplet groups. Unfortunately, the

r

number of droplet groups was rather limited to begin with!
169

The current spray model significantly taxed computer time

resources. The needed computer time would scale directly

with the number of droplet groups. The experimentalist can

widen the measuring probe volume to pick up more droplets

where the spray is less dense. There isn't a numerical

modeling counterpart to this. If one mentally averages the

droplet axial velocities, the simulation seems about at

third of the experimental average for both the 20.0 and 50.0

measurement stations. Intuitively, this should have a

severe effect on the spray vaporization. Dr. Bulzan integrated his spray flux measurements across the flowfield.

Correcting these measurements gives 21.3 % of the metered

flow rate at 2.5 mm, 26.8 % at 5, 50.0 at i0, 23.1 at 20,

and 2.6% at 50.0 mm downstream of the nozzle. It is obvious

that much of the droplet data are not picked up at the 2.5

and 5.0 stations. The predictions presented here used 65%

of the metered fuel flow for droplet input at 2.5 mm. The

•simulation predicted 51.7 % of the metered flow in the form

of spray at 5, 24.0 % at i0, 6.3 % at 20, 0.13 % at 50.0 mm.

Thus, the simulation does predict over-rapid fuel spray

vaporization. However, as velocity predictions near the

injector zone are overly diffusive, it would be inappropriate to blame the spray model. The droplet velocity is about

0.8 of what it should be for the 5.0 and 10.0 location. At

the 20.0 location it is roughly a third and at the 50.0

location it is about a quarter. Assuming the amount of
170

droplet vaporization is directly proportional to time,

applying a correction for the time by factoring in how much

I

the gas phase velocity differs from the data and compounding

this correction at successive axial locations to the above

percentages gives 57.4 % at 5, 33.3 % at i0, 15.5 at 20.0

and i.i % at the 50.0 location. These amounts are still too

low. If one looks at the 20.0 location, this suggests that

the spray is off by a third. This would then imply that all

the nearly all the fuel should be in liquid form at the

calculation start. This does not agree with the data, as

inlet temperature is very high.

The droplet radial velocities are shown in figure 5.61-

64. At the 5.0 mm station, there is surprising little

velocity difference between the gas phase and spray droplets. This was also exhibited by the plots of axial velocity. The predicted gas phase radial velocity peak is only

slightly less than the measured one. The axial velocity

peak was a factor of three off. However, as previously

noted, the profile is twice as wide as it should be, which

causes additional nozzle flow dissipation. The numerical

simulation predicts droplets at larger radii. The cause is

a combination of the predicted axial and radial gas phase

velocities. The calculated radial gas phase velocity profile is much broader than the experimental one, and the predicted axial velocities are much lower than measured. Thus,

the predicted spray spreads out further than measured. At
171

E

O

O

W

>

<

<

n4O

3O

2O

I0

O

OO

I I

°

[]

-10

-0.020 -0.010 0.000 0.010

RADIUS (m)

A

gas

6.9

15.4

23.8

40.8

66.2

97_?.J

0.020

Figure 5.61 Droplet Radial Velocity at 5 mm

40

3O

_o

o

_: 1o

0

%_ o

/

-10 i I I J

-0.04 -0.02 0.00 0.02 0.04

RADIUS (m)

Figure 5.62 Droplet Radial Velocity at i0 mm
172

25

6

[]

0 0

19 O O g []

_" " o

.-. _0 0

-- O v A

O v v vD

ow ;

> [yv

_ 7

__ []

121

<

1

-5 i J i J

-0.050 -0.025 0.000 0.025 0.050

RADIUS (m)

Figure 5.63 Droplet Radial Velocity at 20 mm

I0

8

vE 6

gas

v O 6.9

°_°_ [] 15.4

z_ 23.8

<> 40.8

<> v 66.2

O AAA

A

O A []

4

An

o g_n

0

-2 i l i

-0.070 -0.035 0.000 0.035 0.070

RADIUS (m)

Figure 5.64 Droplet Radial Velocity at 50 mm

Iii1 i!
173

i0.0 mm the predicted gas phase radial velocity peak for the

injector flow is much lower than the measured one. The

experimental data at this location only shows one group with

negative radial velocity near the centerline. The average

diameter of this group is small, which allows it to better

track the local gas phase velocity. As mentioned previously, the reverse flow velocity at this point is decreasing.

As previously discussed for droplet axial velocity, the predicted droplet groups display more radial dispersion than

measured. At 20.0 mm, the predicted radial velocities are

generally much smaller than the measured. There are 5

predicted droplet groups at this axial location. Some of the

groups are in the high temperature flow. All of the measured groups are within the high temperature flow. At 50.0

mm there are only two predicted droplet groups and both of

these are in the lower temperature flow. The measurements

display many droplet groups at the 50.0 mm station. This

suggests that the spray model is overpredicting evaporation.

The experimental data at this location shows much higher

radial velocities, both in gaseous and liquid phases. The

spray is going to provide for gas phase source terms further

downstream as the fuel evaporates and burns. The small

number of predicted droplet groups means that most of the

fuel is already gaseous and thus, already burning. From

this information, it may be inferred that predictions should

not noticeably improve with downstream distance. This was
174

seen in previous figures of velocity and temperature.

Tangential droplet velocities are shown in fig. 5.65-

68. These velocities are also under-predicted, at times,

half the measured values. The calculated droplet velocities

are nearly the same as the calculated gas phase velocity,

while the experimental droplet velocities can be considerably lower than the measured gas velocities. The experimental droplet data shows noticeable stratification by droplet

size. The smallest diameter droplets are closest in tangential velocity to the gas phase, while the largest diameter

droplet groups differ the most. At the 20.0 mm station, the

droplet groups plot nearly horizontally according to droplet

size. The predicted droplet groups show slight pattern

according to size at the 5.0 mm station. At later stations,

the smaller number of predicted droplet groups does not

allow a conclusion. The low number of predicted droplets

suggests over-rapid vaporization and/or combustion. The

numerical calculations do not adequately predict the fuelair nozzle flow. Within a short distance, the injector flow

has been seemingly dissipated.

The experimental data shows higher peak radial velocities than axial velocities for the 2.5 mm station to the 20

mm station. D. Bulzan [108] concluded that the flow was

expanding more in the radial direction due to the strong

central reverse flow caused by the recirculation zone. The

poor velocity predictions for the injector nozzle flow are

IF[li
175

E

v

O

O

W

>

_J

<

_--

Z

LU

0

Z

I--

15

10

5

0

o

[] 15.4 '_>q/ u-"

Z_ 23.8

<> 40.8

V 66.2

U 97.2

I

-0.010

I I )

-5

-0.020 o.ooo O.OLO 0.020

RADIUS (m)

Fig. 5.65 Droplet Tangential Velocity at 5 mm

10.0

" 7.5

LU

> 5.0

_.1

Iz

I..U

L9

Z

< 2.5

I--

I I I

0.0

-0.04 -0.02 U 0.00 0.02 0.04

RADIUS (m)

Fig. 5.66 Droplet Tangential Velocity at i0 mm
176

10.0

7.5

n

0

0

._1

LU

> 5.0

<

FZ

uJ

0.0

gas

O 6.9

[] 15.4

Z_ 23.8

0 40.8

v 66.2

U 97.2

0000

I I 1

-0.050 -0.025 0.000 0.025 0.050

RADIUS (m)

Fig. 5.67 Droplet Tangential Velocity at 20 mm

I0

.. 8

g)

gas

O 6.9

[] 15.4

z_ 23.8

<> 40.8

0 6 _ 66.2

Z

0 I

-0.08 -0.04 0.00 0.04 0.08

RADIUS (m)

Fig. 5.68 Droplet Tangential Velocity at 50 mm
177

probably not the fault of the combustion model. The overly

large radial dissipation of the injector flow must be investigated. This researcher feels the most probable modeling

improvement would be an improved numerical differencing

scheme to adequately calculate the high velocity nozzle flow

which is severely skewed to the current cartesian grid.

Modification of the spray or combustion models will probably

provide little improvement. Indeed, a preliminary calculation was done using a pdf combustion model. While there are

many differences in the modeling than used here, the gas

phase velocities were very similar. The largest error was in

calculating the injector flow, which is simply where the

fuel spray is. If the fuel source terms to the gas phase

are incorrectly calculated, a combustion model for the gas

phase isn't the correct fix.
178

5.3 - Monte Carlo Pdf Model Predictions

Four cases are predicted using the hybrid pdf method.

The stochastic part of the coding solves for enthalpy and

species while velocity and pressure are solved using the

SIMPLE algorithm of Patankar. The hybrid pdf predictions are

compared with predictions made by a more traditional combustion model which is a combination global arrhenius reaction

rate and eddy breakup or mixing model, which was also used

in the previous spray combustion calculations. These predictions were generally optimized to obtain good comparisons

between measurement and predictions.

The first case compares species concentration and temperature for an enclosed diffusion combustor. In this combustor, the fuel and air supply ducts are similar in size.

The second case also involves a diffusion combustor, but the

geometry is much closer to an actual gas turbine combustor.

The fuel enters through a small conical duct which is at a

large angle to the combustor axis. The third case involves

premixed combustion. Premixed combustion is being studied

as a way of lowering emissions and improving fuel efficiency. As simple combustion models may poorly predict premixed combustion, the ability of a combustion model to accurately predict both diffusion and premixed combustion is a

desirable feature. The fourth case involves a swirling

hydrogen jet diffusion flame. This combustor is extremely

Irl! ii
179

interesting as temperature pdfs were measured.

5.3.1 Non-Swirling Diffusion Combustion

The data of M. H. Lewis and L. D. Smoot n9 was chosen

because of the simple geometry. Fuel and air are supplied

by concentric tubes to a large insulated cylindrical combustion chamber. The combustion chamber is 152.4 cm long and

20.32 cm in diameter. The fuel inlet tube has an inside

diameter of 1.6 cm and wall thickness of 0.32 cm. The fuel

inlet tube has a blunt end. The outside radius of the air

annulus is 2.86 cm. This type of research combustor is

known as a dump combustor. The large increase in crosssectional area is easy to model and causes a hot recirculating flow which provides a continuing ignition source for

combustion. Lewis and Smoot did not measure velocities,

although it appears that Smoot and another of his students,

Philip Smith 12° used a CFD code to calculate flow properties.

The fuel is city-gas which is 88.53 % methane, 7.44 %

ethane, 2.55 % nitrogen, 1.39 % carbon dioxide, and 0.09 %

hydrogen. In this work, the fuel is modeled by assuming

pure methane. The experimental inlet velocity and temperature of the fuel was reported as 21.3 m/s and 300 degrees

Kelvin. 20.14 m/s was used in the numerical simulation so

that the overall fuel to air ratio was stoichiometric. The

inlet air velocity and temperature was 34.3 m/s and 589

degrees Kelvin. Temperature was measured at axial distances
180

of 9.5 and 39.5 cm from the combustor inlet. Species were

measured at axial distances of 9.5, 17.5, 24.6, 32.7,47.6,

63.2, 78.5, and 137.5 cm. The calculation modeled the

combustor as being only 1.0 meter long, as initial pdf

simulations exhibited numerical difficulty due to large

aspect ratio grid cells. The grid was stretched in the

axial and radial directions with the finest grid near the

inlet.

The grid has 80 axial and 40 radial gridpoints. The

fuel jet has 5 radial cells and the air has Ii radial cells.

A recirculation zone is predicted. The length of the recirculation zone varies between the hybrid pdf and SIMPLE

calculations. The SIMPLE simulation predicted a reattachment length of 46.5 cm while the pdf simulation predicted

57 8 cm.

Smith and Smoot used an assumed shape pdf model with an

equilibrium combustion model. This model significantly

over-predicted carbon dioxide, especially along the centerline. Oxygen was under-predicted in the same area. This

implies that burning was not occurring as rapidly as predicted by Smith and Smoot's calculation. A finite rate

single step combustion model is used in the simulations

reported here. The model is:

CH 4 + 202 -" CO 2 + 2 H 20

If!II
Where

R e =min[A(P_Ye)"(To)be-_/aT,-- we g" (9)]

A = 8,5 xl013 [ Kg-mole ]

m 3 sec

181

(5.6)

The exponents used for the fuel and oxidant concentrations are 0.2 and 1.3 respectively. These are the third set

recommended by Westbrook and Dryer [25] for Methane. The

first two sets recommended a negative exponent on the fuel

which resulted in some numerical difficulty with large

reaction rates at very low fuel concentrations.

Comparisons of fuel, oxygen, water, and carbon dioxide

mole fraction at an axial distance of 9.5 cm are made with

experimental measurements in figure 5.69. There is little

difference in the calculated species concentrations at this

axial station. The simulations do a very good job of predicting the fuel concentration, except for over-predicting

the peak fuel mole fraction. Up to a radius of 2 cm, the

oxygen concentration is well predicted by both simulations.

Beyond a radius of 2 cm, the oxygen concentration is overpredicted. The pdf simulation gives slightly lower oxygen

concentrations at large radii.

The concentration of CO 2 was found from the concentration of fuel, oxygen, and atom balances. Using a finite

rate reaction model for the prediction of CO 2 resulted in
182

0.10

0.08

._ 0.06

O3

:D

0.04 .......................................................... i............................. i.............................

0.02 _ ..........................

0.00

0.0 0.2 0.4 0.6 0.8

FUEL MOLE FRACTION

0.10

0.08

0.06

n- 0.04

0.02

0.00

0.00

............................ .a...............;.............. i...............................................

0

0

i

...................... _ .................. .:...................."_,_..................... !.......................

............. i............. ;..........

0.05 0.10 0,15 0.20 0.25

OXYGEN MOLE FRACTION

0.10

0.08

_ 0.06

O3

0.04

0.02

0.00

0.00

f 0.10

i

i !

, O i

+ +

+ J

+ J + O

...........i ............................ ! o .....!.......

! I i !

............................. _............ _...................... :-'""+0 .....................................

_o + +

..................... ...................... •........

0.20

0.08

._ 0.176

D

5

<

0.04.

0.02

i

#

: o

.......................... i ......................................................................................

i o

, ! O

................... _ ......................... i-o................................................

,o i++i

O

Q

0.00

0.04 0.08 0.12 0.16 0.00 0.03 0.06 0.09

WATER MOLE FRACTION CARBON DIOXIDE MOLE FRACTION

0.12

Figure 5.69 Species Mole Fractions at 9.5 cm

i +--
183

superfluous product concentrations at high temperatures.

While the predicted concentration of fuel and oxygen can go

to nearly zero, the rest of the arrhenius rate term is very

large at high temperature. Initial computation product

concentrations were 30 per cent above the theoretical maximum concentration. This has serious consequences for multiple reaction step models.

Beyond a radius of 2 cm, product concentrations are

underpredicted by about 50%. This is about the level oxygen

was overpredicted. The simulations underpredict combustion

in the outer recirculation zone. There is some disparity in

product concentration near the combustor axis. The predictions show zero product concentration at this location.

Measurements show significant CO_ concentration, but near

zero water concentration. The temperature profiles for this

axial location are shown in figure 5.70. The temperature is

lowest near the combustor axis. The simulations predict

"temperatures about 200 degrees cooler along the combustor

axis. Predicted temperatures off the combustor axis are very

good, despite the error in predicted species concentrations.

Temperature increases with radius up to 5 cm. Beyond a

radius of 5 cm, temperature is fairly constant. The pdf

simulation which shows a slight increase in temperature

towards the combustor wall. The outer wall was modeled as

adiabatic since the combustor wall was insulated. The

experimental profile shows temperature slightly dropping
184

0.10

0

0.08 ............................................................. O.

E 0.06

<

_ 0.04

0.02

0.00

300

i li

i,

il

!

o :i

_ri i

=

....................................... _........... 0 'i ......................................

i /

............. ........... i.................... l' 'r...........

700 1100 1500 1900 2300

TEMPERATURE

Figure 5.70 Temperature at 9.5 cm

near the wall.

M. Nikjooy and R. So [28] did some calculations for

this same combustor. They used a two-step finite rate

combustion and an equilibrium combustion model. The use of

an equilibrium combustion model implies that infinite time

exists for combustion or that combustion is extremely fast.

The finite rate constant was higher than that employed in

this work and very rapid combustion was predicted. Their

calculations predicted peak temperatures above 2000 degrees

close to the air and fuel inlet ducts, which is not in

agreement with the experimental data.

Comparisons of species at an axial distance of 17.5 cm

for the present simulations are shown in figure 5.71. The

SIMPLE simulation best predicts the fuel concentration. The

pdf simulation slightly over-predicts the centerline peak

fuel mole fraction and the fraction at a radius of 2 cm.

I[I:iii
185

0.10

0.08

,_ 0.06

09

5

0.04

0.02

0.00

0.0

............................. I ° I"per i

i

0.1 0.2 0.3 0.4

FUEL MOLE FRACTION

0.10 '

0.5

0.08 ....................... _....................... 7,_i................ T................... :.......................

i ', i

i ,, i !

_" 0.06 i _, i

_ i i ,,i

....................... .................... o.................... "4,\......................................

re 0.04 i ! i_,,

0.02 .......................... i .........................

# i "

0.00

0.00 0.05 0.10 0.15 0.20 0.25

OXYGEN MOLE FRACTION

0.10

0.08

._ 0.06

(o

:D

<

n- 0.04

0.02

0.00

.

/ o !

........................_...................... '_ ................... 0 ................ .......................

! ! I !

:._ i 0

..................i" ";/ ..............i.......O ............+..............................................

. ": o !

.................... i................. .....................

! :

0.00

0.10

0.08

"_ 0.06

0.04

0.02

j i

J

! / O

i :

............................."r........... ,..............r'"'" O ..............................................

/ ! o

...... i , ........ i 0 ............ :..................

t ! :

:" 0

,/

0

0.00 --.-0

0.04 0.08 0,12 0.16 0.20 0.00 0.03 0,06 0.09

WATER MOLE FRACTION CARBON DIOXIDE MOLE FRACTION

0.12

Figure 5.71 Species Mile Fractions at 17.5 cm
186

The hybrid pdf simulation better predicts oxygen mole fraction. Oxygen concentration is well predicted up to the peak

oxygen concentration at a radius of 2 cm. Beyond the peak,

both simulations overpredict oxygen concentration. The

largest error is at the combustor wall. Product concentrations are underpredicted by the simulations as they were

at 9.5 cm, although the hybrid pdf simulation is improving

compared to the SIMPLE simulation. The anomaly in product

concentrations near the combustor axis still exists.

The species mole fractions at 24.6 cm are shown in figure

5.72. The peak fuel mole fraction is significantly overpredicted. At larger radii, the fuel predictions are excellent. The pdf simulation again predicts higher fuel concentrations. Paradoxically, the pdf simulation also predicts

lower oxygen concentration than the other simulation. Peak

oxygen concentration is well predicted. Oxygen concentration is slightly underpredicted by both simulations near the

combustor axis. Oxygen concentration beyond a radius of 2

cm are overpredicted, but are improving compared to previous

comparisons. Both simulations underpredict product species

\

concentration, with the exception of- water near the combustor axis. The hybrid pdf simulation better predicts oxygen

and product concentration. Carbon dioxide is detected near

the combustor axis, but not predicted. The H20 concentration near the combustor axis is very low for the predictions

and measurements. Predictions and measurements of H20

I[I11
187

0.10

0.08

0.06

0.04

0.02

..... -- _seline 1

............................ I 0 J exper I....................

0.10

:.................... i....................... i.............................................. i........................

o "-.. ............. _................. i............. !..............

io i _. ! !

0,1 0.2 0.3 0.4 0.5

FUEL MOLE FRACTION

0.08

0.06

o')

5

<

tr 0.04

0,02

0.00 0.00

0.0 0.00

i

0

0 ,_

! 0 ,

...................... +....................... ........ 0...,.-.','.'_ ..................... ;........................

: : '//[0

0.05 0.10 0.15 0.20 0.25

OXYGEN MOLE FRACTION

0.08...................... ;.................................................... o......... _...................... o.o_ : ,-,

, i 0

i /' i _ / i

i l !

"_ o.o_ ................i..................... i ° ..... i ................ _ 0.06 ......................... i i ° ....................................

0

oo,I....................... _.. .............. _....................... _........................ / 004 // .................... o..........

0.02 ..................................... !......................002 i o ...... i..........................................................................

I" I '

0.00 0.00

0.00 0.04 0.08 0.12 0.16 0.20 0.00 0.03 0.06 0.09 0.12

WATER MOLE FRACTION CARBON DIOXIDE MOLE FRACTION

Figure 5.72 Species Mole Fractions at 24.6 cm
188

of H20 concentration increase towards the combustor wall. A

peak in CO2 concentration is measured at a radius of 6 cm,

while the simulations show CO 2 concentration increasing

towards the combustor wall.

Predictions of fuel concentration along the combustor

axis worsen at the next axial station. At 32.7 cm, both

simulations over-predict peak fuel concentration by thirty

per-cent. Paradoxically, predictions for oxygen are excellent, much improved from previous axial comparisons. Under

a radius of 4 cm, the pdf simulation only slightly underpredicts oxygen. Much improvement is also seen in product

species predictions. The carbon dioxide concentration is

slightly underpredicted. Near the combustor axis, H20

concentration is overpredicted by the hybrid pdf simulation,

but the pdf simulation best predicts carbon dioxide. As a

single step reaction is being employed in the combustion

models, predicted product concentrations are related. The

experimental data does not show this dependence between

product concentrations near the combustor axis. Thus, evaluating only one product species at this location could lead

to incorrect conclusions about the predictive capability of

combustion models.

Figure 5.74 gives the predicted and corrected temperature profiles at 39.5 cm. The hybrid pdf simulation best

predicts temperature near the combustor axis at this location. However, temperature is significantly overpredicted

IFlI i
189

0.10

0.08

_v0.06

O9

:D

<

c_ 0.04

o.o2

0.00

0.0

• _ _fseline l

............................................... i........... I o I.xp_ i

................... i....................... i....................... !..............................................

o...i i i

•" o _'. ..................!..................... i.......................!..............

o "-< i

_ii, i

0.1 0.2 0.3 0.4 0.5

FUEL MOLE FRACTION

0,10

0.08

_,_ 0.06

O9

_=

0.04

0,02

0.00

0

....................... i...................... i....................... i..............................................

! o

i o, i

i 0 _

...................... 7.................. !........_, "-i ....................... :.......................

................... .i...................... i................ i: ....................

i i,

J

0.00 0.05 0.10 0.15 0.20 0,25

OXYGEN MOLE FRACTION

0.10

0.08

.._ 0.06

(n,

D

5

<

0.0.4

0.02

0.00

0.00

0.10

i , O

! i b

................ i............... ! -, oi .... _........... _ o.o6 f , f )

i i O9

i ' O D

t

q¢

.....................;................... '_......................_.............................................. _ 0.04

," i : 0.02

0.04 0.08 0.12 0.15 0.20

WATER MOLE FRACTION

0.00

o.0o

/i O

............................ i....................................................................................

i io

.............i.........., o .i.................. i.....................

i , i

" i

/ i "

2.-.-.0 i , :,.

0.03 0.06 0.09 0.12

CARBON DIOXIDE MOLE FRACTION

Figure 5.73 Species Mole Fractions at 32.7 cm
190

0.10

i I

i

0

i _ i /,

' " 0 i

_o._ ...... I-P-I _ ..................... i...................

0 exper i

D : Oi

_ _ _ ,

< _ ! i ,

0._ ......................+....................... i'"O ...............!...... ;_...................................

0.02 .................. i.............. . .................. _........................................

• i

x E

zO !

0.00

_0 700 1100 15_ 1_0 230

TEMPE_TURE

Figure 5.74 Temperature at 39.5 cm

beyond a radius of 2 cm. At a radius of 8 cm, the error in

predicted temperatures is 400 degrees. Experimental temperatures near the combustor wall suggest some heat transfer is

taking place. Measurements of species concentration were

not reported at this axial location.

Species mole fractions at an axial distance of 47.6 cm

are shown in figure 5.75. Fuel mole fraction is overpredicted by over 50%. Comparing fuel profiles for previous

axial stations, this is the station with the largest relative error. The peak fuel concentrations have dropped by

about a factor of two from the previous axial station. The

fuel concentration profile is also spreading. Oxygen concentration is overpredicted near the combustor axis. Predicted oxygen concentration is lower than at the previous

station. There is significant combustion near the combustor

axis. Significant H20 and CO 2 concentrations are measured

H!I i
191

0.10

0.08

0.06

D

w

<

n- 0.o4

0.02

0.00

0.0

I--Iba,°,,_/

............................................................... I 0 Jexper Ipdf ...................

i :

O. 1 0.2 0,3 0.4 0.5

FUEL MOLE FRACTION

0.10 i

J tl

o.o81 ............... l°i....................... T......................T...............................

i i i

- i i i

0.04 .............................. ;....................... !...............................................

!o' i

0.02............................... ...................... ;...................

i o° ',

O.O0

0.00 0.05 0.10 0.15 0.20 0.25

OXYGEN MOLE FRACTION

0.10

0.08

.E 0.06

¢,o

0.04

0.02

0.00

O.O0

,j

...................... T....................... ....................... ,............. f.-/.....................

.......................i ....................... i...................... i ..... _ ............

i : i 'o

....................... i....................... i..........................................

i i o i

i _ ,_ _ :

.................. !_.................... :i : ,: O _............ .... :.........................................

} ,o i

i , [

0.04 0,08 0.12 0.16 020

WATER MOLE FRACTION

0.10

0.08 I......................................................... =............................................... Jl

l ! i i°/

_ i i ,

_ 0.04............................ i............................. i-.'o................... i............................

0.02 ....................... i"--'_ 0 i .............. i .....................

o.oo , ; n

0.00 0.03 0.06 0.09 0.12

CARBON DIOXIDE MOLE FRACTION

Figure 5.75 Species Mole Fractions at 47.6 cm
192

and predicted near the combustor axis. The pdf simulation

of product concentrations are very close to the experimental

values. As the pdf simulation overpredicts reactant concentration near the combustor axis, this shows some error

possibly due to the one-step reaction mechanism. By not

allowing intermediate chemical species, product species can

be over-predicted. The conventional combustion simulation

underpredicts product concentration near the combustor axis

where fuel and oxidant are overpredicted.

Reactant concentrations at an axial distance of 63.2

cm, shown in figure 5.76, continue to rapidly decrease.

Both simulations under-predict oxygen near the combustor

wall. The pdf simulation predicts higher reactant concentration and lower product concentration than the conventional simulation. Predicted product concentrations are much

more uniform than the measurements. Both simulations overpredict product concentrations. Product concentration predictions are worst near the combustor axis and the combustor

wall. Except for fuel prediction, hybrid pdf predictions

are closer to the experimental data, especially at larger

radii.

At the 78.5 cm station, measurements show some oxidant

and fuel concentration. This is despite the high temperatures that must exist. The pdf simulation overpredicts

oxygen and fuel. The conventional simulation predicts a

very slight amount of fuel near the combustor axis and some

II] I F
193

0.10 0.10 ,

,,

I

I

0.08 , , _ , i _ i

,...................... _....................... .............................. 0.08 ....,--.o...._............................................... _..............................................

0.06 E 0.06 ,oj !i i _ :

i illiiiiiiill ...... ?.............................

0

0.00 0.00

0.0 0.1 0.2 0.3 0.4 0.5 0.00 0.05 0.10 0.15 0.20 0.25

FUEL MOLE FRACTION OXYGEN MOLE FRACTION

0.10

o

o',

! I

_ ............. a...................... ; .......... : ..................... L

g o.o6 i : o;.......

f

rP 0.04 ....................... _....................... _....................... ' ....................... LO,..], ..........

0.02 ....................................................... 0 .........

!

0,00 0

0.00 0.04 0.08 0.12 0.1( 0.20

WATER MOLE FRACTION

0.10

I

o ,,

I

0.08 .............................................................................. O"'r ........................

I

o:

I

I

_" 0.06 ........................ _ ................. _......... O_ .....................

0

C]

;_ o.o4 i :!

i ,

0.02 ............................................. i 0 _ ..................... .

0 I

0.00 , C_

0.00 0.03 0.06 0.09 0.12

CARBON DIOXIDE MOLE FRACTION

Figure 5.76 Species Mole Fractions at 63.2 cm
194

0.10

0.08

E 0.06

El

<

0.04

0.02

0.00

0.0

I--I b-,,oeI

........................................................ I iexpenment- .... o pc.-- ...............

i

0.1 0.2 0.3 0.4 0.5

FUEL MOLE FRACTION

0.10

0.08

0.04

0.02

0.00

P;

t

i

t

0.00 0.05 0.10 0.15 0.20 0.25

OXYGEN MOLE FRACTION

0.10

0.08

E 0.06

El

<:

E: 0.O4

0.02

0.00

0.00

I i i

1 ................t...................... _...................... !...................o :

/

i ! i i/

i : 0 ,

t

...................... i....................... _...................... _......................... -O ---_.....

7

i _ O ,

i i o 1

0.04 0.08 0.12 0.16 0.20

WATER MOLE FRACTION

i

i

O ,

i

o.o6............................................................................. o :

El

<

_ 0.04

0.02

0.10 ,,

i

o I

_ u

0.0_ ................................................................................ _....',.........................

0

........................... i......................................... o ..........................................

i o

i io

i

0,00 O

0.00 0.03 0.06 0.09 0.12

CARBON DIOXIDE MOLE FRACTION

Figure 5.77 Species Mole Fractions at 78.5 cm

II1 I i:
195

oxygen towards the combustor wall. Predicted product concentrations are nearly uniform across the combustor at this

axial location. The experimental data shows considerable

variability. Carbon dioxide concentrations vary more than

the water. Product concentrations are overpredicted. The

hybrid pdf simulation best predicts product concentration.

Experimental data showed nearly complete reaction at an

axial distance of 137.5 cm.

combustor was not simulated.

50% longer than simulated.

The complete length of the

The actual combustor was about

In the pdf simulations, the much

smaller distances in the radial direction cause the number

of particles exchanged in the radial direction to vastly

outnumber the number of particles exchanged in the axial direction. This hampers information transfer in the axial

direction. Increasing the number of particles increases the

number of particles exchanged axially, but is more computationally expensive. The number of particles used in this

_alculation was 250 particles per cell. Using low numbers of

particles resulted in a gradual prediction of flame or

combustor blow-off. Using fewer particles caused larger

temperature oscillations in the iteration process.

Particles from three different computational cells for

an axial distance of 9.18 cm from the combustor step are

shown in figure 5.78. The particles are plotted as a function of mixture fraction and temperature, where mixture

fraction is the sum of the mass fraction of fuel in all
196

w

n1400

UJ

o.

LUI--

2500

195o .................................... i.............................................................................

x=9.18 cm.

r=1.35 cm.

2500

LU

n"

1950 .................................................................................................................

X=9.18 cm.

r=2.78cm.

.................................................................................................................. _ 1400 ...................................................................................................................

850 ...................................... ;..........................................................................

30O

0.0

LU

85O

300

0.2 0.4 0.6 0.0 0.2 0.4

MIXTURE FRACTION MIXTURE FRACTION

0.6

25OO

1950

uJ

1400

o,.

uJIx=9.18 cm.

i r=5.62 crn.

850 ..................................... _.........................................................................

3O0

0.0 0.2 0.4 0.6

MIXTURE FRACTION

Figure 5.78 Particle plots at an axial distance of 9.18 cm.

IF!!]
197

species. Pure mixing with no combustion is represented by a

line between mixture fraction of zero and temperature of 589

K (air) to a mixture fraction of 1.0 at a temperature of 300

K, which is pure fuel. The line would have some slight

curvature due to variable Cp. Combustion is represented by

an elevation off this line. Complete combustion is represented by a line from a mixture fraction of zero and temperature of 589 K (pure air) to a stoichiometric mixture fraction of 0.055 at 2500 K and hence to a mixture fraction of

1.0 at a temperature of 300 K.

The particle plots show that a range of mixtures is

being simulated in various cells. The particle plot for the

cell closest to the fuel inlet (r=1.35) shows various combinations of air and fuel with very little combustion. The

maximum mixture fraction in this cell is 0.55. No pure fuel

particles are shown. The plot for a radial location of,2.78

cm shows very lean, but complete combustion. The particle

plot for a radial distance of 5.62 cm shows slightly richer

combustion and higher temperatures.

Particle plots for an axial distance of 39.6 cm are

shown in figure 5.79. These particle plots show much higher

temperatures. As before, the level of combustion increases

with distance from the combustor axis. The cell with a

radial distance of 1.35 cm displays the largest number of

unburned and partially burned particles. The cell at a

radial distance of 5.62 cm displays only a few partially
198

2500

1950

LU

rr 1400

W

i--

850

3OO

e i

f.

f

...... ............................ i............................................................................

x=39.6 cm.

a r=1.35 cm.

i

_o t _

"""ii_'_£ *r.................... i.....................................................

_.t_"-

• Ig o'qb _a° __:

0.0 0,2 0.4

MIXTURE FRACTION

2500

• !

|

'1

1950 ......,_

| x---39.6 cm.

w

_: _ r=2.78 cm

14oo-".....................................

(2.

LU

I--

o •

3O0

0.6 0.0 0.2 0,4

MIXTUREFRACTION

0.6

2500

• 1950

IJJ

1400

W

I,--

850

300

t

i " i !

x--39.6 cm.

r=5.62 cm.

.... | ...........................................................................................................

..................... • .............. <............................................................................

0.0 0.2 0.4 0.6

MIXTURE FRACTION

Figure 5.79 Particle plots at an axial distance of 39.6 cm.

I[I I !:
199

burned particles. The reaction model implemented in the

hybrid pdf model gives partially combusted particles for low

temperatures, but most high temperature particles are almost

fully burnt.

The use of a finite-rate reaction scheme here has allowed significant overlap of fuel and oxidant concentrations. The assumed shape pdf, equilibrium chemistry model

used by Smith and Smoot predicted little overlap of oxygen

and fuel concentration profiles. The finite-rate combustion

models used here underpredicted the carbon dioxide concentrations in the recirculation zone, but the temperature

predictions were much improved from their calculation. The

hybrid pdf initially overpredicted fuel concentration, but

predictions using the hybrid pdf were better than the baseline predictions at successive axial stations. Hybrid pdf

predictions of oxygen, carbon dioxide and water were almost

always superior to the baseline simulation.
20O

5.3.2 SWIRLING DIFFUSION FLAME

Attempts were made to simulate the combustor of Jones and

Wilhelmi nl. This is a newer set of experimental data with

LDA measurements of axial, radial, and tangential velocity,

temperature and species data at numerous axial stations

including data very close to the injector face. This combustor flow had a 45 degree gas turbine type swirler, was

I00 mm in diameter and was 300 mm long. The fuel was gaseous propane introduced by a central conical annular injector. The axial velocity profile was completely positive by

100 mm, but the large centerline defect was still apparent

at 289 mm, near the combustor exit. All SIMPLE calculations

showed a nearly uniform axial velocity at this position.

This axial velocity defect may be caused by unsteady effects. Also, there was some ambiguity in the initial axial

velocity measurements. The first measurements showed two

axial velocity peaks with a single radial velocity peak.

This implied that the secondary peak was from the recirculating flow. This recirculating flow was stronger than any

calculated. The second velocity peak may have something to

do with a quarl at the top of the inlet swirler. The simulations did not model the quarl. At the second measurement

station, the axial velocity exhibited strong reverse flow

next to the fuel and air inlet flow, tapering off towards

the centerline. All simulations attempted calculated the

largest negative axial velocity along the centerline. The

i_!1i
201

experimental flow along the combustor centerline was quiescent. The high reverse flow next to the inlet flow and

adjacent quiescent flow seems more characteristic of flow

being entrained into a jet.

It was decided to study a similar combustor with weaker

swirl, as it is known that the k-£ turbulence model poorly

predicts areas of high swirl. A similar combustor is that

of Jones and Tober n2. This combustor has a thirty degree

swirler, with no quarl. The fuel is propane which enters

through a conical annular nozzle with a cone angle of 90

degrees just inside of the air swirler. The internal diameter and length of the combustor are 196 mm and 700 mm,

respectively. The inner and outer diameters of the swirler

are 23.5 mm and 42 mm. The combustor was run at two air to

fuel ratios. The case studied in this work had an air flow

of 26.7 g/s and fuel flow of 1.37 g/s.

A 70 by 55 grid rectilinear grid system was used. Grid

stretching was used to concentrate grid points near the

inlet fuel injector. The axial inlet velocity was calculated by projecting back the velocity data from the x/D=.l data

using the flow angle from the peak axial and radial velocities, and then scaling the profile to give the correct mass

flow. The radial velocity inlet data was found by doing numerous runs until a good match was obtained for the experimental data at the first couple of measurement stations.

The hybrid pdf simulation used a stronger inlet radial velocity.
202

The flowfield is characterized by two hot recirculating

flows. One recirculating flow is a toroidal vortex adjacent

to the combustor walls. This recirculation zone extends to

about x/D=0.5. The second recirculating flow occurs along

the combustor axis and is referred to as a Central Recirculating Zone or CRZ. This recirculation zone is characteristically caused by high inlet swirl. Based on a contour plot

in Jones and Tober, the CRZ appears to extend to x/D=l.6,

and takes up half the diameter of the combustor at x/D=0.8.

The combustor was simulated out to x/D=2.5. The two recirculation zones are separated by higher velocity flow which

is at an angle of about 30 degrees to the combustor axis.

Initially, this jet-like flow is largely composed of air and

fuel from the combustor inlet. This flow is typically fuel

rich at smaller radii and oxygen rich at larger radii.

Predictions and measurements of velocity and temperature at the x/D=.l location are shown in figure 5.80. The

calculated axial velocity peak for the flow separating _he

two recirculation zones is well predicted. This was one of

the criteria for selecting inlet conditions. The width of

the peak is overpredicted. The pdf simulation better predicts the lower portion of the velocity peak. The experimental data exhibits nearly constant, low velocity flow

across the CRZ. Both calculations show the magnitude of the

reverse-flow velocity increasing towards the combustor axis.

The predicted maximum reverse flow velocity in the CRZ is
203

0.10

0.08 ....................... : ...................

_ 0.06

o9

D_

D

0.04

0.02

0.00

-3O

l_I base_"ei

i .............. I..... O Iexperiment I/.....

io !

i o i"-.

-16 - 2 12 26 40

AXIAL VELOCITY (m/s)

0.10

0.08

0.06

u]

5

0.04

0.02

0.00

t / !

I!o i

............. o...-_..=.._.-.'._- - .-.=._.................................. i....................... 1

o i i i I

-5 0 5 10 15 20

RADIAL VELOCn'Y (m/s)

0.10

0.08

._. 0,06

u')

<

C_ 0.04

0.02

0.00

......................... i ............................. i............................. i .............................

................... i°o.............. i................. i.........................

i ° i

..................... i o i

i i",\o o o

;.................. _.0.._.::o si....................................

0 10 15 20

TANGENTIAL VELOCITY (m/s)

0.10 i

00

_ '

0.08

6

b

o,oo_.............. i................ _.... ,,

_ 0.04

0.02

0._

3OO

i 0 s

............................ i ............................ _)'"'_ ....... _"_""" .............................

i i o

ii ii0o \"'_

850 1400 1950 2500

TEMPERATURE K

Figure 5.80 Velocity and Temperature at x/D = 0.i
204

much higher than measured, which was characteristic of trial

calculations of Jones and Wilhelmi's combustor. The outer

recirculation zone has rather low velocity which is well

predicted by the simulations. Near the outer combustor

wall, calculations exhibited slightly larger negative axial

velocity than measured.

Radial velocity predictions for the outer recirculation

zone are very good for both simulations. Both simulations

underpredict the radial velocity peak and overpredict the

width of the higher velocity jet-like flow by almost a

factor of two. Measurements of radial velocity across the

CRZ are nearly zero. The simulations predict slightly

higher radial velocity flow. The radial velocity profiles

near the combustor axis are parabolic, and scale on the

cross-sectional size of the CRZ. The hybrid pdf radial

velocity profile for the jet-like inlet flow is slightly

shifted outward due to a larger inlet radial velocity profile. In this and other simulations of this flowfield, the

pdf simulations exhibited lower radial displacement of the

inlet flow with axial distance.

The predicted tangential velocity peaks are under-predicted and are highly diffusive. The experimental measurements show very low tangential velocity in the CRZ. The

simulations show nearly constant angular rotation across the

CRZ and up to the location of peak tangential velocity.

Tangential velocity in the outer recirculation zone is well

!F[_i
2O5

predicted by both simulations. The simulations do not

predict an experimentally measured depression in tangential

velocity at a radius of 7 cm.

The numerical simulations correctly predict high temperatures in both of the recirculating zones. However, these

temperatures are over-predicted by hundreds of degrees. The

high temperature predictions are the result of assuming i00

% combustion efficiency. Temperature predictions would improve if a lower heat release were used. The pdf module

used here uses thermodynamic data of Gordon and McBride and

assumes 100 % efficiency. This model is more difficult to

manipulate to predict lower heat of reaction. The experimental data shows depressed temperatures near the outer

wall. The combustor wall was treated as adiabatic in the

both of the simulations as the pdf module does not have a

wall heat transfer model in it. Thus, the predictions do not

exhibit a drop in temperature near the wall. The large drop

in temperature between the two recirculation zones is predicted, but the minimum temperature predicted is hundreds of

degrees too high. Both simulations predict the same minimum

temperature.

Fuel, oxygen, carbon dioxide, and mixture fraction

concentrations at x/D=.2 are shown in figure 5.81. The

mixture fraction is a combination of the mass fractions of

hydrogen and carbon in all of the species. Other concentrations are in mole fractions. Mixture fraction, Carbon
206

0.10

0.08

0.04

0.02

0.00

0.00

0.10 I

0 D

baseline b _o

o v°_o.,_,,,1 i _ o.o,............................_._

i E ! _" o.oe _:

! o!

,o o,:, i 0.04............................i_.....o.....-.i..;.; .................

_ -. _ E i . 0 0

i o "-__o i o _ i i "'-,o o

................................ ..........................

_o

i ,i

i 'i

"__ o.oo

0.02 0.04 0.06 0.08 0.10 0.12 0.00 0.02 0.04 0.06 0.08

CO2 MOLE FRACTION MIXTURE FRACTION

o i i ii o

,oZ r

0.08 ..... 0.08

_006,, _ _'_o.o_ i........... i!'_:o..................i_..............._.......

o_"".._' _o_,_,_:_ °__'_° '___.ooo.........

0.02 0.02

o ! !

0.00 0.00

0.00 0.02 0.04 0.06 0.08 0.00 0.04 0.08 0.12 0.16 0.20

FUEL MOLE FRACTION OXYGEN MOLE FRACTION

Figure 5.81 Species concentrations at xlD=0.2

I_I1i
207

dioxide predictions and measurements agree in the outer

recirculation zone. The traditional type simulation overpredicts carbon dioxide concentration at other locations.

The conventional simulation overpredicts the minimum carbon

dioxide concentration by a factor of almost three. The

hybrid pdf simulation does a better job of predicting carbon

dioxide concentration. Both simulations show carbon dioxide

initially increasing with radial distance from centerline.

The experimental concentrations are constant in this area.

The mixture fraction and carbon dioxide profiles match in

the outer recirculation zone as the fuel is fully burned in

this region. The drop in mixture fraction at a radius of 4

cm isn't shown in the simulations. This is due to the

overprediction of carbon dioxide at this radial location.

The peak experimental mixture fraction occurs at a larger

radius than the peak for carbon dioxide. Thus, there must

be significant fuel at a radius of 3 cm. The simulations

predict significant fuel at this location. Unfortunately,

the predicted fuel concentration is not large enough to

predict the mixture fraction peak at the correct radial

location. Instead, the simulations predict mixture fraction

peaks at a lower radius where the fuel and carbon dioxide

concentrations overlap. The conventional simulation predicts the magnitude of peak mixture fraction. Unfortunately, this is caused by an over prediction of carbon dioxide.

Curiously, the conventional simulation predicts peak fuel
2O8

concentration in the same region as peak carbon dioxide.

The oxygen concentration predictions show the same trends as

the experimental data. The predicted concentration of

oxygen in the outer recirculation zone is excellent, as was

other concentrations. At a radius of 4 cm, a peak is shown

due to oxygen from the inlet flow. At 2.5 cm, the oxygen

concentration is at a minimum due to large concentrations of

carbon dioxide and fuel. The conventional simulation underpredicts oxygen, especially across the inlet flow. The pdf

simulation best predicts the oxygen concentration. The

radial location of peak oxygen concentration is well predicted. The hybrid pdf predicted peak lies midway between

the experimentally measured peak and conventionally predicted peak. The low oxygen concentration inside the CRZ is

best predicted by the pdf simulation, but the width of this

part of the profile is predicted by the conventional simulation.

Temperatures corresponding to the above concentrations,

are shown in figure 5.82. Temperature predictions in the

recirculation zones are high. The temperature predictions

across the outer recirculation zone are nearly constant, as

were the species concentrations. The minimum temperature

predicted by the pdf simulation is midway between the experimental minimum and conventionally predicted minimum. At

the previous axial location, the simulations predicted the

same minimum temperatures. The radial location of minimum

If!I !_
209

0.10

0.08

E 0.06

o_

.-j

E3

<

n- 0.04

0.02

0.00

300

0 o i

0 i

• I

io

I

baseline O

pdf O ,'

• I

expedrnent O j

....................................... ii....................................... .:.Oo .......................... ii_//......................................

i 8.-"

_ [

i o -, \

i oo

850 1400 195o 2500

TEMPERATURE K

Figure 5.82 Temperature at x/D = .2

temperature is the same as that for minimum carbon dioxide

and maximum oxygen concentration. The hybrid pdf simulation

best predicts temperature in the CRZ.

Velocity and temperature at x/D=.3 are shown in figure

5.83. Both simulations predicted nearly identical axial

velocity profiles. Predicted positive axial velocity profiles are much wider than the experimental profile. Peak

axial velocity is underpredicted, although slightly beyond

the peak, the predictions are excellent. The size of the

CRZ is underpredicted. The maximum reverse velocity is

over-predicted in the CRZ, although it is much improved from

the x/D=.l location. Again, the experimental measurements
210

0.10

0.08

0.06

5

0.04

0.02

0.10

0.08

0.06

_ 0.04

0.02

0.00

i 0"-

! i ! , 8

i ol ° i --"

i 8 ..--"

{o >Ji

,d I--:- I_ I

/ill o I.xp_n.,-,I

0

-10 0 10 20 3O

AXIAL VELOCITY (m/s)

. [ " . :

i', i

i ".o i

i "o. i

iI _ i;,,_ J° _

...................... i-----.--.-._..-!-....o ....... N.'...'.......................................

................. ! , ........... _......................J.............................................

0.00

-20 -5 0 5 10 15 20

RADIAL VELOCITY (m/s)

o oo.°l°l oI

..................... _......................................... '6 ...................... ....................... ............................ .............................iO ...................... ,': ..........................

i i oi _ '_o ,'

! o i '._ o ,'

: o X i i ,5 '

........................ _.............oO.i ......... , ........... £ 0.08 i ....................... _.............. _" 0.08 , o _ ,

0 i '

-' 80 __ _ q o,'

i i '_ 0 m _ i • Oo:" •,_o.-..._. .......... -.................. _ ............................ ............................. .................. ........................

o -" _ i oo "', !o '-/" i

// ! i _ o

0.02 o i ,,"_:; ..............4....................... -................ 0.02I......................... i................................ o ............_. -/

t '

i i o : //

• | ! : I

0.00 0.00 L 0

0 2 4 6 8 10 300 850 1400 1950 2500

TANGENTIAL VELOCITY (m/s) TEMPERATURE K

Figure 5.83 Velocity and Temperature at x/D=0.3

![I| i:
211

show nearly constant, low axial velocity across the CRZ.

The experimental radial velocity profile shows a narrow high

velocity peak. Predicted peak profiles are blunter and

wider. The general shape of the predicted profiles is very

similar, only differing in magnitude. The hybrid pdf simulation predict a velocity profile about 30 % lower than the

conventional simulation.

The experimental tangential velocity profile shows two

peaks. The first peak is due to the swirling inlet swirling

flow. The second velocity peak is adjacent to the outside

combustor wall. Tangential velocities in this region are 4-

6 m/s. At x/D=0.1, tangential velocities were less than 5

m/s near the outer wall. The increase in tangential velocity

at larger radii is due to spreading of the peak tangential

velocity which has decreased by about a factor of two from

the x/D=0.1 location. The numerical simulations only predict the first velocity peak. Examination of the axial and

radial velocity profiles does not give an apparent reason

for the second outer tangential veiocity peak or a tangential velocity defect in the profile.

The shape of the predicted temperature profiles is

qualitatively correct. Predicted temperatures are too high.

The conventional simulation predicts the minimum temperature

at about the correct radial location. The hybrid pdf profile shows the radial location of minimum temperature about

one cm less than was measured. Both simulations signifi-
212

cantly underpredict the radial position of peak temperature.

The pdf simulation predicts lower temperatures than the

conventional simulation. Experimental temperature gradients

are much larger than predicted.

Velocity and temperature results at x/D=.4 are shown in

figure 5.84. The agreement between the simulations is

better than the agreement between the predictions and measurements. The simulations predict the same magnitude of

peak axial velocity, but the hybrid pdf simulation predicts

a smaller radial displacement of the velocity profiles.

Predicted profile peaks are again wider and blunter than the

experimental peaks. Predicted radial velocities are almost

identical out to a radius of 4 cm. Beyond 4 cm, the hybrid

pdf simulation predicts lower radial velocities. The velocity predictions within the CRZ are closer to the experimental

measurements than at previous stations. The width of the

CRZ is again underpredicted. The maximum reverse flow

velocity in the CRZ is decreasing. The experimental velocity

gradient between the reverse flow in the CRZ and adjoining

flow is greater than predicted. Predicted velocity profiles

show diffusive, gradual changes. The comparison between

predicted and measured tangential velocities is very good

out to a radius of 3 cm. Beyond 4 cm the hybrid pdf simulation predicts higher tangential velocity at a lower radii,

which agrees with the data. This is due to the reduced

radial expansion predicted in the hybrid pdf simulation.

!!II i
213

0.10

0.08

0.06

(/)

o

0.04

0,02

0.00

-I0

-,--"! i I olo

_-.._ i

"_..o i

__ 0.08

0 o

i ! " q_

i i , :0 _ 0.06

"i _ _ p o °° o

............. o ........... !.................... "i............................ _............................ <_

o i /" i _ 0.04

°-= j_,°,ool

/o i ° io_0,o0o, l 00_

0,00

o 10 20 30 -5

AXIAL VELOCITY (m/s)

i ";'. 0

......................_......................._.. .._..._........ :.........0 "0...........................

:: i "i 0 i i , o i _i__ o

! _ ._ o o

i o _ "'i

.......... :i ......... !....................i....... _.............

0 5 10 15

RADIAL VELOCITY (m/s)

2o

0.10

: _ i O0

,,o

_" i /!/°°

i o/"

o.o4 ........................... _--7""_ ............................. ............................

00, '" i........................ ...........................

0.00

0 4 6 8

TANGENTIAL VELOCITY (m/s)

0.10

0.o8 .......... 7............................. i'°° .................. :..........................

i i ° : \

............... i............. i ....°i i i.....

,0.02 ....

0.00 O

300 850 1400 1950 2500

TEMPERATURE K

Figure 5.84 Velocity and Temperature at x/D=0.4
214

Again, experimental data shows a second tangential velocity

peak, which isn't predicted by the calculations. Temperature predictions are high. Temperatures continue to become

more uniform with axial distance. The pdf simulation predicts lower temperatures than the conventional simulation,

but the conventional simulation better predicts the radial

location of peak temperature.

Concentrations for x/D=.4 are shown in figure 5.85.

The simulations predict relatively constant carbon dioxide

concentration across the combustor. The experimental measurements predict a 30 % drop in concentration at a radius

of 7 cm. The conventional simulation predicts peak carbon

dioxide at 5 cm. Carbon dioxide is underpredicted at smaller diameters and overpredicted at larger diameters. The

experimental mixture fraction shows only a slight dip at 7

cm, unlike the carbon dioxide profile. Also, the mixture

fraction shows a significant peak at 5.5 cm while the peak

carbon dioxide is at 3 cm. This implies that there is

significant fuel or partially burned fuel between these

points. The fuel mole fraction should be on the order of

0.01. The conventional simulation predicts almost no fuel

at this axial station. The pdf simulation shows a maximum

mole fraction of 0.0025 around 5 cm.

The conventional simulation does a good job of predicting oxygen concentration out to 6 cm. The peak experimental oxygen concentration at 7 cm isn't predicted in the

IllIi
215

o' ! _'

! i ; : _ I

! i i :

0 !' _ 0.06 i i _ o

v i_ ii i°bi _ " i i ",, o_

] _ _ ]

i _. ,', ol 5_ i i ' o

............................ 4 ............................. i ................ _-.....................................

er" 0.04 ............... n" 0.04 _ _

i ! , ! i

0.02 ......................................... 0.02 ............. i.................................... , o .................

i i i o i o

ooo

iiI2__ ooo o

0.00 0.02 0.04 0.06 0.08 0.I0 0.12 0.00 0,02 0.04 0.06 0,08

002 MOLE FRACTION MIXTURE FRACTION

oo i i ooi _o i i

'0

i

i '0

0.08 ......................................................... _.......................................................... 0.08 ....................... ÷..... "_..0.? ....................... _..................... ;.....................

, O i i

i !

_" 0.o6 ',,........................................................ _ _................ _" o.o_ ........ d oi° ,i................................. ..................

0.04 1.................................................................................................................... o.0_

.......•,-,_..-o _.................. _.............. i.............. i.................. °°'................................................

....................... °°',

000 000 _ :' ! _ ;

0.00 0.02 0.04 0.06 0.08 0.00 0.04 0.08 0.12 0.16

FUEL MOLE FRACTION OXYGEN MOLE FRACTION

0.20

Figure 5.85 Species concentrations at x/D=0.4
216

simulations. The pdf simulation predicts the highest oxygen

concentration at this station. Neither simulation is superior in the prediction of oxygen. Both simulations overpredict combustion. This results in more expansion at smaller

axial and radial displacement from the combustor inlet.

Velocities at x/D=.6 are shown in figure 5.86. Negative axial velocity is measured over 60 % of the combustor

inner diameter at this station. The simulations underpredict

the size of this flow. Thus, the experimental velocity

profiles show higher positive axial flow at greater radii.

The back-flow velocity along the combustor axis is well

predicted. The hybrid pdf simulation shows lower reverse

velocity near the combustor axis, which may be due to increased combustion at smaller radii. The conventional simulation best predicts the axial velocity profile at this

location. Experimental radial velocities are all positive.

Simulations predict negative radial velocity for at least

half of the inner diameter. The pdf simulation shows higher

radial velocity at larger radii, unlike previous stations.

Both simulations predict very similar tangential velocity

profiles at this station. Predicted tangential velocities

are nearly constant from 3 cm to 9 cm. The experimental

data shows a tangential velocity peak at 5 cm, which is

significantly underpredicted in the simulations.

Species concentration and temperature at x/d=.6 are

shown in figure 5.87. The experimental carbon dioxide

IEl11
217

0.10

0.08

0.06

GO

5

0,04

0.02

0.00

-20

0

i ! o i,

i ,i

Oi ,

i ol ,'i

0

i o ,I i

...................... +............... _,.<..-! ..................... i...............................................

j , i

-10 0 10 20

AXIAL VELOCITY (m/s)

0.1o

0.08

"_ 0.06

_ 0.04

0.02

0.00

3O -5

', o i

/o

I ,i' O

Ii " :

.......................... i........O ............................................. i............................

o

......... o ............ .....................!................

o !

i

5 10 15

RADIAL VELOCITY (m/s)

o.lo _

i o

!

008 ............................ _.......... _............_............................. ,.............................

b

if 0

/ o

_ o.os ..................... ; o ..................................

GO

D . i\ 0

0.04 ............................ }...................... 10 ......................... ..............................

E i

_ /

_ z

0.02 ........................................................

0.00

2 4 6 8

TANGENTIAL VELOCITY (m/s)

Figure 5.86 Velocities at x/D=0.6
218

0.10

0.08

_. 0,06

CO

0.04

0.02

0.00

0.00

1

: l

O I expedmen! I

........... i!............... i _................ _.................. _"_.......... ii! tIII0{'

---0 ..........

0

0

0

0

O

•-O ..........

O -

O-----

0.02 0.04 0.06 0.08 0.10 0,12

CO2 MOIP FRACTION

0.10

0.08

_E 0.06

O3

0

.......... 0 .............................

O

O

O

O

i ,

i

: I

t

I

0,02 .................................................... :! ....... ,' 0 ...................

i i 'I

I

i ! ',

0.00 0

0.00 0.02 0.04 0.06 0.08

MIXTURE FRACTION

0.10

0.08

0.06

CO

_=

0.04

0.02

_0j 0%:' ; 0.10

...................... ";' ....................... _....................... -_.............. "" : .................. 0.08

.............. _ii................. i i:................... ::i !ioo!o......... i 'i'_../" __o.o6 Co_

..................................................................... 0 .............................. _ 0.04

........0 ............... "/'_' v................ _......................._....................... i......................

O /i ' i

O i _

O ' i

:1 ; {

..... 0 .................. _................. _ .......................

i io

.................. i.......................................... io ............ ,,........

i o

O

740 1180 1620 2060 2500

TEMPERATURE K

............. O

I

0.02 ] .... 0

i

I

0.00 I (:

0,00

il

i,

i,

i,

.... it..........................................

h

il

iJ

il

il

0.00

300 0.04 0.08 0.12 0.16 0.20

OXYGEN MOLE FRACTION

Figure 5.87 Species Concentration and Temperature at x/D=0.6
219

concentration is constant out to a radius of 8 cm. The

mixture fraction shows a slight peak at 8 cm, which implies

that there is still unreacted or partially reacted fuel at

this location. Insufficient fuel is predicted at this axial

station. The shape of the predicted mixture fraction profile is the same as the carbon dioxide profile. Measurements show a temperature peak and a drop in oxygen concentration at 8 cm, which suggests ongoing combustion at that

location. The simulations show increased oxygen concentration past 6 cm.

Velocity and temperature at x/D=0.8 are shown in figure

5.88. The experimental data show the CRZ is about the same

diameter as the previous station, but the reverse flow

velocity has slightly increased. The simulations show a

smaller diameter CRZ and reduced reverse flow velocity from

the previous station. Predicted radial velocities are all

negative at this station. The experimental data for radial

"velocity is rather sparse, but suggest a different flow

pattern. The discrepancy in prediction of peak tangential

velocity grows to 50 %. The second tangential velocity peak

isn't seen at this location. Measured tangential velocities

at this station are greater than at the previous axial

location at all radii. This implies there is some error in

tangential velocity. Predicted temperatures across the

combustor are nearly constant. The experimental data shows

a slight temperature peak at 9 cm. The pdf simulation shows
220

0.10

0.08

_ 0.06

cO

5

<

n- 0.04

0.02

0.00 _

-20 -10

i o

i _:/,' !

i ° /,,

// L--

........... ol _ ............1..... 1_ L

i I ' 1 0 lexperimentl

0.10

0,08

0.06

cO

C_

0.04

0.02

0.00

....................... r'O ............................ T......................................................

!0

I i

' i 0

10 20 30 -5 0 5 10

AXIAL VELOCITY (m/s) RADIAL VELOCITY (m/s)

15

0.10

0.08

0.06

03

_=

(3

0.04

0.02

0.00

0.10

i ',

I

............................ _............. +............ i.-..O ................... _............................ 0.08 ..................... _....................... ' ..................... ;.........0" _'l........

i _ i 0 0

i ' i o o

........................... i............ , i :, o ......... "_ o.oe i ..................................................... o

i ............ ,,,i..............o _ 0.o4 _ o

!_ i-,i o......................................... ............ _.............................. ........................... o ....................

ii'l

i ' ' 0 _ 0

....................... i.......*"o i............................................ o.o_ ......................................................... o.......... :

0

0.00 _ ; 0

0 2 4 6 300 740 1180 1620 2060

TANGENTIAL VELOCITY (m/s) TEMPERATURE K

2500

Figure 5.88 Velocity and Temperature at x/D=0.8

IIII
221

lowest temperature, closer to experimental data.

Mixture fraction isn't perfectly conserved in these

calculations. If everything were perfectly conserved in the

calculations, the predicted temperatures would be closer at

this location. Summing the conventional mixture fraction

flow across the combustor shows a slight increase in total

mixture fraction at this location. Summing the hybrid pdf

predictions of mixture fraction flow shows a slight decrease. Thus, some of the improvement seen between the

calculations should be attributed to error, not better

modeling.

Some cells on which to perform particle plots were

chosen for examination. Particle plots of mixture fraction

versus temperature for x/D of 0.0034 are shown in figure

5.89. The particle plot for r/D of 0.054 displays a wide

variety in mixture fraction and temperature. This cell is

within the CRZ. Most particles that appear in this cell are

fully reacted. The particle plot for r/D of 0.74 displays

three distinct particle groupings. One grouping is air with

a slight amount of unburned fuel. The second group is a

single particle at 1914 K and a mixture fraction of 0.14.

The third particle grouping is almost pure unburned fuel.

This particle plot is of a cell between the two recirculation cells. The particle plot for r/D of 0.198 is actually

of 250 particles clustered around 1940 K and a mixture fraction of 0.044. The fluid from this cell is within the outer
222

2700

2100

,,/

L,U

,re

'D

1500

UJ

IL

.......... '_........ T...................... i........................................................

\ i x/D=-.oos4 : ! r/D=.054

,i

900 .................................................................................._.....................

i

300 " '

0.00 0.25 0.50 0.75 1.00

MIXTURE FRACTION

ILl

n"

1500

UJ

g,,.

tU

I.-

2700

2100 .......................................................:........................................................

x/D=.O034

r/D=-.074

• " i S "!

300

0,00 0.25 0.50 0.75 1.00

MIXTURE FRACTION

2700

2100 ..............................................................................................................

,z ' x/D=.0034 i

ELl n" r/D=.198

1500 ..............................................................................................................

IJJ

0.

UJ

I--

900 i

3OO

0.00 0.25 0.50 0.75 '1.00

MIXTURE FRACTION

Figure 5.89 Particle plots at x/D=0.0034

I[[I]:
223

recirculation zone, next to the combustor wall. The fluid

is almost completely mixed. The fluid velocity in this cell

is very low. Information transfer in this zone of the

combustor is very slow in the hybrid pdf calculation.

A second set of particle plots for x/D of 0.136 is

shown in figure 5.90. All of these plots show considerable

variability. No particularly fuel rich particles are shown.

The plot for r/D=0.054 is within the CRZ. Most particles

are fully burned except for a smaller number centered around

900 K. The plot for r/D=0.074 is closer to the fuel stream.

The number of partially burned particles centered near 900 K

is larger than for the r/D=0.054 plot. There are also a few

more particles with higher mixture fraction. The plot for

r/D=.198 has the greatest number of particles with moderate

temperature and mixture fraction. These particle plots are

used to show variability within the fluid cells. Particles

can be plotted over one another. Some sort of statistical

analysis would probably help in interpreting this data.

However, much of the easier statistical analysis assumes

normal or Gaussian type distributions. The distributions

are obviously bimodal at low temperature in and near the low

temperature fuel-air stream.

Neither simulation captures the experimentally measured low fluid velocity of the CRZ near the combustor

inlet. The pdf simulation predicts slower combustion and

generally less radial displacement. The inlet radial
224

2700

21oo

%,

w

1500

o.

w

k9OO

300

J " .x/D=-.136

,I rlD=-,0S4

• !

• !

I'' !

' i i

j

0.00

2700

2100

%,

LL]

n"

_ 1500

LU

I--

900

3OO

i °° x,_-. 136

, riD=.074

,;&....................................................................

Oo

_._................ _.......................................................

I; i..

i

0.25 0.50 0.75 1.00 0.0O

MIXTURE FRACTION

0.25 0.50 0.75

MIXTURE FRACTION

1.00

l/

300

0100 0.25 0.50 0.75 1.00

MIXTURE FRACTION

Figure 5.90 Particle plots at x/D=0.136
225

velocities had to be increased so that the hybrid pdf results would compare to the conventional simulation. While

the lower rate of combustion predicted by the pdf simulation

predicted temperatures closer to the experimental data, the

larger radial displacement of fuel and its products by the

conventional simulation produces profiles that were closer

to the data in regards to shape. Both simulations use

hybrid differencing for calculation of coefficients. As the

flow is about 30 degrees to the grid (initially 45 degrees

for fuel), hybrid differencing can be expected to give

larger errors due to numerical diffusion, compared to higher

order schemes. The pdf scheme should be tested with higher

order schemes. Unfortunately these schemes will to have all

positive coefficients, if the same differencing scheme is

used to calculate the pdf particle number exchanges as

transporting a negative number of particles does not have a

physical meaning. A quick fix for this may be to use lower

order differencing for the transport of particles.
226

5.3.3 PREMIXED COMBUSTIONCASE

To explore the further applicability of the hybrid pdf

model, a lean premixed combustor case is considered. In a

premixed combustor, fuel and oxidant are mixed upstream of

the actual combustor. This results in uniform air to fuel

ratios. In normal combustors, there is a wide variation in

air to fuel mixture. Some areas are lean, some are rich,

and some are near stoichiometric. Burning near stoichiometric results in high temperature which markedly increases

nitrous oxide emissions. Lean premixed combustors do not

have this problem. Thus, lean premixed combustion is one

way in which to lower nitrous oxide emissions. Nitrous

oxides are one component of smog. Currently, some European

countries are proposing taxes based on nitrous oxide emissions of the aircraft flying throughtheir airspace. Since

the combustor inlet mixture is premixed and is combustible,

the location where combustion occurs must be fixed or anchored. Combustion can flash back into the inlet duct or

can be blown downstream. This is deleterious to performance

and safety, particularly in aircraft. Combustion is anchored by forming hot recirculation zones. One way to form

a recirculation zone is to put an obstruction, or flame

holder, in the combustor. Some examples of flame holders

are disks, disks with holes, v-gutters, and conical bluff

bodies. In the combustor studied here, a large increase in

flow area is used to form a recirculation zone. This con-

!_I11
227

figuration is also known as a step or dump combustor. The

basic flow configuration is the familiar flow over a step.

The combustor studied here is that of Gould, Stevenson_ and

Thompson In at Purdue University. Premixed fuel and air

goes through a converging nozzle, with an exit diameter of

76.2 mm, into the combustor chamber of 152.4 mm outer diameter. The inlet nozzle was contoured to provide a nearly

uniform velocity flow across the combustion chamber inlet

plane. There is a slight velocity defect due to the boundary layer. The inlet mixture is gaseous propane and air at

an equivalence ratio of 0.5. Velocities were measured using

Laser Doppler Velocimetry. Temperature was measured using

high speed thermometry. Various turbulent correlations were

also measured. A corrective lens was used to take optical

measurements. This lens degraded due to combustion, limiting the number of axial stations where data was taken.

Gould performed computational combustion simulations as

part of his dissertation study. A pressure based solver

with the eddy-breakup combustion model of Magnussen and

Hertajer _4 was used. In this combustion model, combustion

is proportional to turbulent eddy lifetime, and minimum

species concentration of fuel, oxygen or combustion product.

Combustion is initiated in this model by temporarily introducing product species at specified cells in the simulation.

Recently Magnussen _2s proposed an improved combustion model

called the Eddy Dissipation Concept combustion model. This
228

model was used in a calculation of Brum and Samuelson's

combustor 126, which was then compared to the predictions of

Nikjooy [28].

An inlet axial velocity of 22 m/s was specified in the

combustor calculations, except for the boundary layer. Zero

radial velocity was specified. There was no swirl component. The measured inlet turbulent kinetic energy was very

low, 0.0006"U 2, except for the boundary layer, which was

measured at i0.0 m2/s2. Combustion simulations did not

satisfactorily converge using such low inlet turbulence.

Values of 0.0054"U 2 were used, except for the boundary

layer. This allowed the simulations to converge much further. High inlet temperature was initially specified to

initiate combustion simulations. The simulations were then

restarted using the correct inlet temperature of 300 K. The

pdf simulation started from a converged SIMPLE calculation

restart file. As the current pdf model does not allow heat

transfer to the wall, adiabatic wall treatment was also used

in the conventional simulation. The heat of reaction term

used in the conventional simulation was adjusted to give the

same final temperature as the thermodynamic model used in

the pdf module for a single cell calculation.

Experimental values of velocity, temperature, and

various double and triple correlations of these quantities

were measured. The experimental data taken near the combustor axis shows radial velocities around -2.0 m/s. This is
229

physically unrealistic. Evidently, the thermometry apparatus affected radial velocity measurements. Radial velocities measured without the thermometry apparatus for the

isothermal case were reasonable. Predicted radial velocity

profiles do not agree with the measured profiles for the

burning case. There was significant anisotrophy in the

turbulent Reynold's stresses. Axial velocity and temperature at a location of one step height past the combustor

inlet are shown in figure 5.91. The axial velocity is well

predicted by the two different combustion model simulations.

The axial velocity in the main flow is slightly over-predicted. The experimental data shows a slightly sharper

velocity gradient between the cold core flow and the hot

J

recirculation zone. Temperature measurements in the recirculation zone are hundreds of degrees lower than predicted

temperatures. Low temperatures are measured near the combustor wall. The predictions assumed no heat transfer to

the walls. Both the conventional simulation and hybrid pdf

simulation predict complete combustion for much of the

recirculation zone. As velocity is lowest next to the

combustor wall, this portion of the flow has the longest

time for combustion to occur. With the high temperatures

measured at other axial stations, much of the recirculation

zone should be fully reacted. This researcher feels that

the peak experimental temperature at this axial location

should closely match the peak temperatures found at other
230

0.08

•.._ ! .....PC" 1

0.06 .......... '_ .......................... }.......................... ............ I o texperimentl......

ol "'_-. i !

i O i "'--!_. i i

1 i o "---!__ i

.................................................................._ _-_...................

0.02 ....................................................................................................................................... O ................

i o

E _

0.00 C

v

03

Z) 0.04

5

<

rr

-5 0 5 10 15 20

AXIAL VELOCITY (m/s)

25

E

co

0.04

(:3

<

rr

0.08

O

O

O

O

0.06 ............................................................... :........................... 0 .............. i.............................................

O

!-

0.02 ).................. _.......................... _..............................................................................................

)

0.00 >

200 500 800 1100 1400 1700

TEMPERATURE K

Figure 5.91 Axial Velocity and Temperature at x/H=l.0

lil I F
231

axial stations. Allowing some heat transfer in the calculations would improve the temperature predictions along the

wall, but shouldn't produce a 400 degree drop in temperature

across the whole recirculation zone. Experimental velocity

data was only taken out to 80 % of the combustor diameter.

Predictions of axial velocity in the recirculation zone

appear quite good, despite the large temperature discrepancy

in the recirculation zone. This implies that temperature

measurement in the recirculation zone at this station is in

error or, velocity is not as strongly coupled to density as

it is in the core flow.

The results at an axial location of 3 step heights are

shown in figure 5.92. The axial velocity predictions and

measurements are very close. The pdf predictions show a

slightly higher axial velocity outside of the cold central

core. The hybrid pdf simulation shows a slightly faster

combustion progress into the core flow. High temperature is

predicted at slightly lower radii in the pdf simulation.

The pdf simulation predicts a wider blunter temperature peak

than the conventional simulation. The conventional simulation predicts slightly higher temperature than measured in

the hot recirculation zone.

At 5 steps heights past the inlet, the axial velocity

predictions, shown in figure 5.93, are very similar. The

conventional simulation comes closest to the measurements in

the outer recirculation zone, but axial velocity is
232

0.08

0.06

E

09

0.o4

a

<

n0.02

0.00

I

i "b_. ! ! | O experiment

i O 4.

i o'-. i i i

! c

IC

: ! i

....................................................................................................... _.......................... i,,,,,0

(

-5 0 5 10 15 20 25

AXIAL VELOCITY (m/s)

0.08

0

O

i o I

0.06 ............................................................................................... i ............0 0 ........ ,'.......................

E

0.O4

m

<

n"

0.02

O

O

0.00

200 500 800 1100 1400 1700

TEMPERATURE K

Figure 5.92 Axial Velocity and Temperature at x/H=3.0

it!11
233

0.08 i

"- ! i

i\",i ! 1.....Ip_, I......

0.06 .................................. i

.......................... _.................i "-,

0.04 ...................... ' : _ _ ",,

0.02 ................. • i........................ i......................... i........................i i

i i i :, i

ooo_ _o 15 ' o 2_

-5 0 5

AXIAL VELOCITY (m/s)

0.08

0.06

0.02

0.00

2OO 500 8(30 1100 1400 1700

TEMPERATURE K

Figure 5.93 Axial Velocity and Temperature at x/H=5.0
234

bracketed by the predictions in the inner recirculation

zone. The small negative axial velocities shown near the

wall in the conventional simulation mean that the recirculation bubble is almost closed for that simulation. The pdf

simulation shows slightly positive axial velocity near the

wall. The axial velocity near the combustor axis is underpredicted by both the simulations. Previously velocity was

slightly overpredicted. In this same area, predictions and

measured data for temperature are very close. Beyond the

cold core flow, the hybrid pdf simulation noticeably predicts higher velocity and temperature than the both the

conventional simulation and experimental measurements. Up

to a radius of 4.5 cm, the conventional simulation predicts

temperature almost perfectly. At the previous station,

temperature was slightly overpredicted. Beyond a radius of

4.5 cm, temperature measurements are bracketed on the low

side by the conventional simulation and on the high side by

•the hybrid pdf simulation, except near the wall. Qualitatively, the shape of the pdf simulation temperature profile

is closest to the measured temperature profile at this axial

location.

At the last measurement station of 12 step heights past

the inlet, both numerical simulations deteriorate, as shown

in figure 5.94. The difference between the simulations and

experimental data is largest at this station. In most

experimental combustors the last measurement station is

:!rlIIi-
235

0.08

0.06

E

v

Or)

O.O4

<

n0.02

0.00

....... 0 :'_! ......

! o ",

] : \ :

:: 0 X

i o ,, _

................................ ......... x......._.,.;. ...............................

I Ibaselin' ! o,,

.......... I..... o I:'0°_°t............ i....................... I ..................

0 5 10 15. 20 25

AXIAL VELOCITY (m/s)

0.08

0.06

E

v

03

O.04

a

<

m

0.02

0.00

200 500 800 1100 1400 1700

TEMPERATURE K

Figure 5.94 Axial Velocity and Temperature at x/H=12.0
236

usually at the completion of combustion and mixing, and exit

profiles are fairly uniform. In this experiment, a large

fraction of the cold inner flow is still unreacted. Thus,

there are significant gradients in combustor profiles. Both

predictions show a higher axial velocity along the centerline than measured, while temperature is underpredicted.

This is somewhat of an anomaly. Beyond a radius of 2 cm,

the conventional simulation underpredicts axial velocity,

while the pdf simulation continues to overpredict it. The

experimental measurements show an almost linear increase in

temperature off the combustor axis. At a radius of two

centimeters, both temperature predictions are off by 300

degrees K. The hybrid pdf simulation best predicts the high

temperature zone at this temperature. Gould noted similar

deficiencies in his predictions and theorized that it was

the fault of turbulence model. At previous stations, the

predictions in and near the low turbulence level cold core

flow were excellent. At this last station, the cold flow is

hotter and more turbulent. The results at this station

could be more closely predicted by increasing inlet turbulence. Unfortunately, this would work at the detriment of

previous predictions. With the extremely low inlet turbulence, turbulence modeling may be more in the realm of

turbulent transition. The turbulence model used in these

simulations assumes high levels of turbulence, that is, it

is a high-Reynolds number turbulence model. It would be

IF!li
237

interesting to try a low Reynold's number turbulence model

in another simulation.

The hybrid pdf temperature simulation exhibited sharper temperature gradients between the cold inner flow and hot

r

outer flow, and higher temperature, particularly at the last

station. At the last station, the pdf simulation best predicts temperature.
238

5.3.4 SWIRLING HYDROGENDIFFUSION FLAME

Recently, a group of researchers at the University of

Dayton, Research Institute, have taken benchmark data for

model validation. This includes non-swirling and swirling

air jet flows 127, methane jet flames 128, and hydrogen jet

diffusion flames rig. The combustion data for the burning

hydrogen is the most complete. Three components of velocity

were measured along with many turbulence quantities using

LDV equipment. Temperature was measured using Coherent

Anti-Stokes Raman Spectroscopy (CARS). Typically 500 temperature measurements were taken at each location. These

measurements were processed to get rms temperature fluctuations and temperature pdfs. The CARS method is still being

developed and the data is ideal for validating the composition pdf method employed here.

The experimental system consists of a central fuel tube

of 9.45 inner diameter, a concentric annular air tube of

26.92 inner diameter, and external non-swirling air supplied

to a vertical 150 by 150 mm semi-rectangular test section

which is 486 mm long. Various combinations of velocity and

swirl were tested. The case studied here is case no. 3

which had a fuel jet bulk velocity of I00 m/s, annulus air

bulk velocity of 20 m/s with thirty degree swirl, and external air bulk velocity of 4 m/s. Radial data was taken at

1.5, i0, 25, 50, 150, 250 mm from the inlet. Additional

combustor centerline data were also taken. Data at the

IF!11
239

first measurement plane were used for inlet data in the

simulations. The mass" fraction of species were inferred

from the geometry and temperature measurements. The high

temperature combustion interfaGe between the fuel and air

jets provided an ignition source. A one step chemistry

submodel was employed. This mechanism was inferred from

several larger reaction mechanisms for heavier fuels which

included a hydrogen-oxygen to water reaction step. The

activation energy used was 40 Kcal/mole. The exponent used

for the fuel concentration was 0.25, and the oxygen exponent

was 1.5. The frequency number had to be increased to 2.7 x

1013 to simulate continuing combustion in the calculations.

Baseline calculations were done with a coarse grid of

40 by 30 and a fine grid of 70 by 60. These simulations

predicted nearly identical axial velocity and temperature

profiles. The pdf calculations were performed using the

same grids. An adiabatic wall boundary condition was imposed at a radius of 75 mm.

Temperature, turbulent kinetic energy, axial and tangential velocity at an axial distance of 25 mm from the

inlet plane are shown in figure 5.95. Temperature at this

plane is best predicted by the baseline coarse grid simulation. The coarse grid temperature profiles are largely the

same as the fine grid profiles except the coarser grids

truncate the temperature peaks. The peak temperature is
240

21001800f _ ,_" t 5OO

] t ! -- b,lt+lli_I41_I3(I

;_ .... mm,.+_ I

, I-'- _"= /

_+ imo ,z

t< _ 3oo

++'I .... IL) +,+"=

P-. 9oo ......... l

I.-- 10o

_"< i i

3OO 0

0.000 0.005 0.010 0.015 0.020 0.025 0.030

RADIUS (m)

I/ " +/I._,I ............................ + I-- _+_= I

+,/L'i,_ +Io_ I

"_o+,t_ i I [] =r,,*., I

...... ##:ti_ ............ _......................l+ -+.,l...........

,, i i

.................... [................ i................:............

0.000 0.005 0.010 0.01S 0.020 0.025 0.030

RADIUS (m)

150 10

!

8 ...............................I":"_":"_........................... ;...................:..................

+" /-" + i +

/; + ++ +

+,oo + + .....,," +........ +........

£ ,/ i ;

i _ ./j + +

+ ,_ . .+................... +....

.....................i{............. ..................i ...... i

5O

, + + _ _ . -i,._, i

o a \'_

0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.00_000.005 0.010 0.015 0.020 0.025 0.030

RADIUS (m} o CP RADIUS (m)

Figure 5.95 Temperature, turbulent kinetic energy, axial and

tangential velocity at x = 25 mm.

'[!l i
241

slightly overpredicted by the baseline fine grid calculation. The pdf simulations underpredict temperature. Efforts to increase the coarse grid peak temperature by increasing reaction rate and maximum inlet temperature at this

axial location were largely unsuccessful.

There are three sets of experimental measurements.

Measurements of velocity and turbulence quantities are

conditional. That is, the fuel jet, annular jet, and external air were seeded in successive runs. Calculated turbulence kinetic energy and velocities are not conditional. An

average, unconditional velocity is calculated. To calculate

conditional velocities a joint velocity-scalar pdf would

have to be developed.

Predictions of turbulent kinetic energy are very good.

The baseline fine grid calculation best predicts the peak

turbulent kinetic energy. The fine grid pdf simulation

predicts slightly higher peak turbulent kinetic energy. The

two coarse grid solutions predict the highest peaks. The

pdf simulations best predict centerline and near centerline

turbulence kinetic energy. At larger radii, the fine grid

calculations slightly overpredict the turbulence kinetic

energy.

Axial velocity predictions tightly bracket the experimental measurements out to a radius of 5 mm. The baseline

simulations predict higher axial velocity due to the higher

predicted temperatures. The fine grid pdf simulation best
242

predicts the axial velocity deficit between the fuel jet and

annular air flow. At larger radii, predicted axial velocities converge. The pdf simulations predict slightly

higher tangential velocity out to a radius of 12 mm due to

lower temperature predictions and reduced expansion of the

flow. The tangential velocity predictions are quite good.

The coarse grid calculations predict higher tangential

velocity beyond the tangential velocity peak more in agreement with the data. Overall, the baseline coarse grid

calculation best matches tangential velocity data.

Some plots of temperature versus mixture fraction

(total fuel component) at an axial distance of 26 mm are

shown in figure 5.96 for the coarse grid solution. There

are 250 particles in each cell plot. Each computational

grid cell in this region is of the order of one mm in the

radial direction and 3 mm in the axial direction. The cells

are in the mixing layer region where the combustion rate is

highest. The pdf simulation predicts a wide variation of

temperature and mixture fraction in the mixing layer. Each

particle generally falls close to a state function composed

of a straight line connecting mixture fraction of zero,

temperature 300 degrees (pure air) with mixture fraction

0.029, temperature 2550 (stoichiometric or complete burning

of both reactants) and thence a curved line ending at mixture fraction of 1.0, temperature of 300.0 (pure fuel).

IFlI_
243

2500 _"

II

1950 ........................................................................................ ,............................

w

m r= 10.9 mm

<

m 1400 ...... _-_........................... _........................... i......................... W

W

850 ............................................................................. _..........................

300

• i

r= 12.5 mm

• i

I1"

25O0

1950

W

<

1400

W

W

t

85O

• _ r= 8.3 mm

N

i t._.. i

..........................._................_.v.._..._.......................................................

• • i • i°$ •

' " i

% _ r= 9.5 mm

% i

=. i

•

i.

l••

• i

• ° ! 1

300

0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8

MIXTURE FRACTION MIXTURE FRACTION

Figure 5.96 Coarse grid temperature vs mixture fraction at x

= 26 ram.
244

This composite line represents complete burning of fuel for

lean mixtures and complete consumption of oxidant for rich

mixtures, effects of variable Cp, and a fixed ratio of

nitrogen to oxygen. Temperature is highly nonlinear with

mixture fraction. Temperature increases most rapidly near

stoichiometric. This thermodynamic behavior causes very

large temperature gradients in the mixing layer.

The cell plots show that the total fuel component

decreases with increasing radius. Only a few particles are

richer than stoichiometric at a radius of 12.5 mm. There

are at least a dozen partially reacted particles. These

partially reacted particles have temperatures less than 850

degrees. Most particles above 850 degrees are either fully

reacted or very close to being fully reacted. As the majority of particles in these plots are fully reacted, this

explains why efforts to increase reaction rate didn't increase peak temperature predictions. The mixing of fuel and

oxidant species would have to be increased to increase

temperature predictions in the pdf simulation.

Initial coarse grid pdf simulations were done using 100

and 250 particles per cell. Predictions of velocity and

temperature were largely unaffected by varying the particle

numbers. What did seem to affect initial calculations was

the method of treating residual or partial particle transfers. The results shown here were obtained by keeping track

of residual particle transfers and actually exchanging

I!!I I
245

particles when the residual was greater or equal to one,

rather than randomly adding in particles as was done previously. This improved temperature predictions at larger

radii.

Four cell plots of temperature vs mixture fraction for

the fine grid pdf simulation are shown in figure 5.97. There

are only i00 particles in these cells. Most of the particles in these cells plots are fully reacted.

Five coarse grid pdfs of temperature at an axial distance of 26 mm are shown compared to experimental measurements at 25 mm in figure 5.98. The simulated pdf is taken

from the final iteration. The temperatures passed to the

velocity solver are averaged over hundreds of iterations.

Thus, calculated pdfs should be also averaged over hundreds

of iterations. This would require more memory than is

currently needed for the complete calculation, unless a

weighted running average is employed. Temperature measure-

_ents or particle temperature predictions are partitioned

into bins which are multiples of one hundred degrees. The

experimental measurements don't show any temperatures beyond

2500 degrees while the simulations show a few particles just

beyond 2500 degrees. The simulations assume i00 percent

combustion efficiency. The experimental pdfs seem to be

composed of one or two gaussian pdfs. At the largest radial

location the predicted and measured temperature pdf are very

similar, a single peak near 300 degrees. The pdfs taken at
246

2500

i

e "

1950

'" i

rr r: 10.4 mm

rr 1400 ............................................................................................................ LLI

0..

LLI

I-"

850 ............................................... .:................................. i............................

300

i

t

I

I

" i r=l!"_ i ;i PPi

 ̧

• !

25OO

1950

w

_ 1400

o•

!

i

%

i

_. r= 8.4 mm

850 ............................................. -_i ........................................................

300

0.0 0.2 0.4 0.6 0.8

MI_UREFRACTION

-Q

S !

'_. r= 9.7 mm

%

' i

0.0 0.2 0.4 0.6 0.8

MIXTURE FRACTION

Figure 5.97 Fine grid predictions of temperature vs mixture

fraction at x = 26 mm.

_I1i
247

0.8

0.6

i1

£3

a. 0.4

0.2

0.0

0.5

0.4

0.3

n

0.2

0.1

0.0

0.4

0.3

U.

C3

o. 0.2

0.1

0.0

0.4 -

0.3 -

LL

E3

0.2 -

0.1

0.0

0.4 -

0.3 -

LA.

E3

o.. 0.2 -

y=12.8 mm

I I 1 I

y=l 0.4 mm

y---g.6mm

[-],-,___ __

y--8.2 mm

r= 12.5 mm

y=7 mm

A

500 1000 1500 2000 2500

TEMPERATURE (K)

0.8

0.6

0.4 -

0.2 -

0.0 --

0.5 -

0,4 -

0.3 _-

i

0.2 F

0.1

0.0

0.4

0.3

0.2

0.1

0.0

0.4

0.3

0.2

0.1

0.0

0.4

0.3

0.2

i

0.1

0.0

r= 10.9 mm

r= 9.5 mm

l] nnnFIn s-(hn Ooo 

r= 8.3 mm

4- nnnr nFbOOnFh

r= 7.2 mm

0.1 F1nnFnmnnn n 

0.0

0 0 500 1000 1500 2000 2500

TEMPERATURE (K)

i

300O

Figure 5.98 Measured temperature pdfs at x : 25 mm, and

predicted coarse grid temperature pdfs at 26 mm.
248

the next smaller radial distance still show a pdf peak near

300 degrees, but differ largely in the number of particles

at higher temperatures. The simulation correctly predicts a

wide variation in temperature at this location. At the 9.5-

9.6 mm radial cell location, the experimental peak pdf

occurs at higher temperature, while the simulation shows a

pdf peak at much lower temperature. The average temperature

plot showed an experimental peak temperature at 9.6 mm,

while the coarse grid pdf simulation predicted a peak temperature at 8.3 mm. This agrees with the pdf plots shown

here. The pdf coarse grid simulation shows the largest

number of high temperature particles at 8.3 mm. The pdf is

fairly uniform at this radius and this is one of the reasons

why a fine grid pdf solution was done.

A comparison of fine grid pdfs versus the same experimental pdfs is shown in figure 5.99. The fine grid pdf

simulation predicts a peak temperature at 9 mm. This cell

"isn't shown here, but the plot for 8.4 mm shows more particles at higher temperature than does the coarse grid pdf at

8.3mm.

Profiles of temperature, turbulent kinetic energy, L

axial and tangential velocity at 50 mm are shown in figure

5.100. Peak temperature is best predicted in the coarse

grid pdf simulation. At this station, the fine grid pdf

simulation underpredicts peak temperature. The baseline

simulations overpredict peak temperature by over 300

IllT!i_
249

0.8 j"

I

0.6 -

LL

E3

O.. 0.4

0.2

0.0

0.5

0.4

0.3

0.

0.2

0.1

0.0

0.4

0.3

LL

o

a. 0.2

0.1

0.0

0.4 -

0.3 "

IJ.

o

o. 0.2

0._ "

0.0

0.4 [

0.3 ["

0.1 t

I

0.0

0

I I

y=12.8 mm

I I

y=10.4 mm

nOFLn _ n n nF On nn ,

y=9.6 mm

y=8.2 mm

y=7 mm

500 1000 1500 2000 2500

TEMPERATURE (K)

r= 12.7 mm

r= 10.4 mm

r= 7.4 mm

500 1000 1500 2000 2500

TEMPERATURE (K)

i

3000

Figure 5.99 Measured temperature pdfs at x = 25 mm and

predicted fine grid temperature pdfs at x =26 mm.
25O

2100

1800

1500

I,IJ

n-

<_ 12oo

LLI

Q.

ILl

I'- 900

6OO

................ :................... ;1_<.................. t- -- _ _o I..........

i II _ I--_ 1 >- 400

........_ .... ,Y p', ..............................i............ z

i /I /.,d i '"

................ i//_i't _............................................................... i _

_f I V',: i _=

7/-_......... _'....... T.............. -:=

................ i !

3OO

0.000 0.005 0.010 0.015 0.020 0.025 0.030

RADIUS (m)

i

I--I

,'............ i .... ; .......... i................... :.................. _..............

/ i i

0

0._ o._s o.010 0.015 o.o20 0.o25 0.03o

RADIUS (m)

15o

A

lOO

i

I.LI

:>

.-I

,_ 5o

,,<

10

<?A i ! : i

.J

.............. f_

Z

_ 2

0.000 0.005 0.010 0.015 0.020 0.025 0.030

RADIUS (m)

................................ !........i................................... i.......

• / i ".,.. !

.................. _i7 _.................. i...................

............ L_:I i..........................

....................i................

OJ_o 0.005 0.010 0.015 0.020 0.025 0.030

0 RADIUS (m)

Figure 5.100 Temperature, turbulent kinetic energy, axial

and tangential velocity at x = 50 mm.

Iii I !:
251

simulations overpredict peak temperature by over 300 degrees. All numerical simulations predict very sharp temperature peaks while temperature measurements show more rounded

peaks. The pdf simulations best predict turbulent kinetic

energy. The baseline simulations predict larger radial

displacement of the turbulent kinetic energy peak, overpredict the peak, and underpredict turbulence kinetic energy at

and near the combustor axis. Axial velocity is best predicted by the pdf simulations. Baseline simulations overpredict centerline axial velocity by about 20 %. The axial

velocity predictions between the simulations converge with

increasing radius. The baseline coarse grid simulation best

predicts the tangential velocity profile. The coarse grid

simulations are almost identical beyond a radius of 13 mm.

The fine grid tangential velocity predictions are almost

identical beyond a radius of 12 mm. The tangential velocity

predictions are affected by predicted heat release (radial

flow expansion) and also show grid dependence.

Comparisons at an axial distance of 75 mm from the

inlet are shown in figure 5.101. The pdf simulations underpredict the peak temperature by about 100 degrees, while the

baseline simulations overpredict temperature by 300-400

degrees. Predicted temperature peaks are quite sharp, while

the temperature measurements are starting to show a plateau

like area of high temperature just inside the peak temperature. The fine grid pdf simulation shows the bluntest
252

2100

1800

_r

iii

_- 900

6OO

200

. : = :

f

................................ o:_ .........._i;.................. i................... i......

..... !......... I.............. ?............... _...............

................ ',o......i.................. i...................

0.000 0.00S 0.010 0.015 0.020 0.025 0.030

RADIUS (m)

r_ 4O0

rr

LU

Z

LU

_z

k--

-J

I- 100

i : I-- -='_'_ I

......._............ _ ] .... _"_I

i _ T-'- _,._o I

/'X, i I_ .._*,I

./............... :_............ :................... !......L__ .-_,_ ...........

...... o_-.._.. ................ _............... _................_.............

-\\_.\ i ! i

................... :.................... .................. .................. ...................

- _._, i i

0.000 0.005 0.010 0.015 0.020 0.025 0.030

RADIUS (m)

150

100

LU

>.

...I

,<

_5o

i

( "_ i

o:%_._ ! i i

0 '' "

0,000 0.005 0,010 0.015 0.020 0.025 0.030

RADIUS (m)

10

_o6

..J

LU

.-I

UJ

Z

2

...................... /<._:_ _,.............. '_.................. i................

/

/ i !

................. :_.L....... i................... i..................

/'/, i i

,/

0.(3_00 (_0q_j 0.010 0.015 0.020 0.025 0.030

"CO RADIUS (m)

Figure 5.101 Temperature, turbulent kinetic energy, axial

and tangential velocity at x = 75 mm.

I[1 1-i
253

temperature peak. The pdf simulations excellently predict

turbulent kinetic energy and axial velocity, while the

baseline simulations overpredict both quantities out to a

radius of 15 mm. From a radius of 9 to 15 mm, the predictions of axial velocity are surprisingly close to each

other, despite the large difference in temperature predictions at this axial station. The largest difference in

axial velocity predictions occurs along the flame centerline. The fine grid simulations predict the lowest tangential velocity which is also closest to the experimental

measurements. The coarse grid pdf simulation predicts

highest tangential velocity, in part due to lower predicted

heat release. At the previous axial station, the coarse

grid simulations best predicted tangential velocity beyond

the peak. Under a radius of 8 mm at this station, the shape

of the tangential velocity appears to be poorly predicted.

Predicted tangential velocity profiles show almost linear

increases with radial distance out to 10 mm. The experimental data shows near zero, or even negative tangential

velocity out to a radius of 5 mm and then a fairly linear

profile to the tangential velocity peak. This predictive

behavior is also shown in previous plots.

The comparisons at an axial distance of 150 mm are

shown in figure 5.102. Baseline simulations overpredict the

peak by almost 500 degrees. Pdf simulations slightly

overpredict the temperature peak. All simulations show a
254

2100

12O0

t-- -- " p,_ _eo i

, L2__

t,

s t !

.................. +............................ . .....................

I

I-- 900 .............................................................................................................

0.00 0.02 0.04 0.06 0.08

RADIUS (m)

5OO

)- 4o0

r_

W

Z

LLI

z

I--

pelf,Io_3o

i i "°'"

m "''_ ; z

I

= "-'_- % ! i .............

0.00 0.02 0.04 0.06 0.08

RADIUS (m)

B0

Lu

>

<

2o

0.00 0.02 0.04 0.06

RADIUS (m)

0.08

4

_3

LU

:>

,-I

_2

LU

;E

1

(

(

o I

0.0(

/ i

i ....................... :-.................... "....................

i i !

i i ',

,: ................. _............................ i................

o :.,",,, i

0

0.02 0.04 0.06 0.08

RADIUS (m)

Figure 5.102 Temperature, turbulent kinetic energy, axial

and tangential velocity at x = 150 mm.
255

strong drop in temperature towards the combustor axis, while

temperature measurements show a nearly constant profile.

Turbulent kinetic energy is best predicted in the pdf simulations. The baseline simulations overpredict centerline

turbulent kinetic energy and axial velocity. The jet seeded

axial velocity measurements are larger than the annulus

seeded velocity measurements at the same radii. The pdf

simulations predict axial velocity profiles lying between

the sets of conditional measurements. Fine grid predictions

of tangential velocity are very similar to each other, as

are the two coarse grid simulations. The fine grid simulations best predict tangential velocity, while the coarse

grid simulations overpredict tangential velocity by 50

percent.

The largest disagreement in the sets of conditional

velocity measurements occurs for radial velocity, which is

shown in figure 5.103. The experimental data does not

benchmark or validate the predictions. Instead, the data

largely bounds the various velocity predictions. Measurements of fuel seeded flow show peaks at least twice as high

as all predictions. Radial velocity measurements taken at

the 75 and 150 mm axial locations by seeding the fuel jet

are much higher than measurements found by seeding the

annulus and external air flows. To best compare traditionally calculated velocities, velocity measurements should be

unconditional or averaged. All inlet flows should be
256

)-

__ s

8

.J

LU

>

3

<

n-

-I

_D ----- _3o

o_ o,,

annulus

o

[ A ex'ie.-,ot

.............. ! .......................................................................

o

DE]DOE ]

O O

....._..................... ._.................. .i..................

0.000 0.005 0.010 0.015 0.020 0.025 0.030

RADIUS (m)

s

O

...I

UJ

>

5

,<

n'-

:S0mm

7 ............................................................................. _......................................

i '

!oo i i

......... o? Oo.........i.........._.....................

O O_ oC)[3DDO

[3

'"_ .................................. _5"'"_i.................. _i ......................................

....... /, - '_._._ ' ];

-1

0.000 0.0o5 0.010 0.015 0,020 0.025 0.o3o

RADIUS (m)

7

8,

ILl

>

•J 3

<

rr

X = 75 mm

i

i

ioO°_ : §

oo o i '"

o ! i >

............. i................... _"................. _................... i.................. i.................

o i : ! <

0.0o0 o.oos 0._ 0.0_5 0.020 0.025 0.0o0

RADIUS (m)

x= 150ram

7 .................................... '....................................................

_ 5 ............................ ";"............................ _".........................................................

,.(:7_0

1 0 ........ i ...................................

-1

0.00 0.02 0.04 0.0O

RADIUS (m)

0.08

Figure 5.103 Radial velocity at x = 25, 50, 75, and 150 mm.

I!Il i
simultaneously seeded.

257
258

Chapter 6

CONCLUDING REMARKS

A joint pdf solver for species and temperature was

incorporated and further developed. Velocity, pressure,

turbulent kinetic energy and dissipation are solved for

using the SIMPLE algorithm of Patankar. The pdf model uses

Monte Carlo solution technique to solve for the evolution of

the pdf. The effects of transport between fluid cells,

mixing and combustion within fluid cells are treated sequentially. The composition pdf handles the turbulence combustion interaction due to varying species and temperature.

This eliminates the need to correct reaction rate terms due

to turbulence effects, as is typically done in current

combustion models. Turbulence causes large local variation

in species and temperature. As reaction rates are nonlinear functions of species and temperature, reaction rates

determined from mean species and temperature are physically

incorrect. Pdf modeling allows the proper incorporation of

reaction rates without corrections.

The pdf method used in this work employs a constant

number of particles in each fluid cell. An Eulerian framework is to solve for the evolution of the pdf. The solution

of particle transfers, mixing and combustion gives a finite

number representation of the true pdf. This corresponds

nicely to measurement techniques which also use a finite

number of data points to ascertain pdfs.

iF!I;
259

Various improvements were made to the pdf model.

First, variable Cp and proper treatment of thermodynamic

properties was incorporated into the model. The original

model assumed a single specific heat for all species and

temperatures. Temperature rather than enthalpy was used to

keep track of energy. The improved thermodynamic model uses

data in the form of a two temperature range fifth order

polynomial to calculate specific heat and enthalpy. This

way of handling thermodynamic data has been widely used in

equilibrium and finite rate single cell or one dimensional

combustion calculations. The method is accurate over a wide

temperature range and is very highly regarded in the thermodynamic field. An improved Newton-Raphson species solver

was incorporated into the composite pdf for solving finite

chemical reaction rates.

The validation of the improved hybrid pdf model involves the comparison of pdf simulations, baseline calculations, and experimental data. The baseline combustion model

uses the SIMPLE algorithm to solve for velocity, pressure,

species and enthalpy. The baseline model' also includes the

effect of variable specific heat. Originally it was desired

to calculate spray combustion, so the baseline combustion

model was evaluated for these flows. The liquid phase fuel

spray was solved in an Lagrangian framework. The interaction between the liquid phase and gaseous phase was handled

by including drag effects on liquid drops and source terms
260

in the gas phase for the evaporating spray. The baseline

spray predictions for an evaporating methanol case were very

good. Acceptable results were also obtained for the isothermal swirling air in a spray combustor by including some

turbulence modeling changes. The baseline model also gave

qualitatively good results for the same combustor with

heptane spray combustion. Turbulence model changes did not

need to be done for the combusting case, as turbulence

levels generated by combustion were higher. The region

containing spray was not particularly well calculated despite various changes. The largest source of error in the

spray calculation was caused by an overly quick dissipation

of the gaseous velocity components surrounding the spray.

The spray occurs in a high shear region where the gaseous

flow is severely skewed to the computational grid. Velocity

predictions further from the spray atomizer were much better. As this deficiency in calculating spray combustion

•would not be improved by turbulence combustion models,

further testing was confined to single phase flows. It was

decided not to pursue pdf modeling of combusting sprays.

There are a few misconceptions about pdf modeling in

the CFD community. Pdf modeling does not predict emissions

if they are not included in the reaction mechanisms. Also,

pdf modeling will not make up for deficientcies in turbulence modeling. The hybrid pdf model which was used here

used the k-£ turbulence model to predict mixing. The k-£
261

model is known to be deficient for high swirl and turbulence

transition.

The hybrid pdf was first tested for an enclosed diffusion combustor. A set of co-annular jets supplied methane

and air to a dump combustor. Simulations predicted an outer

recirculation zone. The pdf simulation predicted a recirculation zone length almost 25% longer. As the recirculation

zone was hot, and the central fuel core was cool, this

provided an continuing combustion source. Both simulations

predicted substantial overlap of fuel and oxygen concentration profiles in agreement with measurements. Hybrid pdf

predictions were most often superior to results obtained

with the baseline method, especially along the combustor

centerline. The hybrid pdf produced superior predictions of

oxygen, carbon dioxide and water concentration. Fuel was

significantly overpredicted by both hybrid and baseline

simulations at three axial stations. Both simulations did

an excellent job of predicting the first temperature profile. Temperatures on and near the combustor axis were much

better predicted by the hybrid pdf simulation for the second

and final temperature profile. The methane combustion

simulations were beginning to show larger temperature differences with axial distance. More temperature measurements

would have helped to better evaluate the combustion models.

Also, measured velocities were not available to validate the

calculations. Temperature plots of the pdf particles at
262

selected cells showed substantial numbers of partially

reacted particles co-existing with very hot, fully reacted

particles.

The second case involved a swirling combustor flow

characteristic of an actual gas turbine combustor. The high

swirl component causes a central recirculation zone. Inlet

fuel and air flows were at fairly high velocity at a severe

angle to the computational grid for this combustor. This is

where the baseline method was identified as being deficient

for spray calculations. The fuel jet was very coarsely

modeled. Combustor predictions of velocity and species were

more diffuse than measurements. The hybrid pdf results for

the second case did not display as much radial displacement

of fuel as was measured. Thus, the pdf modeling did not

make up for the deficiency identified in the baseline combustor spray calculations. It is likely that a pdf calculation done in a Lagrangian framework would prove to be less

diffusive, and improve flow predictions. Substantial grid

modification to better model the fuel jet would also help

modeling this flow. Velocity and turbulence predictions for

this type of flow would probably be improved with the implementation of higher-order numerical scheme. However, such

fixes are not known for their robustness. The hybrid pdf

simulation did predict a slower rate of combustion more in

agreement with experimental temperature and species data.

The third case involved premixed combustion in a
263

simple dump type combustor. This tested whether pdf modeling could be used in cases were the fuel and oxidant are

intimately mixed. The hybrid pdf results were generally

superior to the baseline simulation. The baseline calculation best predicted velocity and temperature profiles for

the first two axial stations. The baseline calculation

predicted lower temperatures than measured for substantial

portions of the flow at succeeding axial stations. The

baseline simulation of combustion seemed overly limited by

the eddy-breakup reaction rate for latter portions of the

flowfield. Efforts to improve the baseline prediction by

increasing the reaction rates resulted in significant change

in the flowfield and overall solution degradation. The pdf

simulation predicted higher temperature gradients than the

baseline simulation. This particularly improved predictions

for the latter portion of the combustor flowfield. A test

was done to evaluate the molecular mixing submodel in the

pdf model. A second simpler relax-to mean mixing submodel

was incorporated. This change didn't significantly affect

the pdf predictions. The simpler model is more economical

to use.

Very low levels of turbulence were measured in the

experimental apparatus. A slight deficientcy in the turbulence model is suggested. While the pdf model produced

improved temperature profiles, it should not be implied that

the pdf model make up for deficientcies in turbulence model.
264

The fourth case involved a swirling hydrogen-air jet

flame. The swirl mass flow component was small enough that

recirculation was not involved. The hybrid pdf solver produced superior predictions of temperature, axial velocity,

and turbulent kinetic energy. The pdf simulation and baseline coarse grid simulation both predicted higher tangential

velocity than measured later in the flowfield. A fine grid

pdf calculation better predicted tangential velocity decay.

The baseline fine grid calculation better predicted tangential velocity decay.

The pdf simulations correctly predict turbulent phenomena, specifically, large variation in local species concentration and temperature. This represents an improvement in

turbulent-combustion modeling where turbulent combustion has

been either ignored or various ad-hoc corrections have been

applied to chemical reaction rates. Current multi-dimensional CFD modeling generally produces predict too large

temperatures and superfluous product species concentrations.

The inclusion of multiple species concentrations and temperatures ala pdf modeling helps correct this deficiency. The

pdf method usually improves combustor flowfield predictions.

These cases prove that the pdf method can be used for both

diffusion flames and premixed flames provided a proper

reaction model is used.

All simulations here used a single reaction step.

Including multiple reaction rates is known to improve pre-
265

dictions of species and temperature for zero or one dimensional flows, and allows the calculation of additional

species, specifically, pollutants. The inclusion of multiple finite-rate reaction steps in current models, while

possible, involves increasing numerical inaccuracy. Multiple reaction steps can involve the introduction of highly

reactive intermediate species. The solution of reaction

rates at higher temperatures in current models is only accurate for very small time steps. Errors generated by using

too high of reaction rates can produce physically impossible

predictions of temperature and product species concentrations. Use of finite reaction rate models in current multidimensional models involves overwriting negative species

concentration. Such corrections depend a lot on the experience of the CFD modeler. The pdf method can easily incorporate chemical kinetic rate solvers which have been developed

to overcome these problems. These solvers can handle the

largest chemical kinetic reaction schemes that have been

proposed. The drawback to this is a rather large amount of

computer resources and time to perform calculations. Pdf

modeling may be used to improve the solution accuracy of

combustor predictions done with simpler combustion models.
266

References

forest

i. Patankar, S.V., "Numerical Heat Transfer and Fluid Flow",

Hemisphere Publishing Corp., New York, 1980.

2. Chen, K-H, Duncan, B., Fricker, D., Lee, J., and Quealy, A.,

"ALLSPD-3D, Version 1.0", NASA Lewis Research Center, Nov., 1995.

3. Joshi, D.S., and Vanka, S.P., "A Multigrid Calculation Procedure

for Internal Flows in Complex Geometries", AIAA 90-0442, 1990.

4. Tafti, D.K., and Vanka, S.P., "Hot Gas Environment Around STOVL

Aircraft in Ground Proximity: Part 2: Numerical Study", AIAA 90-

2270, 1990.

5. Vanka, S.P., and Krazinski, J.L., "Efficient Computational Tool

for Ramjet Combustor Research", J. Propulsion, Vol. 5, pp. 431-437,

1989.

6. Radhakrishnan, K., and Bittker, D., "LSENS, A General Chemical

Kinetics and Sensitivity Analysis Code for Homogeneous Gas-Phase

Reactions", NASA RP-1328, 1329, and 1330, 1994.

7. Jones, W.P, and Launder, B.E., "The Prediction of Laminarization

with a two-equation model of turbulence", Int. J. Heat Mass

Transfer, 17ol 15, pp. 301-314, 1972.

8. Correa, S.M., and Shyy, W., "Computational Models and Methods

for Continuous Gaseous Turbulent Combustion", Proq. EnerqvCombust.

Sci______., Vol. 13, pp. 249-292, 1987.

9. Leonard, B.P., "A Stable and Accurate Convective Modeling

Procedure Based on Quadratic Upstream Interpolation", Comput. Meth.

Appl. Mech. En_ng., Vol. 19, pp. 59-98, 1979.

i0. Faeth, G.M., "Current Status of Droplet and Liquid Combustion",

Prog. Enercfv Combust. Sci., Vol. 3, pp. 191-224, 1977.

II. Law, C.K., "Recent advances in Droplet Vaporization and

Combustion", Proq. Enerqy Combust. Sci., Vol. 8, pp. 171-201, 1982.

12. Faeth, G.M., "Vaporization and Combustion of Sprays", Proq.

Enercrf Combust. Sci., Vol. 9, pp. 1-76, 1983.

13. Sirignano, W.A., "Fuel Droplet Vaporization and Spray

Combustion Theory", Proq. Enerqy Combust. Sci., Vol. 9, pp. 291-

322, 1983.

14. Onuma, Y., and Ogasawara, M., "Studies on the Structure of a

Spray Combustion Flame", Fifteenth Symposium (International) on

Combustion, The Combustion Institute, Pittsburgh, pp. 453-465,

!_[TII
267

1975.

15. Spaulding, D.B., "The Combustion of Liquid Fuels", Fourth

Symposium (International) on Combustion, The Combustion Institute,

Pittsburgh, pp. 847-864, 1953.

16. Faeth, G.M., "Mixing, Transport and Combustion in Sprays",

Prog. Energy Combust. Sci., Vol. 13, pp. 293-345, 1987.

17. Tsai, J.S., and Sterling, A.M., "The Combustion of Linear

Droplet Arrays", Twenty-Third Symposium (International) on

Combustion, The Combustion Institute, Pittsburgh, pp. 1405-1411,

1990.

18. Odgers, J., Kretschmer, D., Norster, E.R., Lefebvre, A.H., and

Fletcher, R.S., "The Design and Development of Gas Turbine

Combustors", Northern Research and Engineering Corp., Woburn,

Massachusetts, 1980.

19. Tong, A.Y., and Sirignano, W.A., "Multicomponent Transient

Droplet Vaporization with internal Recirculation: Integral

Formulation and Approximate Solution", Numerical Heat Transfer,

Vol. I0, pp. 253-278, 1986.

20. Raju, M.S., and Sirignano, W.A., "Multicomponent Spray

Computations in a Modified Centerbody Combustor", J. of Propulsion,

Vol. 6, pp. 97-105, 1990.

21. Reid, R.C., Prausnitz, J.M., and Sherwood, T.K., "Properties of

Liquids and Gases", McGraw-Hill, 1977.

22. Kumar, R. K., "Ignition of Hydrogen Oxygen_Diluent Mixtures

Adjacent to a Hot, Nonreactive Surface," Combustion and Flame, Vol.

75, pp. 197-215, 1989.

23. Kim, S.-W., "Numerical Investigation of Chemical Reactionturbulence Interaction in Compressible Shear Layers," Combustion

and Flame, Vol. i01, pp. 197-208, 1995.

24. Dryer, F.L., and Glassman, I., "High Temperature Oxidation of

Co and CH4", Fourteenth Symposium (International) on Combustion,

The Combustion Institute, Pittsburgh, pp. 987-1003, 1972.

25. Westbrook, C. K., and Dryer, F. L. "Simplified Reaction

Mechanisms for the Oxidation of Hydrocarbon Fuels in Flames",

Combustion Science and Technoloq%I, Vol. 27, pp.31-43, 1981.

26. Coffee, T. P., "On Simplified Reaction Mechanisms by Oxidation

of Hydrocarbon Fuels in Flames by C. K. Westbrook and F. T. Dryer",

Combust. Sci. and Tech., Vol. 43, pp. 333-339, 1985.
268

27.Coffee, T. P., Kotlar, A. J., Miller, M. S., "The Overall

Reaction Concept in Premixed, Laminar, Steady-State Flames. I.

Stoichiometries",Combustion and Flame, Vol. 54, pp. 155-169, 1983.

28. Nikjooy, M., So, R.M.C., "On the Modelling of Non-Reactive and

Reactive Combustor Flows", NASA CR 4041, 1987.

29.Srivatsa, S. K., "Computations of Soot and Nox Emissions from

Gas Turbine Combustors", NASA CR-167930, 1982.

30. Hautman, D.J., Dryer, F.L., Schug, K.P., and Glassman, I., "A

Multi-step Overall Kinetic Mechanism for the Oxidation of

Hydrocarbons", Combust. Sci. and Tech., Vol. 25, pp. 219-235, 1981.

31. Abdalla, A.Y., Bradley, D., Chin, S.B., and Lam, C., "Global

Reaction Schemes for Hydrocarbon Oxidation", Oxidation

Communications, Vol. 4, pp. 113-130, 1983.

32. Duterque, J., Borghi, R., and Tichtinsky, H., "Study of QUasiGlobal Schemes for Hydrocarbon Combustion", Comb. Sci. and Tech.,

Vol. 26, pp. 1-15, 1981.

33. Westbrook , C.K., and Dryer, F.L., "Chemical Kinetic Modeling

of Hydrocarbon Combustion", Prog. Enerc[y Combust. Sci., Vol. i0,

pp. 1-57, 1984 . _ ....

34. Paczko, G., Lefdal, P.M., and Peters, N., "Reduced Reaction

Schemes for Methane, Methanol, and Propane Flames", Twentv-first

Symposium (International) on Combustion, The Combustion Institute,

Pittsburgh, pp. 739-748, 1986.

35. Seshadri, K., Peters, N., "Asymptotic Structure and Extinction

of Methane-Air Diffusion Flames", Combustion and Flame, Vol. 73,

pp. 23-44, 1988.

36. Kiehne, T.M., Matthews, R.D., and Wilson, D.E., "An Eight-Step

Kinetics Mechanism for High Temperature Propane Flames", Combust.

Sci..and Tech., Vol. 54, pp. 1-23, 1987.

37. Jones, W.P., and Lindstedt, R.P., "Global Reaction Schemes for

Hydrocarbon Combustion", Combustion and Flame, Vol. 73, pp. 233-

249, 1988.

38. Bilger, R.W., and Starner, S.H., "On Reduced Mechanisms for

Methane-Air Combustion in Nonpremixed Flames", Combustion and

Flame, Vol. 80, pp. 135-149, 1990,

39. D. B. Spaulding, "Mixing and Chemical Reaction in Steady,

Confined Turbulent Flames", Thirteenth Symposium (International) on

Combustion, The Combustion Institute, Pittsburgh, pp. 649-657,

1970.

If!11
269

40. D. B. Spaulding, "Concentration Fluctuations in a Round

Turbulent Free Jet," Chem. Enq. Sci, vol 26, pp. 95-107, 1971.

41. Mason, H. B., and Spaulding, D. B., "Prediction of reaction

rates in turbulent pre-mixed boundary layer flows", Combustion

Institute European Symposium, pp. 601-606, 1973.

42. Spaulding, D. B., "Mathematical Models of Turbulent Flames; A

Review," Comb. Sci. and Tech., Vol 13, pp. 3-25, 1976.

43. Magnussen, B. F., and Hjertager, B. H., "On Mathematical

Modeling of Turbulent Combustion With Special Emphasis On Soot

Formation and Combustion," Sixteenth Symposium (International) on

Combustion, The Combustion Institute, Pittsburgh, pp. 719-729,

1976.

44. Srinivasan, R., et. al., "Aerothermal Modeling Program Phase I

Final Report," NASA CR-168243, GARRETT 21-4742-1, 1983.

45. Libby, P. A., and Williams, F. A., " Turbulent Reactinq Flows,"

Springer-Verlag, 1980.

46. Richardson, J.M., Howard, H.C., and Smith, R.W., "The Relation

Between Sampling Tube Measurements and Concentration Fluctuations

in a Turbulent Gas Jet", Fourth Symposium (International) on

Combustion, The Combustion Institute, Pittsburgh, pp. 814-817,

1953.

47. Jones, W. P., "Models for Turbulent Flows with Variable

Density and Combustion," in Prediction Methods for Turbulent Flows

(Edited by Kollman, W.), pp. 379-422, Hemisphere Publishing Corp.,

1980.

48. Tsuji, H. and Yamoka, L. "Structure Analysis of Counterflow

Diffusion Flames in the Forward Stagnation Region of a Porous

Cylinder, " Thirteenth Symposium (International) on Combustion, The

Combustion Institute, Pittsburgh, pp. 723-731, 1971.

v

49. Bilger, R. W., "Reaction Rates in Diffusion Flames", Combustion

and Flame, Vol. 30, pp. 277-284, 1977.

50. Jeng, S.-M., "An Investigation of Structure and Radiation

Properties of Turbulent Buoyant Diffusion Flames, Ph.D. Thesis,

Pennsylvania State University, 1984.

51. Spaulding, D. B., "GENMIX: A General Computer Proqram for Twodimensional Parabolic Phenomena", Pergamon Press, 1977.

52. Liew, S. K., Bray, K. N. C., and Moss, J. B., "A flamelet model

of Turbulent Non-premixed Combustion", Comb. Sci. and Tech., Vol.

27, pp. 69-73, 1981.
270

53. Liew, S.k., and Bray, K.N.C., "A Stretched Laminar Flamelet

Model of Turbulent Nonpremixed Combustion", Combustion and Flame,

Vol. 56, pp. 199-213, 1984.

54. Drake, M.C., and Blint, R.J., "Structure of Laminar Opposedflow Diffusion Flames With CO/H2/N 2 Fuel", Combust. Sci. and Tech.,

Vol. 61, pp. 187-224, 1988.

55. Borghi, R., "Turbulent Combustion Modelling", Proq. Enerqy

Combust. Sci., Vol. 14, pp. 245-292, 1988.

56. Peters, N., "Laminar diffusion flamelet models in nonpremixed

turbulent combustion", Pro_. Enerqy Combust. Sci., Vol. I0, pp.

319-339, 1984.

57. Janika, J., and Kollman, W., "A two-variable formalism for the

treatment of chemical reactions in turbulent H2-air diffusion

flames, Seventeenth Symposium (International) on Combustion, The

Combustion Institute, Pittsburgh, pp. 421-431, 1979.

58. Dixon-Lewis, G., Goldsworthy, F. A., and Greenburgh, J.B.,

"Flame structure and flame reaction kinetics. IX. Calculation of

properties of multi-radical premixed flames", Proc. R. Soc. Lond.,

Vol. 346, 261-278, 1975.

59. Janika, J., and Kollman, W., "The Calculation of Mean Radical

Concentrations in Turbulent Diffusion Flames", Combustion and

Flame, Vol. 44, pp. 319-336, 1982.

60. Kent, J.H., and Bilger, R.W., "Turbulent Diffusion Flames",

Fourteenth Symposium (International) on Combustion, The Combustion

Institute, Pittsburgh, pp. 615-625, 1974.

61. Correa, S.M., Drake, M.C., Pitz, R.W., and Shyy, W.,

"Prediction and Measurement of a non-equilibrium turbulent

diffusion flame", Twentieth Symposium (International) on

Combustion, The Combustion Institute, Pittsburgh, pp. 337-343,

1984.

62. Correa, S.M., "A model for non-premixed turbulent combustion of

CO/H 2 jets", Archvm. Combust., Vol. 5, 223-242, 1985.

63. Pope, S.B., and Correa, S.M., "Joint PDF Calculations of a Nonequilibrium Turbulent Diffusion Flame", Twentv-first Symposium

(International) on Combustion, The Combustion Institute,

Pittsburgh, pp. 1341-1348, 1986.

64. Mongia, H.C., "Combustion Modeling in Design Process:

Applications and Future Direction", AIAA 94-0466, 1994.

65. Kuo, K.K., "Principles of Combustion", John Wiley & Sons, New

York, 1986.

?
271

66. Williams, F.A., "Combustion Theory",

Publishing Company, Inc., Menlo Park, 1985.

Benjamin/Cunu_ings

67. Bartok, W., and Sarofim, A.F., "Fossil Fuel Combustion, A

Source Book", John Wiley & Sons, New York, 1991.

68. Chigier, N., "Enerqv, Combustion and Environment", McGraw-Hill,

New York, 1981.

69. Oran, E.S., and Boris , J.P, "Numerlcal Simulation of Reactive

Flo_____ww", Elsevier, New York, 1987.

70. Anand, M.S., Pope, S.B., and Mongia, H.C., "A PDF Method for

Turbulent Recirculating Flows", in Turbulent Reactive Flows,

Lecture Notes in Engineering, Springer-Verlag, Vol. 40, pp. 672-

693, 1989.

71. Curl, R.L., "Dispersed Phase Mixing, I, Theory and Effects of

Simple Reactors", A.I.Ch.E.J., vol. 9, p. 175-181, 1963.

72. Janika, J., Kolbe, W., and Kollman, W. "Closure of the

Transport Equation for the Probability Density Function of

Turbulent Scalar Fields," J. Non-Equilib. Thermodyn., Vol. 4, pp.

47-66, 1979.

73. Dopazo, C., "Pdf approach for a turbulent axisymmetric heated

jet: Centerline evolution", Physics of Fluids, Vol. 18, 397-404,

1975.

74. Pope, S. B., "Pdf Methods for Turbulent Reactive Flows", Pro_.

Enerqv Combust. Sci., Vol. Ii, pp. 119-192, 1985.

75. Hsu, A. T., and Chen, J. Y., "A Continuous Mixing Model for Pdf

Simulations and its Application to Combusting Shear Flows", _iqhth

Symposium on Turbulent Shear Flows, Technical University of Munich,

p. 22-4, Sept. 9-11, 1991.

76. Chen, J.Y., and Kollman, W., "Segregation Parameters and PairExchange Mixing Models for Turbulent Nonpremixed Flames", TwentyThird Symposium (International) on Combustion, The Combustion

Institute, Pittsburgh, pp. 751-757, 1990.

77. Pope, S. B., "An Improved Turbulent Mixing Model", Comb. Sci.

Tech., Vol. 28, pp. 131-135, 1982.

78. Norris, A. T., and Hsu, A., "Comparison of PDF and Moment

Closure Methods in the Modeling of Turbulent Reacting Flows", NASA

TM 106614, 1994.

79. Masri, A. R., Dibble, R. W., and Barlow, R. S., "The Structure

of turbulent, pilot stabilized Flames of CH4/CO/H2/N 2 fuel

mixtures", Combustion and Flame, Vol. 88, pp. 239-264, 1992.
272

80. Villermaux, J., "Micromixing phenomena in stirred reactors",

Encyclopedia of Fluid Mechanics, pp. 707-768, Gulf Pub. Co., West

Orange, NJ, 1986.

81. Yamazaki, and Ichigawa, A., Intl. Chem. Enq., Vol. i0, pp. 471-

478, 1970.

82. Correa, S.M., and Braaten, M.E., "Parallel Simulations of

Partially Stirred Methane Combustion", Combustion and Flame, Vol.

94, pp. 469-486, 1993.

83. Correa, S.M., "A Direct Comparison of Pair-Exchange and IEM

Models in Premixed Combustion", Combustion and Flame, Vol. 103, pp.

194-206, 1995.

84. J.Y. Chen, and Kollman, W., "Mixing Models for Turbulent Flows

with Exothermic Reactions", Turbulent Shear Flows 7, SpringerVerlag, Berlin Heidelberg, pp. 277-292, 1991.

85. Chen, J. Y., Kollman, W., and Dibble, R. W., "Pdf modeling of

Turbulent Nonpremixed Methane Jet Flames," Combust. Sci. and Tech.,

Vol. 64, pp. 315-346, 1989.

86. Maas, U.A., and Pope, S. B., "Simplifying chemical kinetics:

Intrinsic Low-Dimensional Manifolds in Composition Space",

Combustion and Flame, Vol. 88, pp 239-264, 1992.

87. Chen, J.Y., and Kollman, W., "Chemical Models for Pdf Modeling

of Hydrogen-Air Nonpremixed Turbulent Flames", Combustion and

Flame, Vol. 79, pp 75-99, 1990.

88. Hsieh, K.-C., Shuen, J.-S., Tsai, Y.-L. Peter, and Yu, S.-T.,

"RPLUS2D/3D USER'S MANUAL," Sverdrup Technology, Inc., NASA Lewis

Research Center, Oct. 22, 1990.

89. Yoon, S., and Jameson, A. "An LU-SSOR Scheme for the Euler and

Navier-Stokes Equations," AIAA Paper 87-0600, 1987.

90. Burros, M. C., and Kurkov, A. P., "Analytical and Experimental

Study of Supersonic Combustion of Hydrogen in a Vitiated Air

Stream," TMX-2828, NASA, 1973.

91. Rogers, R. C., and Chinitz, W., "Using a Global Hydrogen-Air

Combustion Model in Turbulent Reacting Flow Calculations", AIAAJ.,

Vol. 21, pp. 586-592, 1983.

92. Strahle, W.C., and Lekoudis, S.G., "Evaluation of Data on

Simple Turbulent Reacting Flows", AFSOR-TR-85-0880, 1985.

93. Kline, S.J., Morkovin, M.J., and Moffat, "Computation of

Turbulent Boundary Layers", !968 AFSOR-IFP-Stanford Conference,

1969.

If|ill_
273

94. Baker, R.J., Hutchinson, P., Khalil, E.E., and Whitelaw, J.H.,

"Measurements of Three Velocity Components in a Model Furnace With

and Without Combustion", Fifteenth Symposium (International) on

Combustion, The Combustion Institute, Pittsburgh, pp. 553-559,

1974.

95. Lockwood, F.C., Ei-Mahallawy, F.M., and Spaulding, D.B., "An

Experimental and Theoretical Investigation of Turbulent Mixing in

a Cylindrical Furnace", Combustion and Flame, Vol. 23, pp. 283-293,

1974.

96. McDannel, M.D., Peterson, P.R., and Samuelson, G.S., "Species

concentration and Temperature Measurements in a Lean, Premixed Flow

Stabilized by a Reverse Jet", Combust. Sci. and Tech., Vol. 28, pp.

211-224, 1982.

97. E1 Banhawy, Y., Sivasegaram, S., and Whitelaw, J.H., "Premixed,

Turbulent Combustion of a Sudden-Expansion Flow", Combustion and

Flame, Vol. 50, pp. 153-165, 1977.

98. Smith, G.D., Giel, T.V., and Catalano, C.G., "Measurements of

Reactive Recirculating Jet Mixing in a Combustor", AIAA J., Vol.21,

pp. 270-276, 1983.

99. Lightman, A.J., Richmond, R.D., Krishnamurthy, L., Magill,

P.D., Roquemore, W.M., Bradley, R.P., Stutrud, J.S., and Reeves,

C.M., "Velocity Measurements in a Bluff-Body Diffusion Flame",

AIAA-80-1544, 1980.

i00. Magill, P.D., Lightman, A.J., Orr, C.E., Bradley, R.P., and

Roquemore, W.M., "Simultaneous Velocity and Emission Measurements

in a Bluff-Body Combustor", AIAA-82-0833, 1982.

101. Schefer, R.W., Namazian, M., and Kelly, J., "Velocity

Measurements in a Turbulent Nonpremixed Bluff-Body Stabilized

Flame", AIAA-87-1349, 1987.

102. Sivasegaram, S. and Whitelaw, J.H., "Oscillations in

Axisymmetric Dump Combustors", Combust. Sci. and Tech., Vol. 52,

pp. 413-426, 1987.

103. Sivasegaram, S., Thompson, B.E., and Whitelaw, J.H., "Acoustic

Characterization Relevant to Gas Turbine Augmentors", J___.

Propulsion, Vol. 5, pp. 109-115, 1989.

104. Charles, R., "Detailed Data Set: Velocity and Temperature

Measurements in the Axisymmetric Can Combustor for a Parametric

Variation in Inlet Conditions", University of California, Irvine,

UCI-ARTR-87-6, 1988.
274

105. Pan, J.C., "Laser-Diagnostic Studies of Confined Turbulent

Premixed Flames Stabilized by Conical Bluff Bodies", Ph.D. Thesis,

University of Dayton, 1991.

106. Heitor, M.V., and Whitelaw, J.H., "Velocity, Temperature, and

Species Characteristics of the Flow in a Gas-Turbine Combustor",

Combustion and Flame, Vol. 64, pp. 1-32, 1986.

107. Bicen, A.F., Tse, D.G.N, and Whitelaw, J.H., "Combushion

Characteristics of a Model Can-Type Combustor", Combustion and

Flame, Vol. 80, pp. 111-125, 1990.

108. Bulzan, D.L. "Structure of a Swirl-Stabilized Combusting

Spray", Journal of Propulsion and Power, Vol. ii, pp. 1093-1102,

1995.

109. Srivasan, R.,"Aerothermal modeling program Phase I final

report", NASA CR-168243, 1983.

Ii0. Sloan, D.G., Smith, P.J., and Smoot, L.D.,"Modeling of swirl

in turbulent flow systems", Proq. Enerqy Combust. Sci., Vol. 12,

163-250, 1986.

Iii. Shih, T-H, Zhu, J., and Lumley, J.L., "Modeling of wallbounded complex flows and free shear flows," NASA TM 106513, 1994.

112. McDonell, V.G., Adachi, M., and Samuelsen, G.S., "Structure of

Reacting and Non-Reacting Swirling Air-Assisted Sprays", Combust.

Sci. and Tech., Vol. 82, pp. 225-248, 1992.

113. E1 Banhawy, Y., and Whitelaw, J.H., "Calculation of the Flow

Properties of a Confined Kerosene-Spray Flame", AIAA. J., Vol. 18,

pp. 1503-1510, 1980.

114. Edwards, C.F., and Rudoff, R.C., "Structure of a SwirlStabilized Spray Flame by Imaging, Laser Doppler Velocimetry, and

Phase Doppler Anemometry", Twenty-Third Symposium (International)

on Combustion, The Combustion Institute, Pittsburgh, pp.1353-1359,

1990.

115. Bachelo, W.D., Brena de la Rosa, A., and Rudoff, R.C., "

Diagnostics Development for Spray Characterization in Complex

Turbulent Flows", AIAA 88-0236, 1988.

116. Zhu, J.Y., Rudoff, R.C., Bachalo, E.J. and Bachalo, W.D.,

"Number Density and Mass Flux Measurements Using the Phase Doppler

Particle Analyzer in Reacting and Non-Reacting Flows", AIAA 93-

0361, 1993.

117. Presser, C., Gupta, A.K., and Semerjian, H.G.,"Aerodynamic

Characteristics of Swirling Spray Flames: Pressure-Jet Atomizer",

Combust. and Flame, Vol. 92, pp. 25-44, 1993.

Q-

?

!FIl!!
"v

275

118. Ghaffarpour, M., and Chehroudi, B., "Experiments on Spray

Combustion in a Gas Turbine Model Combustor", Combust. Sci. and

Tech., Vol., 92, pp. 173-200, 1993.

119. Lewis, M. H., and Smoot, L. D. "Turbulent Gaseous Combustion

Part I: Local Species Concentration Measurements", Combustion and

Flame, Vol. 42, pp. 183-196, 1981.

120. Smith, P. J., and Smoot, L. D. "Turbulent Gaseous Combustion

Part II: Theory and Evaluation for Local Properties", Combustion

and Flame, Vol. 42, pp. 277-285, 1981.

121. Jones, W. P., and Wilhelmi, J., "Velocity, Temperature and

Composition Measurements in a Confined Swirl Driven Recirculating

Flow", Combust. Sci. and Tech., Vol. 63, pp. 13-31, 1989.

122. Jones, W. P., Tober, A., "Velocity, Composition and Temperature Fields in an Axisymmetric Model Combustor", International

Symposium on Applications of Laser Anemometry to Fluid Mechanics,

Lisbon, Portugal, p. 3.9, 1988.

123. Gould, R. D., Stevenson, W. H., and Thompson, H. D., "Turbulence Characteristics of an Axisymmetric Reacting Flow", NASA CR

4110, Feb., 1988.

124. Magnussen, B.F., and Hjertager, B.H., "On Mathematical

Modeling of Turbulent Combustion with Special Emphasis on Soot

Formation and Combustion", Sixteenth Symposium (International) on

Combustion, The Combustion Institute, Pittsburgh, pp. 719-729,

1978.

125. Magnussen, B.F., "The Eddy Dissipation Concept for Turbulent

Combustion Modelling. Its Physical and Practical Implications",

Division of Thermodynamics, University of Trondheim, Norway, 1990.

126. Gran, I. R., Melaaen, M. C., and Magnussen, F., "Numerical

Simulation of Fluid Flow and Combustion in Gas Turbine Combustors",

ASME 93-GT-166, International Gas Turbine and Aeroengine Congress

and Exposition, 1993.

127. Takahashi, F., Vangsness, M.D., Belovich, V.M., "Conditional

LDV Measurements in Swirling and Non-Swirling Coaxial Turbulent Air

Jets for Model Validation," AIAA 92-0580, 1992.

128. Durbin, M.D., Vangsness, M.D., Ballal, D.R., and Katta, "Study

of Flame Stability in a Step Swirl Combustor," ASME 95-GT-III,

1995.

129. Takahashi, F., Vangsness, M.D., Durbin, M.D., and Schmoll,

W.J., "Structure of Turbulent Hydrogen Jet Diffusion Flames," ASME

HTD-Vol. 317-2, Proceedings of the ASME Heat Transfer Division, pp.

183-193, 1995.
Equation

GAS

Appendix 1

PHASE SOURCE TERMS

S_,g

276

Axial momentum OP O I OU_ 1 O iy OV\

Radial momen turn @p 2_effV+ a I aUi

+ I a i aw+2__ -

Tangential momen rum

-

Fuel mass fraction

Oxygen mass fraction -woa_

Energy wto_t

Turbulent kinetic energy

Di ssipa ti on

Concen tra ti on fl uc tua ti on

_IG-pC

C1C.Gkp-C2p¢2/k

|

-/

_=1.44

Turbulent parameters

_=i. 92 _=0.09 Cgi=2.8 Cg2=2 •0

I[I I
Equation

Con tinui ty

SPRAY

Appendix 2

SOURCE TERMS FOR GRIDCELLS

S,,id(Vol)

_nkmk

277

Axial momentum

3 dU,

Radial momentum

3 dVk_

Tangen ti al momen turn

dWk

Fuel mass fraction

Energy nwn k (h_-l_, e::)
REPORT DOCUMENTATION PAG E Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing Instructions. searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson

Davis Highway, Suite 1204. Arlington. VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

:1. AGENCY USE ONLY (Leave blank) [ 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

t., Technical Memorandum November 1998

4. TITLE AND SUBTITLE

Improved Modeling of Finite-Rate Turbulent Combustion Processes

in Research Combustors

6. AUTHOR(S)

Thomas J. VanOverbeke

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Lewis Research Center

Cleveland, Ohio 44135-3191

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

5. FUNDING NUMBERS

WU-523-26--33--00

!8. PERFORMING ORGANIZATION

REPORT NUMBER

E-11430

10. SPONSORING/MONITORING

AGENCY REPORT NUMBER

NASATM--1998-208823

11. SUPPLEMENTARY NO[E5

Responsible person, Thomas J. VanOverbeke, organization code 5830, (216) 433-5867.

12a. DI:_THIBUTION/AVAILAB! LITY STATEMENT

Unclassified - Unlimited

Subject Category: 07 Distribution: Nonstandard

This publication is available from the NASA Center for AeroSpace Information, (301) 621--0390.

12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

The objective of this thesis is to further develop and test a stochastic model of turbulent combustion in recirculating flows. There is _,requirement to increase the accuracy of multi-dimensional combustion predictions. As turbulence affects reaction rates, this interaction must be more

accurately evaluated. In this work a more physically correct way of handling the interaction of turbulence on combustion is further developed

and tested. As turbulence involves randomness, stochastic modeling is used. Averaged values such as temperature and species concentration are

found by integrating the probability density function (pdf) over the range of the scalar. The model in this work does not assume the pd_ type,

but solves for the evolution of the pdf using the Monte Carlo solution technique. The model is further developed by including a more robust

reaction solver, by using accurate thermodynamics and by more accurate transport elements. The stochastic method is used with Semi-Implicit

Method for Pressure--Linked Equations. The SIMPLE method is used to solve for velocity, pressure, turbulent kinetic energy and dissipation.

The pdf solver solves for temperature and species concentration. Thus, the method is partially familiar to combustor engineers. The method is

compared to benchmark experimental data and baseline calculations. The baseline method was tested on isothermal flows, evaporating sprays

and combusting sprays. Pdf and baseline predictions were performed for three diffusion flames and one premixed flame. The pdf method

predicted lower combustion rates than the baseline method in agreement with the data, except for the premixed flame. The baseline and

stochastic predictions bounded the experimental data for the premixed flame. The use of a continuous mixing model or relax to mean mixing

model had little effect on the prediction of average temperature. Two grids were used in a hydrogen diffusion flame simulation. Grid density

didn't effect the predictions except for peak temperature and tangential velocity. The hybrid pdf method did take longer and required more

memory, but has a theoretical basis to extend to many reaction steps which cannot be said of current turbulent combustion models.

14. SUBJECT TERMS

Combustion; Modeling; Stochastic; Numerical

17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION

OF REPORT OF THIS PAGE

Unclassified Unclassified

NSN 7540-01-280-5500

19. SECURITY CLASSIRCATION

OF ABSTRACT

Unclassified

15. NUMBER OF PAGES

28,1.

16. PRICE CODE

AI3

20. LIMITATION OF ABSTRACT

Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

_1__

lr[li