EPA/600/7-87/0
December 1987
FUNDAMENTAL COMBUSTION RESEARCH
APPLIED TO POLLUTION FORMATION
Volume IV: Engineering Analysis
Prepared By:
C. J. Kau, M. P. Heap, W. R. Seeker, T. J. Tyson
Energy and Environmental Research Corporation
18 Mason
Irvine, California 92718
EPA Contract 68-02-2631
EPA Project Officer W. Steven Lanier
Air and Energy Engineering Research Laboratory
Research Triangle Park, North Carolina 27711
AIR AND ENERGY ENGINEERING RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
RESEARCH TRIANGLE PARK, NC 27 711
-------
TECHNICAL REPORT DATA
(Pleate read humicrions on the reverse before completing)
1. REPORT NO. 2.
EPA/600/7-87/—27
3. RECIPIENT'S ACCESSION NO.
PB8S
4. title and subtitle
Fundamental Combustion Research Applied to Pollu-
tion Formation; Volume IV. Engineering Analysis
6. REPORT DATE
December 1987
8. PERFORMING ORGANIZATION CODE
7. AUTHORIS)
C.J. Kau, M. P. Heap, W. R. Seeker, and
T. J. Tvson
B. PERFORMING ORGANIZATION REPORT NO.
0. PERFORMING ORGANIZATION NAME ANO-ADORESS
Energy and Environmental Research Corporation
18 Mason
Irvine, California 92718
10. PROGRAM ELEMENT NO.
11. CONTRACT/GRANT NO.
68-02-2631
12. SPONSORING AGENCY NAME AND ADORESS
EPA, Office of Research and Development
Air and Energy Engineering Research Laboratory
Research Triangle Park, NC 27711
13. TYPE OF REPORT AND PERIOD COVERED
Final; 1/77 - 9/80
14. SPONSORING AGENCY CODE
EPA/600/13
is. supplementary notes aEERL project officer W. Steven Lanier is no longer with the
Agency. For details contact William P. Linak, Mail Drop 65, 919/541-5792. Volume
I of this series is EPA-600/7-85-048.
i6.abstract The report is the fourth and final volume in a series documenting research
activities under the EPA's Fundamental Combustion Research (FCR) program appl-
ied to pollution formation. The FCR program had three major objectives: (1) to gen-
erate an understanding of combustion behavior necessary to help develop control
strategies to minimize NOx emissions from stationary sources, (2) to develop engi-
neering models which would allow effective utilization of a large body of fundamental
information in the development of new NOx control techniques, and (3) to identify
critical information necessary for low-NOx combustor development and to generate
it in a time frame consistent with the needs of EPA's technology development pro-
grams. The report documents FCR program efforts to develop engineering design
models. Six modeling studies were conducted in support of FCR: (1) a computer pro-
gram for general flame analysis, (2) mathematical modeling of microscale combus-
tion of a coal particle, (3) data analysis through inverse techniques. (4) microscale
mixing in turbulent diffusion flames, (5) a numerical model for a two-phase one-
dimensional flow reactor, and (6) generation of elliptic-code test cases.
17. KEY WORDS AND DOCUMENT ANALYSIS
a. DESCRIPTORS
b.IDENTIFIERS/OPEN ENDED TERMS
c. COSATl Field/Croup
Pollution Analyzing
Mathematical Models
Combustion Coal
Research Turbulent Diffusion
Nitrogen Oxides
Flames
Pollution Control
Stationary Sources
13B 14 B
12 A
21B 21B
14F 14G
07B
18. distribution statement
Release to Public
19. SECURITY CLASS (This Report)
Unclassified
21. NO. OF PAGES
210
20 SECURITY CLASS (This page)
Unclassified
22. PRICE
EPA Farm 2220-1 (*-73)
i
-------
NOTICE
This document has been reviewed in accordance with
U.S. Environmental Protection Agency policy and
approved for publication. Mention of trade names
or commercial products does not constitute endorse-
ment or recommendation for use.
ii
-------
ABSTRACT
This document is the final volume in a four volume series describing
research work performed under the Fundamental Combustion Research
-------
TABLE OF CONTENTS*
Section page
ABSTRACT ii
EXECUTIVE SUMMARY iv
PART 1. A COMPUTER PROGRAM FOR GENERAL FLAME ANALYSIS 1-i
PART 2. MATHEMATICAL MODELING OF MICROSCALE COMBUSTION
OF A COAL PARTICLE 2-i
PART 3. DATA ANALYSIS THROUGH INVERSE TECHNIQUES 3-i
PART I*. MICROSCOPIC MIXING IN TURBULENT DIFFUSION FLAMES k-t
PART 5. NUMERICAL MODEL FOR TWO-PHASE ONE-DIMENSIONAL FLOW
REACTOR 5-i
PART 6. GENERATION OF ELLIPTIC-CODE TEST CASES 6-i
FIGURES"
FIGURE 1. SOLID LIQUID FUEL PROCESSING IN TURBULENT DIFFUSION
FLAMES vi
FIGURE 2. FCR PROGRAM STRUCTURE viii
(*)Each port of this volume is accompanied by detailed table of contents, list
of figures, and list of tables.
Preceding page blank
v
-------
EXECUTIVE SUMMARY
The Fundamental Combustion Research Applied to Pollution Control Program
(FCR) had three major objectives. These were:
• To generate an understanding of combustor behavior necessary to aid the
Combustion Research Branch (CRB) in development of control strategies to
minimize NO^ emissions from stationary sources.
• To develop engineering models which would ollow effective utilization of
a large body of fundamental information in the development of new NO^
control techniques.
• To identify critical information necessary for low-NO^ combustor
development and to generate it in a time frame which was consistent with
the needs of the CRB technology development programs.
The overall goal of the CRB, and hence FCR, was to provide the technology for
maximum control of NO^ emissions from stationary sources. The intent of FCR was
to develop a focused program directed at Important Issues relevant to the CRB so
that the results could be applied in a relatively short time period. The program
plan emphasized well defined priority target areas, and the ensuing research
effort was always guided towards engineering solutions to specific problems.
There was a conscious effort to avoid projects that did not provide information
critical to the needs of the CRB. Consequently, it was necessary to isolate
relevant Issues and direct the program away from the studies of physics and
chemistry of combustion which were not relevant to the goals of the CRB.
Fifty-six percent of the annual emission of nitrogen oxides from stationary
sources in this country emanate from the single combustor category of boilers
firing coal or oil. The dominant characteristics of combustors can be described
as follows:
• The flames are large, with energy release zones whose dimensions are on
the order of tens of feet.
e The time mean motion is in steady state.
• Pressure is atmospheric.
• Rodiatian heat transfer dominates exchange from the heat release zone to
the cold confining walls.
• Mass and thermal transport Is turbulent.
-------
• Fuel/air contacting occurs by diffusion and particle penetration.
• Fuels are injected as solids or liquids giving rise to two-phase
transport and both homogeneous and heterogeneous reactions.
• The high fuel-bound nitrogen content of coal and residual fuel generally
requires that NO^ emissions be minimized by the application of staged
combustion.
Thus, the primary initial objective of FCR was to Immediately establish a
subcontractor-oriented program of fundamental investigations that focused on the
control of NO^ from large, confined, one-atmosphere turbulent diffusion flames
burning heavy residual oil and pulverized coal. This section describes the FCR
programming structure and describes the goals of Volume
-------
Burner
I'tfterWhv
Circulating Products
77777
777777
Near Field Dominated Burner
Stahilized flame
Combustion
Ai r
rar Held Domi nan I
Long Flame
External Recirculating
Products
HEAT TRANSFER
T ue 1
Particles/ q q
Contacting
Mixing
Gas Phase Reaction
Pyrolysis
Droplets
°°n°
O O
Thermal
Decomposition
9.
Oxidant
Contacting
I
Products
Heterogeneous Reaction
Char
Figure 1 Solid Liquid Fuel Processing in Turbulent Diffusion names.
-------
liquid droplets of fuels in turbulent diffusion flames is Illustrated in Figure
1. The fuels first decompose, giving a vapor and a solid product. The vapor
may be ejected from the solid as a single jet from a blowhole or may vaporize
uniformly from the particle or droplets. This vapor will then either undergo
pyrolysls, producing soot particles, or undergo gas phase oxidation. The solid
remaining after devolatilization of coal is normally referred to as char. The
char ond soot moy then undergo heterogeneous reaction which includes oxidation
(burnout) ond possible reaction with gaseous nitrogen species. Throughout this
sequence of events the droplets and particles of fuel and their products must
be brought into contact with oxidant and the high temperature combustion
products In order to cause decomposition and subsequent reaction. This
contacting Is important on a macroscale because it dictates the major flame
characteristics (i.e., near-field dominated burner stability or far-field
dominated jet flame), and on a microscale because it ultimately dictates the
production of pollutant species. FCR has been primarily concerned with the
phenomena associated with conversion of nitrogen in droplets and particles
of fuel to the ultimate nitrogen oxide released from the combustion chamber.
FCR was divided into three program areas. These were concerned with (1)
transport processes in reacting systems; (2) gas phase chemistry; and (3) the
physics and chemistry of two-phase systems. Figure 2 presents an overall view of
the program structure, which Is divided into major program areas, specific program
elements and support areas. This program structure was planned to lead to two
major program outputs. These were;
• A description of the chemical limits of NO production in order to
ascertain the lower bounds of both fuel NO and thermal NO production
under a series of process constraints which were not limited in any
way by fuel/air contacting.
• A description of fuel NO formation in turbulent diffusion flames for
gas, liquid and pulverized coal systems.
This overall program was planned as a combined ln-house and subcontract
project effort. EER took responsibility for progrom planning, management and
synthesis of the overall program, and subcontracted separate projects within the
ix
-------
trMjpurt fa
Austin? Syitaw
Fhyifca:! «n* OiwHt*y
at Tmo ftwi«
¦Syitjo*
Major Prcgrw krtu
G4vPti«w. ihtuittry
RMCttr r"I«n
DaconposUlon
ProqrM Cltmnts
Clffvijion fliats
Prsfraa Support Ar»a»
F1 am« todulfs
5
0«icrl3tion of Fu»l
NQ Fflnwtlan
1n Turtulfnt
Olffuilon Him
Chaglcil Limits af
NO Production
NO. Control
Figure 2. FCR Program Structure.
x
-------
program elements to various organizations. Subcontracts accounted for 70 percent
of the program and were used to ensure that FCR had the benefit of the best
scientific talent available. Every effort was made to use subcontractors who had
the necessary experience, and in many cases, equipment. The main synthesis role
was accomplished by EER using analytical tools developed in major program areas.
The reasoning behind the division of the total program into three major
program areas can most easily be described by consideration of the phenomena
occurring to solid and liquid fuels in turbulent diffusion flames. When a pulver-
ized coal particle is heated and decomposes, the fuel-bound nitrogen component is
divided between that which stays with the solid (char) and that which is evolved
with the volatile gases. The processes that lead to the division of fuel nitrogen
are governed by the physics and chemistry of two-phase systems, as are the subse-
quent reactions of the solid phase. Thus, two of the program elements in the
two-phase program area were concerned with the thermal decomposition of the fuel
and heterogeneous NO/char reactions. The fate of the volatiles depends upon the
chemistry of nitrogenous species during combustion of the gas phase pyrolysate
composed primarily of hydrocarbons. This whole area was treated in the gas phase
chemistry program area. The primary goal of the gas phase chemistry program was
to produce a Kinetic mechanism which described the fate of fuel nitrogen in
gaseous mixtures which were llKely to be found in pulverized coal and residual
fuel oil flames. Thus, the gas phase chemistry area was divided into two program
elements: one experimental and one aimed at mechanism development using computer
simulations. A considerable effort was expended in the development of two-phase
reactor experiments. These experiments coupled gas phase chemistry with fuel
decomposition.
All of the above neglected the physics of turbulent transport. In flames the
particles must be mixed with high-temperature gases before they can decompose, and
oxidant must be mixed with the volatile fuel fractions before they can reoct.
These fuel/air contacting processes were treated in the program area defined as
transport processes in reacting systems. Three projects were carried out in this
major program area. Theses were concerned with coupling of Kinetics to turbu-
lence, fuel injection systems and turbulent diffusion flames.
Although each of the projects in the various program areas were planned to
xi
-------
provide information which was of value In itself and was of direct use for data
interpretation, the overall program flow was directed towards the development
and verification of engineering tools. These tools could be used in the
attainment of the two major gools concerning fuel NO formation in turbulent
diffusion flames. It is not possible to show the overall interaction of the
various program areas and program elements because of the additional complexity
that would be added to Figure 2. However, the overall program placed great
importance upon the use of analytical tools to help synthesize the program.
Indeed, all the program elements can be considered to provide input to allow
analytical tools to be developed. The application of these tools then provides
the output of the FOR program.
It was recognized from the outset that within the time frame of FCR the
development of a unified mathematical model for NO^ production in pulverized coal
flames was unattainable. The major modeling effort in FCR concentrated upon the
development of a framework of semi-empirical modular models. Thus, a complex
system could be modeled using a collection of limit-case elements linked together
by means of empirical knowledge of the exchange of heat and mass transfer between
these elements. In its simplest application, modular elements could be used to
analyze data generated in several of the program elements involving reactor
experiments. The modules describing fuel decomposition could be applied to both
reactors and turbulent diffusion flames. Recognizing how complex a
pulverized-coal system is, the limit-case approach was deemed ta be the most
appropriate in analyzing fundamental data for application to real systems. In
limit-case studies certain phenomena are assumed to dominate the process while
other phenomena become suppressed, thus allowing the implications of specific
phenomenological behavior to be examined. For example, limit-case studies could
be used to define the times ond temperatures needed to reach a certain total
fixed-nitrogen content in a fuel-rich combustion system, or to define the extent
of heterogeneous NO reduction in pulverized coal flames.
Although FCR placed considerable emphasis on the provision of information
which was necessary to provide short-term solutions, long-term research was not
entirely neglected. FCR pioneered the use of holography and two-color pyrometry
as tools to study the thermal decomposition and combustion of fuels. Activities
in coupling of kinetics ta turbulence are beginning to show considerable promise
xii
-------
for the future. Effort was initiated to evaluate numerical techniques which could
be used in the ultimate models of furnace performance. Projects such as these
ensured that the CRB had a balanced fundamental combustion research program.
Having established an overall program structure it was necessary to formulate
individual research efforts directed towards meeting the overall program goals.
The following section describes the research efforts that lead to Volume
This volume is concerned with Engineering Analysis and Transport Phenomena.
The research efforts pursued under this Major Program Area focused on the
development of mathematical models, with development efforts closely coupled to
numerous experimental studies. The majority of the modeling work was performed at
EER and directed towards development af the overall madular model discussed
earlier. Models include a description of the following phenomena:
• Particle thermal history.
• Finite rate thermal decomposition.
• Heterogeneous chemical behavior.
• Finite-rate gas-phase chemistry.
• Fluid mechanics and ballistics of turbulent transport in two-phase
systems.
• Complex flow fields associated with the near-field of highly
swirl-stabilized combustors.
• The interaction between turbulent unmixedness and chemistry.
These modeling efforts have resulted in a hierarchy of analysis codes.
Treatment of the fully-coupled finite-rate processes associated with particle
heatup, thermal decomposition, heterogeneous reactions and gas phase chemistry was
emphasized. The codes can handle with ease and efficiency as many as 200 coupled
processes.
In the course of developing an understanding of the dominant mechanisms
governing complex and poorly understood behavior, it is often necessary to start
with a large coupled system and, through sensitivity analysis, reduce It sa it
contains first-order phenomena only. A good example of a situation requiring such
an approach is NO^ generation In p.c. flames. The complexities of HCN, NH^ and
chemistry, in conjunction with finite-rate devolitillzation, defies a simple
xiii
-------
Intuitive understanding and requires complex analysis in order to generate such
intuition.
This volume describes the FCR effort to develop engineering analysis tools.
xiv
-------
PART 1
A
COMPUTER PROGRAM
FOR
GENERAL FLAME ANALYSIS
Prepared by:
C. J. Kau
and
T. J. Tyson
1-1
-------
TABLE OF CONTENTS
Section Paqe
I.0 INTRODUCTION 1-1
2.0 GOVERNING EQUATIONS IN PHYSICAL COORDINATES 1-2
3.0 GOVERNING EQUATIONS IN STREAM FUNCTION COORDINATES 1-5
k.Q GENERAL FORM OF GOVERNING EQUATIONS AND BOUNDARY
CONDITIONS 1-7
<*.1 General Form 1-7
k.2 Boundary Conditions for Symmetric Plane/Axis 1-9
<*.3 Boundary Conditions for Catalytic Surface 1-11
5.0 TRANSPORT PROPERTIES 1-13
5.1 Laminar Transport Properties 1-13
5.2 Turbulent Eddy Viscosity Models 1-16
6.0 CHEMICAL REACTIONS AND THERMOCHEMICAL PROPERTIES 1-19
6.1 Chemical Reaction Rate Equations 1-19
6.2 Thermochemical Data 1-20
7.0 FLAME RADIATION 1-21
8.0 NUMERICAL METHOD OF SOLUTION 1-23
9.0 APPLICATIONS 1-27
9.1 Free Jet Diffusion Flame 1-27
9.2 Confined, Co-flowing LBG Diffusion Flame and
Fuel Nitrogen Conversion 1-33
9.3 Flat Flame Simulation 1-39
9.^ Constant Strain Rate Laminar Diffusion Flame 1-39
10.0 REFERENCES 1-50
II.0 LIST OF NOMENCLATURE 1-52
APPENDIX A. Block Tridiagonal Matrix Equation Solver. . . . 1-55
APPENDIX B. Gas-Phase Reaction Mechanism 1-58
1-11
-------
LIST OF FIGURES
F iqure Page
1 Planck Mean Absorption Coefficients at One Atmosphere
Totol Pressure 1-22
2 Initial Velocity and Initial Temperature Distribution. . . 1-28
3 Isothermal Contours for Laminar H^/Air Diffusion Flame . . 1-29
4 Velocity Components U(r) and V(r) at X/R = 5.0 1-30
5 Distribution of Stable Species and Temperature at
X/R - 5.0 1-31
6 Distributions of Radical Species in on H /Air Diffusion
Flame 1-32
7 EER Diffusion Flame Burner 1-34
0 Influence of CH^ Content on NH^ Conversion in LBG
Diffusion Flames 1-36
9 Influence of CH Content on NO Conversion in LBG
Diffusion Flames 1-37
10 Influence of Temperature on NO Conversion in Diffusion
Flames—10 Percent CH 1-38
11 Flame Profiles of Biordi et al. (1975). Solid Lines Are
Data; Dashed Are Model 1-40
12 Profiles Through the Flame af Biordi et al. (1975) .... 1-41
13 Schematic of Opposed-Jet Diffusion Flame 1-43
14 Velocity Profile of CH /Air Opposed-Jet Diffusion Flame
(a-250 1/sec)
15 CH and 0- Profiles of CH /Air Opposed-Jet Diffusion
Flome (a-250 1/sec) 1-45
16 Flame Temperature vs. Strain Rate for CH /Air Opposed-Jet
Diffusion Flame 1-46
17 Flome Thickness vs. Strain Rate for CH /Air Opposed-Jet
Diffusion Flame 1-48
18 Stagnation Point Displacement vs. Strain Rate for CH^/Air
Opposed-Jet Diffusion Flame 1-49
1-1H
-------
LIST OF TABLES
Toble Poqe
1 Defintion of C, C , C, and C„, and Q and R in Eq. (4-1) . 1-8
0 1 2
2 Coefficient Functions B^ and B^ for Eq. (4-5) 1-10
3 Force Potential Parameters and Diffusion Factors .... 1-17
1 -1 V
-------
1.0 INTRODUCTION
A computer model which is capable of predicting or analyzing various types of
premixed or diffusion flames is reported. Detailed mathematical description of
the problems involved and the numerical technique employed in the solution scheme
are described.
The generality of the model allows computation for systems of various config-
urations; e.g., one-dimensional time-dependent planar/spherical, steady
two-dimensional planar and steady axisymmetric nonrecirculating reacting flaw
systems.
Physical processes such as laminar unequal species diffusion, radiation,
flame holder recombination effects and heat loss are treated extensively. Several
phenomenological turbulent eddy viscosity models based on mixing length theory are
also incorporated into the code. Kinetically, the code is capable of treoting up
to 200 two- or three-body basic reactions and up to 52 species.
A linearized implicit finite difference network is employed. Thus, except
for cross-stream velocity, all dependent variables of all grid points at the same
coordinate line (or at the same integration step) are solved simultaneously in
coupled fashion. Thus, the inversion of a block tridiagonal matrix is required at
each integration step. This numerical scheme is so efficient that the computa-
tional time required for a typical run is only about one-tenth to one-fifth of the
time required for other popular methods such as the implicit/explicit method (Ref.
12) or the operator splitting method (Ref. 13).
For illustrative purposes, detailed analyses of four fundamental types of
laminar flames are presented: flat flame, opposed-Jet diffusion flame, co-flow
diffusion flame, and nonrecirculating confined flame.
1-1
-------
2.0 GOVERNING EQUATION IN PHYSICAL COORDINATES
The set of partial differential equations presented in this section is
capable of describing two-dimensional (or axisyiranetric) laminar (or time averaged
turbulent) boundary layer flow with chemical reaction. These equations are
evolved from full Navier-Stokes equations by incorporating boundary layer assump-
tion, i.e.:
1. Momentum transfer in the main-stream direction is negligible compared
to momentum transfer in cross-stream direction.
2. The momentum equation In the cross-stream direction is negligible
compared to the momentum equation in the main-stream direction.
3. Energy and mass diffusion in the main-stream direction are negligible
compared to diffusion in the cross-stream direction.
The resulting equations are:
Continuity Equation
_3_ (pur*) + (ovra) = 0
9x 6r
Momentum Equotion
„„ ly. + -w 1 3 / a 3u\ .
pu 3x + DV 17 - - 37 + IS 37 (jr 3?) * P9X
Species Equation
9a • 3q . -i a / 3a. \
PJ "ST + <=» — = 75 37 (pD,ra i.
Energy Eouotlon
* Cpi &
i = l
1-2
(2-1)
(2-2)
(2-3)
- (2-4)
-------
Equation of State
P
82.06 T La.
(2-5)
where a is a geometric index. By setting a-0, 1 and 2, the equations can be
obtained in terms of two-dimensional cartesian, polar cylindrical, and spherical
coordinates, respectively. By setting u-1 and ignoring the momentum equation
(Eq. 2-2), the equations become time-dependent and one-dimensional, with the
x-coordinate representing transient time. The nomenclature used in the above
equation is defined at the end of Part 1.
The conservation law which the species equation (Eq. 2-3) must
satisfy is:
which states, by definition, that "sum of mass fractions equals one." For consis-
tency, one of the species equations in Eq. (2-3) shall be replaced by Eq. (2-6).
Another conservation law to be satisfied is that "the chemical reactions
neither create nor destroy mass." Mathematically:
Surrmlng Eq. (2-3) over all species and noting Eq. (2-6) and Eq. (2-7), we obtain
the condition:
Eq. (2-8) Implies that not all diffusion coefficients are independent;
therefore, in actual calculation the dlffusivity of the dominant species is not
specified and Eq. (2-6) is used to replace the species equation of the dominant
species. The fifth term at the right hand side of the energy equation, Eq. (2-if)
(2-6)
(2-7)
(2-8)
1-3
-------
can be rewritten as:
il
Sr
^ "fpi " \ CPR) P°i 3r
(2-9)
which is the consequence of Eq. (2-8). Subscript "R" refers to the dominant
species (usually N^).
1-4
-------
3.0 GOVERNING EQUATIONS IN STREAM FUNCTION COORDINATES
In soma practical situations it is more convenient to write the governing
equations in stream function coordinates via Von Mises Transformation. The stream
function in Von Mises Tronsformation is defined by the following relations
»a 'H) ¦
A
(i)
¦pvr-
(3-1)
(3-2)
which automatically satisfy the continuity equation (Eq. 2-1). The
transformation of differential equations from the physical coordinates (x,r) plane
to the stream function (x , ^) plane requires evaluation of partial derivatives in
I
terms of (x , ip); i.e.,
(w)r * \w).-
r u x
where fh (—| and the streamline angle ;^ is assumed to be small.
Introducing Eqs. (3-1) and (3-2) into Eqs. (2-2), (2-3) and (2-4) and drop-
! I
ping " " from x , we have:
Momentum Equation
3u
3x
p. „ 5U /-L 3 lflr\
dx ^ 3({i I ^ 3'yJ?
pg,
(3-5)
Energy Equation
9T dp . • , ou 3 /, a 4 3T \
DUC n— = U ^ 5 ^ ^ F A -r-r I
D 3x dx Hr .a 3iii \ 3yJ
ouraA /9u\^ , puraA 3T y /c _ r \ oD. -
,,a \3ijJ ,.a 3^^ \ pi Mr pR/ l 3'y
-E h,w(
(3-6)
1-5
-------
Species Equation
da.
PU l_ = pu
3x
1
a o>li
D-rCtA !!il! + **
^ rv. I i
(3-7)
where A
The Von Mises Transformation eliminates the continuity equation and the
cross-stream convective terms in momentum, energy and species equations. These
terms usually cause difficulties in numerical integration processes.
The equation of state in stream function coordinate remains the same as that
in physical coordinate systems.
For a fixed value of x, the r-coordlnate can be related to i> by
numerically integrating Eq. (3-1), giving
1-6
-------
4.0 GENERAL FORM OF GOVERNING EQUATIONS AND BOUNDARY CONDITIONS
4.1 Generol Form
The governing equations in two different coordinates can be combined and
written in one general form:
Cf=CoftCl £ (C2 f) + G + R '4_1)
where F -a^ or T or U. When y = r, the system is in physical coordinates, while
when y = i!/ the system is in stream function coordinates. The coefficients C, C^,
C1 and C^, and functions Q and R for each coordinate system are listed in Table 1.
Equation (4-1) is solved along with the equation of state (in differential form)
to give
dp
- da. - 41 + &
Z_-i i T a
(4-2)
and the species equation of reference species R Is replaced by Eq. (2-6) (In
differential form):
do, = if E da, (4-3)
When the system is in physical coordinates, the continuity equation (2-1) has
to be integrated for cross-stream velocity v. It is noted that Eqs. (4-2) and
(4-3) can also be written in general form (i.e., Eq. 4-1). The coefficient
functions for the general form of Eqs. (4-2) and (4-3) are also listed in Table 1.
The system of equations to be solved is a parabolic system of partial differential
equations which require prescription of initial (x - 0) and boundary conditions
pertinent to the individual flow problem at hand. The simplest and the most
common types of boundary conditions which have been provided in this computer
program are:
1. Zero order boundary condition: The Known value of the dependent
variable is specified at the boundary.
1-7
-------
Table 1. Definition of C, CQ, Cj and C2> and 0 and R in Eq. (4-1)
(a) Physical Coordinate System
F
C
Co
C1
C2
0
R
u
T
ai
aR
P
(b) S
pu
pUCp
pu
1
1
>tream F
- (pv)
- (pvCp)
+ T"pD.(C . -!!i cPB)3ai
Z-r 1 P1 Mr PR 3y
- (pv)
unction Coordinate Systen
i/ya
l/ya
l/ya
is
My"
Aya
pD^y"
- L
i
Wi
1 T m dui
mr rn 1 d*
+ I d£.l
P d*J
dp
p9x - gr
¦fc ~ «r ~ .(|ff
F
C
Co
C1
C2
0
R
u
T
a -
l
aR
P
pu
PUCp
pu
1
1
CPr)^
pu/y01
pu/ya
pu/ya
Arap
AruX
ArapDi
-V>.w.
l l
Wi
-'Li.
\ I^R ' dx
f Mvdai 1 CJT 1 dpi
p L ^ dx ~ T dx + P dxj
dp
pg* " dx
dp . , .2 /3u\2
\ * A " (af)
-------
2. First order boundary condition: The first derivative (with respect
to y-coordinates) or the dependent variable is specified.
3, Second-order boundary conditions: The second-order derivative (with
respect to y-coordinate) of the dependent variable is specified.
In addition to these types of boundary conditions, the program also
provides the choice of the following two types of frequently encountered
boundaries in flame problems:
1. Symmetric plane/axis.
2. Catalytic surface.
k.2 Boundary Conditions for Symmetric Plane/Axis
On the plane/axis of symmetry (generally, y-0) we have:
Due to numerical difficulties one cannot apply this condition directly to the
boundary. The symmetric conditions must be derived from governing equations. By
noting that (L'Hopital's rule):
and exerting some mathematical efforts, the governing equation at symmetric
plane/axis (y-0) can be derived and reduced to the following general form:
C Hf = + B2 77 'or y = 0 (4-5)
where coefficient function C is given in Table 1 and coefficient functions B1 and
Bg are listed in Table 2.
(4-4)
1 i mi t
y o
1-9
-------
Table 2. Coefficient Function B-j and for Eq. (4-5)
F
for (x,r) all a
for (x.ij/) and a = 0
for and a / 0
B1
B2
B1
B2
B1
B2
ai
T
u
*1
^+
-------
To compute p and ct R at the boundary, Eqs. (2-5) and (2-6) are used.
k. 3 Boundary Condition for Catalytic Surface
The surfaces involving chemical reaction and heat and mass transfer (such as
flameholder, burner surface) are commonly encountered in analysis of laboratory or
industrial flames. The boundary conditions on these surfaces can be derived from
the principles of conservation of heat and mass on the surface.
For plane surfaces in physical coordinates, the boundary conditions for
reacting surfaces (say r-0) are:
Species
(».&),
Energy
(* ft), "
In the above, the subscripts "s" and "o" refer to the properties on the surface
and far upstream (where diffusive flux is zero), respectively, w is the
s, 1
production rate of the ith species per unit surface area. hioss denotes enthalpy
loss (such as cooling or heating of reactants before entering the surface) per
unit surface area per unit mass flux m. Energy exchange of any form (such os
surface radiation) is represented by q .
6
For ~ plane surface in stream function coordinates, the boundary conditions
for a reacting surface (say £»0) are:
1 ' ' vcrrr-sec/
"[BbiMVo • ho + hTossJ
+ I>i>s v,-*. (—rH
j \cm -sec /
Species
oDi3r= W("i's - mo] /(pu) (4"8)
1-11
-------
Energy
II
3^
{*[»'
Wo
h + h,
o loss
E(hi
'sWS,
- qs}/(pu)
(4-9)
1-12
-------
5.0 TRANSPORT PROPERTIES
5.1 Laminar Transport Properties
For niulticomponent laminar flow, the mixture viscosity and conductivity are
evaluated from component viscosities and conductivities \ using empirical
formulos of Wilke and Wassiljewa (Ref. 1), respectively. The Wilke formula gives
mixture viscosity as
N
£
1=1
- ~ £ ~,$
j/i 1
gm
cm-sec
(5-1)
Wassiljewa's equation gives mixture conductivity as
N
A= L
i=i
cal
1 + E '-065
J>1 X V
cm-sec- K
(5-2)
where x, is mole fraction and . is defined as
i ij
1J
¦w
Jb (1 +
^Y"
V
(5-3)
The single-component viscosities are computed theoretically from Kinetic theory:
u. = 2.669 x 10
-5
1T
a.2 ^2-2)
l
(5-4)
ond
X. = 1.9891 x 10
(T/M.)Js/oi
2 q{2.2)
(5-5)
(2 2)
where ' is the collision integral.
(2 2)
ft is a function of parameters in
the potential function assumed. The Stockmayer potential function is assumed
1-13
-------
(2 2)
here. Tabulated values of ft , which is a function of T«=KT/e and <5, ore
o
first evaluated by the method of Monchick and Mason (Ref. 2). z /K is minimum
o °
energy of attraction between molecules in K and c is intermolecular distance
when potential energy of interaction is zero and is
in angstroms. 5 is the description of angular dependency of dipole-dipole
interaction energy. K is Boltzmann's constant.
The binary diffusion coefficient between unlike molecules is evaluated from
0.5
D.. = 1.858 x 10
J
-3 .1.5
mi * mj
m.m.
i J
P° ..
IJ 1.1
(TjT
cm
sec
(5-6)
where the binary collision integral ^ is a function of potential parame-
ters e .,/K and £... 3,., e .and b^ . ore evaluated from single component
oij ij ij oij' ij
parameters by following combining laws.
aij " v=(~i + aj»
'01
j "voiVl
(5-7)
(5-8)
a = %(c + o)fc
np n p
1
6
For polar-nonpolar pair of molecules, the combining laws are modified as:
(5-9)
(5-10)
(5-11)
The subscripts p and n refer to the polar and nonpolar constituents.
onp v on op c
For pair of polar molecules, the combining laws for cr^ and c ^ remain the
same as for nonpolar pair. The parameter 6 is evaluated from:
ofJ . (Si4 )*
aij
(5-12)
1-14
-------
The empirical equations for collision integrals, which fit the Monchick and
Mason (Ref. 2) calculation, provided by Neufeld (Ref. 3) et al. are:
„(1.1) _ 1.06036 . 0.19300 . 1.03587
"ij T*0.15610 exp(0.47635T*) exp(1.52996T*)
, 1.76474 + 0-19 6i/ (5 ,3)
exp(3.8941IT*) T* (5_13)
A2.2) _ 1 .16145 , 0.52487 . 2.16178
"ij " T*0.14874 exp(0.77320T*) exp(2.43187T*)
- 6.435 x 10"4 T*0-14874 sin(18.0323 y*-0-7683 . 7.27371 )
2
(5-14)
0.2 6,2
T*
The diffusivity of a single component into a mixture of n components is approxi-
mated by (Ref. 1):
1 -X •
n = ] (5-15)
m)
where x^ is mole fraction.
Estimation of the multicomponent diffusivity requires the calculation of
entire matrix of , which is ane of the most time-consuming computational steps
in the model. Evaluation of diffusivities by the so-called bifurcation approxima-
tion (Ref. 4) provides a less time-consuming but less accurate solution.
The basic assumption of bifurcation method is that the binary diffusion
coefficient can be represented by
D. - « f.f • D (5-16)
"i J "i j rr
where f. and f. are diffusion factors and D is a reference self-diffusion
i J rr
coefficient. For general flame calculation, the reference species, r, is usually
taken to be oxygen. By setting J»r in Eq. (5-16), we obtain:
f. „ llL (5-17)
i Drr
1-15
-------
where f - 1.0 by definition. Using Eq. (5-6), Eq. (5-17) becomes:
f. =
1
a \ 2
rr \
M. + M \ H
i r *
if we take first terms from Neufeld's empirical equation for
ond substitute into above equation, we have:
(1.1)
ij
(5-18)
(E<|. 5-13)
(VK>rr
(Eo/K'ir
0.1561
M, + M
(5-19)
Here, f is found to be independent of temperature and pressure and can be
precalculoted. and ( c 0/K)^r can bo obtained using combining laws given in
Eqs. (5-7) and (5-8). can then be evaluated by Eq. (5-15) when all Dij's are
calculated by Eq. (5-16). The laminar viscosity and conductivity are obtained by
assuming constant Lewi6 and Prandtl numbers.
The force potential parameters for species which usuolly appear in the
combustion products of hydrocorbon fuel are listed in Table 3. The diffusion
factors calculated using Eq. (5-19) are also listed in Table 3. The diffusion
factors, for unstable species for which the potential parameters are not avail-
able, are approximated by:
These value* of f will be improved as soon as their potential parameters become
available.
5.2 Turbulent Eddy Viscosity Models
The attempt to treat detailed chemistry of combustion and turbulence
simultaneously is not only impracticol but also currently unfruitful. For turbu-
lent flow systems, the program provides two phenomenological eddy viscosity
models: Ferri's model (Ref. 5) and Schetz's model (Ref. 6). Both models are
1-16
-------
TABLE 3. FORCE POTENTIAL PARAMETERS AND DIFFUSION FACTORS
i
o-j (A)
e0/k
6
*
+25 3
anxlO cm
*«
"1
***
U
H2
2.827
59.7
.
7.9
2.016
3.69
CO
3.69
91.7
-
19.5
28.011
0.983
ch4
3.78
142.1
-
26.0
16.043
1.09
nh3
3.15
358.0
0.7
22.6
17.032
1.20
Nj
3.798
71.4
-
17.6
28.016
0.973
H20
2.71
506.0
1.23
14.70
18.016
1.31
C02
3.941
195.2
_
26.5
44.011
0.777
0
2.529
218.9
-
-
16.0
1.55
H
2.192
36.8
.
4.0
1.008
6.60
OH
3.147
79.8
4.01
.
17.008
1.35
HO2
4.196
289.3
0.75
20.0
33.008
0.752
ch3
3.79
142.1
-
-
15.035
1.12
ch2
-
.
.
14.027
1.28
CH
-
-
-
-
13.019
1.40
CH2O
4.304
340.4
0.71
-
30.027
0.739
' CHO
-
-
-
-
29.019
1.03
nh2
-
-
-
-
16.024
1 .22
NH
-
-
-
-
15.016
1.25
N
.
.
14.008
1.28
NO
3.489
117.2
0.15
.
30.008
1.00
N02
3.712
404.3
0.028
46.008
0.774
n2o
3.776
248.8
0.01
30.0
44.016
0.797
CHN
3.63
569.1
-
25.9
27.027
0.875
CN
-
-
-
-
26.019
1.06
HNO
3.492
116.7
31.016
0.994
HNCO
-
-
.
43.027
0.934
ch3o
-
-
-
-
31.035
1 .01
H202
4.196
289.3
-
-
34.016
0.746
c2h2
4.033
231.8
-
33.3
26.038
0.849
C2H4
4.163
224.7
-
42.6
28.054
0.806
c?H6
4.443
215.7
44.7
30.070
0.739
h
3.467
106.7
-
16.0
32.0
1.0
NCO
-
-
-
-
42.019
0.938
C2H
-
-
-
-
25.030
1 .07
C2H3
-
-
.
-
27.046
1.04
C2H5
-
-
-
-
29.062
1 .02
*Polan'zabi 11ty of molecules.
**Molecular weight.
• 6H' [SScT'"' CW5)'
f, - ** (for species with no potential parameters
' available).
1-17
-------
locally dependent. Ferri's model is recommended for systems such as jet into
still surroundings or axisymmetric wake. The form of the expression is
p = 0.025bh I 3 u - o u
1 ' n n ' a .
1 u - j u
oo e e
where b1/2 is the value of r at u = 0.5 ( p u + o u ). The subscripts "o" and
' o o e e M
"e" refer to the values at the axis of symmetry (r=0) and the edge of the mixing
layer. Schetz's model has the farm
proportional to the excess mass flow per unit width in the mixing region."
Schetz's model is recommended for coaxial jets with large density gradients across
the mixing zone (e.g., hydrogen Jet into air). Constant turbulent Lewis number
and constant turbulent Prandtl number are assumed for all turbulent flows. Thus,
conductivity is calculated by
and turbulent diffusivities of all species are evaluated by:
where Le and Pr are turbulent Lewis number and turbulent Prandtl number, respec-
tively .
pgue - pu I 2rd r
where r is initial radius of Jet
It Is said that "the turbulent viscosity if
1-18
-------
6.0 CHEMICAL REACTIONS AND THERMOCHEMICAL PROPERTIES
6.1 Chemical Reaction Rate Equations
In generol form, the Kth chemical reaction can be written as:
IX,Kzi ^ ZN;,Kzi <6-')
i i
where N and N'. ore integers representing the stoichiometric coefficients of
11 K 1« K
the ith species in the Kth reaction, Prime denotes products and no prime denotes
reactants. is the chemical symbol (or species name) of the ith species. The
net production rate (in forward direction) of the ith species in the Kth reaction
is given by:
wi,K = (Ni,K ' Ni,K*
ENf K N. .
O i IT i ^
f »K J 1 aj
y" n 1
Kf K f "j.K N.
Jc p n a- 1,J
e,K i J
Ek (6-2)
where the forward reaction rate has Arrhenius form. The equilibrium constant,
K^ , is calculated by:
K ,'- 1 'Ni,K " N,''K)
(6-3)
where AG^ is net production af Gibbs free energy and is defined as:
ENi,KGi - £Ni,KGi ¦ (6-1)
i i
The term E in Eq. (6-2) is the modification of reaction rate due to a third body
K
and is defined as:
Ej( = ^ for three body reoctions (6-5)
= 1 for all others
where E is the third body efficiency of the ith species in the Kth reaction.
K i 1
1-19
-------
Th8 overall production rate of the ith species is the sum of w
i, K
Over all K;
i.e.
w = y w
i ^ i,K
K
(6-6)
6.2 Thermochemicol Data
The thermochemical properties (specific heat, enthalpy and Gibbs free energy)
for each species are taken directly from JANAF Thermochemical Table (Ref. 10) in
tabular form as a function of temperature. Linear interpolation is used to
compute the properties at local temperature.
1-20
-------
7.0 FLAME RADIATION
The major radiating species in a hydrogen flame is water vapor. For a
hydrocarbon flame, carbon dioxide and carbon monoxide are olso omong the major
radiating species. For the optically thin limit ( photon mean free path is much
greater than the characteristic dimension of the flame), the Planck mean
absorptance coefficient of species i can be obtained from gas emittonce data (Ref.
7):
Kpi = 4 ("L1) L+o
where L is the characteristic dimension of the flame, and is the omittance of
91
the ith species. The values of for H^O, CO^, and CO calculated by Kelly and
Kendall (Ref. 8) based on omittance data of Edwards and Balakrishnan (Ref. 9) are
shown in Figure 1. The Planck mean absorptance for a mixture can be estimated by:
K = £P-K •
P V l pi
where P^ is partial pressure of the 1th radiating species.
The volumetric radiative heat transfer term, q^, in governing equation can be
expressed as
$ = -4K c(T4 - T 4) ca1
Hr p w ->
cm->-sec
The second term represents the radiation back from the environment having tempera-
ture Tw.
1-21
-------
s:
u
c
o
CL
i-
o
-------
8.0 NUMERICAL METHOD OF SOLUTION
The governing set of parabolic partial differential equations [Eqs. (4-2),
(<*-2), and C»-3)] are first rewritten in backward finite difference form which
gives implicit formulation. The implicit formulation is necessary due to stiff-
ness of the partial differential equations which describe the flow system involv-
ing fast finite-rate chemistry.
Numerical marching is in the x-direction. That is, all equations at all grid
points in the r-direction will be solved simultaneously at each x station. The
following basic differencing formulations, based on nonuniform numerical grids,
are used.
IE
3x
n+1 ,m
F ^ - F
n+1 ,m n,m
Ax
4y_
m
(8-1)
(2L) = _
\3y/n+l,m (Aym+] + Aym)
Fn+1 ,m+l " Fn+1 ,m
Ay,
m+l
+ Aym+1 Fn+1,m " Fn+1,m-1
(Ay J.1 + Ay )
m+l m' Ay
(8-2)
_2_ /(; IE)
3y I 3y /
n+1 ,m
1
0,5 (AVl + Ay-J
m'
Fn+1,m+1 ~ Fn+1 ,m
'n ,m+J5
Ay
m+l
^ Fn+1,m " Fn+1,m-l
n,m-b Aym
(8-3)
where iy ¦ y - y , and C „ - 0.5 (C , + C ). The subscript n
'm m Tm-1 n,m+1/2 n,m+1 n.m
denotes the nth grid point in the x-direction and m denotes the mth grid point in
the y-direction.
Applying Eqs. (8-1), (8-2), and (8-3) to Eq. and linearizing the
nonlinear terms in Q function via Taylor series expansion (truncated after the
first-order term) gives:
i^(Fn+l ,m " Fn,m)= |d (Fn+l,m+l " Fn+1 ,tn) + b (Fn+1 ,m Fn+1 ,m-l) |
c1(
F - F
n+1,m+l
i \
n+1 ,m/
+ d (Fn+1,m " Fi+1,m-1
+ rn.m + (Fn+l,m " Fn,in)
+R
n,m
1-23
(8-4)
-------
where
Aym
a1 - C1
'°'m (Aym+1 + Aym) Aym+1
k1 r1 Aym+1
~ = I
'°-m (4^l + 4ym)
c1' c1
i 1,m 2,m + h
' °'5 (^m+l + %)
rii . "Cl.m CZ.m-h
~ °'5 (A>Vl * 4ym)
In the above, the superscript '1" was introduced as an equation index. Species,
energy and density equations have to be solved In coupled fashion, while momentum
and continuity equations can be solved subsequently in uncoupled fashion. After
some algebra, Eq. (8-<0 can be rearranged ta form:
A1 F1 , , + Vb..FJ , + C1" F1 . , = (OR)1 (8-5)
n,m n+l,m-l ij n+1 ,m n,m n+1 ,m+l 'n,m
J
where:
a! - - (b1 + d1')
n ,m v '
B. . s 0. . if i ? j
ij ij J
= fBii - £ - a'-c^b'+d1) if 1 = j
(OR)1 = - (
-------
Incorporating proper boundary conditions into Eq. (8-5), the resulting
equation can be rearranged to obtain the following block matrix equation:
M F - Q
(8-6)
where T is a solution vector of elements F , and Q is a vector of elements
n+1, m
(OR)
n, m
M is a tridiagonal block matrix:
B1 C1
A2 B2
M =
^ 1 B ¦. 'C 1
m-l m-l m-l
A S"m
m m
The matrix B is a full matrix of elements B and the matrices A and B
i i
diagonal matrices of elements A and B
3 n, m n, m
(8-7)
are
Block matrix equation (Eq. 8-6) can be solved directly using the band-matrix
factorization method (Appendix A). The solution vector F contains species
distribution, temperature and density of all grid points at the fixed value of x.
Momentum finite difference equation (Eq. 8—^, F* » U) is solved individually via
tridiagonal matrix method for mainstream velocity U at all grid points.
Cross-stream velocity v is obtained via numerical integration of the continuity
equation point-by-point from the boundary where v is known.
The method is basically noniterative. Therefore, the truncation error due to
the linearization of nonlinear terms has to be carefully controlled. The error is
controlled by varying the marching step size Ax, since the finite difference
equations asymptotically approach the differential equations they represent when
A x -»¦ 0. The general strategy used is to limit the step size Ax so that no
component of solution can vary by more than the same small percentage of its value
at the last step; i.e.,
1-25
-------
max
Max
i ,m
F1 - F1
n-t-1,m n ,m
F1
n,m
< Ei
(8-8)
If ^mQX is greater than e ^ the numerical integration is repeated, cutting the
step size in half. The integration is repeated, continually halving the step size
until criterion (8-8) is satisfied. A lower limit (say, e_) for A is also
c max
set. If A is less than £-(< e,). the step size is doubled for next integra-
max 2 i
tion step. Based on past experience In numerical experiments, e = 2% and e =
5# are good error limits for general steady-state or time-dependent problems. For
the class of problems such as flat flame and strain flame in which only asymptotic
solutions are of interest, the error limits and e can be increased as long
as solutions remained stable. Since the calculations seek only asymptotic solu-
tions, the calculations will stop when the time-wise change of all variables
becomes very small, i.e.
Amax (n) ' £3 (8-9)
where e is a small number to be input to the program. The value of vories
with the class of the problems. Generally, e} ¦ 0.001 gives sufficiently
accurate solutions to problems such as flat flames or strain flames.
1-26
-------
9.0 APPLICATIONS
9.1 Free Jet Diffusion Flame
As an example of a typical flame calculation, a hydrogen-jet is considered
with maximum nozzle exit velocity of 1000 cm/sec. The jet issues into still air.
The nozzle diameter is 0.635 cm. Both the air and hydrogen temperatures are 300°K
at atmospheric pressure. This case Is the same as the sample case in Kee and
Miller's work (Ref. 13). The initial velocity distribution at nozzle exit plane
is assigned in the same way as with Kee and Miller (Ref. 13). The initial
ignition source is somewhat arbitrary, has no physical significance and represents
a very small percentage of total energy. Distributions of initial velocity and
temperature are plotted in Figure 2. The calculation was done in a physical
coordinate system with uniform good spacing in the r-direction. Radiative heat
transfer is not considered in this calculation. Initial number of grid points is
25, which includes two free-stream points. The flame shows very rapid expansion
in the radial direction near the nozzle exit plane. Therefore, the code allows
the addition of a grid point (at free-stream condition) whenever the
"next-to-last" radial grid point value of temperature or velocity or species
differs from the corresponding free-6tream value by more than a specified percent-
age of free-stream value (one percent has been selected in this calculation). The
number of grid points is not allowed to grow unbounded, with the number of points
being halved (eliminating every other point) when the maximum number of points
(50) is reached.
The reaction mechanism used is the H^-subset of the basic EER kinetic set
(Appendix B). The H^-subset involves 9 species (H^, O2. 0, H, OH, and
HO^) and 18 reactions.
The numerical results of this calculation are shown in Figures 3-6. Figure
3 shows the isothermal contours of the flame. The flame length, which Is defined
as the peak temperature at the center line, is 5.8 cm. Figure 1* shows the distri-
butions of downstream velocity and cross-stream velocity at a distance of 5 nozzle
radii from exit plane. The distributions of major stable species and temperature
(at X/R-5.0) are plotted in Figure 5. It is noted that the peak temperature is on
the rich side of the flame, which Is consistent with the result of the equilibrium
1-27
-------
(u0)^= 1000 cm/sec
DQ 3 0.635 cm
Tai r 3 300°K
Tu = 300 K
m2
P = 1 atm
2500
2000
1000
1500 _
T (r)
500-
1000
U(r)
500
0.4
0.2
0.3
0.1
r (cm)
Figure 2. Initial Velocity and Initial Temperature Distribution.
1-28
-------
20.0
in I
15.0
a£
x
10.0
5.0
(Uo)^ = 1000 cm/sec
Dq =» 0.535 cm
air
fH2
= 300 K
= 300°K
= 1 atm
4.0
6.0
8.0 10.0
r/R
Figure 3. Isothermal Contours for Laminar h^/Air Diffusion Flame.
1-29
-------
300
300
200
200
U(r)
100
i- 100
V (r)
-40
-40
r (an)
Figure 4. Velocity Components U(r) and V(r) at X/R =
1-30
5.0.
-------
« 0.4
a 0.3
r (an)
Figure 5. Distribution of Stable Species and Temperature at
X/R = 5.0.
1-31
-------
1Q
-1
c
o
-------
calculation.
temperature.
Figure 6.
The peak temperature is very close to the actual adiabatic flame
The distributions of radical species (at X/R ¦ 5.0) are plotted in
The entire calculation for the process, which occurs 13 cm from exit plane,
takes 930 noniterative integration steps in the x-direction. The computer time is
about 200 sec CPU time (in a CDC 7600 computer) which is equivalent to about one
twentieth of the computational time needed by Kee and Miller (RBf. 13) far the
same calculation (using a split-operator method).
The numerical results of the present calculations are not expected to be
identical with those of Kee and Miller, since the formulation of diffusion coeffi-
cients and the reaction mechanism used in present calculations are different from
those used by Kee and Miller.
9.2 Confined. Co-flowing LBQ Diffusion Flame and Fuel Nitrogen Conversion
In order to study the effect of methane level in low 8tu gas diffusion flames
on fuel nitrogen conversion, a series of laboratory scale laminar diffusion flame
experiments has been conducted in EER's laboratory (Ref. 1
-------
Fuel Flow
Figure 7. EER Diffusion Flame Burner.
1-34
-------
Case No.
1
2
3
Percent CH, in Fuel
0
5
10
Fuel Mean Velocity (cm/sec)
685
685
685
Air Mean Velocity (cm/sec)
13.15
17.53
21 . 92
Air Temperature (°K)
367
367
367
, o .
Fuel Temperature ( K)
367
367
367
Visual Flame Length (cm)
35.5
^0.6
50.8
This gives ~ constant overall stoichiometry equal to 150 percent T.a. when
the effects of dopants on stoichiometry are neglected.
In comparison with EER's experimental data, the predicted exhaust NO levels
are shown in Figures 8 and 9 for NH^ doping and NO doping, respectively. EER's
latest kinetic set was used in these calculations, quantitatively predicting the
exhaust NO level for the flame with various types of LBG fuel and dopant.
The exhaust NO levels without dopant, i.e. thermal NO, are also shown in
Figures 8 and 9 in comparison with experimental data. The discrepancy is due to
the sensitivity of thermal NO level to the flame sheet temperature. The predicted
flame temperature for the flame sheet of the flame with 10 percent CH^ (in fuel)
o o
ranges from 1960 K to 1720 K and is, of course, sensitive to the accuracy of the
heat transfer model.
Calculations were made for the NO doping case with 10 percent methane to
establish the effect of flame temperature on fuel nitrogen conversion. Figure 10
shows the peak flame sheet temperature, total NO, and NO with a thermal level
subtracted out, all as a function of increasing heat transfer; i.e., decreasing
flame sheet temperature. Heat transfer increases linearly with the heat transfer
index shown in Figure 10. The index has a value of zero for the adiabatlc case
and a value of one for the heat transfer associated with a simple radlotion model.
Although the peak flame sheet temperature varies by over 200°K, the total NO level
remains approximately constant. However, due to the change in production of
thermal NO, the uncoupled "fuel nitrogen conversion" Is predicted to increase
substantially as flame temperature is depressed.
1-35
-------
Dopant in Fuel
Exp. O NK^ 3059. ppm
Exp. A NH^ 0 ppm
Calc. # NH^ 3050 ppm
LBG
N2 55%
H2/C0 1.0
CH^ Varying
— Calc.
NH3 Q ppm
Total NO
This calculated point has a high degree of
uncertainty due to a heat-loss-1nduced wall
recirculation bubble. Refined calculations
are currently under consideration.
Thermal NO
Percent of CH^ on Fuel
Figure 8. Influence of CH4 Content on NH3 Conversion
in LBG Diffusion Flames.
1-36
-------
1500
Dopant in Fuel
Exp. O NO 3059 ppm
Exp. A NO 0 ppm
Calc.# NO 3059 ppm
Calc.A NO 0 ppm
LBG
N2 55%
h2/co 1.0
CH4 Varying
»>
S-
-a
Q.
Q.
1000 -
CO
3
rO
This calculated point has a high degree of
uncertainty due to a heat-loss-induced «_
recirculation bubble. Refined calculations
are currently under consideration.
otal NO
Thermal NO
Percent of CH^ in Fuel
Figure 9. Influence of CH^ Content on NO Conversion
in LBG Diffusion Flames.
1-37
-------
NO Conversion
(Thermal NO Included)
— 2500
NO Conversion
(Without Thermal NO)
2000
20
EER - LBG Flame
(C0/H2/CH4/N0/N2)/Air
= (17.5/17.5/10/0.3059/55)/A1r
— 1500
10
1000
Radiation Index
Figure 10. Influence of Temperature on NO Conversion
in Diffusion Flames—10 Percent CH^.
1-38
-------
9.3 Flat Flame Simulation
Flat flame calculation is a typical time-dependent problems with an
asymptotic solution. Here the methane flame of Biordi et al. (Ref. 15) has been
chosen as an example for comparison to GFAP's prediction. Biordi's study used a
low-pressure (32 torr) premixed methane flame (10.3 percent CH^, 21.6 percent 0,
68.1 percent AR). Flow velocity at 300°K was 80 cm/sec from a flat-flame burner.
Molecular beam mass spectroscopy was used to determine stable and radical species
concentrations. The flame was simulated as a free (unattached to the flat flame
burner surface) flat flame. The free flame velocity calculated is 76.2 cm/sec,
which agrees very well with the experimental value. The stable species profiles
calculated are plotted in Figure 11 in comparison with Biordi et al.'s measure-
ments. The major deviation appears to be underprediction of CO oxidation rate by
the current EER kinetic set. Figure 12 illustrates the comparison between experi-
mental data and GFAP predictions for a number of radical species. In general, the
agreements are less than satisfactory. The predicted OH shows fair agreement with
data, while both 0 and H radicals are about 50 percent below the experimental
values. The predicted temperature profile is higher than the measured value,
which is probably partially due to gas radiation which was not taken into consid-
eration. The only intention here is to illustrate the application of GFAP;
further discussion of kinetic aspect of this calculation is beyond the scope of
this report.
9•^ Constant Strain-Rate Laminar Diffusion Flame
The constant strain-rate laminar diffusion flame is a planar nonpremixed
diffusion flame with external fuel and oxidizer flow undergoing uniform accelera-
tions (e.g. u-ax) such that the strain rate, du/dx, is independent of position
along the flame. This situation exists, for example, near the stagnation point in
opposed-jet combustion. More importantly, a constant strain-rate laminar diffu-
sion flame characterizes the flamelet behavior of the fuel/oxidizer "eddy" inter-
face in nonpremixed turbulent flame. In these regions, the strain rate is slowly
varying.
It can be shown that the solution to the Navier-Stokes equation describing
the behavior of constant strain-rate laminar diffusion flames is self-similar in
1-39
-------
0.15
0.10 -
(O
ai
"o
0.05 -
Prediction
0.4 0.5 0.6
Distance (cm)
Figure 11. Flame Profiles of Biordi et al. (1975). Solid Lines are Data ;
Dashed Are Model.
-------
Prediction
¦y.
Distance (cm)
Figure 12. Profiles Through the Flame of Biordi et al. (1975).
-------
the streamwise (x) direction and can be represented by the boundary layer equation
under uniformly accelerating flows.
The GFAP model properly analyzes a constant strain-rate laminar diffusion
flame. The computation can start with an arbitrary distribution and march forward
in the streamwise direction until a steady solution is achieved. The pressure
gradient in the momentum equation is obtained from Bernoulli's equation in differ-
ential form:
du dp
ou dx ~ ~ dx
Noting that the fuel side acceleration u equals ax, it can be shown that
dp 2
d£* -pa x •
The sample calculation here is for a CH^/air opposed jet diffusion flame (see
Figure 13). The fuel side strain rate is set to 250 (1/sec). Both fuel and air
are at 300°K and one otm. The velocity (v) profile for the steady solution Is
platted in Figure 14. The displacement 6 Is the distance between imaginary
(nonreacting and inviscid solution) stagnation point of the fuel jet and that of
the oir jet. For this case 5 is 0.15 cm. The CH^ and 0^ profiles are shown in
Figure 15. The temperature profile is also shown in Figure 15. The peak flame
o
temperature is about 1950 K. Peak temperature decreases with increasing strain
rate since increase in strain rate means a thinner flame zone and shorter effec-
tive residence time of reactants In reaction zone. Figure 16 shows the flame
temperature as a function of strain rate. The critical value of strain rote is
shown to be 650 (l/sec), beyond which a flame cannot exist (i.e. the extinction
limit). The critical value is a function of reaction chemistry alone and is
independent of molecular transport coefficients (i.e. independent of Reynolds
number). Thla is due to the fact that the molecular transport coefficients (e.g.
viscosity) influence the solution only through an offine transformation of the
r-coordlnate; I.e.
= r/f
The self-similar solution of temperature and species concentration is found as a
1-42
-------
r FUEL SIDE
Ji
Imaginary Fuel Jet
/ Impingement Plane
» X
FLAME
ZONE
Imaginary Oxidizer Jet
Impingement Plane
OXIDIZER SIDE
Figure 13. Schematic of Opposed-Jet Diffusion Flame.
1-43
-------
12C
100
FUEL
r (cm)
-20
AIR
-40
-60
Figure 14. Velocity Profile of CH./Air Opposed-Jet
Diffusion Flame (a*250 1/sec).
1-44
-------
30
OJ
©
0)
u
i.
dJ
o.
un
20
(Percent CH4)/5
X
o
a;
u
u
QJ
G_
10
0.5 0.4 0.3 0.2
Percent
2000
1500
1000
Figure 15. CH4 and O2 Profiles of Cfy/Air Opposed-Jet
Diffusion Flame (a=250 1/sec).
1-45
-------
CH./Air Stretch Flame—Full Kinetic
2000
1500
1000
700 800
600
400
a (1/sec)
500
300
100
200
Figure 16. Flame Temperature vs. Strain Rate for
CH^/Air Opposed Jet Oiffusion Flame.
1-46
-------
function of r only. Viscosity has on influence on flame thickness and therefore
on total burning rate.
Figure 17 shows the flame thickness, which is defined as the distance between
the intersections of maximum slope of temperature curve and the T = 300°K coordi-
nate line, as a function of strain rate. Figure 18 shows the stagnation point
displacement as a function of strain rate. The computations are carried out using
the basic EER kinetic subset including 71 reactions and 28 species. The computa-
tional time per case is approximately k minutes CPU time in CDC 7600.
1-47
-------
CH./AIR
300 400 500
a (1/sec)
800
Figure 17. Flame Thickness vs. Strain Rate for
CH^/Air Opposed-Jet Diffusion Flame.
1-48
-------
0.4
C
O)
£
O)
o
05
0.3
•5 0-2
CL
T3
C
CD
03
0.1
100 200 300 400 500 600 700 800
Strain Rate a (1/sec)
Figure 18. Stagnation Point Displacement vs. Strain Rate for
CH^/Air Opposed-Jet Diffusion Flame.
1-49
-------
10.0 REFERENCES
1. Reid, R. C., Prensnitz, J. M., and Sherwood, T. K., The Properties of
Gases arid Liquids. 3rd edition ( 1977).
2. Monchick, L. and Mason, E. A., "Transport Properties of Polar Gases," J.
of Chemical Physics, Vol. 35, No. 5, p. 1676 (1961).
3. Neufeld, P. D., Janzen, A. R., and Aziz, R. A., "Empirical Equations to
Calculate 16 of the Transport Collision Integrals (l.S) for
Lennard-Jones (12-6) Potential," J. of Chemical Physics, Vol, 57, No. 3,
p. 1100 (1972).
t*. Bartlett, E. P., Kendall, R. M. , and Rindal, R. A., "A unified
Approximation for Mixture Transport Properties for Multicomponent
Boundary Layer Applications," NASA CR-1063 (1967).
5. Ferri, A. J., "Review of Problems in Application of Supersonic
Combustion," J. Roy. Aeron, Soc. 68 p. 575 (1964).
6. Schetz, J. A., "Turbulent Mixing of a Jet in Co-flowing Stream," AIAA,
Vol. 6, No. 10 (1968).
7. Sparrow, E. M. and Cess, R. D., Radiation Heat Transfer, Brook/Cole
Publishing Company, Belmont, California (1970).
8. Kelly, J. T. and Kendall, R. M. , "Further Development of the Premixed
One-dimensional Flame (PROF) Code," EPA-600/7-78-172a (NTIS PB286243)
(1978).
9. Edwards, D. K. and Balakrishnan, A., "Thermal Radiation by Combustion
Gases," Int. J. Heat and Mass Transfer, Vol. 16, p. 221 (1973).
10. Stull, D. R. and Prophet, H., "JANAF Thermochemical Tables—Second
Edition," Dow Chemical Company, (1970).
11. Peaceman, D. W., "Fundamentals of Numerical Reservoir Simulation,"
Elsevier Scientific Publishing Co., New York (1977).
12. Kau, C. J., Pergament, H. S., and Mikatarian, R. R., "A New Computational
Technique for Chemically Reacting Free-Shear Layer Flows," Symposium on
Application of Computers to Fluid Dynamic Analysis and Design,
Polytechnic Institute of Brooklyn, N.Y., Paper
-------
10.0 REFERENCES (Continued)
Ife. Folsom, 8. A., Courtney, C. W., Corley, T. L.. and Clark, w. D.,
"Advanced Combustion Concepts for Low Btu Gas Combustion," Proceedings
of the Third Stationary Source Combustion Symposium; Volume II.
EPA-600/7-79-090b (NTIS PB292540) (1979). p. 163.
IS. Blordl. J. C.. Lazzara, C. P.. and Papp, J. F., "Flame Structure Studies
of CF Br-Inhlblted Methane Flames, II, Kinetics and Mechanisms,"
Proceedings of the Fifteenth Symposium (International) on Combustion. The
Combustion Institute (1975), p. 917.
1-51
-------
11.0 LIST OF NOMENCLATURE
defined as
,C„,C,
tit
defined in Table 2
defined in Table 1
specific heat of species i
diffusion coefficients (cm^/sec)
third body efficiency defined in Eq. (6-5)
defined as
diffusion factor of species i
Gibbs free energy
enthalpy of species i (cal/mole)
equilibrium constant
forward reaction constant
gas omittance
characteristic dimension of flame
Lewis number
molecular weight of species i
pressure
partial pressure of species i
Prandtl number
defined in Table 1
radiative heat transfer rate (cal/cm3-sec)
defined in Table 1, or gas constant
cross-stream coordinate (cm)
1-52
-------
11.0 LIST OF NOMENCLATURE (Continued)
T Temperature (°K)
u downstream velocity (cm/sec)
v cross-stream velocity (cm/sec)
chemical production rote (mole/cm^-sec)
x downstream coordinate (cm)
species mole fraction
y : cross-stream coordinate
- r for physical coordinate
= for stream function coordinate
Superscripts
1 equation index
Subscripts
£ center line
e edge of mixing layer
i species index
m grid index in r or '¦), direction
n grid Index in x-dlrection
R,r reference species
1-53
-------
11.0 LIST OF NOMENCLATURE (Continued)
Greek Symbols
P
X
Q (1 -1 >
o<2-2>
:1 ' e2 ' e3
geometric index
- D for 2-D cartesian coordinate system
= 1 for polar cylindrical coordinate system
¦ 2 for spherical coordinates
ith species concentration (mole/gm)
density (gm/cm^)
heat conductivity (col/cm-sec-°K)
stream function coordinate
molecular/turbulent viscosity
defined in Eq. (5-3)
collision integrals
minimum energy of attraction between molecules
numerical error controls
1-54
-------
APPENDIX A. Bloch Tridiaqonol Matrix Equation Solver
To solve a matrix, the equation used is
M F - a
(A-1 )
where Q and F are column reactors. Q is given, F is the solution matrix, and M is
o block tridiagonal matrix of the form:
M -
B1C1
A2 B2 ^2
A3 B3 C3
\ \ *
\ *
\ N
A i B i C ¦.
m-1 m-1 m-1
A Bm
m m
(A-2)
here are full matrices, A± and are diagonal matrices,
This matrix equation con be solved by band matrix factorization method (Ref.
11). The method will be described briefly in the following section.
The block tridiagonal matrix M can be factorized into lower and upper band
matrices:
M » L U
(A-3)
where L is the lower band matrix and U is the upper band matrix:
A2 L2
A3 L3,
Am-1 '"m-1
\,L
m
1-55
-------
ond
T U.
I U2
I "3
I Vl
I is an identity matrix ond A^ remains the same as in matrix M. l and u are
full matrices and are defined by:
L1 - B1
-1
and
Li " Bi " Ai L i-1 Ci-1 for 1 " 2'3'
1C^ for i - 1,2, .... m-1
(A-5)
where L. is the inverse of matrix L,.
1 i
Now the matrix equation (A—1) becomes:
L U F - Q
(A-6)
The equation can be solved by the following matrix algorithm:
Let U F =» G, equation (A-6) becomes
L G - Q
or
A2 L2
A3 L3
V' L
1 m-1
A L
m m
in
m
1-56
-------
Now we con solve for by following relations (forward pass):
- l^-1 ((^ - AlGi_1) for i - 2,3 m
After G is obtained, we can solve matrix equation
U F - G
for solution vector F by following relations (backward pass):
F =¦ G
m m
and
F - G - U. F. , for i - m-1, m-2 1
i i i i+1
1-57
-------
Appendix. Basic EER CH^/Air Reaction Mechanism and Rate Constants
HtACTiUiMS btluG CUNSlOtHtO
'
-
KF«AM
£
fa
N)*E*P(
-b/HI)
a
'
1 1 H2
~ 1
u
~
0
3
1
HU ~
1 H
~ 0
i.aoot
11
10
-1 .00,
D
8900.0
)
£ J
t 1
HQ
t -
1
_1 .H
t.. Q -
l.JOQt
g9
-1.30
'i J I H2
t 1
04
~
0
s
1
HO t
1 HU
~ 0
2.51Ot
12
0.
39000.0
¦ '
4 I h
~ 1
H
~
1 M
B
1
H2 ~
1 M
t 0
9.000E
1 /
1 .00
0.
->
3KD bgp *
EFF .
... £•
000/H20
fc 1 H
V 1
HO
t
1 M
s
1
H20 t
1 M
~ 0
2¦200t
22
2.00
0.
3rtD BODK
tFF.
5.
000/H2O
(-1 L » _
~ I
02
~
I M
*
I
H02 ~
1 M
~ 0
1 .040t
15
0.
-1000,0
j 3«i> buor
IFF.
2.
500/H2
i
.800/CH4
15.000
/H20
3.750/CO2
:.! 7 | U
~ 1
u
~
1 M
s
1
02 ~
1 M
~ 0
l.ooot
1H
1 .00
0.
6 1 h
~ 1
02
~
0
s
1
HO ~
1 u
~ 0
4.S00E
10
-1.00
14UU5.0
| V 1 HO
~ i
H(J
~
o
•
1
H20 ~
\ u
~ 0
2.140b
09
-1.10
0.
,=
10 1 H
~ 1
H02
~
0
¦
1
HO t
1 HO
t o
2.Soot
14
0.
1900.0
U l hO
~ i
h02
~
0
«
1
H2U t
1 02
~ 0
5.000E
13
u.
1000.0
/#
12 I 0
~ i
H02
~
0
*
1
HO ~
1 02
~ 0
5.000b
13
u.
1000.0
U 1 H
t i
H02
~
0
a
1
02 ~
1 H2
~ 0
2.500E
13
0.
700.0
>„
14 1 H202
~ i
M
~
0
a
1
HO t
1 HU
~ I M
1 .200t
17
0.
45500.0
15 1 H202
~ i
h
t
0
a
1
H20 t
1 'HO
~ 0
3.000E
14
0.
9000.0
1ft 1 H202
~ i
HU
t
o
a
1
H2U t
1 H02
~ .0
1 .OOOE
13
0.
1000.0
! 17 1 HQ2
~ i
H2
~
¦o
¦
1
H2U2 ~
1 H
~ 0
3.OOOt
11
0 .
1B60O.O
<4
It 1 HU2
~ i
H02
~
0
a
1
H2U2 t
1 02
t 0
2.OOOE
12
0.
0.
<> 19 i co
~ i
0
~
I «
X
1
CU2 ~
1 *
~ 0
3.8UUE
24
3.0U
6170.0
>' 20 1 CO
~ i
02
~
0
a
1
C02 ~
I u
t 0
6.91Ot
07
-1.00
34B10.0
i'i J CO
~ i
HO
~
0
a
1
CO 2 ~
1 H
t 0
1.510b
07
-1 .30
-7*5.0
• a i cu
t i
H02
~
0
a
1
C02 ~
1 HO
~ 0
1.ooot
11
0.
10000.0
il 1 CH4
~ l
M
~
0
a
1
CH3 ~
1 M
~ I M
1f 2b0t
17
0.
88400.0
¦i 2« 1 Ch4
t i
G
~
0
a
1
CH5 t
1 HO
~ 0
i.wot
07
-2.10
7bl0.0
25 1 CH4
t i
h
t
o
a
1
CH3 t
1 H2
~ 0
1,410t
07
-2.00
68440.0
I CHM
~ i
HU
~
0
s
1
CH3 t
1 H2U
t 0
1,550t
06
-2.10
2452.0
' 27 1 CH4
~ i
02
~
0
a
1
CH3 t
1 H02
~ 0
e.5iot
06
-2.00
51970.0
iB 1 CH4
t i
HU2
~
o
a
1
CH3 t
1 H2U2
~ 0
2.ooot
13
0.
18000.0
• 29 1 CH4
~ i
CH
o
a
CH3 ~
1 .CH^
t 0
2!SOOt
1 1
- . 70
6000,0
•! so l chm
~ i
Cn2
~
0
=
1
CH3 ~
1 CH3
t 0
1,260t
1U
-JO
20000.0
11 1 Lh3
~ i
0
t
0
3
1
Crt2(J ~
1 H
t 0
9.OoOt
13
0.
2000.0
12 i CH}
~ i
HU
~
0
a
I
Cirf ~
I nm
t o
6f 30Ot
10
-.70
2000.0
33 1 CH)
t i
U2
t
0
a
1
CMJU ~
1 Q
~ 0
».uoot
09
"1 .00
28000.0
Reproduced from
tx>st available copy.
-------
14
1 CHiU ~
M
~ 0
•
1 CHiU
Ifi
I tMiU ~
U
~ 0
¦
1 Cl
36
1 CHiU ~
H
t 0
X
1 CH20
37
1 CHiU ~
HU
t 0
-
1 Cm
-------
/i
1
1 HL)^ ~
0
¦
NH
74*
1 NH2
1 mO
1 Nh2 ~
0
a
N2U
102
1 NU
1 nh2 t
0
z
N2
-101
1 NU
1 nh2 ~
0
s
n2
104
1 NU
I NH t
0
X
N2U
lot
1 N2U
1 M ~
0
X
N2
10*
1 N2U
1 o t
0
m
N2
107
1 N2U
1 u ~
0
•
NU
t.|Ot
1 N2U
1 H ~
0
s
M2
|09
1 N2U
1 NU t
o
X
N2
M®
1 N£U
1 NH t
0
s
N£
ill
1 N02
1 U ~
0
a
NO
112
1 NU2
1 H t
0
*
no
111
1 N2
1 M ~
0
X
N
H202 ~
0
i .oootfiu
0.
1200.0
U2 ~
0
2.000t*10
-.50
2590.0
NO ~
0
S.OOOE+11
-.50
2740.0
NH ~
0
1 .2U0t»11
-.50
5640.0
H2 ~
1
i,9801~13
0.
12000.0
H ~
0
6.31Otfl1
-.50
0.
HU ~
o
6.31Uttil
-.50
aooo.o
M t
0
2.00Uttlo
.50
0.
H2 ~
0
0•310£f 11
-.50
8000.0
H ~
0
5¦0O0E~ 11
-.50
5620.0
H2U t
0
S.OOUEtl 1
-.50
4450.0
U ~
0
b.UUOEf 11
-.50
43/0.0
hO t
u
6.300fcfl1
-.50
0.
H2 t
o
3.(>0ULtl 1
-.55
1900.0
NU ~
0
5.0006+11
-.50
£740.0
NO ~
0
1.UOOttl 3
0.
2500.0
NU ~
0
1 ,OOOEtl2
-.50
2740.0
M ~
0
3. 000£t16
.50
0.
J t
0
6.400Et09
-1.00
e>260.0
h t
0
6. 3l Qt* 1 1
-.50
0.
NM2 t
0
l.OOOEtl1
-.50
31640.0
NH t
0
2.000L+l1
-.50
22160.0
H ~
0
6.3lOEt11
-.50
0.
U ~
0
3.1OQEt13
0.
334.0
0 t
0
3.OOOEt11
-.50
0.
H ~
0
1,OOOEt21
1 .86
0.
M t
0
2.200Etl7
.50
0.
HU t
0
3.000Etl1
-.50
0.
H2
0
S.OOOEtOS
-.50
o,
H ~
1
HO
1•3B0E~ 1 7
1.85
0.
H2U t
0
2.3/UEtl7
1 .85
0.
H ~
0
b.aOOEt 16
1.40
0.
U t
1
M
1 .<420EtM
0.
51200.0
02 ~
0
6.230E+I3
0.
24540.0
NU ~
0
3,ObOfc*13
0.
21800,0
HO ~
0
7,60 OE~ 1 3
0.
15100.0
NU2 ~
0
2.OooEt m
0.
50000.0
HnO ~
0
2.000ttl0
-1.00
12000,0
U2 ~
0
I,0 OOE~ 1 3
0.
600.0
NU ~
0
2.900E*J«
0.
800.0
N ~
I
M
3.700E*^1
1 .60
225000,0
-------
114
1 l2*c> ~
1 M
0
-
I Cihb
115
1 C2no ~
\ HU
0
a
\ c2H5
116
1 L^ho t
1 Lriit
u
=
1 CH3
117
1 L 2Mb t
1 Chi
0
z
1 C2H5
11«
1 Lh2 t
I Chi
0
3
t C2H2
ll«
1 CU t
1 CM4
0
5
1 C2H<4
120
1 C hi t
i Chi
0
5
1 C2"4
121
1 Ch3 ~
1 CH4
0
s
1 C2M5
\a
I Ch3 t
JLJCitS
1 1*1
-
1 C2Mp
123
1 Ch4 t
1 CH3
u
s
1 C2rtb
124
1 CHN t
1 H
0
3
l Cn
L4S
... 1 ~
1 0 ..
3
i n .
126
1 CHN t
1 HO
0
£
1 CM
127
1 CHN t
1 HU
0
-
1 nnCU
H*
4 tN ~
1 0
0
1 CU
ta«
1 CN t
1 HO
0
x
1 NCO
130
1 CN t
1 HU
0
s
1 CHN
I>1
I CN t
1 02
0
x
1 NCO
132
1 CN ~
1 H2
0
z
1 CHN
133
1 CN t
1 CU2
0
s
1 NCO
u«
1 y ~
I NH2 _
0
5
1 CHN
us
1 CN t
1 nh3
0
a
1 CHN
136
1 NCO ~
1 0
0
B
1 NO
111
I NCO ~
1 H
0
8
1 NH
13*
1 NCO t
1 H2
0
Z
1 hnCo
119
1 NCO t
1 NH2
• 0
X
1 Nrl
140
1 HNCU t
1 HO
0
•
1 NCU
141
I HNCU ~
1 H
0
B
1 nh2
Ml
I CHO t
1 N
0
S
1 CHN
1 CHN t
1 N
0
S
1 CH
140
1 CHN t
1 0
0
X
1 CH
US
1 CH4 t
1 CN
0
X
1 Chi
146
1 CH3 ~
1 CN
0
3
1 Ch2
14?
I CH2 t
1 NO
0
9
1 CH20
146
1 CH2 t
1 N2
0
S
1 CHN
}49
1 CH ~
1 N
0
3
1 CN
ISO
1 CHO t
1 N
0
m
1 CH
1S1
1 C2H t
1 0
0
3
1 CH
1S2
1 C2H t
1 HU
0
Z
1 Cri2
1S3
i can t
1 0 2
0 "
¦
1 CHO
1S4
1 C2«2 ~
1 M
0
•
1 C2H
|S5
1 ci"j
I Q
0
¦
J t
Hi ~ 0 5.370ttu2 -1.50 4200.0
N2U ~ 0 b.3l0fctll 0. 2450.0
C2H5 t 0 1.22i)Etll 0. 15700.0
CH4 ~ 0 5.500E-01 -4.UQ bittO.0
H2 ~ 0 S. 2 0 01 ~ 1i 0. 0.
H ~ 0 1,50(Jt ~ I 0. 0,
H ~ 0 6. OOOEt13 0. 0.
Mi t 0 1 .000E+13 0. £4000.0
* t 0 2,400ttl3 Of o,
H tO B.OOOttli 0. 40000.0
H t 1 M 5 . 700Et1b 0. 1 [7000.0
NiiJJ. . f... 0 l.JOOtfH -.50 6510.0
"20 t 0 2.000Etll -.50 5000.0
H to 2.000ttl0 -.50 5620.0
'« to b.ilOEtJl -,50 0,
H t 0 b.20l)E*13 0. 0.
0 t 0 3) 1bOEt12 0. 3000.0
_0 ~ 0 _ 3,200E*13 0, 1000,0
H tO 6.200Et1 2 0. 5300.0
CU t 0 4.000Etll -.50 7000.0
NH t0 5,QOOEt10 -.70 2000,0
Nh2 t 0 r.OOOE+iO -.70 2000.0
CO t 0 5. OOOEt1 1 -.50 6875.0
CO ~ _0 S.OOOttll -.50 6870,0
H * 0 3.000EM0 -.50 5700.0
HoitO t 0 l.OoOEtll -.50 5000.0
H2J t 0 1 .000EM1 -.50 6290.0
CO ~ 0S.OOOttil -.50 4610.0
0 tO 1.000EM4 0. 0.
N2 t 0 2.SOOEt I 1 0. 16000.0
NO ~ 0 ~ l.OOOttM 0. 72000,0
CHN t 0 3.lbOEtI 1 -.70 5000.0
CHN t 0 l.OOULtl) -.70 3000.0
N ~ ~ 0 1•bOOE 112 0. 7000.0
NH ~ 0 1.000E+14 0. 60000.0
H t 0 fa.loOEt11 -.50 0.
NO t' 0 ~ l.OoOEtia 0. 48600.0
CU t U 5»OOOEt15 0. 0.
Cut 0 fa.OOOEtl2 0. 0.
CUto 1.OOuEtI3 0. 0.
If ~ 1 * 4. 170E116 0. 107000.0
tft. ~ 0 to,7lOEt 11 0, 4000,0
-------
156
i C2M*
1 M
0
8
1 C2M
157
1 L2M2
1 u
0
a
1 LJH
1 SB
i C2"2
1 no
U
¦
1 C2H
159
i 12M2
1 MO
0
a
1 CHi
160
l 12M2
1 02
0
s
1 CnO
,1*1
1
1 CM
g
3
I CH«>
162
1 C2H2
1 Chi
o
s
1 CM3
'{163
1 12M3
1 M
0
3
1 C2H2
- .
1 C«»M3
1 H
0
=
1 C£*d
[l6S
I C2H«
1 M
0
3
1 C2M2
1*6
1 C2H4
1 M
0
a
1 C2M3
1*7
1 C2H«
1 u
V!
s
1 CM2
US
1 C2M«
1 u
0
=
I Crii
119
1 C2h
3
l cdim
iau
1 C 2M5
1 C2rlb
0
a
1 L2M©
i»i
1 12M6
1 u
0
a
1 C2HS
"2
0
7 . 7t>0t*00
-3.20
boo .0
liU
0
3.200tt
b
.60
1 7000.0
Hi (J
0
1 ,«v)0£t
3
0.
7000,0
CO
0
b.bOOEt
3
0.
13700.0
CMU
0
4.000t+
d
0.
20000.0
C^M
0
/. JOOtt
3
0.
1 1920,0
C2M
0
1 . OUOEt
d
0.
0.
ri
1
M
7.9<<0t +
a
0.
31500.0
m
0
1 r000t +
3
0.
0.
M2
1
M
2.600E*
7
0.
79290.0
H
1
M
3.80o£t
7
0.
96170.0
CH2U
0
af060tt0b
-2, 00
307 ,0
CMU
0
5.930E+
0
-.7H
560.0
rid
0
1 . 1 OOEf
<4
0.
8500.U
M
0
1 ,oouE*
3
o.
IbOO.O
CH20
0
4.500E*
2
0.
220.0
M20
0
1 .OOOE*
U
0.
3SOO.O
C2M3
0
l,090fc*
3
0.
12400,0
C2M3
0
b.OOOt*
2
0.
13000.0
CMi
0
S.ldQt*
3
0.
0.
M d
U
1,0bOt~
d
0 .
u.
MO 2
0
1.OOOE*
d
o.
5000.0
C2M
-------
PART 2
MATHEMATICAL MODELING
or
MICROSCALE COMBUSTION OF A COAL PARTICLE
Prepared by:
C. J. Kau
and
T. J. Tyson
2-1
-------
TABLE OF CONTENTS
Section Page
1.0 DEVOLATILIZATION MODELING 2-1
2.0 CHAR COMBUSTION MODEL 2-6
3.0 MICROSCALE COMBUSTION OF A SINGLE COAL PARTICLE 2-8
4.0 SUDDEN IMMERSION OF A COAL PARTICLE INTO HOT
COMBUSTION PRODUCTS 2-10
5.0 MACROSCALE VS. MICROSCALE 2-21
6.0 CHARACTERISTIC TIME SCALES ON COAL PARTICLE COMBUSTIONS . . 2-24
7.0 REFERENCES 2-27
LIST OF FIGURES
F iqure Page
1 Schematic of Microscale Combustion of a Coal
Particle in ~ Confined Space 2-9
2 Mass Flux Rate of 50u and 100- Bituminous Particles
in Higher Heat Loss Environment (T = 1500 K) 2-11
3 wall
3 Particle and Gas Phase Temperature Histories for T
- 1500°K WQl1 .... 2-12
4 Mass Flux Rate of 50 - and 100g Bituminous Particles
in Lower Heat Loss Environment (T - 1800 K) 2-13
wall
5 Particle and Gas Phase Temperature Histories for
T , - 1800°K 2-14
wall
6 Temperature Distribution Around a Coal Particle 2-16
7 Oxygen Distribution Around a Coal Particle 2-17
8 Distribution of NO Production Rate Around a Coal
Particle 2-18
9 Distribution of N Production Rate Around a Coal
Particle 2-20
2-11
-------
LIST OF TABLES
T able Page
1 Rate Coefficients 2-3
2 Asymptotic Yields 2-4
3 PPM Levels of NO and Percentage Conversion of Fuel-N
Calculated by Microscole and Macroscale Models 2-22
i* Physical and Chemical Time Scale 2-25
2-111
-------
-------
1.0 DEVOLATILIZATION MODELING
For the purpose of mothematicol modeling of cool decomposition, the
constituents of raw coal are catalogued into five categories. They are:
1. Light volatiles: including H^, CO, CO^, H^O, moisture and aliphatics
2. Volatile nitrogens
3. Tar: having chemical structure similar to parent coal
Char: basically carbon plus some residual nitrogen
5. Ash: including sulfur, sulfur compounds, and all other mineral
matters which are not treated chemically in this model.
Since the devolatilization (thermal pyrolysis) of coal is a purely chemical
decomposition process, the evolution of light volatiles, volatile nitrogen and tar
can be viewed as a first-order chemical reaction, having rate equation
dY.
1 _ 1, V
where Y is the mass fraction of ith solid species (or functional group, in
Solomon's terminology). is a rate constant of Arrhenius form:
Ki = Ai exp^ Ei/RT) (1/sec)
Currently, the devolatilization model in EER computer codes (OFAP and GFAP)
can handle two types af coals: Montana lignite and Pittsburgh bituminous. The
rate constants, K , were derived from the following experimental works:
1. Suuberg's (1979) measurements: providing the rate data and
information on functional group partition of light volatiles.
2. Kobayashl (1977) measurements: providing data on total weight loss
and rate constant, which can be used to estimate tar yield and
evolution rate.
3. Pohl (1977) measurements: providing nitrogen retention data and
data on the rate of devolatilization of fuel nitrogen.
The aliphatic group, which is the collection of all measurable light
2-1
-------
hydrocarbons such as CH^, C2H2* and other higher HC, is treated as one single
component that is evolved at a single rate. The instantaneous secondary pyrolysis
which follows the initial pyrolysis decomposes aliphatics into kinetically more
manageable gas species such as CH^, ^2^ an^ ^2 or soot" i-n_US8
fractions of CH^, and are determined from MIT's asymptotic yield data
shown in columns 1 and 2 of Table 2. Whether or soot is yielded is determined
by the conservation of H ond C atoms in aliphatics.
The difference between total weight loss and light volatiles yield is
considered to be heavy volatile (or "tar") yield. The resemblance of Infrared
Spectra of Tar and Coal (Solomon, 1979) suggests that the tar consists of
"monmers" released from parent coal. Therefore, it is logical to assume that tar
has a similar structure to parent coal. However, due to difficulties in balancing
O-atoms in MIT data, basic contents of tar are modified to be similar to those of
the parent coal, except the solid CO group is assumed to be the only
oxygen-containing group in tar. There Is only one "tar" group in the EER rate
set. The instantaneous secondary pyrolysis of tar decomposes "tar" into its basic
contents (in gas form).
Volatile nitrogen is treated as a separate category since it is of major
interest. The asymptotic yield of volatile nitrogen is expressed as a function of
temperature by a cubic polynomial
N# (as volatile nitrogen) - A + BT + CT3
The values of A,B, and C used are listed in Table 2. The polynomial agrees
reasonably well with most of the experimental data (Blair, 1977; Pohl, 1977).
Experiments also revealed that the nominal pyrolysis product of volatile nitrogen
is HCN. The small amount of NH^ observed (Axworthy, et al., 197B) is neglected in
the current model. The balance of coal nitrogen, which stays with char, is
eventually oxidized along with the char.
Functional group partitions, coefficients for Arrhenius rate constants and
asymptotic yields currently used in EER flame codes are listed in Tables 1 and 2.
MIT (Suuberg et al.) and UTRC (Solomon et al.) rate constants are also listed in
Table 1. EER rate constants are practically the same as those of MIT, except for
2-2
-------
TABLE 1. RATE COEFFICIENTS
Func-
tional
Group
MIT
(Lignite)
MIT
(Bitjminous)
UTRC
(General)
EER
(Ligni te)
EER
(Bituminous) !
co2 (1)
co2 (2!
C02 (3)
2.138{11)/36.2*
5.128(13)/64.3
5.495(6)/42.0
1.0(135/40
1.0(135/65
75/11.6
2.138(115/36.2
3.65(135/64.3
1.0(13)/40
1.0(135/65
h2o (ir
H20 (2)
7.95(13)/51.4
i
1.0(13)/35 | 370/13.6
7.95(135/51.4
1 .0( 1 35/35.0
CO (1)
CO (2)
CO (3)
1.82(12)/44.4
2.63(12)/59.5
5.89(9)/58.4
1.0(13)/55
1.0(135/65
8.7(5)/24.0
4.0(9)/70.0
1.82(125/44.4
1.85(12)59.5
1.0(13)/55.0
1.0(13)/65.0
HCN
N (S)
-
-
200/15.2
290/25.4
9.3(3)/22.7
Same as Char
Rate
9.3(3)/22.7
Same as Char
Rate
H2
1.588(18)/88.8
1.0(17)/90
3600/25.4
1.585(18)/88.8
1.0(175/90.0
ch4 (1)
CH4 (2)
ch2h4
(1)
c2h4
(2)
C2H6
HC
1.62(14)51.6
4.677(4)/69.4
1.78(20)/74.8
7.08(12)/60.4
1.7(16)/70.1
1.0(13)/55
1.0(13)/65
1.0(13)/55
1.0(13}/65
1.0(13)/55
3.3(5)/24.0
6.48(13)/51.6
5.56(125/55.0
Tar (1)
Tar (2)
7.58(11)/37.4
2.0(17)/75.3
8.7(2)/l3.2
9.3(4)/19.6
1.3(7)/40.0
1.3(75/40.0
*A/E: for A exp (-E/RT), A 1s In (1/sec), E is in (Kcal/mole)
~Moisture
2-3
-------
Table 2. Asymptotic Yields {% Mass)
1 MIT
Ligni te
i
MIT UTRC
Bi tuminous jLignlte
UTRC
Bi tuminous
EER
Lignite
EER
Bituminous
co2 (1)
5.7
0.4
5.0
1.2
5.7
0.4
C02 (2)
2.70
0.9
-
-
3.79
0.9
C02 (3)
1.09
-
-
-
-
h2o (1)+
6.8
1.4
-
-
-
-
H20 (2)
9.7
5.4
6.0
2.4
13.0
6.8
CO (1)
1.77
0.4
9.0
1.3
1.77
0.4
CO (2)
5.35
2.1
19.9
9.8
7.61
2.1
CO (3)
2.25
-
-
-
HCN
f(T)*
f(T)
0.38
0.38
f(T)
f(T)
N (S)
f(T)
f(T)
0.8
1.2
f(T)
f(T)
H2
0.5
1.0
1.3
1.7
0.5
1.0
ch4 (1)
0.34
0.7
CH4 (2)
0.92
1.8
c2h4(D
0.15
0.2
17.0
24.1
2.77
8.6
C2H4(2)
0.41
0.6
C2H6
0.95
0.5
HC
-
4.8
Tar (1)
Tar (2)
2.45
2.93
-
16.0
f(T)
f(T)
A
-92.08
-92.08
B
0.129245
0.129245
c
-9.0290522(-9)
-9.0290522 (-9)
*As a function of particle temperature.
+Moisture
2-4
-------
minor variations in speciation and asymptotic yields. The variation is portly due
to atomic imbalance in data and partly due to the consideration of mathematical
convenience for further improvement of the model. Therefore, EER's model should
agree with MIT data if experimental coal temperature history is known.
Due to the lack of thermochemical information on pyrolysis reactions, the
devolatilization processes are assumed to be neutral thermal. Thus, the
thermochemical properties of functional groups are set equal to those of their
gas-phase counterparts. The heat of formation of char is treated as a variable to
accommodate the uncertainty and is calculated by
Ij^char
Hfcoa1 ~ £ ("ft /"!_)
,chlr (WW
where h_. is heat of formation at 298.15°K, MJ is molecular weight and YJ is
fi i 3 i®
asymptotic yield of group "i."
Heat of formation of coal, H ,, can be derived from the heating value
vcoal
(LHV) of coal by
hra2 * °-5 hm2o' wv
" V*V WStU
The apparent molecular formula is CHx®y"zAS' wliere A stands for ash.
It is understood that error in calculation of energy release might be
incurred locally due to the uncertainty as to thermochemical properties. However,
the model predictions are consistent with conservation of overall thermal-chemical
energy.
2-5
-------
2.0 CHAR COMBUSTION MODEL
Moss rate of oxidation of char is taken from the work of Smith (1971):
(dC) = 20.4 exp(-1900/R7) gm of carbon
02 cm - sec - atm of 02
Rate is calculated per unit external surface area of char particle.
At 1000°K the oxidation and gasification rates of char can be found from the
following relationships (Walker, 1959):
(d£) /(dC) = 3.16(- 3)
ar co2 01 o2
(dC) /(dC) = 1.732
at h2o co2
Thus, the gasification rate of char by C02 is
h r
(tt) = 192.0 exp(-35000/RT)
1 co2
and the gasification rate of char by H20 is
(4r) = 27.3 exp(-30000/RT) •
1 h2o
The elementary rate constants are related to mass rate constants by
6 RT_
Pr
K< ¦ (^> $>, ¦
where gas constont R - 82.06 (cm3-atm/mole-°K). T is particle temperature in °K
P 3
and dp is particle diameter in cm. pg is particle material density in gm/cm .
For maeroscole calculation, the elementary char oxidation and gasification
reactions are expressed as
1. C + 1/2 0 -» CO A - -21.6 Kcal/gmole
2. C + CO - 2C0 AH - +<»-1 .2 Kcal/gmole
2 C
2-6
-------
3. C + H^O ->¦ CO + L Hc ¦ +31.4 Kcal/gmole
In the computer programs, these reactions are built-in and are not listed in
Appendix A of Part 1.
Additional heterogeneous reactions listed in Appendix A of Part 1 are:
(i) Soot combustions (Reoctlons 181 and 182)
(ii) Soot NO reduction (Reaction 184)
(iii) Char nitrogen combustion (Reactions 187 and 188).
2-7
-------
3.0 MICROSCALE COMBUSTION OF A SINGLE COAL PARTICLE
Detailed mathematical formulation of the system of partial differential
equations which describe the physical and chemical processes around a burning coal
particle will not be presented here. Instead, only highlights of physical and
mathematical assumptions are listed:
1. Coal particles are assumed to be nonswelling and spherically
symmetrical.
2. Time-dependent one-dimensional (in spherical coordinates) boundary
layer equations and equations-of-state are employed.
3. Unequal binary diffusion coefficients, based on the bifurcation
approximation, are utilized with Fick's law of species diffusion.
-------
RADIATION
;ADIATION
FLAME ZONE
1 VOLA
TILE
IMPERMEABLE
BOUNDARY
Figure 1. Schematic of Microscale Combustion of a Coal
Particle in a Confined Space.
2-9
-------
<~.0 SUDDEN IMMERSION OF COAL PARTICLE INTO HOT COMBUSTION PRODUCTS
In order to look into the microphysical and chemical processes of cool
particles during their lifetimes in typical coal-fired furnaces, in this section
numerical simulation of a single cool particle suddenly immersed into a puff of
lean combustion products will be examined.
Particle sizes are chosen to be 50y m and 100 um , which covers the range of
average particle sizes in coal combustors. Combustor wall temperatures are chosen
o o
to be 1500 K and 1900 K. The particle is assumed to be initially surrounded by
equilibrium combustion products of CH^ at an equivalence ratio of 1.67. The
initial ambient gas is assumed to be at its equilibrium temperature (i.e. T
o »
1662 K). The thermodynamic state of the gas and the ratio of interparticle
distance to particle size are chosen independent of wall temperature and are
functions of particle size alone. If we choose the interparticle distance to be
57 times the particle diameter, the final combustion products of the system will
maintain approximately two percent oxygen. The initial particle temperature is
300°K.
Figure 2 shows the rate of devolatilization of volatiles and the rote of
combustion of char for a 50 u m bituminous particle at T - 1500°K. The first
wall
rate peak at 2 ms was determined to be from water contents in the coal. The
second peak is the contribution from evolution of CO^. CO and aliphatic "tar"
evolved mainly between tightly bounded light volatiles (e.g. CO, CO^, and H^) and
combustion of char. The lifetime of a 50 u m coal particle is about 95 msec while
the lifetime of a 100um particle is about 350 msec. Doubling the particle size
approximately quadruples the burn-out time.
Figure 3 shows the temperature histories of particle and gas (far field)
under a higher heat loss environment (Twall " 1500°K), Increasing particle size
reduces heating rate which in turn reduces peak devolatilization rate. The
ignition time in gas phase is increased from approximately 5 msec to 15 msec when
particle size is increased from 50 um to 100 um .
Similar results for particles in lower heat loss environment (T - 1800°K)
r wall
are shown in Figure k and Figure 5.
2-10
-------
.0
1500
2% Excess 0
m
Bi tuminous
50 urn
100 ym
C\J
u
1/1
Li-
LO
LO
2
Tar
Char,
O
ro
3
10
10
10
10
t (Sec)
Figure 2. Mass flux rate of 50 ym and 100 urn, bituminous particle
in higher heat loss environment (T 11 = 1500OK).
2-11
-------
2200
2000
T (100 um)
1800
1600
1400
£ 1200
CL
1000
Microscale
wal 1
400
200
140 160
100
120
80
t (ms)
Figure 3. Particle and Gas Phase Temperature Histories for
2-12
T ,, = 1500°K.
wal 1
-------
10'
2% Excess 0,
O
in
a;
to
ca
o
o
O)
-------
2200
(for 50 urn)
Tq(for 100 jjT^
2000
50 urn
100 um
Burned Out
1800
1600
1400
1200
1000
800
Mi croscale
wal 1
600
1662 K
400
200
140 160
100
120
t (ms)
Figure 5. Particle and Gas Phase Temperature Histories
for T ,, = 1800°K.
wa 11
2-14
-------
Temperature distributions for three different stages (for Twall " 1800°K) of
combustion of a 100um particle are plotted in Figure 6. At t - 30 ms (final
stage of energetic light volatile evolution) the temperature peak is very close to
o
the particle surface. The temperature over—shoot is about 130 K. At t = 50 ms
(during the state of "tar" evolution) the temperature peak has moved away from the
o
particle surface. However, the temperature over-shoot increases to 370 K. At t =
75 ms (during char combustion stage) the temperature peak is not obvious since the
char combustion rate is slower than the devolatilization rate, and thus the
reaction time becomes comparable to the diffusion time. The diffusion flame is
not established. The system is controlled by char combustion rate. Oxygen
distributions are also shown in Figure 7.
One of the main objectives of microscale simulation of the particle is to
investigate effects of the microscale process of fuel bounded nitrogen conversion
to fixed nitrogen species (NO^, NH^, and CHN). Figure 8 shows the production rate
of NO around a lOOum particle (T^ - 1800 K) at three different stages. In
general, gas-phase NO reactions occur within 10 particle radii of the particle
surface. Kinetic information can be calculated u6ing equations found in the
computer code. The percentage contribution to the production of each species by
each reaction can be readily obtained.
The reactions which contribute significantly to the direct production of NO
at t • 50 ms are as follows (The overall contribution refers to the spatial
integration of the locol NO production rate):
Reaction No.
(99)
(89)
(112)
(92)
(128)
H20 + NO
HN0 -~ H + NO
HNO + OH ->
-+ OH + NO
+ H + NO
NCO + 0 -»¦ CO + NO
N20 + H
N2 + OH
% Contribution (Vol. Rate)
Overall At Tmax
27.0
20.3
20.2
6.8
6.8
*f*f.O
11.3
<0.5
17.5
<0.5
HNO and NCO are the important intermediates in CHN and NH^ oxidation
processes. In these processes, NO from the overall thermal process [reaction
(92)] contributes less than 7 percent; however, the local (peak temperature zone)
2-15
-------
2300
2300
Microscale
dp = 100 u Bituminous
2200
2200
2100
2100
T = 300 K
po
2% Excess Or, Overall
2000
2000
75 ms
1900
1900
ms
CD
h-
30 ms
1800
1800
1700
1700
1600
1600
1500
1500
0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22
Distance From the Center of the Coal, r (cm)
Figure 6. Temperature Distribution Around a Coal Particle.
-------
N
I
01
e
3
OJ
o
at
o
L.
aj
a.
4.0
t = 50 ins
Micro-Scale
dp = 100 p Bituminous
o
1662 K
T ,, = 1800 K,
-------
.0
Microscale
dp = 100 urn Bituminous
T ,, = 1800°K T
wall po
T = 1662 K
U
50 ms
3.0
CO
30 ms
2.0
75 ms
1.0
0
1.0
0.12
0.10
0.08
0.06
0.04
0.02
Distance From the Center of the Coal Particle (cm)
Figure 8. Distribution of NO Production Rate Around a Coal Particle.
2-18
-------
thermal NO contribution is as high as 17.5 percent.
The production rotes of at the three stages ore shown in Figure 9. The
key reactions which contribute to the direct production of N^ at t = 50 ms ore as
follows:
% Contribution (Vol. Rate)
Reoctlon No. Reactions Overall At Tmax
(183) (N) + (N) ->• N2 72.6 53.1
(96) NO + N + 0 + N2 12.0 M.5
(108) N20 + H ->¦ OH + »2 9.3 1.9
Reaction (183) represents the path by which 90 percent of char nitrogen converts
to N^. The reaction is an ad hoc simulation of char nitrogen conversion and
should be improved later. The hey reaction for N2 production is reaction (96),
which is the well known reverse Zeldovlch mechanism. This reaction is responsible
for 12 percent of the overall N2 production and contributes 41.5 percent at its
peak temperature.
2-19
-------
Microscale
dp = 100 ym Bituminous
Twa11 = 180n°°K Tro
T = 1662 K 6 =
go go
2% Excess 0o Overall
50 ms
IT)
O
«/>
1.5 —
at
75 ms
0.5-
-0.5
0.14
0.12
0.10
0.08
0.04
0.06
0.02
Distance From the Center of the Coal Particle (cm)
Figure 9. Distribution of Production Rate Around a Coal Particle.
2-20
-------
5.0 MACROSCALE VS. MICROSCALE
Data from Pershing's (1981) experimental work suggests that the presence of
char significantly reduces the yield of NO from fuel nitrogen species (NH^ was
used in the experiments).
It was speculated that microscale phenomena surrounding the coal particle
might have a significant impact upon fuel nitrogen species conversion to NO. In
order to examine the impact af microscale phenomena on fuel nitrogen conversion,
numerical experiments were conducted using the microscale model described in
previous sections. The computed results are compared to the macroscale results
obtained from OFAP (One-Dimensional Flame Analysis Program). OFAP codes operate
in plug flow mode where the coal particle has identical chemical characteristics
to those used in microscale calculations. However, the volatlles evolved are
assumed to mix instantaneously with the surroundings. Therefore, no microscale
structures around particles are considered.
Table 3 shows the NO concentration and percent conversion of fuel nitrogen to
NO from microscale and macroscale calculations. For smaller particles (50ym ),
the results from microscale and macroscale calculations are almost the same for
o
the lower heating case (Twau ¦ 1500 K). There is about 2 percent deviation due
to the microscale structure at the higher heating case " 1800 K). For
larger coal particles (100um ) the microscale effects are somewhat more
significant. The deviations are 1.8 percent and 11 percent for Twa - 1500°K and
o
T - 1800 K, respectively,
wall
For heterogeneous NO reduction reaction by char, I.e.,
NO + CO + 1/2N2 + C02
the microscale simulation calculations are also listed in Table 3. The results
show that heterogeneous NO reduction reaction has only a modest effect on overall
NO concentration. For 100 un particles, the effect is about 5 percent. The rote
constant used in these calculations Is derived from DeSoete's work (1977):
"NO s = 1-75 X 10? EXP (_29'600/RT) p 0'7 CN0 CC0 cm2_sec
This reaction is built into the computer code and is not listed in Appendix A.
C. „ and are in (mole/cc), W is the mixture molecular weight, and d is the
NO CO
2-21
-------
Table 3. PPM levels of NO arid Percentage Conversion of Fuel-N*
Calculated by Microscale and Macroscale Models.
dp = 50 urci
Twa)l
Macroscale
Mi croscale
Microscale
(With Heterogeneous
NO Reduction**)
1500
625 ppm/77.5%
621 ppm/77.0%
-
1800
886 ppm/109%
871 ppm/107%
-
dp = 100 ym
Twal 1
Macroscale
Microscale
Mi croscale
CWith Heterogeneous
NO Reduction**)
1500
592 ppm/73.4%
578 ppm/71.5%
540 ppm/66.8%
1800
909 ppm/111%
826 ppm/100%
752 ppm/94.6%
*Pittsburgh Bituminous
**no + co -~ y*2 + co2
2-22
-------
mixture density. Thus, the microscale processes simulated here do not explain the
degree of Impact of char on NO conversion which was observed in Pershing's work.
Pershing's work shows only 25-30 percent conversion of fuel nitrogen to NO for
most of the coal in various ranges of excess oxygen. Our microscale model
predicts 70-100 percent conversion, which is far off. Possible areas of
improvement in the microscale model are:
1. Fuel nitrogen pyrolysis chemistry. The current model assumes
instant conversion of fuel nitrogen to CHN upon devolatilization,
which might partially be responsible for overprediction of NO
conversion. Some intermediates of fuel-n furing pyrolysis might
have an effect on final conversion.
2. Gas-phase CHN/NH^ chemistry. The validity of the current CHN/NH^
chemistry may have to be reexamined, particularly in rich mixtures;
Exxon data suggest that our chemistry set seriously overpredicts NH^
conversion to NO.
3. Heterogeneous NO reduction by soot. Although the mechanism of
NO reduction by soot is included in these computations (Reaction
184), the soot formation mechanism in our model is an ad hoc
mechanism (based only on an atom balance of the aliphatic group).
The impact of soot concentration on fuel-N conversion will be
reexamined when the mechanism of soot formation Is better
understood.
2-23
-------
6.0 CHARACTERISTIC TIME SCALES ON SINGLE PARTICLE COMBUSTION
Several physical and chemical characteristic time scales, which may help to
better understand the processes occurring around a single burning coal particle,
are listed in Table k. The definitions of those characteristic times are
described as follows.
1. Heat-up time: the time interval required for the coal particle to
be heated by the amount which is equal to 95 percent of the
difference between and initial particle temperature.
2. Diffusion time: the time interval required for a gas molecule to
diffuse from the particle surface to the outer boundary of the gas
2 —
puff. This is estimated by L /D where L is the dimension of
g g g
the gas puff. is the mean diffusivity of the gas molecule.
3. Devolatilization time for light volatiles/tar/volatile nitrogen:
the time at which peak evolution rate occurs. Values listed in
Table are taken directly from Figure 2 and Figure (f. For light
volatiles, the last peak (for the tightest bounded species) was
taken.
k. Particle life time: the period during which char is completely
burned.
5. CO combustion time: an estimate of the time required for CO to be
burned. The values (order of magnitude only) are taken to be the
reciprocal of the multiplication of ^i^wall^ mean ^0
concentration within the zone of influence.
6. NO formation time: an approximate estimation of the time required
for NO to be formed from fuel nitrogen (i.e., CHN) destruction. The
values (order of magnitude only) are taken to be the reciprocal of
the multiplication of by the mean OH concentration within
the zone of influence.
7. Nj formation time: a rough estimate of the time required for N^ to
be formed from NO which is presumably coming from fuel nitrogen.
The values (order of magnitude only) are taken to be the reciprocal
of the multiplication of K_.(T ,,) by the mean N atom concentration
96 wo11
within the zone of influence.
2-24
-------
Table 4. Physical and Chemical Time Scale
Time Scale
dp = 50 um
dp = 100 ym
(ms)
Twair'50°O|<
Twair1800°K
i
Twair1500°K
Twair1800°K
Heat-up Time,
7.0
25.0
25.0
60.0
Diffusion Time, t^f
ro
7.2
28.8
28.8
Devolati1i zati on
Time, idv
Light Volatile
7.0
6.0
20.0
18.0
Tar
25.0
20.0
55.0
42.0
Volatile N
50.0
30.0
75.0
60.0
Particle Life
Time, xp
95
62
350
180
CO Combustion
Time, tco
<1.0
<1.0
<1.0
<1.0
NO Formation
Time, t
<0.01
<0.01
<0.01
<0.01
N„ Formation
Time, t,,
N2
-10
-10
-10
-10
CHN Destruction
Time, tchn
<0.1
<0.1
<0.1
<0.1
Ignition Time, ijg
4.7
4.2
15.5
13.5
Heterogeneous NO
Reduction UmpHeter
104
io4
io4
io4
2-25
-------
8.
9.
10.
It is noted that for all cases the particle heat-up time is either comparable to
or longer than light volatile evolution time. The major portion of light
volatiles are evolved during the particle heat-up period. The diffusion time is
either comparable to or shorter than the devolatilization time in most
cases.
For a 100 um particle, the zone of action (r < 0.05 cm) only represents
approximately one-half percent of the total volume (or mass). The characteristic
convective time in this region, which is estimated by
_ Length scale of the zone (Li)
•conv Average convective velocity in the zone (V)
is about 0.7 ms (- 0.05/70).
The characteristic diffusion time in the 2one of action, which is defined by
2 — 2
t /Di is abou't 2.5 ms (- (0.05) /1.0). These times are all
comparatively shorter than the formation time. However, they are either
shorter than or comparable to CO combustion time and much longer than either NO
formation or CHN destruction time. This suggests that although NO will be formed
from fuel nitrogen within the zone of action, it will not be reduced within the
zone of action either.
The heterogeneous reaction times are on the order of minutes, which is too
long to make heterogeneous reaction significant in NO reduction.
2-26
CHN destruction time: an estimate of the time required for CHN to
be destroyed. Values are taken to be the reciprocal of the
multiplication of bY the mean OH concentration within
the zone of influence.
Ignition time: the time at which gas temperature suddenly
increases. The values are taken directly from Figure 3 and Figure
5.
Heterogeneous NO reduction time is defined as
it overall NO in the puff (mole)
NO Heter - 7"^ z ?, rf.—\ \ , '2 T
(particle surface area, cm ) \WN0)S» mole/cm -sec
-------
7.0 REFERENCES
1. Axworthy, A. E., Dayan, V. H.p and Martin, G. 0. "Reactions of
Fuel-Nitrogen Compounds Under Conditions of Inert Pyrolysis." Fuel. Vol.
57, p. 29, 1978.
2. Blair, D. W., Wendt, J. 0. L., and Bartok, W. "Evolution of Nitrogen and
Other Species During Controlled Pyrolysis of Coal." Sixteenth Symposium
(International) on Combustion. The Combustion Institute (1977).
3. DeSoete, G. G., "Mechanisms of Nitric Oxide Reduction on Solid
Particles." Proceedings of the Fifth EPA Fundamental Combustion Research
Contractors Workshop, EPA-600/9-83-001 (NTIS PB83-164483) (1983), p. 174.
4. Kobayashi, H., Howard, J. B., and Sarofim, A. F. "Coal Devolatilization
at High Temperatures." Sixteenth Symposium (International) on
Combustion. The Combustion Institute (1977), p. 411.
5. Chen. S. L.. Pershing, D. W.. Heap, M. P., and Martin, G. B. "Influence
of Coal Composition on the Fate of Volatile and Char Nitrogen During
Combustion." Nineteenth Symposium (International) on Combustion. The
Combustion Institute (1982).
6. Pohl, J. H., and Sarofim, A. F. "Devolatilization and Oxidation of Coal
Nitrogen." Sixteenth Symposium (International) on Combustion. The
Combustion Institute (1977), p. 491.
7. Smith, I. W. "The Kinetics of Combustion of Pulverized Semi-anthracite in
Temperature Range 1400-2200 K." Combustion ond Flame. 17, 421, (1971).
8. Solomon, P. R., and Colket, M. B. "Coal Devolatilization." Seventeenth
Symposium (International) on Combustion. The Combustion Institute (1979),
p.131.
9. Suuberg, E. M., Peters, W. R., and Howard, J. B. "Product Compositions
and Formation Kinetics in Rapid Pyrolysis of Pulverized Coal—Implication
for Combustion." Seventeenth Symposium (International) on Combustion.
The Combustion Institute (1979), p. 117.
10. Walker, D. L., Rusinko, F.. and Austin, L. G. "Gas Reactions of Carbon."
Advances in Cotolysls. 11, 133 (1959).
2-27
-------
i
i
i
i
-------
PART 3
DATA ANALYSIS THROUGH
INVERSE TECHNIQUES
Prepared by:
C. J. Kau
and
T. J. Tyson
3-1
-------
TABLE OF CONTENTS
Section Page
1.0 INTRODUCTION 3-1
2.0 ANALYSIS 3-2
2.1 Basic Equations 3-2
2.2 Inverse Analysis 3-3
2.3 Numerical Consideration 3-k
3.0 SAMPLE CASE 3-6
4.0 CONCLUSION 3-9
5.0 REFERENCES 3"9
6.0 NOMENCLATURE 3~10
LIST OF FIGURES
Figure Page
1 Comparisons of the Distribution of Viscosities 3-7
2 Comparison of the Distribution of NO Production Rates .... 3-8
3-11
-------
1.0 INTRODUCTION
Due to the great volume of field data obtained from experimental combustion
devices, a need has arisen to develop an analytical tool to use these data to
better understand turbulence/chemistry interaction and pollutant formation.
Experimental data usually concern flow field variables such as velocity compo-
nents, temperature and species concentrations. It has been the environmental
engineer's dream to have an analytical tool to calculate distribution of turbulent
transport properties of an entire flow field. This would eventually allow calcu-
lation of the distribution of pollutant formation rates, such as NO production
rate, in a flame zone. Thus, ultimately, an effective control strategy could be
devised. This methodology is called "Inverse Analysis."
Generally speaking there are two types of methods that have been used. The
first method is the control volume method (Ref. 1). This method requires mapping
out of the streamlines to form the field data, and then subdividing the flow field
into small control volumes between streamlines. The turbulent eddy diffusivity is
then derived from a material balance of nonreactive tracer in the control element.
Using the eddy diffusivity, the production rate of any chemical species can be
calculated by balancing turbulent diffusive flux and convective flux. The draw-
backs of this method are:
(i) Requiring tracer
(ii) Mapping of the flow field
The second method is a differential equation method (Ref. 2). One has to write
out an entire set of the governing equotions and then evoluate the evaluable terms
in these differential equations from measured data. The turbulent eddy
diffusivity is the only unknown in the momentum equation to be derived from
measured velocity. The production rate is the only unknown in the species equa-
tion ond can thus be evaluated. This method does not require tracer; however, it
requires very smooth data and cannot handle reclrculatary flows due to the fact
that there are more unknown turbulent stress terms than the number of governing
momentum equations.
The differential equation method for analyzing weakly swirling flow will be
described and tested in the following sections.
3-1
-------
2.0 ANALYSIS
2.1 Basic Equations
The vector forms of the time-averaged continuity equation, Navier-Stokes
equation of motion, and mass conservation equation for a turbulent reacting flow
can be written:
§!+oU-u)=0 (1)
p = -AP + A'T (2)
P ^ = - a ^ (3)
where t is a Reynolds stress tensor, is the mass flux vector of the ith
species, and is the chemical production rate of ith species. Y^ is the mass
fraction of the 1th species. Mechanical dissipation and external body forces are
neglected in the above equations.
For description of a steady axisymmetric flow field, it is convenient to
write the equations in a cylindrical coordinate system; i.e.,
Continuity Equation
£ (cu2) + (rpur) - 0 (4)
Axial (z-direction) Momentum Equation
3 ~!&(«„) -& (5)
Radial (r-direction) Momentum Equation
p (u + u iHr . Hfil) = A (T j + 1 J_ (rT ) . M . IE. (6)
\ Z 32 uf ar r ) 3Z KTzr> r 3z K Trr> r 3r '
Azimuthal (8-direction) Momentum Equation
/ 3ug 3ue Uru0\ 3 . 1 3 ,2 * m
c K IT * "<¦ IF + —) = « 'Tez' sF Trs> <7)
Species Conservation Equation
p (uz TF+ ur ir)= JI (Jiz) + r W (rJir} + "l (8)
3-2
-------
For weakly swirling flow, there are no recirculating zones in the flow field
and a boundary layer type of approximation can be employed to further simplify the
above equations. By invoking a boundary layer approximation, via order of magni-
tude argument, the above equations can be reduced to:
Continuity Equation
to (ouz' + = 0 (9)
Momentum Equations
I 3uz «u z \ 3P 1 3 i \ /-ifiv
p (uzlT + ur —J + S? FT O0)
P^-l§ (ID
r or
+ ueM _ 1
Jz^FTUr^ 1 r9;
Species Equation
/ 3Yi 1 5 , , v •
0 (uz TF + ur ' r Jz ^ ir = ui
/ 9u0 . 3U0 + "eM 13/2 ,
P + — = (r T ) (12)
(13)
Generally, Eqs. (9)—(13) and an equation for conservation of energy are used to
predict the flow field if boundary conditions of the flow region are known and if
Reynolds stresses, diffusive flux, and chemical production term can be adequately
modeled. For chemically reacting turbulent flow, the physics of the interaction
between turbulence and chemistry are largely unclear. It is difficult to model
these terms in an adequate manner; however, when the details of the flow field are
measurable, the known information (i.e. velocities, temperature, pressure, and
species distributions) can be used to compute the Reynolds stresses (or turbulent
viscosity) and chemical production rates by using the same set of equations.
2.2 Inverse Analysis
When the time-mean velocities, temperature, pressure and species distribu-
tions of a flow field are measurable, the ideal gas law
p?¥,-
0 = ] con be used to calculate the density distribution. The terms on
3-3
-------
the left-hand side of equation (10) can readily be evaluated. Numerically inte-
grating equation (10) in the r-direction gives the Trz(z,r) distribution.
Similarly, t „ (z,r) can be calculated by integrating Equation (12). If pres-
re
sure is known only at the outer edge of the flow field [i.e. P(z,°° ) ], Equation
(11) must be used to integrate inwardly in the r-direction for P(z,r), If u (z,r)
are not given, the continuity equation [Equation (9)] can be used to solve for
them.
Apparent turbulent viscosities can be defined by analogy to Newton's law of
friction for laminar flow; i.e.,
T
"rZ
rz
and
ur
rG
(3uz/3r)
(>fl/r)
J? K/r)
Or ^
Similarly, apparent species diffusivity ( aj) of turbulence can be defined by
analogy to Pick's law of diffusion; i.e.,
3V.
i r 1 3r
For engineering purposes, it is generally assumed that the ratio of turbulent
viscosity to turbulent mass diffusivity Is approximately constant. Thus, the
species flux can be evaluated from the known value of u . Schmidt number and
measurable species distribution; thus
u 3Y.
i = -II _L
1 r sci 3r
where Schmidt number it defined as S - -5-?. Then the chemical production rate
ci a i
of species 1 can be evaluated by Equation (13).
2.3 Numerical Consideration
The experimental data are usually very crude and unsmooth. Generally,
smoothing of the raw data before the inverse analysis process is not only desir-
3-4
-------
able but also necessary since we are interested in derivatives of the data.
In the "inverse" analysis code, the experimental data are input tabularly as
functions of two physical coordinates. A subroutine named SLP is used to compute
the directional derivatives for the tabulated function. SLP calculates the end
point derivatives by parabolic interpolation while the derivatives of the interior
points are evaluated by a cubic spline fit procedure. If spatial resolution
higher than the resolution of original input data is required, a subroutine named
CUBIC is used to perform the cubic interpolation and generate an extended data
table. The derivatives required for cubic interpolation are calculated from
original data by subroutine SLP. The derivatives required for integrating the
governing partial differential equation are also computed by SLP using this
extended tabular function.
3-5
-------
3.0 SAMPLE CASE
The first sample case chosen to examine the correctness of the code's arith-
metic is an NHj-doped hydrogen laminar diffusion flame. The flame parameters are:
Fuel Tube: 0.236 cm ID
0.635 cm 0D
Oxidant Tube: 5.08 cm ID
Fuel/Oxidant Ratio: 167 percent T.A.
Percent 0^ in Oxidant: 10.51 (vol)
Percent NH^ in Fuel: '~.Ol (vol)
Fuel Flow: 2178 cm/sec
Oxidant Flow: cm/sec
o
Tube Wall Temperature: 353 K
The fuel nitrogen conversion experiment in the diffusion flame has.been physically
carried out in EER's El Toro laboratory (Ref. 3); however, the flame was not
probed in detail, only stake values were measured. For the purpose of program
checkout, the "experimental data" of the flame is actually "simulated" in detail
using EER's diffusion code (GFAP, in diffusion flame mode). The "numerical
simulated experimental data" is then fed into the "Inverse Analysis" model to back
out viscosity and NO production. It should be noted that the choice of laminar
flame instead of turbulent flow will not affect our program checkout since laminar
and turbulent flows are described by the same set of PDE.
The viscosity distribution computed from "Inverse Analysis" code is shown in
Figure 1 (at two locations downstream from the exit of the fuel tube) in compari-
son with the "viscosity" used to simulate the "experimental data." The NO produc-
tion rates obtained from the "Inverse Analysis" code are shown in Figure 2 in
comparison with the experimental data "simulated." The slight deviations in both
comparisons probably can be attributed to the fundamental differences between
derivative evaluation methods. In inverse analysis code a cubic spline fit
procedure was used to evaluate derivatives while a central differencing scheme was
used in the "experimental data simulation" model.
3-6
-------
Simulated Experimental -
Data (by GFAP)
7.
Predicted by Inverse
Analysis
6.
5.
4.
3.
2.
2.5
2.0
1.0
0
Figure 1. Comparisons of the Distribution of
Viscosities.
3-7
-------
Simulated Experimental
Data (by GFAP)
Predicted by Inverse
Analysis f*
-7 —
U
o
o
2:
Figure 2. Comparison of the Distribution of NO Production
Rates.
3-8
-------
t*. 0 CONCLUSION
So far, the "Inverse Analysis" code has been successfully tested only on very
smooth data (e.g., computer-generated). The tests using raw field data such as
the combustion data from IFRF (Ref. <~) were unsuccessful due to extreme sensitivi-
ty of the solutions to the quality of data, due to the fact that we have evaluated
the divergent terms in the equation from sparsely measured data points. Several
data smoothing techniques have been tried and the results obtained are disappoint-
ing.
5.0 REFERENCES
1. Sadaka, M. and Beer, J. M. "Spatial Distribution of Nitric Oxide
Formation Rate in Swirling Turbulent Methane-Air Flame." Proceedings of
the Sixteenth (International) on Combustion. The Combustion Institute
(1976). p. 93.
2. Lilley, D. G., and Chigier, N. A. "Nonisotropic Exchange Coefficients in
Turbulent Swirling Flames from Mean value Distributions." Combustion and
Flame. J6, p. 177 (1971).
3. Folsom, B. A., Courtney, C. W., Corley, T. L., and Clark, W. D.,
"Advanced Combustion Concepts for Low Btu Gas Combustion,* Proceedings
of the Third Stationary Source Combustion Symposium; Volume II.
EPA-600/7-79-050b (NTIS PB292540) (1979). p. 163.
it. Michels, J., and Payne, R., "Fundamental Combustion Research Applied to
Pollution Formation, Vol. Ila: Physics and Chemistry of Two-Phase
Systems: Flame Combustion Processes, Part II: Detailed Measurement of
Long Pulverized Coal Flames for the Characterization of Pollutant
Formation," International Flame Research Foundation Report, 1982.
3-9
-------
6.0 NOMENCLATURE
ith species flux
Molecular weight of ith species
P Static pressure
R Gas Constant
r Radial coordinates
Sci Schmidt number (ith species)
T Temperature
t T ime
u Velocity
Species mole fraction
Species mass fraction
z Axial coordinates
Greek:
p Density
a Turbulent species diffusivity
9 Azimuthal coordinates
t Reynold stress
Chemical production rate of ith species
w Turbulent eddy viscosity
Subscripts:
i ith species
r Radial direction
z Axial direction
6 Azimuthal direction
3-10
-------
PART 4
MICROSCOPIC MIXING IN TURBULENT
DIFFUSION FLAMES
Prepared by:
C. J. Kau
and
T. J. Tyson
4-1
-------
TABLE OF CONTENTS
Section Page
1.0 INTRODUCTION 4-1
2.0 MODEL DEVELOPMENT 4-2
3.0 NUMERICAL PROCEDURES AND PHYSICAL ASSUMPTIONS 4-10
4.0 TUNING OF THE MODEL 4-11
5.0 REFERENCES 4-14
LIST OF FIGURES
Figure Page
1 Schematic of Turbulent Free Jet 4-3
/mt\ / o.
2 Flame Length vs. Np) 4-5
\ o/ tt V IB
3 Schematic of Reactor Model 4-8
4 Effects of Uncertainty on Flamelet Residual Time 4-12
5 Fine Tuning of a and 6 4-13
4-11
-------
1.0 INTRODUCTION
New information on the rate of microscopic mixing in turbulent diffusion
flames discovered by re-examination of classic jet diffusion flame data obtained
by Hawthorne, Weddel and Hottel (Ref. 1), and the later data by Bilger and Beck
(Ref. 2) has led to the development of a simple turbulent diffusion model. The
simple model postulated here consists of a well-stirred reactor and a plug-flow
reactor to simulate the turbulent flamelet and flame-core in a turbulent diffu-
sion flame. The model is used to predict NO formation in H^/air flames. Compari-
sons are made against experimental data in terms of NO concentrations. Calculated
-1/2
peak NO concentrations for various flames show Re dependency, which has also
been suggested in the work of Bilger and Beck (Ref. 2).
4-1
-------
2.0 MODEL DEVELOPMENT
Considering a circular Jet of fuel issuing into stationary ambient air (see
Figure 1), the momentum flux based on macroscopic entrainment is:
2 2 2 2
J = Hp u; R„ = nPn uf FT
hi m m ooo
The mean local velocity can be expressed in terms of initial density p and
velocity Uq as
/p. U R
u„-\Ar -§^ (1»
" »P, Rm
The mass flux per unit Initial jet mass flux is
\ I-pmumRm _ / pm ^
Hp U R2
K0 0 0
(2)
This expression is obtained by elimination of U through Equation (1).
m
Since the flame spreading angle is approximately independent of fuel type and
jst conditions, we have:
R - k2X (3)
where is the tangent of half Jet angle. The virtue origin X - 0 is a distance
(- upstream of the nozzle exit plane. This distance is very small in
comparison with the physical dimension of the flame and can be neglected. Substi-
tuting Eq. (3) into Eq. (2) gives:
The turbulent entrainment model low-derived by Ricou and Spalding rscornnended -
0.16. The axial distance with a turbulent-jet-entrained stoichiometric amount of
air is evaluated from
-------
Figure 1. Schematic of Turbulent Free Jet.
-------
Eq. (5) Is plotted in Figure 2, where macroscopic mean mixture density p is
rtl
token to be p /2.
a
Recently, Broadwell (Ref. k) hos re-examined the classic turbulent-jet
diffusion flame length data obtained by Hawthorne, Weddel and Hottel (Ref. 1) and
has discovered that the flame length, which is a legitimate representation of the
microscopic mixing distance of the stoichiometric amount of fuel and air, is also
linearly proportional to 5t • Data from Hawthorne's (Ref. 1) paper and a
later paper by Bilger and Beck (Ref. 2) are plotted in Figure 2. These data
represent various types of fuel over small range of Richardson numbers. It is
recognized that the buoyancy forces do effect flame length. The flame length here
is the analytical flame length, which is the axial distance of the point of 99
percent completion of combustion (based on sample gas). The slope of the straight
line best fit to the data point is 9.2 percent. Figure 2 suggests a constant
ratio of microscopic mixing rate to macroscopic mixing rate; i.e.,
m . - k rft
micro 1 macro »
Figure 2 implies that k - 0.339. Thus, microscopic mixing rate is approximately
one-third of macroscopic mixing rate
L
„c , . „ „ (6)
A logical extension of the above conclusion is that the microscopic mixing rote is
linearly proportional to
'iTk
micro
(7)
where t<3 - 2h.,k2 - 0.1085.
The definition of velocity and Eq. (1) give
Um E 3* =. U0Ko
dx _ n%r
dt
k2X
p
Integrating the above equation with respect to t by assuming a constant mean
results in the following space-time relation:
X
D„ ~ (8)
4-4
-------
300
200 —
100
O CH/
c3h8
^2h2 \ Hawthorne, et al^ /
/ 9.217
/
_ ~ H2
x CO
(3) /
H2 Bilger and Beck* ' /
/
J
A
/
/
/
/
/
/
%
/
/
/
/
3,25
20
30
VI
Figure 2.
0/ st
Flame Length vs.
40
50
4-5
-------
Substituting Eq. (8) into Eq. (7) results in an expression for microscopic mixing
in terms of time:
0 / micro
(9)
Taking the time derivative of the above equation gives microscopic air mixing
rate:
where AF is air/fuel mass ratio. The overall unburned raw fuel left In the flame
per unit initial raw fuel flux is evaluated by
Experimental evidence of turbulent-Jet mixing suggests that the macroscopically
entrained air is first engulfed into a large coherent structure and is subsequent-
ly cascaded down to the Kolmogoroff scale and burned completely. The flame
surface per unit volume of the turbulent structure is considered to be small until
the Kolmogoroff scale is reached; therefore, a simple turbulent diffusion flame
model is postulated which contains two major zones: a flamelet zone and a
flame-core zone. The flamelet Is the Kolmogoroff scale structure which contains
stoichiometric (or near stoichiometric) amounts of microscopically mixed fuel and
olr. The flamelet structure is subsequently destroyed by the molecular diffusion
process and conglomerated into the flame core structure which contains the micro-
scopic mixture of combustion products and yet unreacted fuel. It is recommended
that the flamelet structure be simulated by a well-stirred reactor with proper
stolchiometry and residence time and that the flame-core zone be simulated by a
plug flow reactor. The plug flow reactor is continuously fed by combustion
products from the flamelet stirred reactor. A schematic of the reactors is shown
(10)
The microscopic fuel mixing rate Is
(11)
(12)
4-6
-------
in Figure 3. The flamelet fresh air mixing rate is calculated from Eq. (10),
while raw fuel mixing rate is evaluated from Eq. (11). The diluted (by combustion
products) fuel mixing rats is calculated by:
w
The proper residence time for the flamelet reactor should correspond to the
molecular viscosity and the Kolmogoroff scale:
\2
T = C, ±-
I v
where is proportionality constant and \> is kinematic viscosity. The
Kolmogoroff scale Is evaluated by
where is another proportionality constant. The rate of dissipation of turbu-
lence is estimated by
U 3
e ¦ c, 4-
3 D
C1C22 vV»
7^ Um3/2
Eliminating U using the expression following Eq. (7), the residence time expres-
fTl
sion becomes:
Q.\W/*v\ Dq ' \ /x_Y
(14)
2
where B = CiC2 ^C3 iS 0 maJ°r Physical uncertainty of the model. Expressing the
residence time in terms of clock time using Eq. (8) gives
/ Q-\h t
' rn %
T"B4k'U) ^7 051
-------
micro
_d_ / mfd\
\ mo / micro
FLAME CORE
P.F.R.
Figure 3. Schematics of Reactor Model.
4-8
-------
where Reynolds number ReQ is defined as
U D«
Ren - -V1
o v
v here will be evaluated at some typical flame temperature instead of that of
the exit plane.
4-9
-------
3.0 NUMERICAL PROCEDURES AND PHYSICAL ASSUMPTIONS
A computer model has been developed for the interacting well-stirred reoctor
ond plug flow reoctor described in the previous section. The highlights of the
physical and mathematical aspects of this code ore:
1. Implicit finite-difference method is used to solve time-dependent
governing equations (continuity, species and energy conservation
equations and equation of state) implicitly.
2. Fully detailed kinetic mechanism with 17 species and 60 fundamental
reactions far H^/alr combustion is used.
3. JANAF Table: Thermochemical properties for each species are
calculated by interpolating the JANAF table.
k. Optical-thin approximation is made, which considers only
gas-to-ambient radiative heat transfer.
One major uncertainty of the model is 0 in Eq. (15), which defines the flamelet
residence time. The other major uncertainty which may effect prediction of NO Is
the optical-thin approximation. It is recognized that the optical-thin approxima-
tion overpredicts the rate of heat transfer, which in turn underpredicts the flame
temperature, and thus underpredicts the NO formation rate.
4-10
-------
^.0 TUNING OF THE MODEL
As a first example and check of the computational procedure, the H^/air flame
measured by Bllger and Beck (Ref. 2) is modeled. The flame is a hydrogen Jet with
0.635 cm diameter and a mean exit velocity of 19323 cm/sec, which is issuing into
o
still, 300 K air at one atmosphere. In order to compare the calculated NO concen-
tration (mass mean of both reactors) to experimental data, it has to be multiplied
by k^. This is because the experimental data are macroscopic and the reactor data
are microscopic.
The computed results, for various g values and a - 1.0, are shown in Figure
4. a is a multiplicative factor for the gas radiative heat transfer based on
the optical-thin approximation. For instance, a= 0.5 means that the heat
transfer rate considered is 50 percent of the optical heat transfer rate. When
comparing the results against data in Figure 1, it can be seen that 0 is on the
order of 10. However, the peak NO concentration point, which generally coincides
with peak flame temperature, is shifting away from measured value as B increases.
This is due to fact that with longer the flamelet residence time heat release
occurs earlier.
If we take 6 ¦ 10, the computed results for various a values are shown in
Figure 5 in comparison to measurements. When a - 0.7 the computed NO concentra-
tions are best fitted to the experimental data.
4-11
-------
D = 6.35 mm
Uj = 19'323^
Re. = 3,067
150
Varying
° 100
l/)
IA
z:
¦s»"S 8 = 10
20
40
60
80
0
100
120
140
(X/D)
Figure 4. Effects of Uncertainty on Flamelet Residual Time 8.
4-12
-------
200
6.35 mm
19,323
sec
Varying
150
E
o.
o.
o
o
^ 100
u
to
s_
>»
V
Li-
en
-------
5.0 REFERENCES
1. Havrthorne, W.R., Weddell, D.S., and Hottel, H.C., Third Symposium of
Combustion. Flame and Explosion Phenomena, p. 266, Williams and Wilkins,
19*»9.
2. Bilger, R.W., and Beck, R.W., "Further Experiments on Turbulent Jet Diffusion
Flames," Fifteenth Symposium (International) on Combustion, p.<*51, The
Combustion Institute, 1975.
3. Ricou, F.P., and Spolding, D.B., J. Fluid Mech., 1_1_, 21, (1961).
i*-. Broadwell, J.E., "A Model of Turbulent Diffusion Flames and Nitric Oxide
Generation. Part I," TRW Document No. 38515-6001-UT-00, EERC Final Report,
PO No. 18889 (1982).
4-14
-------
PART 5
NUMERICAL MODEL
FOR
TWO-PHASE ONE-DIMENSIONAL FLOW REACTOR
Prepared by:
C. J. Kail
and
T. J. Tyson
5-1
-------
TABLE OF CONTENTS
Section Page
1.0 INTRODUCTION 5-1
2.0 FORMULATION 5-2
3.0 CHEMICAL REACTIONS AND THERMOCHEMICAL PROPERTIES 5-7
3.1 Chemical Reaction Rate Equations 5-7
3.2 Thermochemical Data 5-8
<~. 0 HEAT TRANSFER 5-9
*+.1 Gas Radiation 5-9
h.2 Particle Radiative and Convective Heat Transfer 5-9
5.0 WELL-STIRRED REACTOR 5-12
6.0 NUMERICAL METHOD OF SOLUTION 5-1
-------
LIST OF FIGURES
F iqure Poqe
1 Heat and Mass Fluxes Balancing in a Control
Volume (Flow Reactor) 5-3
2 Percentage Fuel-N Conversion for a Premixed
CH^/air Flame at t ¦ 300 ms and 500 ms 5-17
3 HCN and 0 Concentration in a Premixed CH,/air
If'
Flame as a Function of Time 5-18
k Comparisons of Predicted NO Concentration
Against Measurement (Wyoming Subbituminous) 5-21
5 Comparisons of Predicted Exit Oxygen Concentrations
Against Measurements (Wyoming Subbituminous) 5-22
LIST OF TABLES
Table Poqe
1 Definition of A^ and R^ for Eq. (2-8) 5-6
5-111
-------
-------
1.0 INTRODUCTION
A computer model capable of analyzing various types of two-phose
one-dimensional nondiffusive reacting flow6 (e.g., plug flows) is presented in
this report. Detailed formulotions and the numerical method of solution are
described. The model also is applicable to zero-dimensional time-dependent flow
problems and zero-dimensional steady problems with asymptotic solutions (such os
well-stirred reactors).
For illustrative purposes, the detailed analyses of a gas-phase reactor
(i.e., well-stirred reactor followed by a plug-flow reactor) and a caal
well-stirred reactor are presented. Comparisons of calculated results to experi-
mental data are also made.
5-1
-------
2.0 FORMULATION
A control volume Ax " * is shown in Figure 1 through which fluid flows steadi-
ly. The main stream direction is the x-direction. A is the surface area af the
x
reactor (or streamtube) perpendicular ta the flow direction. A^ is assumed to be
varying very gently and only in the x-direction. A^ is the external surface area
af the reactor or streomtube through which externol energy or mass is added or
removed. p is the overall density (i.e., P= Pp+^g). u is the convective veloci-
ty in the x-direction. W. and Iare volumetric chemical production rate of
ig ip
gaseous and condensed species, respectively. 0^ is the volumetric rate of heat
transfer to the particles. is the volumetric heat transfer rate to the gas
that is not associated with external mass flux, qQ is the energy flux to the gas
that is associated with external mass flux. For exomple, qg could be the heating
flux of external mass before the mass enters the control volume. a is the
species concentration af the ith species (either gaseous or condensed species) and
has units of mole per gm of total mass, h is enthalpy. The subscripts "g" and
"p" refer to gaseous phase ond condensed phase, respectively, while the subscript
"e" stands for external flux.
The conservation equations of ma66 and energy can be derived directly from
balancing convective fluxes, as shown in Figure 1; i.e.,
outgoing flux - incoming flux + production.
In the reactors considered here, convective fluxes are assumed to be dominating
and thus diffusive fluxes are neglected.
The resulting conservation equations from flux balancing are shown below:
Conservation of Mass
"(suAx)
dx
Ax ¦ PeveAy ; me (2-1)
Species Conservation Equations
dt"f A. ,
dx " Axflx t'Ji ,e
•l) * (£&) U"2)
5-2
-------
puAx
u
h„
(Condensed-Phase
Irradiation)
(Volumetric Heat Addition Rate)
•- puAx
dx
puAx
(p ~i£Ax)
Ax
q (Heat flux due to heating/
cooling of external mass flux)
Figure 1. Heat and Mass Fluxes Balancing in a Control Volume (Flow Reactor).
-------
Momentum Conservation Equation (Dynamic Equilibrium)
du m , , ijn
QU li OJ (ue - ul ' C1 hT (2-3)
Gas-Phase Energy Conservation Equation
dTn K
~uC pg 3 - -j—^
dT'-JT^ (hg,e ? " Ehi tT9) V £ W "1
<-» fc) (sic) <»>
Solid-Phase Energy Conservation Equation
puCPP^ =I^T (hp,e £Vl.e) +°p (2'5)
Equation of State
PM
0
rng
RTn
or 1n differential form
dc = = dp = dTg . pMg £ dflki (2-6)
y g
Auxiliary Species Conservation Equation
EV, ¦10
p+g
or in differential form
daR3 a:SMidai
H ijn
(.2-7)
where and are conversion factors, Is molecular weight of 1th species and
the subscript "R" refers to a reference species (usually N2). Furthermore, it is
assumed that the condensed and gaseous phases are In dynamic equilibrium (i.e., u
- Up - u) except In well-stirred reaction. However, the condensed and gaseous
phases are not necessarily in thermal equilibrium (T ji Tp^' since heat
transfer mechanism (radiation) for the solid phase in partlclB applications is
very much different from the mechanism for the gaseous phase.
5-4
-------
Depending on the problem at hand, the set of ordinary differential equations
can be solved In several ways;
• If the pressure or pressure gradient and the cross-sectional
area distribution, A^, along the reactor or streamtube are
known. Eqs. (2-2), (2-^f), (2-5), (2-6), and (2-7) ore
integrated simultaneously to solve for , T , T^, and ,
while Eq. (2-3) is integrated individually for u.
• If the cross-sectional area distribution, A^, and the main
convective velocity distribution, u, are given, Eqs. (2-2),
(2-*>), (2-5), (2-6), and (2-7) are integrated simultaneously to
solve for a. , T , T , ond p while Eq. (2-3) is used to solve
i g p
for p.
e If the pressure or pressure gradient and the convective
velocity, u. are specified, Eqs. (2-2), (2-<»). (2-5), (2-6),
and (2-7) are solved simultaneously for a., T , T , ond P ,
i g p
while Eq. (2-1) is integrated for A^.
• If the main convective velocity, u, is set to unity, ond the
pressure or pressure gradient of the reactor is known (i.e.,
eliminating the momentum equotion), the remaining set of
equations can be used to describe time-dependent problems with
the x-direction considered the time domain (t - x/u). Further-
more, this set of equations can be used to solve zero-dimensional
problems with asymptotic solutions via a time-dependent Iterative
method. For example, well-stirred reactor simulation requires
a zero-dimensional asymptotic solution.
In summary, Eqs. (2-2), (2-<0, (2-5), (2-6), and (2-7), which are to be
integrated simultaneously, can be rewritten in the following general form:
C. -j— * A.F. +
-------
Table 1. Definition of Ai, <3^ and R.. for Eqn. (2-8)
F.
Ci
Ai
ri
me
me
a.
l
pu
A Ax
X
AxAx "i,e
Tg
puC
pg
- ^ h .w. + Q
p+g
•
m / >
AxAx ^hg,e 21 hiaie^
. . - dp
+ AV i u dx
T
P
ouC
PP
--
me \
AxAx /hp,e " £hiai,e^
* P
P
1.0
" T ST " P" VdC1'
£ d£
p dx
aR
1.0
"
j-yii,da--
^ d*
--
Note: c] = 0.0242
c3 = 1.01325 x 10+5
5-6
-------
3.0 CHEMICAL REACTIONS AND THERMOCHEMICAL PROPERTIES
3.1 Chemical Reaction Rata Equations
In general form, the Kth chemical reaction can be written as:
ENi,KZi ^ £Ni.KZi (3",)
i i
where N and N" are integers representing the stoichiometric coefficients of
X ) IV 1
the ith species in the Kth reaction. Prime denotes products and no prime denotes
reactonts. is the chemical symbol (or species name) of the 1th species. The
net production rate (in forward direction) of the ith species in the Kth reaction
Is given by:
*i,K - (N!,K " Ni,K>
f ,K J 1 aj
*¦, , rNj.»< n. .
& „ J „ - '.J
i/ P J Ha.
*e.K 1 J
E, (3-2)
where the forward reaction rate has Arrhenius form. The equilibrium constant,
K is calculated by:
e, K *
/-AGK\ -H/N! - N. )
exp I-—-/ (RT) f I T.Kj
e,K \ RT / ' ' ' (3-3)
where AG is net production of Gibbs free energy and is defined as:
K
iGK = ZNi.KGi -ZK1.KS1 ' (3'4)
i i
The term E^ in Eq. (3-2) is the modification of reaction rate due to a third body
and is defined as:
E » a, for three body reactions (3-5)
K K, 1 1
¦ 1 for oil others
where E 19 the third body efficiency of the ith species in the Kth reaction.
K, 1
-------
The overall production rate of the 1th species is the sum of ^ over all K;
i.e.,
Wf = Zw. ,
- 1,J( (3-6)
<
3.2 Thermochemicol Doto
The thermochemical properties (specific heat, enthalpy and Gibbs free energy)
for each species are taken directly from JANAF Thermochemical Table (Ref. 8) in
tabular form as a function of temperature. Linear interpolation is used to
compute the properties at local temperature.
5-8
-------
k.O HEAT TRANSFER
^•1 Gas Radiotion
The major radiating species in a hydrogen flame is water vapor. For ~
hydrocarbon flame, carbon dioxide and carbon monoxide are also among major radiat-
ing species. For the optically this limit ( photon mean free path is much greater
than the characteristic dimension of the flame) the Planck mean absorptance
coefficient of species i can be obtained from gas omittance data (Ref. 1); i.e.,
K - % M
pi * \ L / L-o
where L is the characteristic dimension of the flame and e . is the omittance of
gi
the ith species. The values of K for HgO, CO^, and CO were calculated by Kelly
and Kendall (Ref. 2) based on omittance data of Edwards and Balakrishnan (Ref.3).
The Planck mean absorptance for a mixture can be estimated by:
where P^ is partial pressure of the ith radiating species.
The volumetric radiative heat transfer term, Q , in governing equation can be
K
expressed as
0R - -V
-------
very high. Treatment of interparticle radiation is beyond the scope of this
report.
The radiative energy exchange between a reactor wall and a spherical particle
surface can be expressed by:
p i4-n
where dp is diameter of particle in cm and e is particle omittance.
The convective heat transfer to the particle is:
ipg-Vhc (VTp) (He) <"-2>
where h is convective heat transfer coefficient defined as
c
Nu A(
p sec-cm"
u _ 9 cal
hc " 3 0l
where \ ^ is thermal conductivity of the surrounding gas and Nu is the Nusselt
number. For a rigid spherical particle in completely quiescent surroundings or
when dynamic equilibrium between the particle and the ambient fluid prevails, the
Nusselt number can be shown to equal two. If particle number density is defined
as the number of particles per unit volume, it can be shown that
n = 6p
fi-Y (#/cm3) (4-3)
p W
where c is particle material density and o is mass of condensed phase per
s p
volume. Then, the total particle surface area can be calculated by:
n x TtdJ" = - jP (4-4)
P P Psdp
and the total radiative heat transfer per unit volume can be expressed by:
%,-Vrp "SHC-V) fe) (4"51
5-10
-------
Similarly, the total convective heat transfer per unit volume is
% - VM " -¥t !T9 " Tp) ("rM (4"6>
^ pS p \cm -sec/
The overall particle heat transfer is the sum of the convective and radiative
components; i.e.,
0 = On + 0 (4-7)
"p ^Rp pg
This allows calculation of the last term in Eq. (2-5).
For a coal particle undergoing devolatilization or combustion, the material
density varies with time. If the fact that the coal particle might swell is
ignored and it is assumed that the particle maintains the same size throughout the
process, the instantaneous particle material density for a monosize system can be
estimated by:
ps = np Tidp^ (4-8)
where Y is the total mass fraction of particle phase. Ash is used as a tracer.
P
If we further assume that the ash content in a cool particle stays constant during
the process, the particle number density can be estimated from the ash content of
the system; i.e.,
Instantaneous Ash Content/Volume _ Yash
n _ in-jtial Ash Content/Particle ,3 o „ o
P Tid P Y
K P s ash
(4-9)
Q
where Y and Y . are the instantaneous and initial ash mass fractions in the
ash ash
coal particle, respectively. cg is the initial particle material density of
coal. Thus, from Eq. C»-S) and Eq. (4-9), the instantaneous material density can
be found fro«:
P. = («-10)
ash
5-11
-------
b.O WELL-STIRRED REACTOR
A well-stlrred reactor, also referred to as a perfectly stirred reactor, is
an idealized reactor consisting of a reaction chamber, and inlet for the entry of
premixed reactants and an outlet for the withdrawal of products at such a rate as
to maintain steady overall mass flow through the reactor. Theoretically, a
well-stirred reactor will provide an infinite mixing rote which allows the reac-
tion to be purely chemically controlled without physical mixing limitation. Thus,
the contents of the reactor are chemically and thermodynomically homogeneous
throughout the reactor and are the same as those of the products leaving the
reactor.
Mathematically, simulation of a well-stirred reactor is a zero-dimensional
kinetic problem with a steady solution. The key parameters are reactor residence
time and reactor inlet conditions. As was discussed briefly in previous sections,
if u - 1,0, the set of governing equations described in Section 2 can be applied
for time-dependent problems with the x-dlrectlon considered the time domain. For
example, Eq. (2-2) becomes:
• ¦
da., m w.
dt pA^Ax (ai,e " ai) + p 15 1)
Asymptotically ^i.e. ^jr ¦+ o) , Eqn. [5.1) becomes:
^ (3i " 3i,e) ° "o" ,5_2)
This is the steady state well-stlrred reactor equation, where a is input
species distribution, a is exit species distribution and w is Instantaneous
OAgr At
chemical production rate, t = is defined as the residence time of a
£
well-stirred reactor, since A A x is reactor volume, and m fo is volumetric flow
x e
rote. The residence time is defined as:
_ Volume
Volume Flow Rate '
In summary, we can apply the set of equations described in Section 2 to a
well-stlrred reactor based on the following conditions:
Pa v A V
• Letting " ¦ ^ be the residence time of the well-stirred reactor.
e Eliminating t?q. (2-3) by setting U - 1.0 and dp/dx- 0.0.
5-12
-------
• Integrating time-dependent equations, Eqs. (2-2), (2-<0. (2-5),
(2-6) arid (2-7), until a steady solution is achieved.
Practically, a well-stirred reactor model can be used to simulate a flame holder,
multiburner or recirculation zone in a furnace.
5-13
-------
6.0 NUMERICAL METHOD OF SOLUTION
Since the system involves fast finite-rate chemistry, the governing differen-
tial equations to be integrated are numerically stiff. The implicit finite
difference technique is considered to be the most effective numerical method. The
implicit finite difference scheme for the governing ordinary differential equa-
tion, Eq. (2-8), is expressed as follows:
AX
n+1
- F
n
= A.F.
1 i
n+1
+ 0-n +
n
n+1
-F/)
~ Ri
n (6-1)
where nonlinear terms in Q are linearized via Taylor series expansion and are
truncated after the first order term.
n
Rearranging the above equation gives:
/afiiV
\aa./+ A.
Lx
F n + 1 + F 0+1
1 j
all jV J'
(6-2)
All terms on the right hand side of the equation are either Known or can be
evaluated at xn station. Mathematically, Eq. (6-2) is a set of simultaneous
linear algebraic equations. A Gauss-Jordan reduction algorithm (Ref. 4) with
diagonal pivot strategy is used to solve these equations for all F simultaneously
at xn+^ (- x" + fix) station.
This numerical procedure is basically noniterative. Therefore, the trunca-
tion error due to the linearization of nonlinear terms has to be carefully con-
trolled. The error controlling strategy used is to control the error via the
integration step size, A x, since the finite difference equations asymptotically
approach the ordinary differential equation they represent when 0. The
method used here Is to limit the step size so that no component of solution may
vary by more than a small percentage ( e ) of its value at the last step; i.e.,
ax
= Max
P n+1 p n
i " i
F.
n
(6-3)
If this criterion was violated, the current numerical integration step would be
5-14
-------
repeated, cutting the step size in half. The integration procedure would be
repeated, continuing to halve the integration step size until criterion (8-3) was
satisfied by all solution components. A lower limit ( £„) is also set for a
2 max
If ii was less than eA
-------
7.0 APPLICATIONS
The computer model can operate in the following three modes:
• Plug-flow reactor (gas-phase or two-phase)
• Well-stirred reactor (gas-phase or two-phase)
• Well-stirred reactor followed by plug-flow reactor (gas-phase only)
For illustrative purposes, here we present two sample computations:
• A gas-phase well-stirred reactor followed by an isothermal plug-
flow reactor for gas-phase fuel nitrogen study
• A coal well-stirred reactor simulating the KVB coal stirred reactor
(Ref. 5). Computational results are compared to experimental data.
Details of the computations are described in the following sections.
7.1 Gas-Phase Fuel Nitrogen Conversion Simulation
Conversion of fuel nitrogen in fuel-rich environments has been one of the
most interesting topics in current N0^ control research. As the first example
calculation, the computer model is used to predict fuel nitrogen (HCN) conversion
in a fuel-rich CH^ flame under controlled conditions (i.e., constant temperature
and equivalence ratio). The reactor is an isothermal gas-phase, 10-ms,
well-stirred reactor (as flame holder) followed by an isothermal plug-flow reac-
tor . The methane is doped with 0.5 percent (by volume) HCN. Computations are
made at three different reactor temperature levels and five reactor equivalence
ratios:
T - 1850°K, 1950°K, and 2050°K
$ - 1.3, 1.5, 1.7, 1.9, and 2.1
The calculated results are plotted in Figure 2. Figure 2 shows percentage conver-
sion of HCN to total fixed nitrogen (TFN - [NO + NH^ + HCN]) at two different PFR
exits (300 ms and 500 ms). It is noted that the minima in TFN are observed under
very rich conditions and the minima are enhanced by increasing temperature.
Figure 3 shows HCN and 0 concentrations versus PFR time for two typical
calculations; i.e.,
• o ¦ 1.7 and T ¦ 105O°K
• $=1.9 and T - 1950°K
As can be seen in Figure 3, the HCN concentrations for both cases decay very
5-16
-------
h% HCN in CH./air flame
•SR
1850 K
50
30-
20
10
1.1 1.3 1.5 1.7 1.9 2.1 2.3 2.5
Equivalence Ratio
-------
1000
1/2% of HCN in CH./air flame
100
300
0
100
200
400
500
Plug Flow Time (ms)
Figure 3. HCN and 0 Concentration in a Premixed CH./air
Flame as a Function of Time.
5-18
-------
rapidly in the early stages (near t = 0) and decay linearly at a later stage (t «
300-500 ms).
For the case" <~> - 1.7 and T ¦ 1850°K, reaction rates have been screened. The
screening result shows that the controlling reaction for HCN destruction in the
early stage is the reverse-Zeldovich reaction; i.e.,
N + NO - N2 + 0 (96)
This reaction is responsible for about 67.2 percent of production. The net HCN
destruction rate, from the screening result, at the same stage is approximately an
order of magnitude greater than the production rate. Thus it can be concluded
that reaction (96) is controlling and is second order. In a later stage ( t = 500
ms), the HCN destruction rate is approximately one order of magnitude less than
the Ng production rate. According to the screening result, the reaction:
HCN + 0 ¦+ CH + NO (116)
Is responsible for 66.7 percent of HCN destruction; thus, It can also be concluded
that this reaction is controlling at this stage. Reaction (116) is virtually
first order since the level of O-atom is practically constant at later stages of
plug flow as can be seen in Figure 3.
It is of interest to note that the conclusions derived above are qualitative-
ly in agreement with findings of the experimental investigation recently conducted
in Exxon Research and Engineering Company by Song et al. (Ref. 7).
7.2 Coal Well-Stirred Reactor
The devolatilization of coal can be considered a first-order chemical reac-
tion (Kobayashi, et al., Ref. 7);
dm.
"3T" Ki (m1 * mt,«) l7-"
where m is the mass of ith component (or functional group) af coal and m is the
1 1 | 00
asymptotic yield of the component. The devolatilization rate constant K is found
by the Arrhenius equation (K^ - A^ exp( i/RT). If it is further assumed that the
coal particle temperature is constant during devolatilization in a well-stirred
reactor, the asymptotic yield can be considered constant. Integrating Eq. (7-1)
5-19
-------
gives:
m. — m- . .
(7-2)
i,®
Thus the time rate of devolatilization can be obtained by differentiating the
abave equation:
¦ ' '"l\ = K.e'kit f-U l7"3>
dt \m1' (
(ife)
Due to the lack of experimental information on coal particle exit age distribu-
tion, an analytical exit age distribution is assumed (Ref. 9), For a reactor
having mean residence time TSR> the °9e distribution is:
E (t) = ~ e"t/TSR
SR
E(t) dt is the mass fraction of coal at reactor exit having age between t and
t+dt. Using this analytical age distribution expression, we can calculate the
rate of devolatilization of the ith functional group as:
"i 7~ e"t/TsR * M'V dt
LSR
= Ki / M , ,
KiTSR + ^ Uec/ (7-4)
The coal stirred reactor model is then used to simulate KVB's experiment (Ref. 5).
KVB's well-stirred reactor fired pulverized Wyoming subbituminous coal. The mean
particle size was approximately 50 pa and the reactor temperature reported was
1511°K. The mean reactor residence time varied from 10 ms to 50 ms. The reac-
tor's exit NO and oxygen concentrations were measured. The computed results of
the coal reactor model are plotted against KVB's measurements in Figures k and 5.
Figure
-------
2000r
Measure-
ment
Prediction
2.63
1000
800
600
400
200
10
30
0
20
40
50
Mearv Reactor Residence Time, msec
Figure 4. Comparisons of Predicted NO Concentration Against
Measurement (Wyoming Subbituminous).
5-21
-------
\
o
\
Measure-
ment
A
O
o
$ Prediction
2.63
1.25 .
\
^ o
1.05
X
O
x
10 20 30
Reactor Mean Residence Time, msec
40
50
Figure 5. Comparisons of Predicted Exit Oxygen Concentration
Against Measurements (Wyoming Subhituminous).
5-22
-------
this report. The fundamental reaction rates for char oxidation and char gasifica-
tion (for 50 ub particles) used in these computations are:
Oxidation
C(s) + 02 -
C02« + C(s)
Gasification
C(s) + h2o
C(s) + CO-
CO.
CO + CO
CO + h2
CO + CO
K • 1.08 x 10
Kf - 1.0 x 1016
10
exp (-20.000/RT)
10
11
exp (-30,000/RT)
Kf = 9.01 x 10
Kf - 6.44 x 10"' exp ( -35,000/RT)
5-23
-------
8.0 REFERENCES
1. Sparrow, E. M. and Cess, R. D., Radiation Heat Transfer. Brook/Cole
Publishing Company, Belmont, California, p. 218 (1970).
2. Kelly, J. T. and Kendall, R. M., "Further Development of the Premixed
One-dimensional Flame (PROF) Code," EPA-600/7-78-172a (NTIS PB2862<»3)
(1978).
3. Edwards, D. K. and Balakrishnan, A. , "Thermal Radiation by Combustion
Gases," Int. J. Heat and Mass Transfer, Vol. 16, p. 221 (1973).
k. McCracken, D. W. and Dorn, W, S., "Numerical Methods and Fortran
Programming," Wiley, New York (l96
-------
APPENDIX
Program Input Data
The input data cards for One-Dimensional Flame Analysis Program are explained
in this appendix. The input cards can be divided into the following groups:
I Title Card
II Namelist Cards (SINPUT)
The namelist cards are divided into the following sub-groups:
• General Input
• Well-Stirred Reactor Option
• Solid Phase Option
• Option to Specify Temperature Distribution
• Option to Specify Pressure Distribution
• Option to Specify Velocity Distribution
• Option to Specify area Distribution (Confined Flame)
• Screen Option
• Optional Table Input for External Flux
• Optional Automatic Air-Staging
• Option to Specify Oxidant/Fuel Ratio for external Flux
• Sensitivity Analysis
• Output Control
• Numerical Control
III Species Cards
IV Reaction Cards
V Namelist Cards (5C0ALRK)
The program is coded for CDC 7600 computer using FORTRAN IV language.
5-25
-------
INPUT FOR OFAP
The following are the descriptions of input data for the One-Dimensional
Flame Analysis Program.
I.
TITLE CARD
Columns 1-70 are used for run ID
Columns 71-75 are used for continuation index which is a right justified
integer indicating the number of files in the continuation tape.
II. NAMELIST (SINPUT)
The following parameters are general input.
Symbols
Descriptions
Default
Value
KASE
NS
NR
ICRU
ITIME
X
XMAX
DXMAX
DXMIN
DX
E3RDB(I,K)
T
U
ALPHA(I)
Case No.
Total number of species (Max. 52)
Total number of reactions (Max. 200)
Index for rate constant units
- 0 molecule-dm-sec units
- 1 mole-cm-sec units
Maximum CPU time allowed for the run (sec)
Current streanrwise distance
Maximum X to be integrated to
Maximum streamwise integration step size allowed
Minimum streamwise integration step size allowed
Initial streamwise step size
Third body efficiency
I - species index
K - reaction index
Initial temperature (°K)
Initial axial velocity (cm/sec)
Initial species distribution (mole fraction)
0
0
0
1
DXMIN
1 . 0
Default
5-26
-------
Symbols
Descriptions
Value
TKINET
P
ITHEQM
IDRADG
IDRADP
TWALL
AREA
TAUSR
HSTATE
ALPHAF(I)
TE
VE
Kinetic cut of temperature ( K)
-------
Symbols
Descriptions
Value
YSOLDE Solid mass fraction in external flux
ICOLSR Coal S.R. Index
j* 0 if run is coal S.R.
PHE S.R. equivalence ratio (for coal S.R. only)
ICONPF P.F.R. continued from S.R.
PHEF P.R.F. equivalent ratio
IDCOAL Coal ID
» 0 not coal particle
¦ 1 bituminous
- 2 lignite
Option to Specify Temperature Distribution
0.0
0
0
0
0
0
If the temperature of reactor is to be specified, the energy equation is
by-passed. The following parameters for the T(X) table should be provided in
$INPUT namelist.
Default
Symbols
Descriptions
Value
NXTRTB No. of entries for T(X)
Table (Maximum 50)
XTRTB(I) X Coordinate of T(X) table
TRTB(I) Temperature (°K) at XTRTB(I)
5-28
-------
Option to Specify Pressure Distribution
If the pressure of reactor is to be specified, P(X) table has to be
provided (under $INPUT) as follows:
Default
Symbols Descriptions Value
NXPRTB No. of entries for P(X) table (maximum 50) 0
XPRTB(I) X Coordinate of P(X) table
PRPTB(2) Pressure (ATM.) at XPRTB(2)
Note:
i) dp/dx will be calculated by simple differencing
ii) If U(x) is constant, reactor area AREA will be
calculated
ill) If reactor area is defined (i.e., NXRWTB 0) U(x)
will be calculated.
Option to Specify Velocity
If reactor velocity is to be specified, U(x) table Mas to be provided
(under (INPUT) as follows:
Default
Symbols Descriptions Value
NXURTB No. of entries for U(x) table (maximum 50) 0
XURTB(I) X coordinate of U(x) table
URTB(I) Velocity (cm/sec) at XURTB(I)
5-29
-------
CONFINED FLAME (VARIABLE AREA)
For confined flame with variable area, stream function coordinate system
must be used (i.e.", IPSICD » 1). The following parameters defined the area:
Default
Symbols Descriptions Value
NSRWTB Number of entries for area specifying table 0
(maximum 50)
XRWTB(I) x coordinate of area table (cm)
RWALTB(I) r coordinate of the wall(cm)
SCREEN
Screen option gives the reaction number and its percentage contribution
(integrated over r at constant x) of production/destruction, of up to five
most active reactions, for each specified species.
Default
Symbols Descriptions Value
NSPSCN Number of species to be screened (maximum 10) 0
ISPSCN(I) Index of the species to be screened 0
TABLE INPUT FOR EXTERNAL FLUX
Default
Symbols Descriptions Value
NXETB
XETB(I)
XMDETB(I)
AETB(I)
TETB(I)
QDETB(I)
Number of entries in external flux table
(maximum 30)
Axial coordinate for external flux table
Mass flux table (gm/sec-cm)
Species mole fraction table
Temperature table (°K)
Heat flux tube (cal/gm)
5-30
-------
OPTIONAL AUTOMATIC AIR-STAGING
Symbols
Descriptions
Default
Value
IAUSTG
XAUSTG
TAUSTG
Automatic air-staging index
? 0 if automatic air-staging (cm)
Duration of automatic olr-staglng (cm)
o
Staged air temperature ( K)
Note: If IAUSTG is not 0 program will set
NXETB - 5 and automatically calculate
external flux tables (e.g.. XETB, SMDETB,
AETB and TETB) based on abave information
and PHE/PHEF
SPECIFYING OXIDANT/FUEL RATIO FOR INFLUX
Default
Symbols Descriptions Value
IOFRSP o/F specification index
f If o/F ratio is specified for influx
ALPHAA(I) Oxidant species distribution (same unit as ALPHA)
ALPHAF(2) Fuel species distribution (same unit as ALPHA)
Symbols Description
Default
Value
YSOLIA
YSOLIF
OFRE
Solid mass fraction in oxidant
Solid mass fraction in fuel
o/F for influx (Mass)
0.0
0.0
0.0
5-31
-------
OUTPUT CONTROL
The following parameters are output control parameters.
Symbols
Description
Default
Value
NPRT
XPRINT(I)
I0UT1
I0UT2
I0UT3
ZPUNCH
IMFCHK
IREWIN
Number of print stations in X direction (max. 30)
X coordinate for output (cm)
If y 0 print net production rates for each species
If + 0 print
If j1 0 print forward/backward rates of each reaction
If )l 0 store the necessary variables in tape for
continuation
Will print mass flux (integrated over r at fix x)
if f 0
Will rewind continuation tape before written if jl 0
NUMERICAL CONTROL
The following variables are designed to gain some internal control of the
accuracy of the numerical calculations.
Default
Symbols
Description
Value
EPSIL
EPSI
Maximum tolerable negative value (EPSIL) for
species concentration which will be set to a
very small positive value (EPSI)
-1 .0 E-9
+1.0 E-16
Default
5-32
-------
Symbols Descriptions Value
EPSTOP Convergent criteria (°K/ms) for time-dependent 1.0
calculation
EPSIPC Minimum solid mass fraction to be considered 1.0 E-6
DELLO Lower limit for step size control 0.02
DELUP Upper limit for step size control 0.05
MINNER Maximum number of inner Iterations allowed 1
NC0UPL Number of coupling eqns. NS+2
SPECIES CARDS
The following group of cards are the species name, molecular weight, heat
of formation and parameters for laminar diffusion coefficients.
Card
No. Column
3.1 1-5
6-10
11-20
21-30
31-40
M-50
51-60
61-70
71-80
3.2 1-5
REACTION CARPS
The chemlcol reaction mechanism for a particular problem is input on the
lost set of cards, one card for each reaction. No particular order is re-
quired. The forward reaction rate constant is expressed as AT N exp(-B/RT).
Description Format
Name of the first species A5
Blank 5X
Molecular weight E10.3
Heat of formation at 298,15°K (kcal/mole) E10.3
Simplified diffusion parameter (bifurcation) E10.3
Stockmayer force constant £0/K (°K) E10.3
Stockmayer force constant a(A) E10.3
Stockmayer for constant 6 (dimensionless) El0.3
Polarizability a (10^ cm^) E10.3
Name of second species A5
5-33
-------
Ten possible reaction types included in the program are:
Reaction Type
(1)
A
+ b + c
(2)
A
+ B + M t;
(3)
A
+ B X I
(<0
A
+ B
(5)
A
+ M % c
(6)
A
+ B - C
(7)
A
+ B + M -*
(8)
A
+ B -+C
(9)
A
+ B - C
(10)
A
+ M -+ C
(11)
S« + B -~ i
(12)
A
s«
+ Q -> <
(13)
«
C
1
+ B I
(1<0
A
C«
+ B ¦+ 1
C * M
C + D + E
C + M
C + D Surface reaction (soot)
C + D Surface catalytic reaction (soot)
C + D Surface reaction (char)
D + D Surface catalytic reaction (char)
Card
No.
Column
Description
Format
4.1 1-2 Stoichiometric coefficient of Species A 12
3-7 Species A A5
8 + sign
9-10 Stoichiometric coefficient of Species B 12
11-15 Species B (or M) A5
16 -i- sign
17-18 Stoichiometric coefficient of Species M (or 0) 12
19-23 Blank or M A5
2k - sign
25-26 Stoichiometric coefficient of Species C 12
27-31 Species C A5
32 + sign (if needed)
33-3
-------
Card
No.
Column
Description
Format
k. 1
M-42
Stoichiometric coefficient of Species
E
12
kZ-H7
Species E (or M/blank)
A5
«-50
Reaction type, 1 to H
12
51
Blank
52-59
A, pre-exponential factor cm-mole-sec
units
E8.2
60-63
N, temperature exponent
Fk. 1
6k-72
B, activation energy kcal/mole
F9.1
73-80
Comments
A0
k.2
Next reaction
5-35
-------
If IDCOAL t 0, the namelist of $COALRK, which defined cool properties,
has to be provided.
SCOALRK
YC
YH
YO
YN
YS
YMOIST
YASH
Mass fraction of coal composition (C, H, 0, N, S, moisture
and ash) from ultimate analysis
HHVWIF High heating value of coal (ash and mineral free) (cal/gm)
AEO(I) ) Parameters for equilibrium constants which define asymptotic
ANEQ(I)f yield of devolatilization
EEQ(0]
(
)) '
ACD(I)
ECD(I)
ATAR
ETAR
AOLE
AACE
ASSO
ACH(f
ECHff
Arrhenius parameters of devolatilization rate for Ith
species, tar, olefins, acetylene, soot and methane
YLGTVO(I) Mass fraction of light volatiles in coal
SEND
5-36
-------
PART 6
GENERATION OF ELLIPTIC-CODE TEST CASES
Prepared by:
C. J. Kau
and
T. J. Tyson
6-1
-------
TABLE OF CONTENTS
Section Page
1.0 INTRODUCTION 6-1
2.0 OPPOSED-JET DIFFUSION FLAME 6-3
2.1 Coordinate Transformation 6-8
2.2 Limitation of the Angle of Rotation - a 6-10
2.3 Test Cases No. 1 and No. 2 6-12
3.0 COFLOWING-JET DIFFUSION FLAME 6-13
3.1 Test Cose No. 3 6-13
k.O REFERENCES 6-10
APPENDIX 6-19
LIST OF FIGURES
Figure Page
1 Schematic of Opposed-Jet Diffusion Flame 6-k
2 Velocity of CO/O^ Opposed-Jet Diffusion 6-5
3 Temperature Profile of Opposed-Jet Diffusion Flame . . 6-6
it Species Distribution of CO/O^ Opposed-Jet Diffusion Flame . . 6-7
5 Schematic of Coordinate Transformation of Opposed-Jet
Diffusion Flame 6-9
6 Limitation of the Angle of Rotation B-11
7 Planar Co-Flowing Diffusion Flame of CO/O^ 6-1
-------
1.0 INTRODUCTION
The EPA's Fundamental Combustion Research (FCR) Program has undertaken a
project to assess "current capability to numerically solve the Navier-Stokes
elliptic field equations, which govern a wide class of combustion behavior. The
intent is to establish whether computational capability exists which could be used
to develop low pollutant emission combustion systems. Numerical techniques which
appear promising will be candidates for future development through FCR sponsor-
ship.
The assessment focuses on the numerical capability of solving the differen-
tial equations, and not on the analysis of the physics and chemistry imbedded
within the equations. To assess this capability, the values of all physical and
chemical variables have been prescribed. The computational results will be
compared with the solutions used to generate the test cases.
This section will describe how the first group of test cases are generated.
The characteristics of the first group are:
• Non-Premixed Reactants
• Steady State
• One Atmosphere Nominal Pressure
• Simple Slow Chemistry
• Simple Flow Field Behavior
• Simple Diffusive Transport Model (Equivalent molecular diffusion
with large "effective" coefficients).
• Boundary Conditions on Rectangular Coordinate Lines
The first group of test cases consists of 3 cases which are generated from
the fallowing two simple parabolic types of flow fields:
• Planar Opposed Jets of Fuel and Oxidizer
e Planar Co-Flowing Jets of Fuel and Oxidizer
For slow chemistry, carbon monoxide and pure oxygen are selected as the fuel
and oxidizer. A simple three step reaction mechanism which includes only two
carbon monoxide oxidation reactions and one oxygen dissociation step is assigned,
6-1
-------
i.e.
Reaction Rote (mole/cc~sec system)
>2 * C°2 + ° kf1
(1) CO + 0 + CO +• 0 k » 6.905 x 107 T exp(-34810/RT)
(2) CO + 0 +' M ; C02 + M kf2 - 3.8 x 102if T~3 exp(-6170/RT)
(3) 0 + 0 + M £ 02 + M k - 1.4 x 1O10 T_1 exp(-340/RT)
Rate constants for the above reactions are taken from EER's basic reaction set.
Thermochemical properties for all four species (C0,02, C02 and 0) are taken
directly from JANAF Tables (Ref. 1). Turbulent eddy diffusivity and turbulent
Prandtl and Lewis numbers are assumed to be constant. EER's General Flame Analy-
sis Program (GFAP) (Ref. 2) is used for preliminary computations. The results
obtained from these computations are then used to extract essential flow field
parameters such as flame thickness, impingement plane displacement, etc. These
parameters are then used to help define the boundary conditions analytically for
the test cases.
Obtoining the exact Navier-Stokes solution for the opposed jet case is
possible because the solution behavior is self-similar, leading to ordinary
differential equations. After solving the equations, the solution is then skewed
in the coordinate system to destroy the self-similarity, and diffusive fluxes are
introduced in both coordinate directions. A problem posed in this manner is a
true test of the computational capability of a Navier-Stokes code.
The second coflowing case is close to being exact. Its accuracy rests on the
validity of the Prandtl boundary layer approximation for coflowing, slowly-varying
mixing layers. Having generated a solution under this valid assumption, the
solution is once again skewed in the coordinate system to provide a flame field
where diffusive fluxes are of the same order in both coordinate directions.
Although a problem posed in this fashion requires the full Navier-Stokes equations
for solution, it is a considerably Ibbs severe test than the opposed Jet case
because of the lack of a strong upstream influence. With this type of problem it
is felt that a calibration of computational capability can be obtained.
The details of this test-case generating process are described in the follow-
ing sections.
6-2
-------
2.0 OPPOSED-JET DIFFUSION FLAME
It is known that there are similitude solutions for opposed-jet flames. The
self-similar solutions for temperature and species concentrations involve the
similarity variable n _where a is the strain rate and \j is kinematic viscos-
ity. Figure 1 is a schematic of an opposed-Jet diffusion flame. Due to the
density change caused by heat release, the impingement plane of the oxidizer jet
shifts away from the plane of the fuel jet by a displacement 5 The
Impingement plane is an imaginary plane where the V component of potential flow
goes to zero (by extrapolation).
The test case selected consists of a CO/O^ flame with a fuel-side strain rate
0^ - 30(1/sec). Prandtl and Lewis numbers both equal 0.9. The equivalent molecu-
lar viscosity, y , equals 0.01 (gm/cm-sec). The similitude solution obtained
from GFAP computations is plotted in Figures 2 through Figure 2 shows the
velocity profile of V versus r. The displacement 6 , which is measured directly
from Figure 2, is 2.29 cm. Figure 3 shows the temperature profile of the same
flame. The flame thickness, tf, measures the distance between points where
o
temperature deviates by approximately 1 K from the temperature of the unreocted
jet. In this case, tf » 6.5 cm. Figure 4 presents the species concentration
distributions for 0^, CO, C0^ and 0.
The potential flow field on the fuel side is defined by: . .
(Z-1 a}
V » -a y
u - v (2-">)
Thus, the displacement <5 will help define the potential flow field on the oxygen
side,
(2-2a)
V - -ao(y-5 )
U - a X (2-2b)
0
where a is the oxidizer side strain rate and is related to the fuel-side strain
o
rate by
fit
" if V p0 (2-3>
5 p
since there is no pressure drop cross the flame (i.e. — -0 ).
Due to the similarity of the solutions, the conditions
ST <7> = 0
£ = 0 6_3 (2-4b)
-------
y
V = -aly - 5)
Oxidizer Side
*T"*VH|'
Imaginary Oxidizer Jet Impingement
Plane
Imaginary Fuel Jet Impingement Plane
Fuel Side
V = -afy
U = a^x
Figure 1. Schematic of Opposed-Jet Diffusion Flame.
6-4
-------
120
100
80
60
40
20
0
-20
-40
-60
-80
•100
— 5=2.29 cm
u = 0.01 (grn/cm-sec)
-Cn
P = = 0.9
= 0.9
Op-side
J I I I I
Figure 2. V Velocity of CO/O2 Opposed Jet Diffusion Flame.
6-5
-------
3000
0.01 (^m/cm-sec)
2500
sec
° 2000
« 1500
1000
500
0L_
3.0
r (cm)
Figure 3. Temperature Profile of CO/Og Opposed-Jet Diffusion
Flame.
6-6
-------
c
o
(J
ro
S_
O)
1.0
0.9
0.8
0.7
0.6
0.5
r 0.4
0.3
0.2
0.1
0.01 (9 /cm-sec)
0.9
0.9
30
-3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 5.0
r (cm)
Flqure 4. Species Distribution of C0/0£ Opposed-Jet
Diffusion Flame.
6-7
-------
(2-4c)
(2-4d)
are satisfied for the entire flow field.
2 .1 Coordinate Transformation
Two elliptic test cases are generated by coordinate transformation of the
opposed-Jet diffusion flame flow field in the following manner (see Figure 5):
(a) Rotating the coordinate axes x-y by a positive (counterclockwise) angle,
to the axes x'-y' (see Figure 5a). The transformation can be expressed mathe-
matically by the equations:
/ x » x' cos a + y' sin a
\ y - -x' sin a + y' cos a
and
/ U » U1 cos a + V sin a
\ V « -U' sin a + V cos a
(b) Carving out a rectangular domain (x' , y' ) to define the
m m
boundaries of the flow field (see Figure 5b). y is taken to be 2y.
ffl I
(i.e. symmetric with respect to y'-axis).
(c) Translating the coordinate axes x'-y' to the axes x"-y" which have
origin at the left-lower corner of the rectangular domain, the axis yn on left
boundary and axis x" on bottom boundary (see Figure 5c), or mathematically:
x' - x" -(x -x.)
m 1
v' " v" "(vm " v1)
and
U" - U
V - v
6-8
-------
Figure 5.
Schematic of Coordinate Transformation of an Opposed-Jet Diffusion Flame.
-------
Boundary conditions for the rectangular domain ore transformed according to the
above equations.
2.2 Limitation of the Angle of Rotation - a
As is shown in Figure 6, inflowing conditions ore physically meaningful on
boundaries AO and BC (i.e. all independent variables are specified), while
outflowing conditions are meaningful on boundaries AB and DC (i.e. only deriva-
tives of the variables are specified). In order to maintain outflow conditions on
boundaries AB and DC, two criteria that limit the angle of rotation (i.e. a )
must be satisfied:
• There is net outflow at points A and C (i.e. S^^_a )• 02 tf19
angle between net flow direction and AB.
• The flame zone does not appear on boundories AB and BC.
Geometrically (see Figure 6), the distance between point D and the intersection of
DC and the x-axis is:
5 - y - x tan a
or
S
_ = tan a = cot(o + 0,) - tan ct (2-5)
*1 X1 1
Assuming 6 is small in comparison with the dimension of the rectangular domain,
02 is approximately equal to s1 since
tan Ql = *
and
tan 0, - U
'2 -V. y-S
In the limiting situation in which velocity at point A is tangential to boundary
AB:
02 — G-| = ot
and the flame is barely touching the upper boundary X6.
— = cot (2a) - tan a a 0
X1
6-10
-------
Figure 6. Limitation of the Angle of Rotation a.
6-11
-------
Solving for gives:
a » 30°
which is the maximum angle of rotation. The same conclusion con also be drown by
analyzing the flow direction at point C.
2.3 Test Coses No. 1 and No. 2
The boundary conditions for test cases No. 1 and No. 2 are derived from
analysis of the same opposed-jet diffusion flame which is rotated by an angle
o
a ¦ 25 . The only difference between case No. 1 and case No. 2 is the size of
the rectangular domain. The parameters for these two cases are specified in Table
1.
The physical flame thickness is obviously the same for both cases No. 1 and
2. However, relative to the y computational domain, case No. 1 represents a thin
flame and case No. 2 is a thicker flame. Numerically, a thick flame is easier to
analyze since more of the available grid points may be positioned within the
region of rapid change. Boundary conditions for these two cases are listed in the
Appendix at the end of Part 6.
Table 1.
^arametars far Cases No. 1 and No. 2
Case No. 1 2
main (*m x yj
52 x 38
30 x 22.5
y, (cm)
19.794603
11.996887
x1 (cm)
26.0
15.0
a (deg.)
25^3
25.0
«f lik1
30.00
30.00
•o life)
28.0685
28.0635
Tf (°k;
300.00
300.00
Ta (°K)
300.00
300.CO
fuel type
CO
CO
oxidizer type
°2
u( ^/en-sec)
0.31
0.01
Pr
0.9
0.9
le
0.9
0.9
6-12
-------
3.0 C0FL0WING-JET DIFFUSION FLAME
The coflowing planar jet-boundary diffusion flome calculation is performed
using GFAP. As shown in Figure 7, the upper flow is from an oxidizer jet which
contains pure at 300°K. The lower flow is from a fuel jet which contains pure
CO at 300°K. The mainstream velocity of oxidizer Jet is 1500 cm/sec and the fuel
jet velocity is 750 cm/sec. Numerically, the cross-stream velocity on the fuel
side is forced to zero at yiO in the coordinates shown in Figure 7, while the
cross-stream velocity on the oxidizer side is left to be calculated. Figure 8
shows the edge cross-stream velocity on the upper flame boundary. Prandtl and
Lewis numbers are assumed to be 0.9. Figure 7 also presents the upper and lower
flame zone boundaries and the locus of peak temperature in the flame zone.
3.1 Test Case No. 3
To produce elliptic test case No. 3 the following coordinate transformations
are made:
(a) Translating the origin of the coordinate system x-y to x-12.433, y-
13.0, which is an arbitrarily selected point on the locus of peok
temperature. The new axes are denoted by x' - y'.
(b) Rotating the coordinate axes x'-y' by a positive angle
(counterclockwise) a - 25°, to the new axes x"-y".
(c) Translating the origin of the coordinate axes x"-y" downward to x" =
0, y" - -15.0. The new coordinate axes are represented by x"'-y"'.
Figure 9 shows the flame after the transformation. Boundary conditions on
all four sides of the domain (250 cm x 180 cm) are obtained in the following
manner:
• The bottom (fuel-side) boundary condition can be determined
analytically by
U"1 - U cos a
V"' - U sin a
since V - 0 is assumed on the fuel-6lde.
e The upper (oxidizer side) boundary condition can be obtained by
U"' - U cos a - V sin a
V"1 » U sin a + V sin a
6-13
-------
Oxidizer Side
Out-Going B.C.
100
U =1500 cm/sec
50
Upper Flame Boundary
Peak Flame Temperature Locus
0
Fuel Side
Lower Flame Boundary
co
-50
250
150
200
100
50
x(cm)
Figure 7. Planar Co-Flowing Jet Diffusion Flame of CO/O^
(in original coordinates).
6-14
-------
300
200
100
0 50 100 150 200 250
x(cm)
Figure 8. Edge Velocity on Upper Flame Boundary (in original
coordinates).
6-15
-------
150
100
E
u
>>
250
200
100
150
x1" (cm)
Figure 9. Planar Co-flowing Jet Diffusion Flame of CO/COg
(in transformed coordinates).
6-16
-------
where V is obtained by curve fitting the numerical result (Figure 8).
• The left-hand side (x"' =* 0) boundary condition can be divided into
three portions:
1. The region from the coordinate origin to the lower flame
boundary (y"' - 10 cm): The boundary condition in this portion
is the same as that of the bottom boundary.
2. Flame zone (h"' « 10 cm to 20 cm): The boundary condition is
provided by the fine plots from the numerical results of GFAP
code.
3. The portion above the flame zone (y"' >20 cm): The boundary
condition in this portion is the same as that of the upper
boundary.
• The right-hand side (x"1 - 250) boundary condition is the
transformed outgoing boundary condition.
The complete details of the boundary conditions can be found in the Appendix.
6-17
-------
4.0 REFERENCES
1. Stull, D. R., and Prophet, H., "JANAF Thermochemical Tables - Second
Edition," NSRDS-NBS37, Dow Chemical Co. (1979).
2. Kau, C. J., and Tyson, T. J., "Fundamental Combustion Research Applied to
Pollution Control, Vol IV: Engineering Analysis, Part 1: A Computer
Program for General Flame Analysis," EERC Final Report, Contract No.
68-02-2631 (1907).
6-18
-------
APPENDIX
The boundary conditions for all three coses of the first group of
elliptic-code test cases are presented in this appendix. For simplicity, the
superscripts (i.e. " ' " " ", and " ") which were used to differentiate
variables in different coordinate systems during coordinate transformation ore all
dropped in this section.
6-19
-------
TEST CASE NO. 1
1. Geometry: Planar with boundary conditions on rectangular
coordinate lines
2.
3.
Species Present: CO, CO^, 0^, 0 (All thermodynamic properties from
JANNAF tables)
Reaction Mechanism:
Reaction
CO + o2 z co2 + 0
CO + 0 + M Z C02 + M
0 + 0 + H : 02 + M
Rate (mole-cc-sec system)
kf = 6.905 x 107 T exp (-34810/RT)
kf = 3.8 x 1024 T"3 exp (-6170/RT)
k_ = 1.4 x 1018 T"1 exp (-340/RT)
4. Transport Properties:
Constant Equivalent Molecular Viscosity: u = 0.01 gm/cm-sec
Constant Prandtl No.: Pr = 0.9
Constant Lewis No.: Le = 0.9
5. Pressure: P = 1.0 atjn at (x = 0, y - 0)
6. Boundary Conditions:
(I) at y = 0, (0 1 x j 52 cm)
T = 300°K
|xCo - 1.0
X. = 0 if i f CO
U = 19.2836x - 919.75789 (cm/sec)
V = 22.98133x - 246.44847 (Cm/sec)
(II) at y = 38 cm, (Ol x 1 52 cm)
T = 300°K
Xn =1.0
°2
X. = 0 if i f 02
U = 18.04084x - 70.64078 (CT/sec)
V = 21.501718x - 857.92S9 (^/sec)
6-20
-------
(Ill) at x = Q, (Q £ y £ 38 cm)
ax. ax.
5—!- + 0.4663 = 0
3 x 3 y
3U
3x
+ 0.4663 + (J32til55_ y)u + (73.9625 - y)
U ~ 0.4663 ~ (73.96k —)u ^(73,09626563. y)>
(IV) at x = 52 cm. (0 1 y 1 38 cm)
3T . ri acc? n
tt— * 0.4663 -r— = 0
3x 3y
3X • 3X¦
1 + 0.4663 1 - "
3 x - 3 y
3U
3x
^•"«f-(^^7)u-(305|rTT)v
3x + °'4663 f
Steady State
8V + n *1 ( 1 V. ( 0.4663 \
3y " 137.5517 + y/ " \37.5517 + y)
Solve for Ignited Solution
6-21
-------
TEST CASE NO. 2
Geometry: Planar with boundary conditions on rectangular
coordinate lines
Species Present: CO, CO^, 0^, 0 (All thermodynamic properties from
JANNAF tables'
Reaction Mechanism:
Reaction
CO + o2 t co2 + 0
CO + 0 + M Z C02 + M
0 + 0 + M J 02 + M
Rate (mole-cc-sec system)
kf = 6.905 x 107 T exp (-34810/RT)
kf = 3.8 x 1024 T"3 exp (-6170/RT)
k, = 1.4 x 1018 T"1 exp (-340/RT)
Transport Properties:
Constant Equivalent Molecular Viscosity:
Constant Prandtl No.: Pr = 0.9
Constant Lewis No.: Le = 0.9
u = 0.01 gm/cm-sec
Pressure: P = 1.0 atm at (x = 0, y = 0)
Boundary Conditions:
(I) at y s 0, (0 1 x 1 30 cm)
T = 300°K
xco= 1-0
X. = 0 if i f CO
U - 19.2836x - 530.6299 (cm/sec)
V = 22.98133x - 142.1818 (^/sec)
(II) at y = 22.5 cm, (0 - x - 30 cm)
T = 300°K
Xn = 1.0
°2
Xi = 0 if i f 02
U = 18.04084x - 39.8421 {""/sec)
V = 21.5017l8x - 480.6995 (cm/sec)
6-22
-------
(Ill) at x = Q (0 1 y 1 22.5 cm)
§7 * t7 * 0
3X • 3X.
_L+ 0.4663 ^=0
3U + 0.4663 3U U °-4653 V
3x u-HOOJ 9y ' (0.4663y - 19.8976) " (0.4663y - 19.8976)
9V A n „cc, 3V 0.4663 U 0.2174 V
3x u-'fODJ 3y " (0.4663y - 19.8976) " (0.4663y - 19.8976)
(IV) at x = 30 cm (0 1 y i 22.5 cm)
0.4663 |I-0
3X. ax.
+ 0.4663
3x 3y
3U n arm y 0.4663 V
3x 3y " (0.4663y + 10.1023) " (0.4663y + 10.1023)
3V , „ 3V 0.4663 U 0.2174 V
31 + °'4663 37 " (0.4663y + 10.1023) " (o.4663y + 10.1023}
Steady State
Solve for Ignited Solution
6-23
-------
TEST CASE NO. 3
1. Geometry: Planar with boundary conditions on rectangular
coordinate lines
2. Species Present: CO, C02> C^. 0 (All thermodynamic properties from
JANNAF tables)
3. Reaction Mechanism:
Reaction Rate (mole-cc-sec system)
CO + 02 J C02 +0 kf = 6.905 x 107 T exp (-34810/RT)
CO + 0 + M J CO, + M L = 3.8 x 1024 T"3 exp (-6170/RT)
1 ft -1
0 + 0-»-M^02 + M kf = 1.4 x lO10 T 1 exp (-340/RT)
4. Transport Properties:
Constant Equivalent Molecular Viscosity: y = 0.01 gm/cm-sec
Constant Prandtl No.: Pr = 0.9
Constant Lewis No.: Le = 0.9
5. Pressure: P = 1.0 atm at (x = 0, y = 0)
6. Boundary Conditions:
(I) y 0. (Olxl 250 cm)
U = 679.73 (^/sec)
V = 316.90 (^/sec)
T = 300°K
XC0 = 1'°
Xi = 0 if i f CO
(II) y = 180 cm (0 < x 1 250 cm)
U - 1323.09 + 0.1472x - 2.5301 x 10'4 x2 (""/sec)
V = 711.72 - 0.3126x + 5.3508 x 10"4 x2 (^/sec)
T = 300°K
X02 = l'°
Xi = 0 if 1 * 02
6-24
-------
(III) x = 250 cm (Oi y i 180 cm)
fl + 0.4964 fl . 0
ax ax
—¦ + 0.4964 -r-J- = 0 for all i
o x d y
— + 0 4964 — - -7 r - ^ r = 0
3x 3y (x + 0.4964y) (2.0145* + y) u
M * n At** *1 {0.4964)U (0.4964)V _ n
3x °-4964 a> " (x + 0.4964y) " (2.0l45x + y) "
(IV) x = 0
(i) 0 — y _ 10 cm
U = 679.73
V = 316.90
T = 300°K
XC0=
X. - 0.0 if 1 f CO
(ii) 10 cm 1 y 5. 20 cm
See Figures A1-A6 for U, V, T, XCQ, Xg , XQ and XCQ
2 2
(i i i) 20 cm 1 y 5 180 cm
U = 1254.39 + 0.7622y - 2.114 x TO"3 y2
V = 859.38 - 1.642y + 4.565 x 10~3 y2
T = 300°K
X02 = 1.0
Xi = 0.0 if i t 02
7. Steady State
6-25
-------
Reproduced from
bed available copy
Figure Al. U and V Distributions (X-0; 10 cm
-------
Figure A2. Temperature Distribution (X-0; 10 cm
-------
o
in
10
p;
5°
U"1
U.2
a
0«
o -
*i
IAJ
4
2C
Figure A3. XCq Distribution (X-0.0; 10 cm
-------
J?
li
=1
*3
= 3
1AI
9
-------
9
.-•-- ""-r1-- [•— L-,. 1 1 f 1;
J 1 1
% -
¦-+ . __ : : 1 —,
1
—^ 1
-1 ¦
«
U : r ' ' |
' r
¦ ¦ ' ! r
I
f • • ;
1
4
: H-
;
. . ... ¦ ; : - —
- —
-—1 . — 1- -
*
. .
1
t
- ; j r T -
•.
( ' i
\ i
; 1 I
?
- -
r- :-z . -
i /
\ \
! . 1
- —-r r :
......
I
- ¦
da.' 2..
: - _ ;
"..:
: / j
V j
I ; :
1
..
.. . ... /
i : f
.
B&.
H —i
1 , ; ! /
h V r i t
a
7
1 - . i—f
—l— 1——
I /
1
5
" "
•
\ t ; :
i
——
.
L_|. . u 1_
—j
— r-1
r-TT^ •"
^7 fr~ ¦ r"
—
: / • : l :
• tp" ; =.U:: : . , : ; •
1
I • "" - •:
- ;• - U - - F -
W1
_£
I
1 ¦ -"
--
\
- i V ~ ~
¦ t
- ¦
¦ f '
r^: 1
!
—V *—¦ —
t
*
l_—_j—
-A- T"-—
!
!¦
1
.
1 *
|—~ -
: V
| . j
r
: i
• - v •
• : : . J 1:
i "r :
. ' ¦ ¦ ' j : .
' - r:
:
to 4
v" : ¦
». ¦ J •
¦ -
- -
: h ¦ :
~ -f •"
-1
• H
-t
"= 1
1 ¦ -
—J-——I
K
- r 1
" f V
I =
4
in"5-
\ | • ¦ :
! 1 ' ¦- ¦ : ' i j \ 1 - !
i
"=•- - - !~
H—-— f- 1 | j \ ! »
; j |
r
"
1
Y—i
i
i
i • • j ¦: : ; : .
\
:
'
¦ i ¦
¦ ¦. L : -- :
' ' r
t —• -
' 1
• V
L.
—
-
- i \ 1 -.
ri \ H
MJ J
.....
.1
- •;
1=1 : ..
--r— ~P^"
• • u.
4
" •' =¦ •
- : .. . • . ;
: l - -
. ¦ -
.
. . _
. - f • : — -
"... - :
J
— -
:'":"
, ¦
__r:~
•.:£f=r--
;--
^ • ,
_ .
— .
•_
. ..
—
;
10 n 1213 14 15 16 17 18 19 20 21 22
y (cm)
Figure A5. XQ Distribution (X-0.0; 10 cm
-------
o
W
2C
h
Mgure A6
Distribution (X-0.0; 10 cm
------- |