id: :34! 5.. . . .
., finfiwfia. ,._ , ”$5
Er... 1...... .v. my
>4 V .
,m¥.vm..af #‘ «MY. 5.
mufififlwzmwfiw x
.2». 9, . “a,
, , a”. ”TE.” 3
r: z

7

n: .
Jr“...
r,
i
.
... y
1“

£45.
p.
. c
77“
M4 PW... .

.flP. .

.
a: ‘
1

finagx .. h .
fifinwumwumm
. awwafififiwn.

a. x
.u . .
.‘fl.
1. h.
.-

.

:2. I: .v. .
.

 

. 5.1.: .1. .

:1 4. PM 3:.

 

 

 

 

Wvflw

: .9

 

 

!.\-..:¥V

"$813

2601

This is to certify that the

dissertation entitled

Studies of Microgravity Diffusion Flames

presented by

Robert W. Vance

has been accepted towards fulfillment
of the requirements for

PhD. . Mechanical Engineering
degree In

 

 

    

Major professor

4/12/01
Date

 

MS U is an Affirmative Action/Equal Opportunity Institution 0-12771

 

 

LIBRARY

Michigan State
University

 

 

 

PLACE IN RETURN Box to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 C‘JCIRC/DateDue.p65-p.15

 

STUDIES OF MICROGRAVITY DIFFUSION FLAMES

BY

ROBERT W. VANCE

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements

for the degree of

DOCTOR OF PHILOSPHY

Department of Mechanical Engineering

2001

ABSTRACT
STUDIES OF MICROGRAVITY DIFFUSION FLAMES
By

Robert W. Vance

Various topics relating to flame stability are discussed. A linear stability analysis
is performed on a set of simplified equations to discern regions of both cellular and
oscillatory flame instabilities of diffusion flames. This is performed for a one-
dimensional model and also for a two-dimensional model representing a simple version
of an “edge flame”. It is shown that the one-dimensional model agrees well with the two-
dimensional model in a qualitative sense. A one~dimensional model of laminar premixed
flame annihilation is examined to identify the effect of non-unity Lewis number on flame
front extinction. Conclusions are drawn regarding the implications for turbulent flame
stability. Additionally, a heat transfer analysis is performed for a slot-bumer
configuration of a diffusion flame. Here, characteristic profiles for the conductive heat
loss from a lifted tri-brachial flame are examined. Correlations are presented and

compared with experimental data. -

ACKNOWLEDGEMNETS

I would like to acknowledge the members of my committee, Dr. I. Wichman, Dr.
C. Somerton, Dr. M. Miklavcic, and Dr. M. Zhuang, who provided considerable support
throughout. This dissertation is a product of their willingness to give their valuable time

and a reflection of their commitment to excellent engineering education and research.

iii

TABLE OF CONTENTS

LIST OF TABLES .................................................................................... v
LIST OF FIGURES .................................................................................. vi
CHAPTER 1
INTRODUCTION .................................................................................... 1
CHAPTER 2
ONE-DIMENSION AL DIFFUSION FLAME STABILITY ................................... 5
CHAPTER 3
TWO-DINIENSIONAL DIFFUSION FLAME STABILITY ................................. 39
CHAPTER 4
ONE-DIMENSIOANL LAMINAR PREMIXED FLAME ANNIHILATION ............ 55
CHAPTER 5
HEAT LOSS ANALYSIS OF A DIFFUSION FLAME LEADING EDGE .............. 100
CHAPTER 6
NUMERICAL TECHNIQUES .................................................................. 140

iv

LIST OF TABLES

Table 5.1 Correction-factor coefficients .................................................. 133

Table 5.2 Reaction information for Methane and Propane ............................. 124

List of Captions for Figues

Figure 2.1. Schematic Diagram of convective film diffusion flame ........................ 27

Figure 2.2. Qualitative plot of 1-D non-premixed flame Temperature vs
Damkohler number. Three distinct branches are traditionally represented.
The upper branch (above A) represents steady burning, the middle branch
is unstable, and the lower branch corresponds to the ignition branch. In
actual plots the lower turn appears at values of D many orders of
magnitude higher than D at the upper turn ...................................................... 28

Figure 2.3. The response curve in the vicinity of the upper turning point.
Here Le=1.0, Pe=0.0, Ta=4, ¢=1.0 ............................................................... 29
Figure 2.4. Leading two eigenvalues for the case of Pe=0, Le=1, Ta=4.0,
To=0.05. Here 02 has multiplicity 2 and is independent of D. Note that 0’1
is positive on the middle branch and negative on the upper branch (dotted
curve) .................................................................................................. 30

Figure 2.5. l<Le<Lec. Real eigenvalues on and 62 merge at B to form a
complex conjugate pair until they split again at E. B moves up and E
moves to the right as Le is increased. Here, Pe=0.0, Le=l.01, Ta=4.0,
To=0.05 ............................................................................................... 31

Figure 2.6. Leading two eigenvalues for Le=2.0 Ta=4.0, Pe=O, To=0.05.
Upper branch instabilities are seen up to point C’. Oscillations are
unstable between B and C’ and stable beyond C’. Pure perturbation
growth exists between A and B ................................................................... 32

Figure 2.7. When Le < 1, the third eigenvalue diverges from the constant
value obtained when Le = 1 by moving lower. As Le —> 1 0'3 merges with
0’2. Pe=0, Le=0.95, Ta=4, To=0.05 ................................................................ 33

Figure 2.8. The thermal-diffusive properties of the reactants plays a
significant role in determining the flamelet size as well as the location of
the onset of instability. Here Pe=0.5, Ta=5.0, Tor-0.05 ........................................ 34

Figure 2.9. The critical wave number for Le=1 for solutions along the
unstable middle branch are relatively unaffected by the crossflow velocity.
However, the increased cross flow results in a smaller flamelet after flame
front breakup. Le=1, Ta=5.0, To=0.05 ........................................................... 35

Figure 2.10 Experimental observations show a sinusoidal waveform prior
to flame front breakup. (Colors inverted) ........................................................ 36

vi

Figure 2.11 After the flame breaks up, the flame pockets that form are
dynamic in nature. The flame pockets do not behave uniformly in that
different sizes are observed. (colors inverted) .................................................. 37

Figure 2.12. Transient solution for Le=2 Pe=O Ta=4 To=0.05 at values of
D on either side of point C’. In one case the oscillations grow with time
and lead to instability. In the other, the oscillations decay leading to a
stable solution ........................................................................................ 38

Figure 3.1 Schematic of 2-D edge-flame model with boundary conditions. . . . . . . . . . .......47

Figure 3.2 Eigenvalue distribution for Le=1 Pe=0.5 Ta=4 ¢=l.0. Similar
distributions are seen for the 1-D model when Le=l ............................................ 48

Figure 3.3 Eigenvalue distribution for Le=0.9 Pe=0.3 Ta=4 ¢=l.0. 01
turns positive at the upper turning point, indicating instability along the
middle branch. The third eigenvalue approaches 0'2 asymptotically along
the middle branch .................................................................................... 49

Figure 3.4 Le=0.5 Pe=0.5 Ta=4.0 ¢=1.0. When Le is decreased below a
certain value, 03 crosses 02. This behavior is not seen in the one-
dimensional model. The Leading eigenvalue is not shown ................................... 50

Figure 3.5 Critical wave numbers for several different configurations. It
is seen that 1) Kc," is not greatly affected by changes in Le or Pe and 2)
The order of magnitude of Kc... is the same between the l-D and 2-D
models ................................................................................................. 51

Figure 3.6 Le=1.l Ta=4.0 Pe=0.0 ¢=1.0. Stable oscillations exist on the
upper branch while the entire middle branch is unstable. The regions
within the circle represent eigenvalues along the middle branch ............................. 52

Figure 3.7 Le=2.0 Pe=0.0 Ta=4.0 ¢=1.0. The gm branch eigenvalues
indicate regions of pure perturbation growth as well as oscillation growth
leading to instability. The oscillations cease to exist when 03 becomes the
leading eigenvalue ................................................................................... 53

Figure 3.8 Le=1.2 Ta=4.0 ¢=l.0. The three leading eigenvalues for
different Pe. Changes in Pe have only a small effect on the flame behavior

which is consistent with Le < 1 ................................................................... 54
Figure 4.1 Upstream extremum for H and location versus Le ................................ 85
Figure 4.2 Le=2, free stream ...................................................................... 86

Figure 4.3 Le=0.8, free stream .................................................................... 87

vii

Figure 4.4 Le=l.O, free stream .................................................................... 88

Figure 4.5 Le=2, T and Y profiles through reaction zone ..................................... 89
Figure 4.6 Le=0.8, T and Y profiles through reaction zone ................................... 90
Figure 4.7 Free stream laminar flame speed as a function of Lewis number ................ 91
Figure 4.8 Le=l.0 Reaction profile through annihilation ...................................... 92
Figure 4.9 Reaction profile for Le=2.0 through annihilation ................................. 93
Figure 4.10 Reaction profiles for Le=1.25 through annihilation ............................. 94
Figure 4.11 Reaction profiles for Le=0.8 through annihilation .............................. 95
Figure 4.12 Le=0.5 reaction profiles during Flame propagation towards

Junction ............................................................................................... 96
Figure 4.13 Temperature profiles for Le=2.0 through annihilation .......................... 97
Figure 4.14 Temperature profiles for Le=0.8 through annihilation .......................... 98
Figure 4.15 Annihilation time versus Lewis number .......................................... 99

Figure 5.1 Model configuration and boundary conditions. The side and
lower walls are porous to reactants. The vertical walls admit only
diffusive fluxes. In the far field the diffusion flame is one-dimensional.
The near field triple flame structure is shown, with (I) as the diffusion
flame, (II) and (III) as rich and lean premixed flame arcs, and (IV) as the
triple point or leading flame edge ............................................................... 128

Figure 5.2 Heat transfer model for infinitesimally thin flame sheet and
isothermal boundary conditions. This model is the heat transfer
simplification of the flame in figure 5.1 because by far the greatest heat
release is along the diffusion flame arc ......................................................... 129

Figure 5.3 Model configuration of mixture fraction Z after conformal
mapping. Note that the straight diffusion flame along Z=Zf occurs only in
the case of zero convection. In addition, lines of constant Z are radially-
directed arcs, Z=9l1t ............................................................................... 130

Figure 5.4 (a) Heat flux profiles comparison for straight flames rq=0.27.
Note the inaccuracy, generally, of the Burke-Schumann calculation of
Equation(12.a,b) and the comparative similarity of the Finite-rate and Heat
Transfer model results. A scaling factor renders the latter two visually

viii

indistinguishable.Heat flux comparisons for straight flame rq=0.9. (b) Heat

flux profiles comparisons for straight flames rq=O.9. The comments to (a)

apply here as well .................................................................................. 131
Figure 5.5 (a) Heat flux profile comparisons for fuel lean flames phi=0.5,

rq=0.38. For this non-symmetric case the flux distribution is asymmetric.

The Burke-Schumann model is once again a poor approximation to the

heat flux. (b) Heat flux profile comparisons for fuel rich flames phi=1.5,

rq=0.28. See comments in figure 5.4a .......................................................... 132

Figure 5.6 Heat flux profile comparison for straight flame with convection
Pe=1,rq=0.395 ..................................................................................... 134

Figure 5.7 Centerline temperature distributions for a straight flame

rq=0.375. The broadness of the actual flame moderates the temperature

variation, compared with the sharp heat transfer profile. The r values

along the diffusion flame (n>o.5) agree to within 043") .................................... 135

Figure 5.8 (a) Centerline temperature profile for a straight flame rq=0.525.

(b) Centerline temperature distribution for a straight flame rq=1.225. For

large quench distances the pure heat transfer temperature profile differs

significantly from the finite-rate chemistry profile ............................................ 136

Figure 5.9 Heat flux distribution for Methane/Oxygen combustion at
various stoichiometries ........................................................................... 137

Figure 5.10 Heat flux distribution for Methane/Oxygen combustion at
various stoichiometries. The second reaction mechanism given in table
5.2 is used ........................................................................................... 138

Figure 5.11 Heat flux profile for Propane/oxygen configuration. Larger
heat losses than Methane are predicted ......................................................... 139

ix

INTRODUCTION

Let the reader be forewarned that the material herein does not consist, as is
typically the case, of a thorough and complete analysis of a single problem. Rather, the
topics, though comprising four chapters, are actually three non-related problems that can
be considered as pertaining, either directly or indirectly, to the broad field of flame
stability. Due to the inhomogeneity of this work, each chapter, save for perhaps the third
for in truth it is an extension of the second, will be treated as a separate entity wherein
pertinent literature will be discussed, results presented and conclusions drawn. Hence,
this introduction shall serve only to provide the reader with a brief description of each
problem and indications of its importance.

As noted above, the topics of flame stability are very broad. This is only to be
expected given the complicated nature of the combustion processes. Instabilities can be
caused by imbalances in thermal and mass diffusion (Chapters 2 and 3), concentrations of
oxidizer or fuel that are insufficient to maintain steady burning, since combustion is a
kinetically controlled process excessive heat losses can cause flame extinction (Chapter
5). Additionally, hydrodynamic concerns can cause flame instabilities. In turbulent
flames large strain rates can effect flame behavior, flame front curvature (wrinkled
flames) is dependant on not only the flow field but thermal-diffusive effects as well.
Furthermore, the complex nature of turbulent flames show that flame fronts are often
dynamic and interact leading to local flame front breakup and healing (Chapter 4). The
above list of different types of instabilities is not complete but gives the reader an idea of
the complicated behavior of reacting flow problems. It is due to the aforementioned

complexity that any analysis seeking to understand the role of one particular element

must seek to systematically eliminate the effects of the others without ignoring the key
characteristics of the problem that cause such instabilities to occur in the first place. In
this spirit, each chapter consists of what may be termed a “minimalist model” in which it
is believed only the barest terms are retained in order to facilitate examination of the
phenomena in question. Such work, though not rigorously accurate, is simple enough
that the conclusions drawn are not obfuscated by unneeded complexity.

In Chapter 2 a linear stability analysis is performed on a one-dimensional model
of a diffusion flame in which the fuel is convected into the reaction zone while the
oxidizer is diffused. Of particular interest is the effect on non-unity Lewis number (a
ratio of thermal diffusion to mass diffusion) as well as the role of convection. The model
produced allowed the examination of the eigenvalues for several parameters: Lewis
number, Damkohler number (a ratio of characteristic flow time to a characteristic
chemical reaction time), and the Peclet number (non-dimensional mass flux). Several
types of instabilities were observed. The first consisted of pure flame extinction, the
second involved regions of flame oscillations (Le>l) that would either lead to extinction
or dampen out to steady burning, and the third resulted in a cellular flame pattern (LeSl).

The one-dimensional model can be considered a “classic” model having been
examined extensively in the past using Activation Energy Asymptotics (AEA). The
question that begs to be answered is, how well do one-dimensional models predict flame
behavior? It is not expected that these simplified models will predict quantitative
behavior but rather give an overview of flame characteristics. To begin to answer the
above question a two-dimensional version of the model in Chapter 2 is examined. In this

model the flame is examined in the presence of a cold, chemically inert wall. This

produces and “edge-flame” that better models real flame behavior. The linear stability
analysis is carried out along similar lines as in Chapter 2. Such an analysis is very
computationally costly but in the end offers reasonable support for the one-dimensional
model.

Chapter 4 investigates the thermal-diffusive effects on laminar, premixed flame
annihilation. Such a model can be thought of as representing the local flame interaction
in turbulent flow fields. How these flames interact under different Lewis number
conditions may offer some insight into how well turbulent flames can survive. The work
in this chapter is a combination of theoretical approximations and numerical analysis,
with the prior used as support for the latter.

The usefulness of one-dimensional models cannot be ignored. However, these
models lack important features that real flames incorporate. Of particular interest is the
interaction with a nearby surface. Such a feature was incorporated in Chapter 3 but was
absent in Chapters 2 and 4. A two-dimensional model is presented in Chapter 5 in which
a tri-brachial flame (“triple flame”) quenched a small distance above a “cold” wall. The
model represents a slotted burner configuration in which the reactants are separated by an
infinitely thin divider plate. Of particular interest is the conductive heat flux distribution
along the cold adjacent boundary. By evaluating this heat flux profile, deductions can be
made about the characteristic flame thickness; provide correlations that may be used in
theoretical work; and when compared to experimental results, provide insight into which
heat transfer mechanisms are important.

It is important to note, that although all the analysis performed in the above

chapters has a theoretical bent, in all chapters some ties to existing experimental results

are made. In Chapter 2 (and subsequently Chapter 3) comparisons are made between
predicted cellular flamelet sizes with cellular flames observed in micro-gravity drop
tower flame spread experiments. The ties to turbulent flame stability are more anecdotal
and compare only the observations made by several experimentalists about observed
Lewis number effects on stability. The correlations obtained for the heat flux distribution
are compared with two experiments burning methane gas. The maximum heat flux from
the bumer-attached flames is compared with the predicted values in Chapter 5.

A brief exposition on the numerical methods used to solve the various partial
differential equations (PDEs) is provided in Chapter 6. The numerical work in this
document does not represent in any way a novel approach to solving non-linear sets of

PDEs. The description is included for sake of completeness.

ONE-DIMENSIONAL DIFFUSION FLAME STABILITY

2.1 Chapter Nomenclature

A -— Pre-exponential factor for Arrhenius reaction rate, k=Aexp(E/RT)
A,B,C,E — Generic points along the “S” curve

D - Damkohler number, D=Alel(YFoY00)/0t

E — Activation Energy

K — Wave number

Le - Lewis number

I-—- Half of separation distance between porous plates
m — Mass flux

Pe — Peclet number, Pe=ml/2por

Q — Heat release

R - Ideal gas constant

T — Non-dimensional temperature, T: T (CpLe)/QY0o
t - Non-dimensional time

w - Non-dimensional reactivity = Dyoypexp(-Ta/T)
x,y - Non-dimensional spatial coordinates

Y - Species mass fraction

y - Non-dimensional species mass fraction, y=Yi/Y,o
Z — Mixture fraction

Greek

[3 - Zeldovich number

(I) - Global stoichiometric coefficient

(b - Generic scalar variable
A - Wavelength
1,5,1“ - Arbitrary spatial function
Subscripts
a — Activation reference
c,crit — Critical value
F — Fuel
f — Flame
0 — Oxidizer
o - Reference
max - Maximum
new — New
Superscripts
/ - Perturbation term
- - Steady state term
2.2 Introduction
In the study of diffusion flames the question of stability naturally arises when the
words “ignition” and “quenching” are used. These are “near-limit” phenomena in the
sense that either the temperature or the reactant concentrations (or both) can barely
support combustion, making the transition to either a burning state for the former and a
non-burning state for the latter a distinct possibility. These concepts can be broadened
when other physical, geometric, and parametric influences are examined, such as nearby

cold surfaces in the case of spreading flames over solid and liquid fuel, or multi-step

chemistry, or disparate Lewis numbers of the reactants, or radiant losses from the gas and
from soot particulates near the flame. Before flame extinction concepts are broadened,
however, the skeletal phenomenon should be understood in the clearest possible terms.
The present article attempts to create such an understanding for a physically idealized
model of diffusion flame stability.

In addition, a second goal is sought, namely a clear exposition of diffusion flame

stability for pedagogical purposes. The authors believe this is absent in the literature.
Our interest in diffusion flame stability lies in applications to flame spread. We found,
however, that there was no complete exposition of the stability calculation for even the
simplest model problems. The original study of Kirkby and Schmitz [1] comes closest to
a full description but is hampered by the use of difficult terminology.

The mathematical analysis of the stability of partial differential equations has a
long history. It has evolved in the last few decades into a theory based on the functional
analysis of PDES in terms of their eigenvalue spectra. The analysis of stability from a
mathematical standpoint amounts to the examination of the eigenvalue spectrum of the
linearized system of equations in relation to the original spectrum of the nonlinear system
of equations. For infinitesimal disturbances the confluence of the two systems can be
rigorously demonstrated, thereby rendering the analysis of the linearized problem
mathematically representative of the original (nonlinear) problem.

The eigenvalue spectra can be examined in numerous ways. The study by Kim et
al. [2] has, for the case Le<1, postulated simplified forms of the conservation equations
on either side of the flame sheet and then used asymptotic methods to derive “jump”

conditions across the flame sheet with the goal of obtaining formulas for the dependence

of the eigenvalues on the parameters D, Le, and Pe. Extending this work with
asymptotics, Kim [3] included higher order terms in the expansion of the inner zone
solution in a more thorough examination of effects of non-unity Lewis number. The
subsequent work of Kim et al. [4] employs no asymptotic methods in Le>1 studies of
diffusion flame stability. With regard to the latter work, which demonstrates a more

complicated response1

than the Le<1 case, stability is examined by integrating the
conservation equations numerically. Transient evolution of a monotonically growing
solution indicates instability.

The overall nature of the spectrum of the eigenvalues for the diffusion flame can
be understood both in a “physical” way and from experimental observations. Consequent
to Linan’s analysis of diffusion flame structure [5], Peters [6] determined by physical
arguments that the middle branch instability must be a “fast time” instability. The “fast
time” stability discussion is the result of a simple scaling analysis: if the reaction sheet is
thin, of order [activation energy]", then in order to retain the transient term in the
diffusion equation the temporal scale factor must be of order [activation energy]1 times
smaller than the spatial scale. Laboratory observations have also produced a physical
understanding of the nature and types of instabilities (cellular, oscillatory, etc.). Among
these experimental works are ref. [7,8], which deal with flame instability in many
configurations. In addition, many spectacular photographs can be found of ceiling flame
instabilities, see especially reference [9].

Our analysis will demonstrate that the eigenvalue spectrum is vitally important,

especially for the case Le > I examined recently by Matalon and Cheatham [10] and Kim

 

' The case for Le < l is less complicated in that no upper branch instabilities are seen and the middle
branch is unstable. Furthermore, no oscillations are present for this case as are seen when Le > 1.

et a1 [4]. We shall demonstrate that the tracking of the leading eigenvalue is insufficient
for understanding the instability process. The behavior of the second and third
eigenvalues often signals imminent transitional behavior.

2.3 Model

The proposed model is a one-dimensional version of the film-diffusion flame
examined by Kim et al. [2]. There, the authors were interested in the formation of
striped, instability patterns in a one-dimensional flame sheet. A minimum of two spatial
dimensions was required to describe the pattern of instability: one coordinate lay along
the flame sheet, while the other lay perpendicular to the flame sheet. The coordinate
along the flame sheet described the striping pattern. The fragmented flame sheets have
been called “flame tubes” [11]. In this study, we are interested in the fundamental
mathematical description of the instability mechanism. This does not require the
formation of different sorts of flame patterns. For this reason, following Kirkby and
Schmitz [1], we retain as little detail as required in order to produce a mathematical
description of diffusion flame stability. The in-flame coordinate is thus deemed
inessential for the stability analysis.

The proposed one-dimensional model consists of flowing fuel and diffusing
oxidizer species on Opposite sides of the reaction sheet. Radiative heat losses have been
ignored, as have both Soret and Dufour diffusion. The flow is assumed uniform. The
chemistry is modeled with a single-step, irreversible, Arrhenius reaction. Figure 2.1
shows a schematic of the one-dimensional model that has been used in previous stability
analyses by Kim and colleagues [2,4] and Matalon and co-workers [10,12]. The latter

investigators also describe a realistic physical configuration from which this model could,

in principle, be realized in a laboratory experiment. The problem formulation proceeds
along similar lines as that of [2], producing the following set of non-dimensional

equations.

2
gl+pe£=_a'£ (2.1.a)
at 8x 8x
83's BYI BZYi -
Le?+ PCLC—a—x— = ax—z - W ,1=O,F. (2..1b)

The equations are coupled through the non-linear reaction term, w=Dyoypexp(-
Ta/T). A single Lewis number is employed for both reactants. All terms in Eqs. (2.1) are
defined in the Nomenclature.

The boundary conditions for Eqs. (2.1) are

T=To, yp=1, yo=0, at x = -1 and T=To, yp=0, yozl, at x = +1. (2.2)

In this article we assume, for simplicity, that the fuel and oxidizer inflow
temperatures are both equal to a reference temperature, To. The symmetry of Eqs.(2.1)
and the boundary conditions implies that yo and yp are interchangeable; hence, from a
stability perspective, this formulation describes both blowing fuel and blowing oxidizer.

The configuration of figure 2.1 is a generic, simplified version of many
combustion problems. These include droplet combustion (blowing fuel from the droplet,
inward oxygen diffusion form the ambient), spray flames, jet flames, flame spread, and
diffusion flames in channels, among others. More complicated configurations require
additional geometry dependent terms; and more complicated model formulations may

require the inclusion of buoyancy or variable properties or pressure variations, or even

10

magnetic fields; but the qualitative nature of the solutions generated herein shall be
representative, as a model, for all of these configurations.
2.4 Steady State Solutions
2.4.1 Analytical Burke-Schumann Solutions

The character of the steady-state solution provides insight into the influences of
disturbances on the flame’s ultimate stability or instability. For clarification of
subsequent discussions, we present the steady Burke-Schumann solutions and briefly
describe their structure. The solution for the mixture fraction is

(1+x)
2

 

z = [exp(LePe(1+ x)) — 1]/[exp(2LePe) — 1] s- [1 — LePe(1 — x)/2 + O(LePe)2]. (2.3)

When Pe=0, we have Z=(l+x)/2. It is relatively straightforward to show that when Pe is
non-zero, Z is smaller everywhere (expect at xdl) than the Pe=0 Z-distribution, i.e.,
Zp¢>o < chzo. The flame is shifted to the right by the blowing fuel, and this shift
increases as Fe increases. We will not examine the large Pe case in this article.

At the flame sheet Z=1/2=Z[=[exp(LePe(1+xf)-l]/[exp(2LePe)-1], which can be
solved to give xr=ln[cosh(LePe)] for the flame sheet location. Burke-Schumann
solutions for oxidizder, fuel, and temperature distributions are obtained by solving
Eqs.(2.1) with w replaced with a delta function:

Oxidizer side (‘1<X<Xf):
Y0=0 (2.4.a)

= _ exp(LePe(x /2)) sinh(LePe(1 + x) / 2)
F exp(LePe(x, /2))sinh(LePe(1+ xf )/2)

 

(2.4.b)

T _ T ___ (T, _ T.) CXP(Pe(x / 2)) srnh(Pe(1 + x) / 2) .
exp(Pe(xf /2))srnh(Pe(1+ xf )/2)

O

(2.4.c)

 

ll

Fuel side (x,<x<l):

_ _ exp(LePe(x/ 2))sinh(LePe(1-x)/ 2)

 

 

O — . (2.4.d)
exp(LePe(xf / 2)) srnh(LePe(1 — xf ) / 2)

YF=0 (2.4.c)

T _ To = (T, _ To\ exp(Pe(x / 2)) sinh(Pe(l — x) / 2) (2.4.0

lexp(Pe(xf /2)) sinh(Pe(1 — xf )/2) '

These distributions are applicable in the regions on either side away from the
flame location. From these distributions it is possible to deduce some important facts (see
[2]). First, the outer-zone formula for Tf produces a Le-dependence in which Tf increases
as Le decreases and Tf decreases as Le increases. The roughly Le’l dependence of Tr is
well known and has been exploited in more detailed formulations of the non-dimensional
combustion equations [2]. Second, the general shape of the profiles in response to the
flow is ascertained, as is the response of the flame position to the fuel inflow.
2.4.2 Numerical Solutions

The steady state Eqs.(2.1.a,b) can be solved numerically by standard methods.
The numerical solution produces the characteristic “S” curve which plots Tmam vs D,
where Tmax is the maximum temperature in the domain. This “S” curve, to the knowledge
of the authors, was first produced and thoroughly examined by Fendell [13]. The top
branch of the curve corresponds to steady burning (under most but not all conditions), the
middle branch is unstable; and the lower branch is known as the ignition branch, see
Figure 2.2. The real “S” curve does not look like an actual S since the lower turn occurs
at a value of D that is usually many orders of magnitude larger than its value at the upper
turn. A representation of the upper turning point for Le=1, Pe=0, Ta=4 and ¢=l is shown

in Figure 2.3. At the lower turning point the Arrhenius exponential is O(exp(-E/RTO))

12

whereas at the upper turning point the Arrhenius exponential is O(CXP('E/RT{)). The
“separation” in D between these two turning points is thus O(exp(E/R)(1/'l‘o-1/'I‘ f)), at
least in the lowest approximation (i.e., without accounting for differences in reactant
mass fractions at the two turning points).An analysis of the “S” curve using large
activation energy asymptotics was carried out by Linan [5]. Combustion regimes in the
various parts of the “S” curve were elucidated and canonical lowest order equations
describing the mathematical character of each region were produced with accuracy that
improved as the “Zeldovich number”, [3, increased. The expansions, which proceed for
both inner (flame zone) and outer (chemical equilibrium zones) in inverse powers of B,
can be continued, at least in principle, indefinitely. Inner and outer solutions are matched
in the characteristic fashioned outlined in [14]. The expansions in B" are valid
everywhere along the “S” curve but the relative accuracy of the solutions diminishes near
the turning point where additional (higher-order) series terms are needed to maintain
accuracy. The leakage of reactants was shown to be a cause of the transition from
vigorous to non—vigorous burning. Numerous features of this S curve have been clarified
since the work of Linan [5], but the fundamental distinctions made therein (and in the
original work of Fendell[l3]) remain valid. Though some aspects of the preceding
description, such as the use of one-step chemistry, or the predominant leakage of the
“wrong” reactant, etc. [15] are unsuitable for a more comprehensive approach, the appeal
of the “S” curve remains strong. Discussions of diffusion flame stability are referred
almost exclusively to “S” curve solutions [2,4,10]. Extending or broadening them to

more realistic flame processes, geometries, and parametric ranges is desirable.

l3

2.5 Stability Analysis

Equations (2.1.a,b) form a set of coupled semi-linear parabolic equations [16,17].
Hence the evolution of small disturbances is governed by linearized forms of these
equations. These linearized equations are derived by writing the solution of each

dependent variable as a combination of a steady solution and a perturbation
¢<x.t) =6<x>+¢’(x.t>. (2.5)
where (I) can be T, yo, or yp, Substituting Eq. (2.5) into Eq. (2.1) and retaining only the

terms that are linear in ¢' in the Taylor series expansion of w yields the following linear

 

 

 

system:

ar’ an 321" ,
+ Pe = + w , 2.6.a
at 3x 8x2 ( )

I I 2 I
Lnfl + PnLnQYi = 9124 — w', i=O,F , (2.6.b)
at 8x 8x
I — — Ta I — I I — Ta

W =D YOYF fiT +YOYF+YOYF XCXP “'I-‘I.’ , (2-6-9)

which is solved subject to the following boundary conditions (no disturbances at the
boundaries):
T'(il,t) = y'F(il,t) = y'0(i1,t) = 0. (2.6.d)
Since the spatial part is a linear operator with compact resolvent [17], the evolution of
small disturbances is determined by the eigenvalues of Eqs.(2.6.a-d). The eigenvalues are
the complex numbers, 0', for which Eqs.(2.6) possess a nontrivial solution in the form
T'(x, t) = t(x)exp(ot) (2.7.a)

yn(x.t) = 5(X)6Xp(ot) (2.7.b)

l4

y'F(x, t) = F(X)exP(6t). (2.7.c)
If all eigenvalues have negative real parts, then all solutions of the nonlinear Eqs. .( 1) that
originate as small disturbances of the steady state solution will decay exponentially to the
steady solution [16,17]. In this case the steady state solution is said to be stable. On the
other hand, if any of the eigenvalues has a positive real part then the steady state solution
is unstable, meaning that there exists a threshold disturbance size such that some
arbitrarily small initial disturbances of the steady state solution will evolve according to
the nonlinear Eqs.(2.1.a,b) to grow eventually past the threshold [16,17].

With fixed parameters the steady state solution was found first, using
Mathematica’s NDSolve in the shooting algorithm for the boundary value problem.
Error controls in the software enable the successive elimination of subsequent higher
order error terms; the steady solution is for all practical purposes exact. In other words,
the errors can be made as small as desired because it is known that the computed errors
depend quadratically on grid size. Three different meshed were used (80, 160, 320 grid
points) to reduce the computational errors.

This steady solution was used in Eqs.(2.6.a-d) along with Eqs.(2.7.a-c) for the

disturbance quantities. The resulting disturbance equations are given below as equations

2.8.a.d.
2
to+Pe§1=fi+wC (2.8.a)
3X 3x2
2
Le50+PeLe§ =9—5- w' , (2.8.b)
3" 3x2
2
LeF6+ PeLe-a—I: = 8—: - w' , (2.8.c)
3X 3x2

15

, _ _ T _ T
w =D[YOYF[:—r-_Z'T]+ yOF+FyF:JXexp(-—%—). (2.8.d)

The equations were then discretized using a second order central difference scheme. The
resulting matrix eigenvalue problem was also solved by using Mathematica. The
computational errors of the eigenvalue also depend quadratically on the mesh size, hence
we were able to reduce eigenvalue errors below a specified tolerance. For example, when
Pe = 0.5, Le = 0.9, To = 0.05, Ta = 5, and D = 21190000 we found that the leading
eigenvalue on the middle branch was equal to 0.1612, and on the upper branch it was
equal to —0.1522. Many spot checks were made by using different numerical methods.
In particular, a Fortran based finite difference scheme was used (with no error controls)
with 300 nodes to produce the steady solution, and the eigenvalues were calculated using
standard Lapack subroutines.

The stability in the neighborhood of the upper turn, point A on figure 2.2, on the
“S” curve has been investigated for a large set of parameters:

OSPe $1, 0.1 SLe 54,1.5 STaS6.5, 0.02 STo $0.05.
A brief description of the findings appears below. In later sections lengthier discussions
of their physical significance are presented. For larger values of Pe the eigenvalue
behavior is more complicated; however, and more detailed analysis is needed before
concrete generalizations are made.

Identifying the Le effect on the eigenvalue behavior is essential in understanding
the various regimes of flame behavior. When Le S 1, the leading eigenvalue is negative
on the upper branch, implying stability, and positive along the middle branch, implying
instability. The leading eigenvalue is equal to zero at the turn, point A (see figures 2.3

and 2.6). When Le = 1 a linear combination of Eqs.(2.1) commonly called the Schvab-

l6

Zeldovich procedure produces two stable linear, homogeneous heat equations with a
convection term, and a nonlinear equation. The two heat equations are responsible for an
eigenvalue of multiplicity two that does not change with D (see figure 2.4 curve denoted
“62”). This constant eigenvalue is given by
-(1t/2)2/Le —(Pe/2)2Le. (2.9)
This result is obtained for yo-yp derived from Eq.(2.6.b) which is a linear heat equation
with a convective term, and can be solved explicitly. When Le=1, the solution has
multiplicity two because 2T+yo+yp satisfies the same equation. The remaining nonlinear
equation produces the eigenvalue denoted “0’1” in figure 2.4. This eigenvalue is always
positive on the middle branch, indicating instability as stated above. (Note that in all the
figures pertaining to eigenvalue distributions, the dotted line indicates a position along
the gage; branch.) It is the character of the remaining two eigenvalues on this middle
branch, however, that changes as Le is varied; the nature of this change dictates the
stability of points along the sector ABC of the “S” curve of Figure 2.2.
When Le is slightly larger than unity, the eigenvalue that has multiplicity 2 when
Le = l splits into 2 eigenvalues. One, 0'3, is still given by Eq. (2.9) and the other one, 62,
interacts with 0'] as Shown on Figure 2.5. 0'] and 62 are negative real numbers between
the turning point A and point B on the upper branch (see figure 2.2). At point B they
merge to form a complex conjugate pair with negative real parts (indicating decaying
oscillations) until point B where they again split. At point C, 0'3 becomes the leading
eigenvalue. As Le approaches unity the region between B and E shrinks to point C where
61 meets 62 on figure 2.4. As Le increases the value of 0’] at B increases, B approaches

A, and B reaches A at some critical value Lec > 1.

l7

When Le increases past the critical value Lee, point B moves into the region of
positive eigenvalues, Re(6) > 0, and the characteristic behavior in region ABC of Figure
2.2 is changed completely (see figure 2.6). Notice that figure 2.6 (Le>Lec) is
qualitatively the same as figure 2.5 (Ledec) with the point E (not shown) moving rapidly
towards higher D as Le is increased. When Le>Lec, the region between points A and B
on the upper branch has eigenvalues that are real and positive, indicating that perturbation
growth will occur and lead to flame extinction. Between B and C the two leading
eigenvalues form a complex conjugate pair, whose real part is positive at B and decreases
as D increases. The real parts become zero at a point C’ between B and C indicating
damped oscillations until point C (not shown) where the leading eigenvalue (0'3) is
independent of D.

On the other hand when Le is decreased below unity, 63 moves under the constant
eigenvalue as shown in figure 2.7. The constant eigenvalue is still given by Eq. (2.9).
Notice that along the upper branch all the eigenvalues are negative, while along the
middle branch 0'] is positive and 623 are negative. Thus, for Le<l the upper branch is
stable and the middle branch is unstable. The remaining constant eigenvalue disappears
in the case Leo It Le]: #1.

It has been shown that the system behavior changes when Le is increased passed
unity. Damped oscillations appear when l < Le < Lec, but the entire upper branch
remains stable. The beginning of the upper branch is unstable when Le > Lee.

2.5.1 Cellular Flames
“Cellular flames” have been observed in both freely propagating premixed flames

[18], burner attached diffusion flames [19,20], and micro-gravity flame spread. All these

18

systems are fairly complex, involving multidimensional processes that generally include
flame-surface interactions. The qualitative nature of the cellular flame instability can be
described by examination of the proposed one-dimensional problem.

In order to allow for spatial variation along a plane of unit depth, the Eiz/ax2
operator in Eqs.(l.a,b) needs to be replaced with the Laplacian aZ/axz + 32/3)!2 and
perturbations take the form ¢’(x,y,t) = f (x)exp(6t + iKy) in place of Eqs.(2.7.a-c).
Here, K is a non-dimensional wave number in the new coordinate direction, y, which lies
in the plane ((y,z) plane, where z is the unit-depth coordinate) of the flame sheet. This
coordinate dependence will show the development of undulations in the flame, whereas
the one-dimensional formulation merely produced fluctuations about the mean flame
position, Xf. When Le = 1, the new eigenvalues, 6mm, are simply equal to the old
eigenvalues o shifted by K2, i.e. on... = o - K2. When Le e 1, the relation between o...w
and 6 is far more complex but qualitatively identical. In particular, the introduction of K
has a stabilizing influence on the flame. Two conditions need to be met in order to see the
flame stripes in the y-coordinate plane. First, the boundary control requires that K be a
multiple of 21t/L where L is the dimensionless length (y-coordinate) of the burner. If, for
example, the ratio between the length and the width of the burner is 10, then K can take
values 0.31, 0.63, Second, the new leading eigenvalue has to be very close to zero to
ensure the temporal changes of the flame structure do not rapidly vary and destroy the
striped flame pattern. The first condition, pertaining more to experimental aspects, will be
ignored in the rest of the discussion. For each unstable steady state solution near point A
on the “S” curve a K value is computed, called Kcm, such that the new leading eigenvalue

on the middle branch is identically zero. When 0 < Le < Lec on the middle branch away

19

from the upper turn on the “S” curve, Km. gradually increases from zero at point A as D
is increased. When Le > Lee, the leading eigenvalue is positive near the turn A; and,
hence, a finite K, sufficiently greater than zero, is needed to stabilize the unstable steady
solution, see figure 2.8. Therefore, it is clear from figure 2.8 that cellular flames will be
more easily witnessed for conditions satisfying Le < Lee.

Of particular interest is how thermal-diffusive effects contribute to the
determination of cellular flame size. Additionally, the convection of fuel as described by
Pe is examined to identify the hydrodynamic effects. Thermal-diffusive effects are
examined by varying the Lewis number and holding Pe constant at 0.5. The results for
this numerical experiment are illustrated in figure 2.8. It is apparent that the influences of
varying Le are relatively small. The critical K values are only moderately changed
between the cases Le = 0.9 and Le = 0.3. This suggests that thermal-diffusive effects
play only a minor role in determining flame cell size. The distance, measured in D, away
from the instability point A is the main determinant of the size of the flame cells

In order to determine the effects of fuel flow on the cellular flame size several
different values of Pe were examined. The range of Pe was varied between zero and
0.75, and Le was set equal to unity so that the diffusion term and the convective term
were approximately of the same order of magnitude. In this range of Pe, the “8” curves
retain their characteristics and are simply shifted towards higher D ranges. Figure 2.9
shows the Kc... values for the three different flow rates, Pe=0, 0.5, and 0.75. There is
little difference between the Km, values for a zero convection diffusion flame (Pe = 0)
and one for which convection is appreciable, i.e., Pe=0.75. This leads to the conclusion

that for order—unity fuel crossflow the size of the cellular flames is little affected.

20

2.5.la Experimental Results

Calculating the wavelength in the direction normal to flame is difficult because no
physical length scale exists. There is, however, a diffusion length scale that appears in
the Damkohler number. This length is used in making comparisons with experimental
observations. Furthermore, as noted above, a smaller wave number does not necessarily
mean a larger wavelength, as can be seen in equations (2.10-2.11). The following
relationships were used to calculate a wavelength based on a diffusion length scale.

I
2
L: 4D. ‘L (2.10)
YOOLeA

(2.11)

Nit“

Here D is the Damkohler number, Otis the thermal diffusivity, Yoo is the mass fraction of
oxygen in the oxidizer stream, A is the Arrhenius pre-exponential factor, and A is the
disturbance wavelength. For comparisons with experimental results, where a thin
cellulose sheet was burned, the parameters for A and or were approximated by using
propane as the fuel. The values for K and Da were taken very near the turning point, i.e.,
(D-Dmm)~1000, where Dmm was found to be 1.7984e7 for the case where Le=1. It is
noted that in this region, near the stable limit, all the curves from figure 2.9 seem to be in
very good agreement. Based on these numbers, the estimated flamelet wavelength was
found to be on the order of four centimeters. This result is compared with experiments
that show the size of the flame pockets to be on the order of several centimeters. Figure
2.10 illustrates a sinusoidal waveform prior to breakup with the characteristic wavelength

shown. The distance between the peaks is 3.5 cm. Similarly, Figure 2.11 shows several

21

flame pockets after breakup. It is clear from this photograph that not all the flame
pockets behave identically. In this instance the flamelet size varies by a large percentage,
but the relative order of magnitude is still several centimeters. These results Show that
that the simplified one-dimensional model was capable of predicting the order of
magnitude of the size of the flame pockets.

2.5.2 Oscillatory behavior (Le>Lec)

Oscillatory behavior in flames has been observed in microgravity candle

experiments aboard the Mir space station, as well as during droplet burning experiments
[21]. Several numerical investigations have examined oscillatory flames [4,10]. The
qualitative nature of this phenomenon can also be described with a linear stability
analysis of a one-dimensional diffusion flame.
When Le > Lee, the leading eigenvalues form a complex conjugate pair beyond the point
B on the upper curve (see figure 2.2 and 2.5). The real part is positive at B and decreases
as D increases. At point C’ the real part reduces to zero although the imaginary part is
non-zero. So, between point C’ and C on the upper curve, small perturbations will result
in damped oscillations. Between B and C’, small perturbations will oscillate with
growing amplitude leading to flame extinction. A stable periodic orbit was not found
numerically near point C’. Near C’, long lasting oscillations are expected with circular
frequency equal to the imaginary part of the leading eigenvalue pair.

The oscillatory behavior predicted by the linear stability analysis was supported
by solving Eqs.(2.1.a,b) directly using a second order correct finite-difference method. In
doing so, a graphical representation of the transient nature of the flame was achieved as

well as a search for possible steady oscillatory behavior. The results in Figure 2.12

22

represent transient solutions for two points along the upper branch. The unstable solution
is representative of a D value below point C’ (Dc~=l .961E+6), while the stable solution is
for a D larger then that at point C’. The numerical investigation of the oscillatory
behavior is similar to that performed by Kim et al. [4]. Although the method is able to
identify all the characteristics of this simplified problemthat have been identified by the
preceding stability analysis, it is time consuming and does not allow Lec to be readily
evaluated.
2.6 Discussion and Conclusions

The mathematical examination of Eqs.(2.1.a,b) and the linearized small
perturbation Eqs.(2.6.a-d) in terms of the first three eigenvalues leads to a clear picture of
the overall phenomenon. The principal virtue of the model considered here is its
simplicity and the absence in it of complications that might difficulties of interpretation.
It has been stated that convection is “essential” for instabilities to occur in the
configuration of Figure 2.1 [2]. Our examination produces results for Pe=0 that are
qualitatively identical with those for PeabO. It appears that a nonzero value of Fe does not
alter the nature of the stability calculation, even though complications enter as Pe
becomes large. The stability results for Le < Lec were qualitatively simple, as the results
of section 2.5.1 demonstrate. The situation for Le > Lee, however, was more complicated.
Here, several possibilities existed, whose realization depended primarily on the
movement of the locus of the first two eigenvalues as a function of D and Le. Since a
complex conjugate pair is formed on the upper branch after point B, oscillatory behavior
there is the norm. Between A and B there are no oscillations. The critical value Lec is the

value of Le > 1 for which point B is located at the turning point A of Figure 2.2.

23

It is possible that asymptotic analyses can provide approximate formulas for the
information obtained here by numerical methods of eigenvalue computation. The
asymptotic calculations, however, are tedious and lengthy (see e.g. [2] for the Le < 1
case), and provide results in terms of parameters that are, in practice, difficult to
calculate. An example of a difficult parameter to calculate is the Damkohler number
derivative [2]. The calculation of the first few eigenvalues as outlined here, however, in
complete accord with the mathematical theory of systems of equations [16,17]. The
eigenvalue computation provides an immediate and complete assessment of the diffusion
flame response, not only in the vicinity of specially chosen points of analysis (such as A)
but everywhere along the entire “S” curve.

The Nyquist-plot analysis of Kirkby and Schmitz [1] also provides
mathematically rigorous results. Nevertheless, considerably greater effort is required to
interpret and employ the results presented therein. The difference between this work and
that of Kirkby and Schmitz [1] is not mathematical rigor but ease of interpretability and
consistency with the function-theoretic analysis of partial differential equations.

For the reasons given in the previous two paragraphs, we believe that our
examination provides a comprehensive and readily understood description of diffusion
flame instability. Additional analysis of the eigenvalue spectrum is suggested by the
present analysis for the case of large Pe (where the behavior of the eigenvalues gets more
complicated) and when the values of Le differ for the two reactants. The permutations in
the latter case are numerous, as there is a blowing reactant and a diffusing reactant; Le for
both reactants may be small, large, one large the other small; and so on. In addition, the

behavior of the solution on the upper branch for Le > LeC near point B may warrant

24

theoretical analysis because a bifurcation exists there (real parts of two complex
conjugate eigenvalues cross from positive to negative values in the direction of increasing
D).

2.7 References

1. LL. Kirkby and RA. Schmitz Combustion and Flame 10 205-219 (1966)

2. J .S. Kim, F.A. Williams, and PD. Ronney Journal of Fluid Mechanics. 327 273-
301(1996)

3. J .S. KimCombustion Theory and Modelling 1 1340(1997)

4. C. H. Sohn, et Al. Combustion and Flame 117 404-412(l999)

5. A. Linan Acta Astronautica 1 1007-1039(l974)

6. N. Peters Combustion and Flame 33 315-318 (1978)

7. Y. Zhang et al. Combustion and Flame 90 71-83 (1992)

8. RH. Chen, G.B. Mitchell, P.D. Ronney Twenty-fourth Symposium (International) on
Combustion 213-221 (1992)

9. M. Kokkala VTI' Tutkimuksia 586 (1989)

10. S. Cheatham and M. Matalon Twenty-sixth Symposium (International ) on Combustion
1063-1070 (1996)

11. J. Liu, P.D. Ronney Combustion Science and Technology 144 21-45 (1999)

12. M., Matalon 37th AIAA Aerospace Sciences Meeting, AIAA 99-0584 (Reno, NV)
(1999)

13. FE. Fendell Journal of Fluid Mechanics 21 281-303 (1965)

14. J. Cole Perturbation Methods in Applied Mgthemaflcg, Blaisdell Publishing(1968)

25

15. A. Linan, and F. A., Williams Fundamental Aspects of Combustion, Oxford
University Press (1993)

16. D. Henry Geometric Theory of Semilinear Parabolic Equations, Lecture Notes in
Mathematics, No. 840, Springer-Verlag (1981)

17. M Miklavcic Applied Functional Analysis and Partial Differential Equations, World
Scientific (1998)

18. C. M. Dunsky Twenty-fourth Symposium (International) on Combustion 177-187
(1992)

19. Class et. Al SIAM Journal of Applied Mathematics 59 n.3 942-964 (1999)

20. R. Kuske and BI. Matkowsky Quarterly of Applied Mathematics 52 n.4 665-88
(1994)

21. S.L. Olson Buoyant Low Stretch Stagnation Point Diffusion Flgmes over a

Solid Fuel, Ph.D. Dissertation, CW RU (1997)

26

 

 

 

 

 

X=-1 X=1

g 7 u 3

IL 5 :

Porous wall Porous wall
with convection Reaction Sheet (OXIdIZBI' Dlflsulon)

Figure 2.1 Schematic Diagram of convective film diffusion flame

27

Oxidizer ‘

 

 

D

Figure 2.2 Qualitative plot of 1-D non-premixed flame Temperature vs
Damkohler number. Three distinct branches are traditionally represented.
The upper branch (above A) represents steady burning, the middle branch is
unstable, and the lower branch corresponds to the ignition branch. In actual
plots the lower turn appears at values of D many orders of magnitude higher
than D at the upper turn.

 

0.55 -

0.5 -

0.45 ~

Tmax

0.4 2

 

0.35 a

 

 

0.3 i i i i
1 .5E+07 2.3E+07 3. 1E+07 3.9E+07 4.7E+07

D

Figure 2.3. The response curve in the vicinity of the upper turning point. Here Le=l.0,
Pe=0.0, Ta=4, ¢=l.0.

29

 

 

 

 

 

 

 

I I ‘I 1 I I I

1.5 I- -

1 P 4

01
0.5 - -
Re(0)0

{1.5 - -

i- d

-1 L ..

4.5L ~. ~ . -

-2 - i ' ~ .

25 - .......... _

G: """"""

-3 l- .

_35 L 1 r I; r r I
1.715 1.72 1.725 1.73 1.735 1.74

x 11]6

Figure 2.4 Leading two eigenvalues for the case of Pe=0, Le=1, Ta=4.0,
To=0.05. Here 62 has multiplicity 2 and is independent of D. Note that 0'1
is positive on the middle branch and negative on the upper branch (dotted

curve).

30

 

 

 

 

l I

_4 l l l l l 1
1.715 1.72 1.725 1.73 1.735 1.74 1.745 1.75 1.755 1.76

D x108

Figure 2.5. 1<Le<Lec. Real eigenvalues 0'1 and 62 merge at B to form a
complex conjugate pair until they split again at E. B moves up and E
moves to the right as Le is increased. Here, Pe=0.0, Le=l.01, Ta=4.0,
To=0.05.

 

31

 

14 I I I l I I I

 

 

 

12 - 61 -

10 - -
. A

8 '- ll '1
Re(6) .

5 ~ -

4 r a 4

IE
2 b I" ~~~~~~~~~ 'l
in." ~~~~~~~~~~ x_“ Cl
0 NEE ---~e,._q__m_
G: M5.”

 

_2 r l r
1.7 1.75 1.8 1.85 1.9 D 1.95 2 2.05 2.1

x106

Figure 2.6. Leading two eigenvalues for Le=2.0 Ta=4.0, Pe=0, To=0.05.
Upper branch instabilities are seen up to point C’. Oscillations are
unstable between B and C’ and stable beyond C’. Pure perturbation
growth exists between A and B.

32

 

 

 

4 _ .
Cl
2 _ -
Re(6)
0 ~ :A -
2 . .................................. 02 ‘
Kr 6
.4-A 3 d
-8 _ ........ - ' ...... .

 

 

 

1.71 1.72 1.73 1.74 1.75 1.76 1.77 1.78 1.79 1.8 1.81

x106

Figure 2.7. When Le < 1, the third eigenvalue diverges from the constant
value obtained when Le = 1 by moving lower. As Le —> 1 0'3 merges with
62, Pe=0, Le=0.95, Ta=4, To=0.05.

33

2.5

 

 

 

 

Le=1.5
2 _/ Le=0.9
Le=0.5
Le=0.3
1.5 ..
I(crit
1 l
0.5 1
0 . i . i
0 50000 1 00000 1 50000 200000 250000

D-DA

Figure 2.8. The thermal-diffusive properties of the reactants plays a
significant role in determining the flamelet size as well as the location of the
onset of instability. Here Pe=0.5, Ta=5.0, To=0.05.

34

 

2.5

 

 

 

2 ' Pe=0.5
Pe=0.7

1.5 -

1 "‘ //

’/

0.5

0 i i i .

0 50000 100000 1 50000 200000

D-DA

250000

Figure 2.9. The critical wave number for Le=1 for solutions along the unstable
middle branch are relatively unaffected by the crossflow velocity. However,
the increased cross flow results in a smaller flamelet after flame front breakup.

Le=1, Ta=5.0, To=0.05.

35

     

2.10 Experimental observations show a sinusoidal wav orm prio
to flame front breakup. (Colors inverted)

36

 

 

Figure 2.11 After the flame breaks up, the flame pockets that form are
dynamic in nature. The flame pockets do not behave uniformly in that
different sizes are observed. (colors inverted)

37

 

0.49
0.485 -
0.48 ‘
0.475 (\f‘ __ =2.05e6
- 0.47 a
0.465 r
0.46 ‘
0.455 ‘
0.45 ~
0445 q D=1.955e6
0.44 r .
0 5 10 15
time
Figure 2.12. Transient solution for Le=2 Pe=0 Ta=4 To=0.05 at values of D on

either side of point C’. In one case the oscillations grow with time and lead to
instability. In the other, the oscillations decay leading to a stable solution.

 

 

Tmax

 

 

 

38

TWO-DIMENSIONAL DIFFUSION FLAME STABILITY

3.1 Chapter Nomenclature

D — Damkohler number

K — Wave number

Le — Lewis number

Pe — Peclet number

T — Non-dimensional temperature

t — Non-dimensional time

w - Non-dimensional reactivity = Dyoypexp(-Tfl)
x,y - Non-dimensional spatial coordinates
Y - Species mass fraction

y - Non-dimensional species mass fraction
Greek

(1) - Global stoichiometric coefficient

2. - Wavelength

Subscripts

a - Activation reference

c,crit - Critical value

F - Fuel

f — Flame

0 — Oxidizer

o - Reference

39

3.2 Introduction

The one-dimensional flame examined in Chapter 2 is a highly idealized model.
However, with only one-dimension there is no possibility of the flame interacting with a
nearby surface; and hence, some robustness may be lost. As noted in Chapter 2, there has
been a large amount of work focused on various configurations of one-dimensional
flames. There has not been, to the knowledge of the author, any work done involving
multiple dimension, stability analysis to help support the validity of the one-dimensional
model. Some work, in similar fashion to [1], states to be a multidimensional analysis
solely because a dimension perpendicular to the base state is added, as in Chapter 2, to
identify critical wave numbers even though the base state used is still one dimensional.
Some interesting work has been done by Lingens et al. [2] in which experimental data is
used to model the flow field of a bumer-attached flame. The stability analysis that is
carried out is done using a one-dimensional model in a radial coordinate but is applied at
various locations above the burner, so in some sense the work took on a
multidimensionality. Some numerical work has been done with two-dimensional models
to identify cellular flames [3] for Le<1, or flame oscillations [4] for Le>l.

The question this chapter seeks to address is, “How does the flame behavior
change with the addition of a second dimension?” Since real flames interact with
surfaces it is important to know if the one-dimensional models are missing some
important characteristics of this interaction.

3.3 Model
The governing equations are similar to those found in the previous chapter, with

the exception of a second coordinate.

40

ar P or a’T azr

—a—t—+ e—a-;;=—a—x—2—+—a;§-+W (3..la)

aYi aYi BZYi BZYI -
LeE+PeLe 8x = 8x2 + ayz —w,1=O,F. (3.1.b)

 

The equations are coupled through the non-linear reaction term, w=Dyoypexp(-T,fl‘). A
single Lewis number is employed for both reactants. Figure 3.1 offers a graphical
representation of the model.

The boundary conditions for equations (3.1) are

T=To, yrs-=0, yo=1, at x = -1 and T=To, yp=1, yo=0, at x = +1. (3.2.a)
30)
T=To, yp=0, yo=0, at y = 0 and 3;— = 0 at y = 3 (3.2.b)

It is assumed that the far field flame, in the “in flame” coordinate, does not change with y
and the flame is essentially one dimensional at this point. This assumption would be
more valid if the “in-flame” boundary could be moved farther away from the “cold”
lower wall, but computational limitations restrict the number of nodes available.
Furthermore, we assume, for simplicity, that the fuel and oxidizer inflow temperatures
are both equal to a reference temperature, To<<T f.

The present model offers a simplistic representation of a flarne-spread
configuration with blowing oxidizer and diffusing fuel. Additionally, the flame develops
a “flame-edge” [3,5], which will offer generally qualitative characteristics of flame-edge
behavior. It is possible to obtain a representative model for slot burners by changing the
convective term to model flowing reactants in the in flame coordinate (y-direction), this
configuration produces another edge flame namely a triple-flame structure (see Chapter

5). The present model, however, is closest to the one-dimensional model studied in

41

Chapter 2 with essentially the addition of a cold inert surface to add the multi-
dimensionality.
3.4 Stability analysis

The stability analysis conducted is similar to that performed in chapter 2 and will
not be expounded on here. There are, however, some limitations that arise when
examining a two-dimensional base state. Namely, the eigenvalue calculations must be
done over the whole domain, which makes the calculation very computationally
inefficient. For this reason a relatively coarse grid was needed and, hence, some loss of
accuracy near the static turning point (see figure 2.2 point A) is unavoidable. A finite
difference algorithm was used employing a 92 x 92 uniformly spaced mesh. The
eigenvalues were obtained by employing standard Lapack subroutines designed for large
banded matrices. The coarse nature of the grid does not produce smooth, continuous
eigenvalue curves, as does the one—dimensional model. However, it is believed that the
results represent a qualitative picture of the flame behavior near the upper turning point
of the “S” curve and will allow conclusions to be drawn about the usefulness of a one-
dimensional model.
3.4.1 Eigenvalue results Le S. 1

As discussed in the previous chapter, the governing equations can be reduced to a
non-linear equation and two “heat equations”, which produce a constant eigenvalue of
multiplicity two when Le=1. Because of this, the eigenvalue behavior for unity Le is not
expected to vary greatly from the one-dimensional model. This is in fact the case as can
be seen in figure 3.2. Here 01 represents the leading eigenvalue which is negative along

the entire upper branch and turns positive (indicating instability) at the turning point. The

42

second two eigenvalues are independent of D and are always negative. This is
qualitatively the same picture as the one-dimensional model (see figure 2.3). It is seen
that the required coarse grid does not produce smooth curves in the vicinity of the upper
turning point on the “S” curve. The solutions along the middle branch are difficult to
obtain, and produce “ragged” curves.

The picture becomes slightly more complicated when Le is reduced below unity.
As shown in the previous chapter (figure 2.6) there existed, among the three leading
eigenvalues, only one constant eigenvalue. It is the behavior of what is labeled the third
eigenvalue, 63, that changes as Le is deviated from unity. In the one-dimensional model
63 departed from the constant value, 62, by bending lower as the solution traveled along
the upper branch and asymptotically approached the constant 62 as the middle branch
was traversed. This behavior is also seen when Le is slightly less than unity in the two-
dimensional model, see figure 3.3. The first qualitative difference appears for the two-
dimensional model when Le is further reduced. If the Lewis number is too small 63 will
not asymptotically approach 62 but will increase past this value as the middle branch is
traversed, see figure 3.4. Although this difference does not change the physical behavior
of the system, which still exhibits stable solutions along the entire upper branch, and
unstable solutions along the middle branch with no oscillations present, one begins to see
the effect of the flame edge and the cold boundary.

As noted in the previous chapter when Le S 1 the instability manifests itself in the
form of cellular flames. The relative size of these “flamelets” is determined in part by the
critical wave number for which the “new” leading eigenvalues become zero. This is

discussed in detail in chapter 2. For the one-dimensional model it was found that Km.

43

was not greatly affected by Le or Pe. This is seen also for the two-dimensional model as
seen in figure 3.5. Here various Le and Pe configurations are presented along with the
one-dimensional result for Le=1, Pe=0, and Ta=4.0. It is clear that the two-dimensional
model produces critical wave numbers that are of the same order of magnitude as
predicted with the one-dimensional model. Figure 3.5 also indicates that the critical
wave number is only slightly affected by thermal-diffusive or hydrodynamic
considerations.

It is therefore concluded, based on the above results, that for Le S l, the one-
dimensional model can accurately predict both the flame behavior, as described by the
eigenvalue spectra, and the relative order of magnitude of the cellular flame instability.
3.4.2 Eigenvalue Results Le > 1

In the previous section it was shown that flame instabilities for Le < 1 result in
cellular flames. When Le > 1 cellular flames are not observed, rather flame oscillations
can occur that either decay to stable burning or oscillations that continue to grow
resulting in flame quenching. This behavior was laid out in Chapter 2 for the one-
dimensional model. Here it was seen that if Le is large enough the upper branch
solutions exhibited pure perturbation growth or oscillation growth or decay. This
behavior is described by the eigenvalue response examined below for the two-
dimensional model. Note that in Chapter 2 the analysis was restricted to low Pe because
the eigenvalue behavior became more complicated and wasn’t the focus of the work. In
this chapter, the restriction for low Pe is used because of the coarse nature of the grid

needed to perform these calculations.

Small deviations from Le = 1 are first studied and are Shown in figure 3.6. Small
deviations from unity Le, as well as all the cases presented below, were discussed in
detail in chapter 2 (see figure 2.4) and will not be expounded here again. If one compares
the nature of eigenvalue behavior for figure 3.6 to that of figure 2.4, many similarities
arise. Primarily, the entire upper branch is stable with decaying oscillations present over
a small portion. Note figure 3.6 does not extend as far as figure 2.4 along the upper
branch, but the behavior is similar. That is, the oscillations will persist until the complex
conjugate pair 612 crosses 63 and 63 becomes the leading eigenvalue.

There exists a critical value of Le in the one-dimensional model for which larger
values produce upper branch instabilities. This is the case also for the edge-flame model.
Figure 3.7 shows the upper branch results for Le = 2.0 and as can be seen the complex
conjugate pair extends into the positive region and then breaks up. This behavior
indicates that not only do unstable oscillations exist but also pure perturbation growth.
This is the same behavior seen in figure 2.5. Both figure 3.6 and 3.7 represent eigenvalue
responses for the case of zero convection. It was seen that the characteristics of the flame
behavior didn’t change appreciably with changes in Pe when Le < 1. To identify any
effect of Pe for Le >1 the eigenvalue behavior was obtained for Le = 1.2 with Pe=0.0 and
Pe=0.5. These results appear in figure 3.8. It is seen that the changes in Pe do not alter
the behavior of the eigenvalues, with the exception that the curve is shifted towards
higher D. This is consistent with the results seen for Le < 1.

3.5 Conclusions
A two-dimensional model of an edge-flame was examined to identify the validity

of a simpler one-dimensional model. This work examined both the sub-unity Lewis

45

number regime whose instability is exhibited as cellular flame behavior and the Le > 1
regime which has the characteristic of producing oscillatory behavior. It was
demonstrated that both of these regimes behave similarly to the one-dimensional model.
Although a more exacting analysis, that is a more refined mesh, will be needed to more
accurately predict the eigenvalue response. This will have to wait until more powerful
computers are available. However, this analysis should provide strong support to all the
previous, and future, work done with one-dimensional models. It is clear the basic
characteristics and qualitative behavior of flames is predicted with only the barest of
terms and dimensions.

3.6 References

1. Mukunda, HS. and Drummond, J .P, Applied Scientific Research, 51, 687-711, (1993)

2. Lingens, A. et al.,, Twenty-sixth Symposium (International) on Combustion, 1052-
1061, (1996)

3. Thatcher, R.W. and Dold, J.W, , Combustion Theory and Modelling, 4, 435-458,
(2000)

4. Buckmaster, J. and Zhang, Y., Conference proceedings First Joint Meeting of the
US. Sections of the Combustion Institute, 677-680, (1999)

5. Buckmaster, J ., Journal of Engineering Mathematics , 31, 269-284, (1997)

46

 

 

 

 

 

 

 

 

' d( )/dy =0
Y = 3 All?
_...
g; I Flame zone
T: O ‘2; T: o
yF= yF=1
Edge-Flame
——-—>
0
:-I _ _ — X =1
x _To yo ‘ O yF ‘ 0

Figure 3.1 Schematic of 2-D edge-flame model with boundary conditions

47

 

 

 

 

 

0 5 G i
1
o I
-O.5 oooooooooooo q
A eeeeeeeeeeeeee
b ooooooooooooooooo
v .1 OOOOOOOOOOOOOOOOOO 4
g oooooooooooooooooooooooooo
-1,5 ’ I
-2 r I
-2.5 ” G q
2,3
-3 l l I L r I r i

 

5 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6

D x 106

Figure 3.2 Eigenvalue distribution for Le=1 Pe=0.5 Ta=4 ¢=l.0. Similar
distributions are seen for the 1-D model when Le=1.

48

 

 

 

 

 

-4~ O' -

3
4.5 l r i 1 1
5.285 5.29 5.295 5.3 5.305 5.31 5.315

D x106

Figure 3.3 Eigenvalue distribution for Le=0.9 Pe=0.3 Ta=4 ¢=l.0. 0'1 turns
positive at the upper turning point, indicating instability along the middle
branch. The third eigenvalue approaches 62 asymptotically along the middle
branch.

 

 

 

49

 

«L5
«L7 —
«L9 —
411 ~
-513 - 4~—_’////////////"“33
415 -

-5.7 ~

Re(6)

-5.9 -

 

 

 

-6J . . .
5877000 5882000 5887000 5892000 5897000

D

Figure 3.4 Le=0.5 Pe=0.5 Ta=4.0 ¢=l.0. When Le is decreased below a certain
value, 63 crosses 62. This behavior is not seen in the one-dimensional model.
The Leading eigenvalue is not shown.

50

 

 
 

‘ X 2-D Le=1 Pe=0

 
 

 

 

 

g 5 o o o <
0'8 938862515555” ° o2-0Le=1Pe=o.5
£50 6 0 2-D Le=0.9 Pe=0
A 2-0 Le=0.5 Pe=0.5

0.4 - —1-D Le=1 Pe=0
0.2 l

o I I I i I

1.0E+3 6.0E+3 1.1E+4 1.6E+4 2.1E+4 2.6E+4

0.0,

Figure 3.5 Critical wave numbers for several different configurations. It
is seen that 1) Kc," is not greatly affected by changes in Le or Pe and 2)
The order of magnitude of Kc... is the same between the l-D and 2-D
models.

51

1.5

 

 

 

 

 

 

I -
I. .’ G .4
L“ 2
G -
3
5 5.05 5.1 5.15 5.2 5.25 5.3 5.35 5.4 5.45
6
D x10

Figure 3.6 Le=1.l Ta=4.0 Pe=0.0 ¢=1.0. Stable oscillations exist on the upper
branch while the entire middle branch is unstable. The regions within the
circle represent eigenvalues along the middle branch.

52

10 i

 

 

 

 

 

I I
8— i .
i
I
1
er ; -
‘O'
1
. 1
A |
I
b 4.- ‘ d
v “
Q} .
l
M x
2- ~ -
\
\
"A
If \
’0' ~
0»- l 2 \\\ .4
,2_ 63 ,
1 4 m l r I L

 

Figure 3.7 Le=2.0 Pe=0.0 Ta=4.0 ¢=l.0. The any branch eigenvalues
indicate regions of pure perturbation growth as well as oscillation growth

leading to instability. The oscillations cease to exist when 63 becomes the
leading eigenvalue.

53

 

E;
"U
t‘D
II

P
In
I

 

 

 

 

 

 

_2.5 I l i l l
5 5.2 5.4 5.6 5.8 6 6.2

D x 106

Figure 3.8 Le=l.2 Ta=4.0 ¢=l.0. The three leading eigenvalues for different
Pe. Changes in Pe have only a small effect on the flame behavior which is
consistent with Le < 1.

54

ONE-DIMENSIONAL LAMINAR PREMIXED FLAME
ANNIHILATION

4.1 Chapter Nomenclature

Greek

(1

Pre-exponential factor

Specific heat

Mass diffusivity of species i

Activation energy

Excess enthalpy, H = t + y —1

Integral; see equation (4.19)

Characteristic length

Lewis number, Le = NpCpD = org/D]:u

Reaction order

Chemical heat release

Non-dimensional reaction rate, R = Qy"exp{—B(l-t)/[1-or(l-t)]}
Time

Flow velocity

Flow direction (physical spatial coordinate)

Reduced mass fraction of limiting reactant, y = Yp/Yt:u

Mass fraction of limiting reactant

X
Mass coordinate, z = jpdx

—OO

Enthalpy ratio, or = 1 - Tu/Tb (also (1.. = hu/puCp)

55

Zeldovich number, B = or(E/RuTb)

8 Deviation of Le from unity, 8 = Ie — 1
1] Stretched mass fraction, 11 = By

1 Thermal Conductivity

v Stoichiometric coefficient

p Density

6 Non-dimensional time, 6 = t/tu

t Non-dimensional temperature, 1: = (T-Tu)/(Tf-Tu)
0) Reaction rate (units mass/vol-s)

Q Damkohler number

E Transformed spatial mass coordinate
C Normalized mass coordinate, C = z/zu
Subscripts

b Burned

f Flame or reaction zone

F Fuel — limiting reactant

i Species

m Mass

t Thermal

u Unbumed

Superscripts

a:

At location of extremal H, i.e., where dH/dC = 0

56

4.2 Introduction

Recent advances in the speed and accuracy of computational methods have made
possible the numerical simulation of detailed flow and combustion problems involving
interacting flames, complicated chemistry, turbulence, and, real property variations. A
particularly difficult feature of highly turbulent flames is the destruction and subsequent
healing of flame fronts; a fraction of the flame surface per unit volume is, in all
likelihood, continually changing through this complicated process. And though some
aspects of flame destruction have been examined, like destruction through stretching, for
instance [1,2], the difficult problem of mutual flame annihilation is virtually untouched.

Our goal in this study is to examine a simplified model for the mutual annihilation
of two identical, oppositely propagating premixed laminar flames. Our intent is to
develop a sound understanding so that the appearance of such flames in larger, complex
calculations as a sub-problem does not cause unnecessary hardship in their interpretation
and eventual use.

The direct inspiration for this study is the previous work of Echekki et al. [3] on
mutual premixed flame annihilation. This work reviews the relevant, and quite scant,
previous research and provides a rather extensive DNS examination with detailed
chemistry. One of the main results was a rather clear demonstration that the dominant
simple free flame balance between reaction and diffusion is completely upended during
annihilation. Also, it was found that the light, fast hydrogen containing species had
begun a mass diffusivity- based interaction before the thermal interaction took place: in
other words, the Lewis number of the reactant was important. In particular, they

determined that the flame acceleration during the final stage of the interaction is largely

57

produced by the diffusion of the hydrogen molecule (H2) into the reaction zone. This
alters the balance between reaction and diffusion in the reaction zone, and even for a
stoichiometric CH4-air flame with a Lewis number of unity based on the deficient
reactant, accelerations to as great as a factor of 20 over the unperturbed flame speeds
were obtained. We shall see that for Lewis numbers very close to unity these sorts of
accelerations do not occur in our simple model; in order to obtain such behaviors a model
employing at least two-step chemistry is needed, with one species being light and highly
diffusive. Another study with great relevance to our work is the recent study of Chen and
Sohrab [4] on substantially the same problem. Their study, which is purely numerical,
focuses on fundamental features of the interaction, using a chemical kinetic model that
simulates real flame situations (asymptotically reduced four step methane chemistry).
Influences of Lewis number were studied and flame accelerations and locally
homogeneous explosions were observed. Accelerations to flame speeds seven times
greater than free flame values were obtained. The general behaviors we shall describe
here were explained on the basis of numerical, not analytical, calculations.

There has also been much research on the related problem of premixed flame
destruction by cold walls [5-8]. This research is driven mostly by the study of
combustion and pollutant formulation in engine cylinders, wherein the quenching layer
and the various crevices are presumed to generate the bulk of the unburned hydrocarbon
emissions. The early work was simplified, qualitative and inconclusive because the
models were not examined thoroughly [5 ,6]. Subsequent work approached the opposite
extreme: every possible ingredient was include in the numerical simulations, like real

chemistry, thermal properties, etc.[7]. Still, these models, by their complexity, were

58

exceedingly unwieldy and the problem remained unresolved: the origins of the
hydrocarbons from the quenching layer still could not be explained. Later work [8]
emphasized the thorough examination of simple models. These have the virtue of
producing describability and understanding without interjecting cluttered detail.
Nevertheless, some kind of realistic chemistry is needed to make further advances with
this approach.

The specific aim of this study is to understand and describe the details of the
interaction and ultimate annihilation of two oppositely directed planar, laminar premixed
flames. We intend both to shed light on the extinction and to produce correlations that
can describe the dominant features of the extinction. We limit ourselves to simple one
step chemistry and simple property variations (constant pk, pzDi, etc.). A good
description of the equations and methodology of our work is given in the article of Peters
[9], where he derives and explains at length the model we shall use here.

In this article we first derive the conservation equations and discuss their
limitations (sec 2.). In sec. 3 we describe the various regimes of burning, with special
emphasis on the final and most dramatic regime, the “strong interaction” stage or
“homogenous explosion” stage [4]. In sec. 4 we present our numerical results. Finally, in
sec. 5 we examine the results of sec. 4 in light of the theoretical work of sec. 3. Possible
extensions and comparisons with the more detailed models of previous investigations are
made in sec. 5.

4.3 The Model
We consider a perfect gas undergoing a one-step irreversible Arrhenius reaction at

the flame front. The limiting component is denoted as Yp, in which the reaction is of

59

arbitrary order, n. The transient one-dimensional equations for conservation of mass,

limiting species, and energy are respectively

p. +(Pu)x =0. (4.1a)
pcprrt + uTx) = (213,), + gm, (4. lb)
p(YFl + UYFX ) = (pDFYFX )X — VFO) . (4.1C)

We assume that pDz, p71. and Cp are constant and we transform to a mass-based

X
coordinate, z = I pdx where the symbol -oo denotes a distant, but not necessarily

infinitely far, boundary. We define the non-dimensional variables C=zlzu (zu=p,,L.,), t=(T-
Tu)/(Tb-Tu), y=YplYpu, 6=t/tu (t=Lu2/or.,), and we put VF = 1 without loss of generality, and
define the flame temperature as Tb=Tu + quu/Cp. We note that the characteristic length
L0 is carried over from the steady flow propagation problem in which 3(0)/ at = 0 in the
above equations, although an important clarification must be made. In the steady
problem the characteristic length contains the eigenvalue mass burning rate m=puSL
which must be derived as part of the solution. In our case we are concerned only with the
unsteady evolution and an explicit eigenvalue calculation is not necessary. The non-

dimensional equations reduce to

2
i=3—3+ (4.2a)
ao agz
3y lazy
_=___R 4.2b
ao Leag2 ( )

where

R =Qy"exp[—B(l—T)/{l—a(l—t)}] (4.3a)

60

-E

Q = AYFueR_Tb (4.3b)
B = (E/RuTb)or (4.30)
or = 1-Tu/Tb (43d)
Le = (Xu/Dfu (4.3e)

We note from equations (2) that the convective term has disappeared in the mass-
coordinate formulation.
The initial and boundary conditions for equations (2) are now specified. At the

symmetry line, §=0, the gradients of species and temperature vanish,

8y 8‘:

— =— =0 at =0 4.4
In the far field the chemical reaction has run to undisturbed completion with zero
chemical potential. Hence,

y=0, t=1 at 5: co. (4.5)

Of course, in our numerical solution, we cannot impose conditions (4.5) at C: co; we must
choose a large downstream value instead. Our initial condition amounts to specifying the
upstream distribution of species and temperature at the instant numerical integration
begins. These distributions are derived in the following section. We note that the
upstream distribution is somewhat idealized, since it obtained from the solution of the
steady problem without the reaction term. This “preheat zone” solution gives y=0,1:=l at
the flame, resulting in a discontinuous gradient across the flame, or better, flame sheet. A
smoothing of this postulated distribution occurs in the first few instants of the numerical

integration so the y and 1: profiles soon resemble their undisturbed free stream values.

61

We observe that most of our discussion will center on the first order reaction, n=l.

For this case, the asymptotic estimate for the steady-state eigenvalue is [9],

 

2 —
Q:_[i__[1+2(3.ot 2.344+Le)

—2
2L6 B +0(B )]. (4.6a)

We shall assign to Q a magnitude consistent with this estimate. The flame velocity is
proportional to (I’m, giving

l/2

{ 2(30t-2.344+Le)]“2
1+ 13

 

 

suggesting that the flame speed increases with Le approximately as Le”2 for Le of the
order of unity and large B. We shall use this correlation to test our numerical integrations
in section 4.5.
4.4 Theoretical Discussion1

In this section we shall discuss the theoretical basis for some of the behaviors we
observe in the numerical solutions of section 4.5. In section 4.4.1, we discuss the
influence of Le ¢1on the expected flame behavior. In section 4.4.2, we discuss the
expected flame behavior during annihilation.
4.4.1 The Le =1 Asymmetry

We begin our discussion of the influences of Le =1, which we shall show are
profound, by examining the combination H=t+y-1, called the excess enthalpy function,
through the combination of equations (2). The physical significance of H is that it

essentially defines the global thermal plus chemical enthalpy of the gas, normalized so

 

I The majority of the theoretical derivations in this section is the work of the author’s major professor Dr.
I.S. Wichman, and is included in this document for sake of completeness.

62

that an unchanged field produces H=0 everywhere. From the equations and boundary
conditions given below, a number of conservation properties can be easily derived (e.g.,
integrate over the spatial coordinate Q from C=O to C: «>0 , for instance). We write the

equations and boundary conditions as

Ho = HQC +(fi—l)y§€ , (4.7a)
Hg = 0 at C= 0 (4.7b)
H=0 at C: 00 (4.7c)
H = 0 initially . (4.7d)

The initial condition (4.7d) was obtained by requiring y=l, 1:=0 everywhere in the gas
prior to igniting the flame. We note that in equation (4.7a) the variables H and y are
obviously interdependent; therefore, the third term cannot be considered as a source or
sink term in the standard sense. However, if the third term can be shown to take a certain
Sign, then it can indeed be considered as a source or sink term and the solutions for H can
be correspondingly interpreted.

When Le=1, the source term vanishes in the differential equation, so that H=0 is
the solution everywhere for all time: there is no excess or defect of enthalpy. When
Le at 1, however, the influences of reactant diffusion produce an upset of the H field.

This upset is modified by the pre-factor Le’l-l. When we let Le=1-+8, where 5 can be

positive, zero, or negative we obtain for small 6,

Le<l: Li— 1 = |5|(1+|5|) + 0(|5|3), (4.8a)
e

Le>l: fi- — 1 = —|0|(1—|5|)+ O(|6|3) , (4.8b)

63

indicating that the Le<1 case is more sensitive to a fixed change of magnitude I6] than the
Le>l case. In the extreme case Le—>0, the amplification Le'l-l becomes infinite,
whereas on the other side the limit Le —-> 00 produces merely Le'l-1—>-l. We conclude
that there is a great asymmetry between Le<l and Le>1, the former being a greatly
amplified influence, particularly as Le approaches zero.

To complicate matters further, we observe that y;; in equation (4.7a) changes Sign
at the inflection point. Upstream, ygg is negative, downstream it is positive. Hence, for
Le<1, the source term is negative upstream of the inflection point, and positive
downstream of the inflection point: we expect H<0 in front of the flame. When Le>l, we
correspondingly expect a region of enthalpy excess H>0 to precede the flame front.

We shall understand the unsteady annihilation process much more thoroughly
when we examine the enthalpy function for the steady flame. The governing equation is
equation (4.7a) with the transient term replaced by Hg. Integration subject to upstream
bounded constraints and vanishing downstream gradients gives, after a straightforward

calculation,

_ I ”Y(S) -(s-i;)
H — 1-—— l- — d 4.9
(C) ( I )y(C)[ iy(C)e 8] ( )

which carries the sign of l-Le’l since the factor in square brackets is always greater than
or equal to zero. For Le<1, HS 0 everywhere; for Le>1, H20 everywhere. As the
simplest possible test case we specify y(§) as the flame sheet solution given by

y =1—exp(Le§) upstream of C=0 and y=0 downstream of §=0. Substitution into

equation (4.9) yields

H(§)=eC—eLeC, §<0 (4.10a)

H(§)=0, §>0 (4.10b)

The upstream vicinity of the flame sheet is a region of enthalpy defect when Le<1,

enthalpy excess when Le>1, whereas the downstream side of the flame is a zero-excess

enthalpy zone throughout. The gradient of H is discontinuous across the origin, with

values l-Le at §=0' and zero at (50+. This is clearly unphysical and should be resolved

easily by a standard asymptotic analysis of the flame zone. Nevertheless, we expect that

the upstream solutions for H will be extremely accurate and, for this reason, we press
ahead with our examination of the steady flame sheet.

Function H attains a local upstream extremum at the location C* where dH/dC

vanishes. From equation (4.10) we find

1

 

*=— lnLe 4.11
C Le-l ( ) ( )
When Le is close to unity we let Le=1+6 to find
Le<1: §*=— iii (4.12a)
|5|
Le>l: §*=-1+—2— (4.12b)

Consequently, when Le<1 the upstream influence of non-zero H extends further, for a
given increment I5] , than when Le>1. A plot of -§* versus Le is given in figure 4.1. We
obtain the extremal of H(§) from equation (4.10) by substituting for C from euqation

(4.11), viz.,

Le
-(——)
H*=Le Le-l X(Le—1) (4.13)

65

The curve of H* versus Le is also shown in figure 4.1. We note that H*(Le)=-H*(l/Le),
so that H* is anti-symmetric about Le=1 under the transformation Le-—>Le". In the
vicinity of Le=1 we let Le=1+8 in equation(4. 13) to obtain.

H*(1+8) z 6(1+5)“(1+8)‘“8 (4.14)

As 8 becomes small the factor (l+8)'”8 approaches e‘l. Hence,

Le<1: H*(1—|8|)=:|—5i(l+A:—l) (4.15a)
e

Le>l: wanting-(pg) (4.15b)
e

For a given lo] , the trough of H for values of Le smaller than 1 is deeper than the hill of H
is high when Le is larger than 1.

The implications of this simplified steady flame analysis are that both the
streamwise extent of the enthalpy change, and the magnitude of the change itself will be
greater for an incremental decrease of Le to 1-8 than for an incremental increase to 1+8.
For our flame annihilation problem, the vicinity of the symmetry plane should be much
better informed of the arrival of a Le<1 flame than a Le>1 flame. For this reason, we
expect gradual and undramatic extinctions when Le<1, but sharp, well defined and
dramatic extinctions when Le>1.

4.4.2 Flame Annihilation Regimes

There are essentially three regimes of flame interaction, or non-interaction as the
case may be. One of these is the strong interaction when two premixed flames (PF) are
actively annihilating one another. A second is the weak interaction where thermo-

diffusive contact has just recently begun. The third, and final, regime is the free flame

66

stage wherein the interaction is negligible and the two PFs behave as though their
separation were infinite.

We shall analyze these separate regimes by transforming equations (4.2) into a
more suitable coordinate system, namely one in which the peak of the reactivity R is

always located at a fixed position, 5:1. The transformation is I';=§/Cf, which yields

 

1:6 — (gig-)1: = ('33)ng 'I' R (4.1621)
Cf Ct
Cr: 1
yo Cr Ye Leg? Yer;

where {I = dCf/do is the propagation velocity, which is negative because Cf decreases as

time increases. The reactivity peak, Rm”, is located at §=1. As the annihilation process
proceeds from start to finish, the gas between the symmetry plane at =0 and the
reactivity maximum at I;=1 undergoes great changes. At first, in the free flame stage, the
temperature and species mass fraction are essentially undisturbed because R occupies an
infinitesimal volume fraction of the upstream gas. Here we use 12:0, y=1 as the respective
distributions. In the weak interaction stage, the flame is close enough to the symmetry
plane that thermodiffusive communication begins to alter the undisturbed gas. In the
final stage of strong interaction, a significant segment of the reactivity has entered the
upstream gas; here the reaction term in equations (4.16) can no longer be ignored.

We begin by analyzing the strong interaction. We postulate that the spatial
distribution of R has substantially eliminated the gradients t§,y§ and their derivatives
T§§,y§§ leaving the lowest order balance between transient and reaction terms, as in

homogeneous thermal explosions. This is reasonable because the upstream preheating,

67

confined to a finite pre-flame zone, eventually heats the gas to a nearly uniform condition,
at least when the Lewis number is larger than one and the thermal layer is thicker than the
mass diffusion layer. In the latter case of Le<1 the two flames communicate via mass
diffusivity which, in effect, chokes the approaching flames by rapidly depriving them of

the needed reactant. In either case, the relevant equations are to = R , ya = -R , which
can be combined to give (1' + y)a = 0. This integrates to t+y=1+H=constant, where the
number H may be positive (Le>1 or excess enthalpy case) or negative (Le<1 or enthalpy

defect case). We substitute t=1+H—y into the equation yo=-R to obtain, in the limit

,6 —) oo , the equation (recall we use n=l presently)
yo = -r‘2ye‘By => no = 42116“ = —r‘2R(n), (4.17)
where S2 = Qexp(Hfl) and R(n)=nexp(-n) is the lowest order asymptotic representation

of the reaction rate, with n=By being the stretched species mass fraction. In the
isenthalpic case we have H=0. Qualitatively, the difference between the isenthalpic case

and the non-isenthalpic cases can be explained from the ratio

(n6)isentropic = -Qne'" ’C—HB

. _ — (4.18)
(n6)non—isentropic -Q1‘|e ‘1

 

When His positive, the rate of change of n, in the non-isentropic case, is greater; when H

is negative, the isentropic rate of change is greater. Equation (15) integrates to

.. Tlo .,
o=o“. i R(x)"dx arr-11mm) (4.19)
Ti

where no is the initial value of n. Integral I is always positive. Consequently equation

(4.19) predicts a monotonic decay of n that is hastened when Oincreases. This, of

68

course, is the excess enthalpy case, H>0. When Odecreases, the decay of y is slowed.
We observe that these changes in the rate of decay of T] are extremely sensitive to small
changes in the non-isenthalpicity, H, since E2 is proportional to exp(HB) and B is large.
Once again, small changes in the excess enthalpy easily produce O(1) changes in the
reactant consumption rate.

In the weak interaction stage the reaction term R is essentially negligible in the

vicinity of the symmetry plane §<<O(1), whereby

 

To - (SC—Ur: =(i2)13§§ (4.208)
Cr Ct

_ $51 = 1 420b

yo ( f )Yg (11%,?ng ( - )

When properly normalized, i.e. with 6~O(§22§2), the transient and diffusive terms are of

the order unity and the convective term is of the order {I {f 52, which for small E is

negligible. We are left with

g(6)‘to = tfig (4.21a)
~ 1 421b
g(G)yo ~EY§§ ( . )

where g(6)=Cf2. The initial and boundary conditions for equations (4.21) are zero
gradients at i=0, burned reactant at §=oo , and given initial distributions when 6:0. The
decrease of g(6) with increasing 6 will raise the effective thermal diffusivity g(6)’l and
the effective mass diffusivity (Leg(6))", thereby accelerating the rate of increase of t and

the decrease of y near the centerline. Depending on whether Le is greater than or smaller

69

than unity the change in one will outweigh the other: when Le<1 the y field should begin
its’ change before the 1 field; vice-versa for Le>1.

We note that a discussion based on length scales would determine only that the
strong and weak interactions occur over the same spatial location. In the former case,
spatial gradients of t and y had become vanishingly small; in the latter case, the reaction
R was negligibly small. The spatial scales, therefore, are the same as for the ordinary
premixed flames: the standard conduction length based on flame speed; and the reaction-
zone scale, which is the former divided by B. We observe that in this problem, as
formulated, there is no characteristic physical length scale. Hence, the preceding length

scales are of a secondary, or derivative, kind.

In the free stream stage we examine equations (4.16). We begin by integrating from §=0

to §=oo to find

7844+ (13.44

St =‘Cr 0 (4.22)

 

Iydé
0

Since the flame is, in real terms, infinitely distant from the symmetry axis at §=0, we
make negligible error by asserting that R=0 between §=0 and §=1. The integration of

equations (4.16), with R=0, between §=0 and §=1 leads directly to

1 . l
CigYodéi’CrCigydfi: y§(l) (423)

which, by comparison to equation (4.22) yields y§(1)= —§% ]Rd§.
o

70

Thus, equation (4.22) with the upper limits changed from 00 to 1 becomes

l
I—yg(1)+ci(i)yad&l

 

C, = - (4.24)

l
Criydg
0
1
But Iyd§=l and y;(1)=xf dy/dx ; also since y=l—exp(x—xf) in the upstream
0

region, we obtain finally -y§(1)=x{. These results, coupled with the observation

1
[yodfi = 0 , yield from equation (4.24) the simple result
0

C, =—1 (4.25)
Since fl = (tu/zu)dz/dt = (tup(Xf)/Zu)de/dt = -(p(xf)/pu)(tu/Lu)SL, this constraint requires

that the non-dimensionalization values tn, L,l be chosen so that their ratio is given by
t./I.=p(xi)st/p. .

There are a sufficient number of degrees of freedom to allow this condition to be
fulfilled. Hence with a suitable normalization of the characteristic length and time scales,
the flame front propagates with unit negatively directed velocity. This undisturbed
propagation is altered once the two PFs begin their thermo-diffusive interaction.

One important observation has not yet been mentioned, namely the absence in our
problem of a single characteristic physical length. Our flame is laminar, which alleviates
the need for Taylor, Kolmogrov, etc. length scales. Indeed, the only length scales are
derived scales, like the characteristic thermal diffusion length (L), mass diffusion length
(14,.) and the reaction zone thickness (14‘). The former two are combined in the Lewis

number, le=L,/I.,,.; the latter appears only when the reaction zone is subjected to a careful

71

asymptotic analysis, where we find Lr-O(Lt/B). Our analysis suggests the two PFS are
aware of one another when their separation is O(L,) or 0%), depending on whether Le is
greater than or less than 1. Subsequently, the extinction is determined by the amount 8 of

excess enthalpy in the region between the flames. From Eq.(16) we obtain

o = 5‘1" 1n(19-) + (no -11) +—1——(n§ —n2) +... z 0“1n(9£)+... (4.26)
11 2-2! 1)
since
xes x2
J—ds=ln(x)+x+—+.... (4.27)
0 s 2-2!

When n/n,~0(10), we see that o is still 0(62'l ), indicating that the time to extinction is
governed by parameter S2 , which is proportional to exp(HB). For large B, the time to
ignition can decrease exponentially, for positive H. We suppose that the rate of change of
H is much smaller than the extinction rate, hence H is expected to be nearly constant.
When Le>1 and H is positive, the elapsed time between the substantially undisturbed

initial state and complete extinction is

o z 0(0“) = wig-96”“) (4.28)

which can become very small as B increases. The constant of proportionality should be
understood as a characteristic magnitude of the square brackets in equation (4.26); this
requires choosing the initial state as represented by the rescaled mass fraction 110 and the
initial to final state ratio, no/n. Our subsequent numerical examination will discuss the

choice of these “parameters”.

72

The influence of non-unity reaction order is easily assessed. When mt 1, the

equivalents of Eq.(15) in the strong interaction stage became ya = —Oy"e’”’ and
770 = —(O,B"")7]"e‘”. Thus, (2 in our explosion equations (4.26) and (4.28) is replaced
by (2 "". We obtain, for 6, equation (4.26) with B2 replaced by B“. Hence when n<l,

the extinction is made more abrupt by the algebraic factor B’“'"’, and when n>1 the
extinction is prolonged by the algebraic factor B“.

Some additional speculations on the nature of the annihilation sequence are
possible, especially when Le>1. Here we observe that prior to the strong interaction the
spatial distribution of R resembles a strongly peaked bell curve. As the centerline is
approached, the upstream segment of the bell begins to rise. During the weak interaction,
the qualitative nature of the bell remains unchanged: one maximum, two local rninima,
and two inflection points. During the strong interaction, a moment arrives when the
reactivity upstream of the previous maximum (at the top of the bell) exceeds the bell-top
maximum, meaning the upstream inflection point is lost. When this occurs, a distinct
qualitative change has occurred and mutual annihilation is now assured.

4.5 Numerical Results

The primary emphasis of the numerical work is directed towards discerning the Lewis
number effects on the flame propagation and annihilation. Cases of Lewis numbers
greater and less than unity, as well as equal to unity, were studied to find the effects on
the reaction term profiles and how it relates to the diffusive and convective terms, the

temperature and species profiles, and the excess enthalpy. Furthermore, relationships

73

between the numerical flame speed, as well as the annihilation time, and the Lewis
number were derived and compared to theoretical results.
4.5.1 Numerical Technique

Given the non-linear nature of equations (4.2 and 4.3) and the coupling between
the energy and species equations through the reaction term, the solution of these
equations is very critical. The numerical technique used was a finite difference model
which incorporated a second order correct central differencing scheme. The non-linear
reaction term was handled by performing a quasi-linearization around the previous
iteration. This involved iterating at each time step until the new values of t and y
converged with respect to an arbitrary small positive value. This allowed larger time
steps to be taken without sacrificing stability or numerical accuracy. Furthermore, the
governing equations were solved simultaneously using Gaussian elimination and back
substitution. To test the code the elemental lengths was reduced by one half with
negligible change in the results. Hence, we are confident that the solution method used
produces accurate results.
4.5.2 The Free Stream

To examine the free stream propagation, the three terms of equation (4.2a) were
plotted along with the excess enthalpy for three cases of Lewis number equal to 1, 2, and
0.8. These plots appear as figures 4.2, 4.3, and 4.4. As can be seen the free stream term
balances look very similar. For this reason it is concluded that, although the magnitudes
of the terms are different, the primary mechanisms for flame propagation, i.e. a reactive-
diffusive balance, are unaffected by the Lewis number. The excess enthalpy, however,

has a different profile for all three cases. After the flame, the excess enthalpy is zero for

74

all three cases, but does not follow the same trend before the flame. For a Lewis number
of unity H remains zero before the flame, but for a Lewis number greater than 1 the
excess enthalpy becomes positive, whereas for Lewis less than 1, H is negative. These
results support the theoretical conclusions drawn in section 4.4. We also note that the
premixed flame enthalpy defect region is larger for Le<1 than the enthalpy excess region
is for Le>1, in agreement with predictions.

To study the effect of diffusion and pre-heating, i.e. Lewis number effects, species
and temperature profiles were plotted along with the reaction term for Lewis numbers of
2 and 0.8. These appear as figures 4.5 and 4.6. Comparing figures 4.5 and 4.6, it is
evident that the species profiles are very different before the flame. For Lewis number of
2 there is very little diffusion towards the flame, whereas for Lewis number of 0.8, there
is a significant reduction of species upstream of the flame, indicating that the upstream
flame region is aware of the approaching flame front. This supports the results that the
excess enthalpy for Lewis number greater than 1 is positive and for Lewis number less
than one H is negative.

The final area of interest in the free stream is the laminar flame speed. As
discussed in section 4.4, the flame propagation speed should have a nominally square root
power dependence on the Lewis number (see the final equation of section 4.3). To
examine how well the numerical flame speed corresponded to this, the various velocities
were plotted against Le and then curve fitted. Figure 4.7 contains the plot and equation

(4.29) shows the functional relationship,

SL = 0.8926Leo'5025 (4.29)

75

As can be seen, the free steam flame speed has the functional relationship expected for
free stream flow. It should be noted that the theoretical relationship is based on
propagation an infinite distance away from any boundary. This condition is impossible to
obtain in a numerical code, however, it was assumed that a non-dimensional distance of
10 would provide an adequate approximation. Furthermore, the functional relationship
becomes stronger as the flames become closer; more will be discussed on this in the next
section dealing with the annihilation process.
4.5.3 Annihilation Zone

The examination of the “strong interaction” zone is the primary concern of this
study. Examining the Lewis number effects on the reaction profiles, the temperature
profiles; and correlating the annihilation time of the two flames is the aim of this section.

To begin the analysis of the reaction profiles through extinction, we will look at
the Le=1 case. From this we will get a sense of some of the characteristics associated
with flame annihilation. Figure 4.8 illustrates the reaction term profile as the two flames
approach the symmetry line. As can be seen, the interaction of the two flames does not
occur suddenly for this case, i.e. there is a small, low-grade reaction ahead of the flame,
which is caused by the thermal upstream diffusion. For this reason, to get true free stream
propagation it is necessary to start very far from the symmetry line. Further examination
of figure 4.8 shows that there is an acceleration into the flame junction; this is indeed
expected for the Le.>_1 case. Furthermore, it is important to note that the most prominent
parts of the flames have not yet come into contact. This occurs when the concavity of the
reaction profile changes from positive to negative. After having come into contact we

0

surmise that extinction is not very far away. It is this time, from first joining to when the

76

reaction term is essentially zero (0(10'2)), that we consider as the numerical annihilation
time. Obviously, the theoretical time cannot be so nicely delineated hence our approach
in comparing theory and numerics will be to scrutinize the functional dependencies.
Reaction profiles for Lewis numbers of 2, 1.25, 0.8, and 0.5 will be compared to the Le=1
case to either support or refute the analytical work of the previous section. These
additional plots appear as figures 4.9 through 4.12.

Primary examination of the various reaction profiles shows that as the Lewis
number increases, the width of the reaction zone decreases; i.e. the reactions are much
sharper and quicker with less upstream thermal and species diffusion. Furthermore, the
acceleration of the flame for Le>1 becomes greater. For the case of Le=1 the acceleration
is almost negligible, whereas for Le=2 it is clear that the reaction rate increases
tremendously. AS mentioned in the introduction, the previous work of Echekki et al.[3]
can in fact produce accelerations when the value of Le of the deficient reactant is close to
unity. However, we have not included multiple - and potentially very light - species into
our one-step chemical model. We also observe that for Le<1 there is no acceleration, but
rather a deceleration as the flame approaches the symmetry line. If one follows the time-
labeled reaction profiles of figure 4.11, it is clear that when the two flames come into
contact there is no jump in the reaction term, but conversely a slow, inevitable extinction.
This effect is more pronounced as Le—> 0. It is for this reason that we expect the flame
speed dependence on the Lewis number to become stronger as the flame becomes aware
of its counterpart. The reaction profile in figure 4.11 shows that for Le=0.5 the reaction
essentially dies out before the flames come into contact. These results support the

theoretical assertions made in section 4.4 dealing with the Lewis number asymmetry.

77

Examination for the corresponding temperature profiles for Le=2 and 0.8 through
annihilation will further shed some light on this problem. These can be found in figures
4.13 and 4.14. i

The time-labeled temperature profiles for these two cases appear to have different
characteristics. The case of Le=2 shows a change of concavity in the temperature profile,
as discussed at the end of section 4.4.2. This is indeed the case for Lewis numbers greater
than 1. Additional analysis shows that this change in inflection is consistent with the
maximum reaction rate; that is, this profile change occurs immediately after the reaction
rate reaches its maximum. This change in the temperature profile conducts heat away
from the reaction and marks the beginning of the annihilation process. It should be noted
that the flame does not extinguish itself because of the heat loss but rather by a defect of
species through the burning process. However, this heat loss may add to the rate at which
the flames annihilate one another. For the Le=0.8 case, the temperature profiles retain
their positive concavity through annihilation. These profiles are consistent with all Lewis
numbers less than 1, although as the Lewis number decreases the change in each profile
becomes less dramatic. For the Le<1 case the thermal energy is being conducted into the
flame junction, as blood is pumped into a dying patient, in hopes to maintain the flame. It
is this additional energy input that may contribute to the slower annihilation times for
Lewis number less than 1.

The final area of interest is finding a correlation between the annihilation time and
the Lewis number. It is observed that the time to extinction from maximum reaction at
the junction has a strong dependence on the Lewis number. We have noted that the

criteria used in evaluating numerical annihilation times are not usable in theoretical

78

evaluations. Therefore, we expect some discrepancies between the theoretical and the
numerical results. Theoretical flame quenching times can be shown, using equation
(4.28), to be proportional to Le{exp[-B(te-1)/Iew‘”"’]}, where B is the Zeldovich
number. The numerical results can be found graphically in figure 4.15, and give a Lewis
number dependence, based on the theoretical result, described by equation (4.30),

cam, = 0.98Le-exp[—0.416B(Le- 1)] (4.30)

If the average value of Lek/‘1‘") is taken from Le=0.l to Le=2.0, its value is
found to be 0.404. This value is very close to the factor of 0.416 found in the correlation
of equation (4.30). We can therefore surmise that the theoretical result corresponds with
the numerical result within an average error of 16 percent. That is, the 16 percent is the
average error between the correlation and the numerical results. The numerical results
can be correlated to an average error of about three percent by using the following

correlation:
om = 1.26071" (4.31)

The more accurate correlation does not remove any of the validity from the one found in
equation (4.30). Given the good agreement with the theoretical result of equation (4.30)
and the accuracy of the correlation given by the above equation, neither is to be preferred
over the other, although equaiton (4.30) does provide further support that the numerical
results are accurate.
4.6 Comparison with experimental observations

It has been suggested that this analysis could offer insight to turbulent flame

stability. Examination of the numerical and theoretical results leads to the following

79

conclusions. For cases when Le>1 the annihilation times are quite small, that is the event
occurs rapidly with an associated increase in the excess enthalpy of the system. When
Le<1 the effect is quite the opposite, that is the annihilation occurs over a relatively long
period of time and is associated with an enthalpy defect. Based on these observations re-
ignition seems more possible when Le>1. Hence, it is postulated that turbulent flames in
which Le>1 will be more stable than those when Le<1.

Experimental observations for Methane/air turbulent premixed flames [10]
concluded that local flame extinction for cases when Le>1 is caused by excessive flame
stretching. That is, the local flame does not propagate in a planar nature but rather
becomes curved and essentially “stretched”. It was observed that when Le<1 flame
stretch and excessive heat losses play equally important roles. This is consistent with the
excess enthalpy defect observed in the present study.

Along similar lines as the above study, and experimental examination of
Propane/air mixtures [11] concludes that flame wrinkling is increased by as much as
thirty percent when Le<1 as compared to Le>1. This result can be explained by the rapid
annihilation and presumed quick re-ignition of Le>1 flames. An additional observation
made by the researchers is that for nearly freely propagating wave fronts, instabilities
were observed when Le<1.

A clear description of thermal-diffusive effects of propagating flame fronts is
given in reference 12. This numerical model of turbulent premixed flames uses a two-
step chemical kinetics scheme in which the Lewis number of the reactants is allowed to
vary. The conclusions reached about flame curvature support those seen experimentally,

though no comments of flame extinction were presented.

80

The above references do by no means represent a complete list of available data.
However, all the results support in some way the results found by the one-dimensional
laminar flame annihilation studied in this chapter.

4.7 Conclusions

We have demonstrated that a simple model of a premixed flame can describe many, if
not most, of the principal features of mutual flame annihilation. These features can be
anticipated theoretically, without requiring a numerical analysis. However, the numerical
analysis provides a seamless connection between the disconnected theoretical submodels
while also providing direct indication of the magnitudes involved. We see, consistent
with reference 3, that the standard free flame reactive-diffusive balance is entirely upset
during annihilation; we also see that Le has perhaps the strongest influence on the
qualitative character of extinction, since Le<1 extinctions are completely different from
the Le>1 extinctions. The former involve significant pre-extinction communication,
hence the extinction is gradual, inexorable, and undramatic. By contrast, the Le>1
extinctions occur with furious rapidity, partly because the two flames are largely ignorant
of their mutual existence until they very close. Our correlation of the extinction times
with the asymptotic formula predicts a dependence on exp[-B*constant*(Le-l)], which
we in fact observe numerically (see equation [4.30]). Thus, the sensitivity of extinction to
Le is closely coupled to the magnitude of the Zeldovich number. When this is small,
extinctions will be tame regardless of Le.

We can only speculate on what occurs in the post annihilation stage when this
model is included in a larger problem of the kind discussed at the outset of this article. In

this case we have many mutually interacting flames and “flamelets” that are in continual

81

relative motion due to turbulence and the other flow requirements. When the Lewis
number of the principle limiting reactant is larger than unity, we expect dramatic
annihilations accompanied by surges in the thermal enthalpy, as shown in figure 4.13.
These local pockets of excess thermal enthalpy could seemingly serve as local hot spots
for subsequent reignition of locally quenched mixtures. If the characteristic turbulent
flow time is substantially larger than the characteristic annihilation time, the subsequent
healing of the flame may occur quite rapidly; this would suggest that this flame might
remain significantly intact with only minimal permanent destruction of flame surface. If,
however, the Lewis number is small, the flames can annihilate one another while
simultaneously depressing the thermal enthalpy (see figure 4.14). In this case, the
prospects for reignition and flame healing are more remote. This is consistent with
observed experimental results. We have not conducted a detailed comparison of our
calculated annihilation times and turbulent flow time scales. Such a comparison would
be interesting indeed. In real flames, however, there are many constituents, each with its’
own characteristic Lewis number. It remains a challenge to determine the reactants
whose influences through Le, are greatest: some composite average of these Lei’s may be
used in a correlation such as equation (4.30) in order to attempt to describe these more
complicated and more realistic flames. It would be interesting to repeat our analysis with
a more realistic two- or even four-step mechanism, in order to determine, with greater
precision, the relation between annihilation time and the relevant Lewis numbers. The
previous work of Echekki et al.[3] suggests that with light, highly diffusive species like
H2, the value of Le for the deficient reactant does not control the flame behavior as

strongly as it does here. The qualitative nature of such a system should be evaluated with

82

a model 2-step reaction scheme, in which one of the intermediate species is much lighter
than the others. Our findings therefore appear to be representative only for flames in
which (1) The overall chemistry can be safely reduced to a global single step; (2) The
lewis numbers of all species are close to one another and to unity; and (3) The Lewis
number of the limiting reactant (deficient species) is used in our formulae. We should
caution the reader, however, because this study was not conducted with such a practical
aim in mind. After more study, such practical questions might be completely addressed,
but given the current paucity of hard evidence, we leave such questions to future work.

In a pedagogical sense, this simple model affords us yet another chance to
examine non-trivial premixed flame propagation. That is, we are able to examine and
appreciate a premixed flame process that does not conform to the restrictions obeyed by
the steadily propagating one-dimensional laminar flame. The latter, though centrally
important, is shown in this and other works [8] to be a singular special case. Future work
may Show that the entire concept of a “flame speed” in complicated transient flows,
though useful for a preliminary understanding, is void of meaningful content.

4.8 References

1. Candel, S., and Poinsot, T., Combust. Sci. Tech., 70: l-15 (1990).

2. Poinsot, T., Veynante, D., and Candel, S., J. Fluid Mech., 228: 561-606 (1991).

3. Echekki, T., Chen, J.H., and Gran, 1., 26m Symposium (International) on Combustion,
The Combustion Institute, 323-329 (1997).

4. Chen, CL, and Sohrab, S.H., Combust. Flame, 101: 360-370 (1995)

5. Daniel, W.A., Sixth Symposium (International) on Combustion, Reinhold, New York,

(1957).

83

6. Adamczyk, A.A., and Lavoie, G.A., SAE Transactiona 87, SAE paper # 780969

(1978).

N

Westbrook, C.K., Adamczyk, A.A., and Lavoie, G.A., Combust. Flame, 40:81-99

(1981).

8. Wichman, 1.5., and Bruneaux, G., Combust. Flame, 103: 296-310 (1995).

9. Peters, N ., Numerical Methods in Laminar Flame Propagation, a GAMM Workshop,
Braunschweig : vieweg, (1982).

10. Yahagi,Y., Ueda, T, and Mizomoto, M., 24% Symposium (International) on
Combustion, The Combustion Institute, 537-542 (1992)

11. Lee, T., North, G, and Santavicca, D, Combustion and Flame, 93, 445-456 (1993)

12. Rutland, C., and Trouve, A., Combustion and Flame, 94, 41-57 (1993)

84

0.2 0.6 A LA

4 0
.1.”
0
0

.5
C
.s
.b
.s
a:

 

 

Figure 4.1 Upstream extremum for H and location versus Le.

85

 

 

+H*
+§e

 

 

 

 

 

 

 

 

Figure 4.2 Le=2, free stream

86

 

 

 

 

 

 

 

+ R
+ dt/dt
‘ . ‘ d2t/dC2
m4— H
-1 ~-
-1.5 ‘L
C
Figure 4.3 Le=0.8, free stream

87

A

01
I
I

0.5 “

II
[Il’Ill/I

 

III A
umIII/nunIIIN’N’H’H’I "
1Io'.r”.;.‘.*.'~1:Z~'2~‘v.'€¢'4§’4'~3'e}‘3"?"’I“.’("r?:¢-I'¥3-?3755553'1‘rid-3'] 5 a; ' «it; 83:53

_0.04 0.68 1.32

-0.5'
-1 ~-

 

-1.5-L

Figure 4.4 Le=l.0, free stream

88

WW : : item‘ .' ‘. .r‘
I L

1.96 2.6 . 3.245.388

 

 

+d‘E/dt
+R

‘fi— d2‘t/dC2

 

 

4.5 "

3.5 ‘”
:3 u-
2.5 “

1.54“

0.5 “

0
0.05 0.65

 

1 .25

1.85 2.45
C

IJAIJIAIJI‘,.,""A,"
AA.
0‘.
‘a
‘3
A
A
-A
A
h?
u.a‘ D
. uh‘p'o' A \
‘ ‘ -' “,‘/DI0I'DIDIAIII‘IOI)IAII(A.

3.05 3.65 4.25 4.85

Figure 4.5 Le=2, T and Y profiles through reaction zone

89

 

 

--‘.--'FR

fi' V

 

 

 

+R

 

‘5' Y

 

 

 

 

0.05 0.75 1.45 2.15 2.85 3.55 4.25 4.95
C

Figure 4.6 Le=0.8, T and Y profiles through reaction zone

90

 

1.4 —~

0.8 —~

Flame Speed

0.4 +

0.2 ~»

 

 

Le

 

 

 

Figure 4.7 Free stream laminar flame speed as a function of Lewis number

91

2--
1.8-

1.6 _. 2-x” .‘\
-r-“’/ \.
1.4 \

1.2 4-... . \
1 4- ‘ \
R 0.8 «.-
0'6 1‘ \
0.4 “

0.2 4~

0 4 -—-——

 

0.02 0.08 0.14 0.26 0.32 0.38 0.44 0.56 0.62 0.68 0.74 0.86 0.92
C

Figure 4.8 Le=1.0 Reaction profile through annihilation

92

 

 

0.02 0.06 0.18 0.22 0.26 0.34 0.38 0.42 0.46
C

Figure 4.9 Reaction profile for Le=2.0 through annihilation

93

 

 

 

94

 

0.37 1
0.35 .
0.33 J
0.31
0.29
0.27
0.25

 

0.39 ._

 

 

0.05 -

    

—+——1
LO
0;
o

 

 

—o— R(4.92)
+ Fl(4.94)
R(4.96)
-—-><-— R(4.98)
+ R(5.0)

 

 

Figure 4.11 Reaction profiles for Le=0.8 through annihilation

95

 

1.2 -~

0.8 4r

0.4 «-

  

 

 

 

0 AW A I: _/’.r
0.05 0.35 0.95 1 .25 1.85 2.45 3.05 3.65 4.25 4.85

C

Figure 4.12 Le=0.5 reaction profiles during Flame propagation towards junction

 

96

 

 

 

0.02 0.06 0.14 0.18 0.22 0.26 0.34 0.38 0.42 0.46
C

Figure 4.13 Temperature profiles for Le=2.0 through annihilation

97

 

 

+ 1:(O.751)
—I— T(0.752)
+ 1:(0.753)
“4*“ T(0.754)
+ T(0.755)
+ t(0.756)

 

 

0.9 1
0.88 .. ,, 9v
0.86» ,,.-xv

0.84 4 1,...
0.82 .9..-

0.8 J-, ”

.m"

 

. arr” .. if.

ear" p

0.78

0.76

0.7411111111111111111411
0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95

C

 

Figure 4.14 Temperature profiles for Le=0.8 through annihilation

98

 

 

+ 1:(4.92)
+ t(4.93)

t(4.94)
+ «4.95)
+ 1:(4.96)

 

 

 

0.45 1
0.4 .- ~
0.35 ~~ \
0.3 + \
0.25 -- \
0.2 -» x
0.15 1- 0‘

t annlhllatlon

0.1 «L \
\0\

 

 

__ \
0.05 \

 

Figure 4.15 Annihilation time versus Lewis number

99

HEAT LOSS ANALYSIS OF A DIFFUSION

FLAME LEADING EDGE

5.1 Chapter Nomenclature

b
C?
D

Da

U,V

u,v

Reduced Damkohler number

Specific heat

Binary diffusion coefficient

Damkohler number

Excess enthalpy

Half channel width

Non—dimensional spatial coordinate

Peclet number

Heat flux

Non-dimensional spatial coordinate

Quenching distance, (§2+112)"2 in gm coordinates of figure 5.1
Quenching distance, (u2+v2)”2 in u,v coordinates of figure 5 .2
Scaling factor

Temperature

Non-dimensional coordinates (u/rq,v/rq)

Dimensional conformal mapping coordinates

Velocity

Reaction term

Non-dimensional reaction term

100

x,y Dimensional Coordinates

Y, Mass fraction of species i

yi Non-dimensional mass fraction of species i
Z Mixture fraction

Greek

0L Thermal diffusivity

a l-To/Tf

B Zel’dovich number, 01(E/RTf)

[3i Curve fit parameter

(1) Global stoichiometric index

1'] Non-dimensional spatial coordinate

0 Heat transfer parameter detailed in Equation(22.a)
9» Thermal conductivity

v Ratio of oxidizer mass to fuel mass

6 Non-dimensional spatial coordinate

p Density

T Non-dimensional temperature, (T 'To)/(Tf'To)
E,§ Non-dimensional spatial coordinate
Subscripts

F Fuel

f Flame

0 Oxidizer

101

0 Reference boundary value
q Quench
5.2 Introduction

Real flames interact with surfaces. These interactions may be primary, in which
case the interaction of the flame with the surface cannot be ignored in an analysis of
flame behavior, or secondary, in which case the solid surface is largely incidental to the
combustion process. The flame-surface interaction includes many physical processes:
flame chemistry; fluid flow dynamics; and conductive, convective, and radiant heat
transfer. For this reason a comprehensive theoretical model is exceedingly difficult,
perhaps impossible, to construct. On the other hand, a pure numerical study without
theoretical guidance is likely to make few novel observations. A similar statement
applies of course to experimentation, thus suggesting that a “democratic” investigation, in
which every conceivable influence is equally weighted and seriously considered will
surely fail for not making distinctions between what is important and what is
unimportant. Consequently, the development of theories for simplified limiting cases is a
necessary part of a rational examination. Although theories almost by definition involve
the suppression of some physics, this suppression ideally should have a positive, not a
negative, aspect. This is adumbrated in the following quotation: “Idealization does not
consist, as is commonly believed, in a subtracting or deducting of the petty and
secondary. A tremendous expulsion of the principal features rather is the decisive thing,
so that thereupon the others too disappear” [1]. Thus, it is the prominence of various
model (theory) features that by necessity force the remaining features into the

background. Finally, we note from practical experience that large-scale numerical

102

simulations often employ simplistic “sub-models” in order to describe processes that
cannot be fully resolved by the global theory. It is imperative that such “sub-models”
faithfully represent the “sub—problem” they purport to describe. They must also be
sufficiently integrated with ease into the larger computation.

The purpose of this study is to first examine the leading edge of a flame in terms
of the heat transfer during the quenching of a flame leading edge by an adjacent, cold,
solid surface. The conductive heat transfer from a flame to a nearby surface is important
for ensuring a continued supply of fuel, as in flame spread, or for ensuring flame survival,
as in bumer-attached flames. In the former case, the transfer of heat from the spreading
flame gasifies the solid fuel beneath it, producing the fuel vapors that feed it with
reactant. In the latter case, the heating of the burner rim by the flame minimizes transient
heat losses, thereby maximizing reactant consumption while preventing possible flame
extinction. Studies have recently been completed on various features of flame-wall
interaction, in configurations reminiscent of both flame spread [2] and burner-attached
flames [2,3].

A second feature of flame-surface interaction that we wish to investigate is the
flame behavior as predicted by the relatively simple analysis of the conserved scalar
mixture fraction combined with infinitely fast-rate chemical kinetics. As will be shown
in comparisons to numerical solutions of the non-linear, finite-rate chemistry problem
(model problem “A”), the examination of the mixture fraction field can provide valuable
insight into the flame quenching mechanism. In addition, the study of model problem

“A” provides the necessary impetus for examining the quenching problem discussed

103

previously. Nevertheless, certain analytical solutions obtained by these means shall
exhibit almost no correspondence with our detailed examination.

From a practical viewpoint, and in terms of material response, the most important
quantity to examine is the conductive heat flux, both to the surface and at various
locations in the gas. Radiative heat losses are ignored for two primary reasons: in small
scale combustion, conductive losses are dominant [4] and under micro-gravity conditions,
soot and gas radiation are generally diminished. Currently, we know very little about
such conductive fluxes, including their characteristic orders of magnitude or their
functional shapes. Although some progress has recently been made in this area [2], we
are far from a satisfactory understanding that might enable the use of such estimates in
any engineering capacity. Our analysis herein has features in common with reference 5,
where the nature of the opposed-flow flame spread over solid fuels was examined in
detail. The simplifications provided by employing infinite rate chemistry and an
idealized geometry for gas flow and flame spread enabled numerous deductions to be
made, which might otherwise have remained hidden in excessive detail. One of these
was the seemingly paradoxical result, long since confirmed by numerical examination [6]
that the streamwise conductivity through the solid fuel bed plays a minor, almost
negligible, role in the rate of flame spread. In addition to the heat flux into walls and
surfaces, it is also worth examining heat fluxes across various planes in the gas [5].

Of particular interest in our model problem is the flux directed from the flame leading
edge towards the nearby surface. As discussed in [2,7], this flux can change significantly
from very large values near the leading edge, where the flame quenches, to much smaller

values near the surface, where the gas loses thermal energy to the wall.

104

1

 

5.3 Background

In this chapter we shall examine a greatly simplified heat transfer model problem
(model “B”) in the light of a fuller (but still simplified) model “A”. The latter includes
finite-rate chemistry. It has been shown in past work that model “A” is an often
quantitatively accurate representation of the actual, physical problem [8]. This
conclusion was reached by comparing flame behavior in model “A” with a variable
property Navier-Stokes simulation.

In this section we shall restate and examine various previous results from the analysis
of model “A”. We shall also present various analytical results from this model for the
infinite-rate chemistry limit. This examination will be conducted with the intent of
providing motivation for a detailed examination of the pure heat transfer model “B” in
sections 5.4—5.6 following.

Model “A” describes qualitatively a fuel injector problem and more rudimentarily a
general lifted triple flame resembling but different in some features than a flame-spread
problem. The model configuration is shown in figure 5 .1. The fuel flows past one side of
an impermeable, perfectly conducting divider and mixes with the oxidizer stream that
flows past its other side. Permeable walls through which fuel and oxidizer separately
diffuse are aligned with the bulk flow direction, which is vertical. Depending on the
global stoichiometry, the diffusion flame inclines either to the left, to the right, or lies
directly downstream of the divider. This model produces a triple flame configuration
with a premixed flame arc anchoring the diffusion flame to the lower boundary. The
diffusion flame extends downstream from the anchor point. Figure 5.1 shows a

schematic of the triple flame configuration, including the fuel-rich and fuel-lean premixed

105

flames and the trailing downstream diffusion flame arc. The system is assumed to be
steady. All boundaries except the divider plate are porous, non-reactive walls, whose
temperature is To. The velocities, v, of the two streams are identical and specified at the
inlet boundary [9,10]. The chemical reaction between the fuel and oxidizer occurs
through a single irreversible step. It has been shown [11,12] that, although not capturing
all the quantitative characteristics of the flame, single step chemistry captures many
qualitative characteristics such as flame structure and temperature. Furthermore, even in
multi-step kinetics models, only a few reactions are responsible for a majority of the heat
released during combustion. The molecular weights of the reactants are assumed
identical, the therrnophysical properties p, 2., cp, D are assumed constant, the Lewis
numbers of the reactants are equal to unity, the influences are gravity and radiation (from
the surface and the gas) are neglected, and Soret and Dufour terms are neglected. The
diffusion velocities are given by Fick’s law, and heat conduction is described by Fourier’s

law. Under these restrictions, the equations for conservation of species and energy

become
BY
WW0: pDVzYo —vW (5.1.a)
BYF 2
pv—a—y—=pDV YF-W (5.1.b)
pvcp g} = MN 2T + qW (510)

where W=pAY0YFexp(-E/RT) and the boundary conditions are
T=To, Y0=Yoo, YF=0 at (x=l/2, yZO),(0<x<l/2, y=0) (5.2.a)

T=To. Yo=o, mayFF at (x=-l/2, y20),(-1/2<x<0, y=0) (5.2.b)

106

aT/8y=BYO lay=aYFIBy =0 at y=°o (5.2.c)
Equations (5.1) and (5.2) can be written more compactly in non—dimensional form. We
rescale the independent variables with 1/2 to obtain E,=x/(l/2), n=y/(l/2) so that i=1- 1 in
the two corners of figure 5.1. The species mass fractions are normalized as yp=Yp/YFF,
yo=YolYob. Hence, both yp and yo are bounded by zero and unity. Similarly, the non-
dimensional temperature is T=(T-To)/(TrTo), whereby 1::0 on all boundaries and 2' —> 1 at
the flame sheet. Finally, the non-dimensional mass flux past the divider is Pe=pchU7t.
Additional physical quantities, which are important to our subsequent analysis, are the
global stoichiometric index

= VYFF (5.3)

 

cP Yoo

where v is the ratio of mass of oxidizer to mass of fuel in the one-step reaction; the
Zel’dovich number [3, which is a measure of the temperature sensitivity of the reaction

B=a£.a=1—-T£; (5.4)

the Damkohler number, Da, as a ratio of the diffusion time to the reaction time

= LZ/Ot/pcp)
[M0O exp(-E / RTf )1“

 

Da (5.5)

With the adiabatic flame temperature Tf defined as TFTo'l'QYFF/Cp(l+¢), the non-

dimensional reaction term becomes

_ 1-
W = Yoy F CXP[1——[:(_l%] (5.6)

The non-dimensional conservation equations (5.1) can now be written as

107

PeQE- =szO —(pDaw (5.7.a)

8r]

Peay—F=V2yF—Daw (5.7.b)
an

Pe§=Vzt+(l+tp)Daw (5.7.c)

with V2 = 82 / 8&2 +82 / 8:12. The dimensionless equations (5.2) become

t=O,yo=O,yp=1 at (§=-1,n 2 O),(-1<§<0,n=0) (5.8.a)
1:030: 1 ,cho at (g: 1 ,n 2 0),(0 <C<1.n=0) (5.8.b)
atlan=8yolan=ayF/&n=Oas n—)oo (5.8.c)

Equations (5.7) and (5.8) can be formulated in terms of the excess enthalpy function
H=t+yo+yr1, the mixture fraction z=(l'Zf)YF1'ZF(l-YO), and 1:. Here Zf-=(1+(1))'l is the
value of Z along the stoichiometric contour. The minimum value of Z is zero at the
oxidizer surface, where yo=1 and yF=O and its maximum value is unity at the fuel surface,
where yo=0, yp=l. For our problem, the maximum possible value of H is zero and the
minimum possible value is -1. Functions Z and H satisfy the homogeneous forms of
equations (5.7.a) and (5.7.b), which amounts to a balance between convection and

diffusion. The system of equations (5.7) become

Peg: = V2't+ (1+(p)Daw (5.9.a)
3n
82 2

Pe— = V Z (5.9.b)
81)

Note that the solution for H is H=0 everywhere; the total enthalpy in the gas is constant.

The boundary conditions for equations (5.9) are

108

1:0,b1 at (§=-1,n 20),(-1<§<0,n=0) (5.10.8)

t=O,Z=O at (§=1,n 2 O),(0<§<1,n=0) (5.10.b)

aZ/8n=at/8n=0at 77—)oo (5.10.c)

The simplest case is Pe=0, which produces a purely diffusive problem. The equation for

Z is easily solved to give Zzn—l ln/2—tan—l {tan(1t§/ 2)/ tanh(1m/ 2)}J. The

stoichiometric contour is the arc along which Z=Zf, see figure 5.1. The trailing diffusion
flame lies very near this are, as shown.

These solutions for H and Z can be utilized to produce estimates for the heat flux

to the bounding surfaces. This becomes especially straightforward in the case of infinite

rate chemistry because the flame is attached to the divider, see figure 5.1. On the fuel

side, VZ = (1— ZF)V yF , whereas on the oxidizer side VZ = -ZF Vyo . From the

constancy of H we have either VyF=-Vt or VyO=-Vt, giving

V t = —(1 — Z,:)"1 V Z on the fuel side and V1: = 21:"1 V Z on the oxidizer side. From the
preceding solution for Z we find

~

A - sinh(1m/2)cosh(1m/2) A sin nfi/ 200s 1t§/ 2
Z: 2 2 + 2 2 (5.11)
~ ~ 2(sin n§/2+sinh Inn/2) 2(sin 1t§l2sinh 11:11/2)

Along the lower surface -l<§<l, 71:0 we find

A

—11
Vt: " ,
~ 2(1—ZF)tan1t§/2

 

—1<§<0 (5.12.a)

109

A

11

VI: " , O< <1 5.12.b
~ 2ZFtan1t§/2 C (, )

 

Along the vertical walls at §=-1 and §=+1, we find

A

C
VT: ~ , g:—
~ 2(1—ZF)tan1m/2

 

1 (5.13.a)

-C
Vt: " , §=+1 (5.13.b)
~ 22,: tan 101/2

 

From Equation (5.12.a) the heat flux to the divider is infinite because the flame is
attached to the divider. However, the fluxes to the nearby surface may possibly be
accurately represented by these equations. Part of our goal in this article is to ascertain
the accuracy of Eqs. (5.12.a) and (5.13.a) away from the divider, and thus to determine
the practical usefulness of infinite-rate Burke-Schumann heat flux calculations.

Equations (5.9) and (5.10) will serve as a benchmark for comparing the Burke-
Schumann results and the heat-transfer model results of Section 5.5. These equations will
be solved numerically using an Alternating Direction Implicit (ADD method. Because
the use of the mixture fraction and excess enthalpy functions de-couples the species
equations from the energy equation, equation (5.9.b) is solved first. Then, using the
definitions of Z and H, equation (5.9.a) is SOIVed. The non-linear reaction term in
equation (5.9.a) is handled by using the results of the previous iteration, thus linearizing

the numerical system.

110

_' L. A"

Nil-in

5.4 Simplified Heat Transfer Model “B”
In model “B” the flame arc becomes a constant temperature line, t=1, extending
downstream from the point rq, see figure 5.2. The boundary conditions on the surfaces

are 1:0 and at n —> 00; we have the zero gradient condition 81:an =0 as before. The

governing equation is V21 = O in the case of zero convection; V21 = Pea'c/ an in the case

with uniform and identical convection from the fuel and oxidizer reservoirs adjacent to

the divider. The introduction of scaled coordinates E=§Irq and N=1]/rq places the flame

leading edge at the intersection of the circular arc R=\/L'.-€2 + N 2 =1 and the line Z——Zf.
Here rq is the quenching distance defined as the distance from the divider plate to the
point of maximum reactivity. The flame usually quenches a very small distance upstream
of the point of maximum reactivity, hence there is no inconsistency in defining quenching
distance in this manner. In the general case this transformation offers no advantages for
solution because the scaling merely changes the separation distance between the two
vertical walls from two to 2/rq. When rq becomes small, however, the vertical walls
effectively disappear, and the flame sheet lies along the radially-directed are 13:12:,
1<R<°°. There are no length scales remaining in this problem, leaving a universal heat
flux distribution to the lower surface given by 31(E,O)/ 0N = g(E; Zf ). Here Zr is the

sole remaining parameter in the analysis. In terms of the original coordinates we have

@4831.) (5.14)
3!] rq rq

The maximum flux occurs at the point where the derivative of 8216.0) / 377 vanishes.

111

 

 

The above expression in equation (5.14) is universal to the extent that the flame
curvature can be neglected. If rq is not sufficiently small, the curvature of the flame
towards the vertical axis (see figure 5.1) becomes important and the equation loses
accuracy, producing an over-estimate for the surface heat flux. From equation (5.14) the

integrated flux along the surface is a function only of Z, i.e.,
J (ax/811mg: Ig(E;Zf)dE=Q(Zf). The function Q(Zf) is not symmetric about

Z1=1/2 because the quench radius is asymmetric with respect to Zf. It has been shown in
reference 2 that Q increases as Zf decreases, largely because of the increase in the flame
temperature. This feature, we observe, is not built into this simplified model.

As shown in previous works [2,9,10] the model configuration of figure 5.1 can be
mapped via conformal mapping methods into the configuration shown in figure 5.3. In
this latter configuration the flame is in fact a radial arc'. The heat flux to the surface is

amt/3r]: (31/ 3v)(8v/ 377). Analogous to the previous derivation, we define the scaled

coordinates U=uqu, V=v/rq, where rq=‘/u(2l + v: and U2 + V2 =1 at the semicircular
arc containing the flame leading edge. Then 81/ 811=(Bt/ 8V)(8V/ 8v)(8v/8n), with

81(U,0)/8V-=-g(u/rq;Zf), aVlavzlqu and 8v/8n=(1t/2)costt§/2 at the surface

n=0. Thus
at sin1t l2 ncosn /2
—= g(___§ :Zf)——-"—'g

(5.15.a)
an rCl 2rq

 

rq=,/u§ + v3, =\/sin2(1t§q /2) “1111120111q /2) (5.15.b)

112

For small fig and 11,, we find rq=(1t/2)rq in which case equation (5.15.a) corresponds with
equation(5. 14).

Numerical solutions of this heat transfer model are generated in Section 5.5 for
the heat flux to the surface.
5.5 Results and Discussion

To begin the analysis, comparisons are made between the finite-rate chemistry
model “A”, the Burke-Schumann, and the conduction heat transfer model “B” predictions
of the heat flux distribution to the lower cold boundary. Here, it is assumed that the
diffusing reactants at 11:0 have the constant temperature, To. In figures 5.4a-b the results
from each model are presented for a vertical flame with Zf=ll2 at two different quenching
distances. Comparisons are also made for a fuel-lean and a fuel-rich flame in figures 5.5a
and 5.5b respectively. The results presented are for the Burke-Schumann (B-S) solution
for the heat flux, the heat flux obtained by solving the finite rate chemistry equations
(model “A”), and the heat flux from the simplified heat transfer model “”.B

Examination of the stoichiometric flames in figures 5.4a and 5.4b reveal two
important features. First, the Burke-Schumann solution given by equations (5.12.a,b)
offers very different results for both large and small quench distances when compared to
either the heat conduction and the finite chemistry models. Second, the simple heat
conduction model captures most of the characteristics of the finite chemistry model,
although the heat flux magnitude is slightly diminished. If a scaling factor is employed,
better agreement is achieved. These results are for a pure diffusion flame with no

convection (Pe=0).

 

‘ This is only strictly true in the zero convection case, Pe=0

113

.11

T “v5 r 1.3x sang-q.
.

 

The reaction zone characteristics for model “A” differ from those used in model
“”.B The conduction model is based on a flame sheet consisting of an infinitesimally thin
reaction zone, which extinguishes a finite distance rq from the divider plate. By contrast,
the finite—rate chemistry case produces a more broadly spread reaction zone. Given this
condition,‘it is expected that the heat transfer model will under-predict the heat loss for
the pure diffusion case, as well as have a narrower peak. These characteristics are indeed
observed in figures 5.4-5. Consistent with this, the conduction model produces better
agreement with the finite rate case as the quenching distance, rq, decreases.

For the non-symmetric case Zn: 1/2 (see figure 5.1), the flame shape will be

defined by the following relation between the flame-shape coordinates E; and 11:

tanh(fl)
2 (5.16)

 

It tan(12t— '" TLZf )

Figures 5.5a and 5.5b show the results for two non-stoichiometric cases. The simplified
model follows well the trends of the finite rate case. The model produces a maximum
heat flux at a slightly different location along the lower surface than predicted by the
finite rate chemistry case. The conduction model shifts the maximum away from the
divider. This can be attributed to the wider flame tip for the finite-rate chemistry case.
Even with this shift the heat transfer model is closer to the finite rate chemistry case than
the Burke-Schumann solution. There also exists a strong asymmetry for the non-
stoichiometric flames. Departures towards the fuel rich side produce closer quench

distances, and consequently, larger heat losses to the lower boundary. The maximum heat

flux location is also closer to the divider.

114

Some generalizations about the shape of the heat flux profile can be made. We
wish to provide a better approximation for the heat loss profiles than those given by the
readily available Burke-Schumann solution. We also seek to ascertain where, if at all, the
Burke-Schumann solution accurately represents the heat flux to the lower boundary.

The heat flux profiles for small quenching distances (rq<0.15) obtained from
model “B” were found to collapse to a single curve. This collapse of the heat flux
profiles is predicted theoretically in section 5 .4 by equation (5.14). The correlations were
obtained for both straight flames and non—stoichiometric flames. These curve fits were
obtained by using an ordinary least squares method. The correlations are non-linear so a
Gauss iteration scheme was used.

For symmetric (ZFl/Z) flames, the heat flux profile can be described by the three

parameter relation given below

$1

rq

0.572

=—e

r<1

 

 

( 5 +132)“ (5.17)

rq

Parameters B), [32, and [33 are 1.07, 1.01 and 0.8 respectively. Note that this result does
not include any scaling factor; it is solely the shape of the heat flux profile obtained from
the conduction model. Although this result is applicable only for small quenching
distances, the same shape can be applied to larger quenching distances by changing the
numerical parameters [31, [32, and 83.

It is more challenging to model the non-symmetric flames because several
quantities must be taken into account. The maximum heat loss is shifted away from the

divider, and the spatial decay rates are different on each side of the maximum. The

115

inflection points of the heat flux profile are located at several different distances from the
maximum. To account for these characteristics a five parameter correction factor, C(fi), is

introduced, viz.

 

 

 

 

a a a B“
C(C)=[134(—-l35)+136(——135)2:lexp('137 ‘135 ) (5-18)
rq rq rq
With this correction factor the heat flux profiles take the form
1.07
= 0572 e ’1 ( C +1.01)“8 —c(§) (5.19)
rq rq

 

The five-parameter values are dependent on the stoichiometry parameter 0. Table 5.1
shows values of B4 to Bg for selected 0 values.

The above results can be used in conjunction with a scaling factor (see next
section) to obtain a good approximation for the heat flux profile to the lower boundary.

The case with convection is not considered in the numerical model or in the
analytical model. However, as the discussion of section 5.4 makes clear; when the
quenching distance is small, the heat flux profiles should closely resemble the zero-
convection case. This occurs because a rescaling of the independent variables in the
governing equations with rq will eliminate the convective terms as rq—9 0. Figure 5.6
shows the results when Pe=1. The quenching distance is computed for the finite rate
chemistry case used to determine the heat flux profile from the simple conduction model.

From figure 5.6 we observe first that the model predicts well the heat flux profile.

The qualitative shape of the heat flux profile is identical to the zero convection case.

116

Once again, it appears all that is needed for complete agreement whenever Pe is changed
is the introduction of a different numerical scaling factor in equation (5.19).
5.5.1 Scaling Factor
In order to account for the difference in reaction zone thickness between the model and
the finite rate chemistry case a scaling factor S was introduced. Factor S scales the
maximum heat flux from the heat transfer model to the value obtained in the finite rate
chemistry case. Values of S were obtained for stoichiometric flames and then tested with
non-stoichiometric flames to determine whether these values were still valid. Factor S
was obtained as functions of the quenching distance, not the chemical parameters such as
the Zel’dovich number, for which the dependence is much weaker than quenching
distance. It is noted that the quenching distance depends on the chemical parameters.
Hence, the scaling factor was correlated with the quenching distance. The relatively
weaker B and Damkohler number dependence manifests itself implicitly through the
effects of the reaction zone thickness on the quenching distance.

It was found that for the stoichiometric flame, S had a nearly linear dependence on
rq. The relationship is shown below

S=1.72rq + 0.66 (5.20)

Equation (20) can be used as a direct multiplication factor in equation (5.19) for small
deviations from stoichiometric flames, 0.33<¢<3.0, with S values within 10 percent of the
numerical results.
5.5.2 Comparison of centerline heat losses

Examination of the temperature profile along the stoichiometric contour Zp=ll2

for vertical flames is used to verify the scaling factor as well as to identify any

117

shortcomings of the proposed heat transfer model. The preliminary comparison was done
for a moderately small quenching distance, rq=0.375, and a pure diffusion flame. The
quenching distance was calculated by using the finite-rate chemistry code to solve for the
quenching location, as defined by the maximum reactivity, as well as the vertical
temperature distribution at the centerline. The resulting quench distance was employed in
the heat conduction model and the temperature profile was calculated. A graphical
comparison is shown in figure 5.7.

From figure 5.7 we observe that the heat conduction model underpredicts the heat
loss to the divider plate. This is expected in view of the relative magnitude of the heat
loss profiles to the lower boundary. Furthermore, as expected, if a ratio of the finite rate
gradient to the heat transfer model gradient is examined it will produce the scaling factor
S already examined. With this information, and noting that the scaling factor had a small
dependence on the Zel’dovich and Damkohler numbers which will be ignored here, the
gradient ratio should behave very similarly to the scaling factor with respect to rq. We
recall that S varied nearly linearly with the quenching distance. Examination of. that
relationship in equation (5.20) suggests that the heat conduction model gradient will
approach the finite rate model gradient as the quenching distance becomes small.

Another feature shown in figure 5.7 is the different flame temperatures of the two codes.
The heat conduction model assumes a flame temperature of 1.0 for the flame and uses
that elevated value as its quasi-heat source. The finite rate model temperature is
calculated from the set of partial differential equations shown in section 5.3. The finite
rate model produces a lower flame temperature, as expected, since Tflamcz l-O(B'1). This

contributes to the more gradual decline in temperature towards the lower wall, and

118

consequentially the more linear gradient. This linear gradient supports the simplified
notion of assuming a linear profile from the quench location to the wall [2,10]. However,
this assumption is typically based on an adiabatic flame temperature close to unity (i.e.,
B—) 001th), and does not take into account any heat losses through the reaction zone.
The heat loss from the reaction zone in the finite-rate model “A” is evident beginning
near 11:05 and continuing to the location of the heat-transfer model flame tip at 11:0.375.
5.5.3 Comparison to Finite-Chemistry Model

A comparison of our heat transfer results may be made to the finite-rate chemistry

model examined in [2,10]. The heat transfer cold wall beneath the flame leading edge is

 

 

given by
_ 19 Da 4 1/2 . 2
q—‘W’ EEO-2f) Sll'l (TCZf) (5.21)
where
8 =1+{(1- z. )[m 11(2, _ 50) -1]+ zf [5m ”.(Z‘ + 8F) - 1]} (5.22.1)
srn 1th srn 1th
8q = 3.752{1 -exp(-7.366Z}'71)}Z?'36 (5.22.b)

and 50=Zg(l—tq), 8F=(l-Zf)(1-tq), where tq, the flame leading edge quenching temperature,
is given by l-K/B. It is shown in section 5.5 .2 of [2] that a value of K of between 1.5 and 2
is reasonable. We use 1c=1.75. For the symmetric case Zr=1l2 we find 50=6p=1cl28,
whereby 0=cos1n<l28 and bq=2.62, giving q=0.21¢os(1tK/28)(Da/B3)m.

The use of 1c=1.75, and [#80 gives the estimate qEO.l98(Da/B3)”2. In order to

relate this result to the quench distance, we employ equation (8) of [2], which gives

119

 

,q = 2 q 3/2 /5_ (5.23)
n 42f (1-zf) Da

where we have used the factor 182 to convert the rq of equation (8) of [2] from the u-v

system of coordinates to the é—n coordinates used here, as outlined in section 5.4, see

equation (5.15.b) et sequation. The preceding numbers give rq_=_3.82([33/Da)“2

, whereby,
q50.76/rq , which compares favorably with equation (5.17), q50.57/rq, especially when we
consider that the latter must be multiplied by a scaling factor greater than unity given by
equation (5.20). Consequently, the heat transfer model prediction is functionally
consistent with the more detailed finite chemistry model. In addition, and in contrast with
reference [2], the heat transfer mode] produces an expression for the distribution of q
along the entire surface, not only its maximum value. For this reason, it is possible to
replace the factor 0.572/rq in equation (5.17) with equation (5.21), using equation
(5.22a,b) for 13 and bq respectively. The pure heat transfer formula now includes the
influences of finite-rate one-step chemistry. A similar replacement, of course, may be
made in equation (5.19) for the asymmetric case.

5.5.4 Zel’dovich number effects

As stated previously, the temperature gradients at the lower boundary produced by the
two models show better agreement as rq decreases. That is not to say that the profiles
become more similar under this criterion. The presented profiles can be divided into
three general regions: the flame tip region, the wall region, and the transition region
which connects the two aforementioned regions. As will become evident from figures

5.8a and 5.8b below, the transition region between the two models have distinctly

different shapes. The conduction model produces a much steeper gradient immediately

120

after the flame tip and a much smaller gradient at the wall. Hence the transition region is
much more curved than the nearly linear temperature profile produced by finite rate
chemistry model. Having noted that the profiles agree well in the wall region when the
flame quenches very close to the lower boundary and that the transition regions will have
their own characteristics, it is of interest to find the conditions where the flame tip regions
exhibit agreement.

The finite rate model produces a flame temperature that is lower than the adiabatic
flame temperature. This suggests that for the flame tip regions to agree, the finite rate
model should produce a flame temperature closer to the adiabatic temperature assumed
for the conduction model. This can be achieved by increasing the Zel’dovich number, [3,
which also has the effect of increasing the quench distance. This produces a higher flame
temperature and lower heat losses through the flame tip, thus improving agreement in the
flame tip region of the temperature profiles. However, the gain in agreement in one
region reduces it in another, as can be seen in figures 5.8a and 5.8b. Agreement may also
require an unrealistically large value of B.

The Zel’dovich number also depends on the wall temperature T0 through the
relation B=(E/RT2)(TrTo). When all quantities other than T0 are fixed, it is possible to
examine the influences of wall temperature on heat flux and flame quenching distance.
Since q.~.0.57/1q and rqocmro)”, we find qocmroym, drq/dTooc-(TrTof’z, and
dq/dTooc(Tg-To)'5’2. Thus as To increases, rq vanishes at a successively decreasing rate
while q increases at a successively increasing rate. Under practical conditions, we expect

that the near-wall flame tip can heat the wall, raising TO and decreasing rq, while

121

increasing q. It is not clear how long this incremental process can continue, for there are
upper bounds on the chemical production of heat that leads to q.

5.5.5 Characteristic Flame-Tip Width

For small rq, a practical estimate of the width of the triple-flame structure, shown
as the lowermost arc of figure 5.1, is obtained by exploiting the relationship for q(§)
given by equation (5.17). We note the closeness of [31,132, and B3 to unity and write
q(§)=Cexp(-E/rq)(l+§/rq), where C=0.572/rq. This function possesses an inflection point
at §=rq, which we interpret as a measure of the triple-flame tip structure half-width, 1,.

Thus, we write
1r = 2rq x113” exp(E/ZRTf) (5.24)

This result agrees with a previously derived estimate [3]. As the flame temperature
increases, all else remaining constant, both [3 and the exponential factor decrease. Hence,
lr decreases exponentially and algebraically as flame temperature increases. Conversely,
as the non-dimensional activation energy increases, Ir increases exponentially.
5.5.6 The q«-rq Relationship

In our study, the relation (5.17) was derived for infinite-rate chemistry and the
modifications of section 5.5.3 were appended under the approximation that the infinite-
rate chemistry limit was a reasonably accurate approximation. In this case, the essential
prediction is qrq~O(1) and constant. We expect this relation to apply in a broad middle
range but not at the extremes q—>oo (rq-—)0) and q—)0 (rq—->oo). The latter limit may be
interpreted as flame liftoff, whereas the former as flame extinction. Liftoff, followed

eventually by blowoff, occurs when heat losses to solid boundaries cannot be sustained.

122

;TR9

'( '5'!)

 

The separation distance between surface and flame increases until the flame effectively
no longer interacts with the surface. At blowoff, even gas-phase flaming becomes
impossible. Extinction, on the other hand, occurs when the flame loses too much heat to
the surface to survive. It produces enough heat that lifting is not necessary, but the losses
to the surface become excessive. Our numerical integrations have in fact produced both
limiting cases, although the detailed description is beyond the scope of this article.
Nevertheless, it appears that the extinction limit is attainable without radiant losses from
the flame tip, without thermal expansion of the gas, or without buoyancy. All that is
required for a minimalist description is chemical heat production and conductive surface
losses.
5.6 Comparison with Experimental Results

The relationships derived in section 5.5.3 allow calculations to be performed
when all pertinent information about the problem is known, such as the overall heat
release, stoichiometry, chemical rate data for a one-step global reaction, etc... Such data
data readily available for many hydrocarbon fuels such as methane (CH4) or Propane
(C3H3), thought the accuracy of one-step reaction data is suspect as will be seen below.
These two fuels are used to estimate the heat flux profiles along the nearby cold boundary
and the methane estimates for the maximum heat flux are compared to experimental data.

The information used in conjunction with the equations in 5.5.3 was found in
reference 13, and appears in table 5.2. The constants A, m, and n correspond to the rate
equation given below.

H
£11(—::lxt—y-J=—Aexp(- Ea/RuT)[C,‘Hy]m [02]" (5-25)

123

where A is a pre-exponential factor, m and 11 represent independent reaction orders, and

Ea is an activation energy. All units are consistent with gmol-cm-sec.

Table 5.2 Reaction information for Methane and Propane

 

 

 

 

 

 

Fuel A m n EalRu (°K) HHV(KJ/Kg)
CH.I 8.3x105 -0.3 1.3 15098 55,528
(:11.2 1.1111018 1 2 15100 55,528
C3Hg 8.6x10” 0.1 1.65 15098 50368

 

 

 

 

 

 

These values were used to obtain estimates for the flame temperature and quenching
distance for various stoichiometries. Note that there are two different global reaction
schemes used for Methane/Oxygen combustion. The two are included to show the large
discrepancies in these models and how they effect the final estimates. The results appear
in figure 5.9 and 5.10 for Methane and figure 5.11 for Propane.

Examination of figure 5 .9 shows that the heat flux distribution is fairly flat across
the lower boundary for fuel lean cases (¢<1) and become more defined when the fuel rich
case (¢>l) is examined. The broadness of the heat flux profile is due to a large quenching
distance of the flame. As 11) increases, the flame tends to be able to survive closer to the
lower boundary and hence produces a larger heat flux. This is a characteristic that is
witnessed for both fuels and is primarily caused by the increase in the adiabatic flame
temperature for fuel rich mixtures. The profiles for Methane combustion using the
second global reaction scheme are much more defined. This is due to the increased rate
of reaction predicted by the second scheme and hence the flame can survive closer to the
cold boundary.

Two experimental studies were used to serve as benchmarks for our model. The

first involved a bumer-attached flame burning Methane fuel [14]. In this study, the fuel

124

. ".‘-I ‘Vfi‘F—w

-‘."—_fi

was flowing into a quiescent oxygen environment so the present model, which assumes
zero convection, may not be entirely appropriate. Additionally, the second experimental
study involves both flowing fuel and oxidizer, in this case air [15]. The heat flux
distribution is not examined in either of the above cases but the maximum heat flux
information can be obtained. With all the above differences it is expected that the
comparisons will be able to yield an order of magnitude estimate, and help validate the
results from our study.

For both the studies examined a stoichiometric mixture is assumed, that is 11) is
taken to be unity. From reference 14 it is found that the maximum heat flux is on the
order of 1 W/cm2 while for the co-flowing reactant case [15] finds q”max to be roughly 2
W/cmz. Comparisons with figure 5.9 show that the maximum heat flux for ¢=l is on the
order of 1 W/cmz, while figure 5.10 shows a maximum heat flux of order 5 W/cmz. Both
these results show the maximum heat flux to be on the same order and agree well with
experimental results thus adding credibility to the proposed heat flux distribution for
burner attached flames.

5.7 Conclusions

The conductive heat transfer sub-model reproduces many of the global features of
the finite rate chemistry case. The profiles at the lower boundary are similar for both
models with only a scaling factor needed to produce accurate agreement. However, when
examining the temperature profiles along the centerline of the lifted flame the two models
can produce significantly different results for the transition region and the flame tip
region, especially as the quenching distance increases. This can be partially attributed to

the assumption of a flame sheet in the heat transfer model, which differs from the non-

125

zero flame thickness produced by finite rate chemistry. The conduction heat transfer
mode] is unable to accurately predict the heat loss from the flame tip, even for pure
diffusion flames. For this reason the conduction model is not useful for describing
thermal phenomena close to the flame tip, but is very useful in examining events away
from this region as is generally true for other heat transfer models of detailed flame
process [5]. This limitation might be overcome by introducing a volumetric heat
generation term into the conduction model that has temperature dependence, but this is a
subject for future study.

For modeling finite rate chemistry influences on the heat flux, we recommend a
combination of equations (5.17) [or (5.19)] and (5.21), with the later being substituted for
the factor 0.572/rq appearing in the formula. This combined equation accurately predicts
the flux beneath the flame leading edge and also its distribution along the surface.

It is also noted that the results obtained herein may be limited to configurations
closely resembling the slot configuration used in this work. The heat flux profiles for
other lifted flames with different physical and flow orientations may produce different
correlations. It is a limitation of the work that the methods used here will need to be
followed for different problems. However, the usefulness of the results is not diminished
for this lack of adaptability. The general result of our research is to show that beneath the
flame tip q~0(1/rq), where rq~O([B3/Da]”2), to within a multiplicative function of global
stoichiometry. Although this result may be deduced by simple dimensional analysis and
scaling arguments, we have in addition quantified the heat flux distribution beneath the

flame tip.

126

I.
31

 

5.8 References

1.

9.

10.

11.

12.

13.

14.

15.

Nietzsche, F., Twilight of the Idols, R.J. Hollingdale (translation), Penguin, London,
(1968)

Wichman, I.S., Pavlova, Z., Ramadan, B. and Qin, G, Combustion and Flame, 118,
651-668 (1999)

Wichman, LS. and Ramadan, B., Phys. of Fluids, 10, 3145-3154 (1998)

Williams, F.A., Combustion Theory, 2“d Ed., Benjamin-Cummings, Menlo Park, CA.
(1985)

Wichman, 1.8., and Williams, F.A., Combustion Science and Technology, 32, 91-123,
(1983)

DiBlasi, C. and Wichman, I.S., Combustion and Flame, 102, 229-240 (1995)
Wichman, I.S., Lakkaraju, N. and Ramadan, B., Combustion Science and Technology,
127, 141-165 (1997)

Ramadan, B., Personal Communications, (1997-1998)

Wichman, I.S., Combustion Science and Technology, 64, 295-313 (1989)

Wichman, I.S., Pavlova, Z., Ramadan, 3., “Attachment of Diffusion Flames Near
Cold, Chemically Inert Surfaces With and Without Reactant Flow”, MSU CRL 09-
24-96 (1996)

Rightly, M. L., Williams, F. A., Combustion and Flame, 101, 287-301 (1995)

Hsu, P., Matthews, R.D., Combustion and Flame, 93 n. 4, 457-466 (1993)

Turns, S.R., An Introduction to Combustion second edition, McGraw-Hill, (2000)
Robson, K. and Wilson, M., Combustion and Flame, 13, 626-634 (1969)

Takahashi, F. et al., 1988. ,Journal of Heat Transfer, 110, 182-189 (1988)

127

 

 

 

 

 

1 7"
; .—
l
: .—
T=T0 ‘- T=To
Yo=0 ‘_ Y0=Yoo
YFYFF Yp=0
T=0 ‘— r=0
Yo=0 —’ -— yo=l
1 _
YF‘: _> 1“” yF—O
(III) *‘
_. ”11,: L
3 1 7‘ ; ’ 1‘ T i
=-l/2 ‘ =1/2
53:1 T=To 1:0 T=To 1:0 :1
Yo=0 Yo=0 Yo=Yoo Yo=1
YF=YFF yp=l Yp=0 yF=O

Figure 5.1 Model configuration and boundary conditions. The side and lower walls
are porous to reactants. The vertical walls admit only diffusive fluxes. In the far
field the diffusion flame is one-dimensional. The near field triple flame structure is
shown, with (I) as the diffusion flame, (II) and (III) as rich and lean premixed flame
arcs, and (IV) as the triple point or leading flame edge

128

 

 

 

 

T-To T”To
_ inn
T on:
s r
C I a 5'1

Figure 5.2 Heat transfer model for infinitesimally thin flame sheet and
isothermal boundary conditions. This model is the heat transfer
simplification of the flame in figure 5.1 because by far the greatest heat
release is along the diffusion flame arc

129

V—)°°

1“ . ‘."- ‘. farm.“ .-
4 1- hull.

Z=Zf

. V

Z=0 u—>+oo

 

‘
“—"°° 2:1 u=v=0

Figure 5.3 Model configuration of mixture fraction Z after conformal
mapping. Note that the straight diffusion flame along Z=Zg occurs only in
the case of zero convection. In addition, lines of constant Z are radially-

directed arcs, Z=0ltt

 

130

 

 

 

 

 

 

 

 

 

 

 

 

Finite Rate
q — — — — B-S
------- Model
,5
. \
I 4 __ 1
I \
.. ,’ 3 \ Finite Rate
q \\
2

 

 

 

 

 

 

 

 

Figure 5.4 (a) Heat flux profiles comparison for straight flames rq=0.27. Note
the inaccuracy, generally, of the Burke-Schumann calculation of
Equation(12.a,b) and the comparative similarity of the Finite-rate and Heat
Transfer model results. A scaling factor renders the latter two visually
indistinguishable. Heat flux comparisons for straight flame rq=0.9. (b) Heat
flux profiles comparisons for straight flames rq=0.9. The comments to (a)
apply here as well

131

 

 

Finite Rate
— — — — B-S
------- Model

 

 

 

 

 

 

 

 

qb

 

 

Finite Rate
qu ————B-S

------- Model

 

 

 

 

 

 

 

 

Figure 5.5 (a) Heat flux profile comparisons for fuel lean flames phi=0.5,
rq=0.38. For this non-symmetric case the flux distribution is asymmetric.
The Burke-Schumann model is once again a poor approximation to the heat
flux. (b) Heat flux profile comparisons for fuel rich flames phi=1.5, rq=0.28.
See comments in figure 5.4a

132

Table 5.1 Correction Factor Coefficients

 

 

 

 

 

 

 

 

 

 

 

 

9 B4 135 96 137 BB
0.333 1.045 -0. 179 -0.3 l 3 1.074 0.977
0.429 0.746 -0.08 -0.210 1.039 1.018
0.538 0.563 -0.096 -0.127 1.093 1.004
0.667 0.392 -0.194 -0.046 1.170 0.974
0.818 0.204 -0.084 -0.024 1.227 0.964
1.222 -0.204 0.084 ~0.024 1.227 0.964

1.5 -0.392 0.194 -0.046 1.170 0.974
1.857 -0.563 0.096 -0. 127 1.093 1.004
2.333 -0.746 0.08 -0.210 1.039 1.018

3.0 - l .045 0.179 -0.313 1.074 0.977

 

 

 

 

 

 

133

 

 

qll

 

 

Finite Rate

——-——— Model

 

 

 

 

 

 

 

Figure 5.6 Heat flux profile comparison for straight flame with convection

Pe=1,rq=0.395

134

" -‘ o
“M n“mmfl
. ..
I

 

1.2

- - . mew-1‘1,

 

I+,+++++++++++-
.1.

 

 

— —-I-— — Finite Rate
Heat Trans.

 

 

 

 

 

 

 

Figure 5.7 Centerline temperature distributions for a straight flame rq=0.375.
The broadness of the actual flame moderates the temperature variation,
compared with the sharp heat transfer profile. The 1: values along the
diffusion flame (n>0.5) agree to within 0(8'1)

 

135

1.2

 

 

 

 

 

 

 

- -+- - FR beta=15.0

Tau Heat Trans.

 

 

 

 

 

 

1.2

 

 

[+++++ ++ ++'-
+

 

 

 

 

- -l-- — FFl beta=10.0
Tau Heat Trans.

 

 

 

Figure 5.8 (a) Centerline temperature profile for a straight
flame rq=0.525. (b) Centerline temperature distribution for
a straight flame rq=l.225. For large quench distances the

pure heat transfer temperature profile differs significantly

from the finite-rate chemistry profile

136

2.5 1

 

 

"A 15 +¢=0.667
5 ‘ +¢=0818
g +¢=1.0
:3 l ‘ +¢=l.5
+¢=3.0

 

0.5

 

0 I I I I I I I I l

-0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02

x(m)

Figure 5.9 Heat flux distribution for Methane/Oxygen combustion at various
stoichiometries

137

161

CI‘IJ'FOZ

q" (W/cmz)

 

 

 

-0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02

x(m)

Figure 5.10 Heat flux distribution for Methane/Oxygen combustion at various
stoichiometries. The second reaction mechanism given in table 5.2 is used

138

10 ‘ C3H8 + 02

    

 

1 ~ 1 ' ' '
0 ‘ . _ _ ‘ _ . p . _
-0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02
x(m)

Figure 5.11 Heat flux profile for Propane/oxygen configuration. Larger heat losses
than Methane are predicted

139

NUMERICAL TECHNIQUES

6.1 Numerical Techniques
The equations examined throughout this document represent models of more complex
systems. The problems are greatly simplified by assuming the flow to be known and
uniform and the fluid properties assumed constant. This eliminates the need to solve for
the flow field through the Navier—Stokes equations. This model of uniform flow and
constant density is a simplifying assumption to be sure; but it allows the thermo-diffusive
effects, i.e. the Lewis number effects, to be examined without the complicating additions
of an unknown flow field. However, even with these simplifying assumptions, the
problem is still non-linear because of the chemical reaction that eliminates or creates
constituents within the domain and adds thermal energy to the system to drive this
transformation. Because of this non-linearity, analytical solutions become quite difficult
and thrust one into the realm of asymptotic analysis. Asymptotic solutions, although
insightful, are generally limited to “simpler” problems, i.e. those with only a single value
for Le, simplified kinetic expressions, and one dimension. To study more complex
combinations of parameters a computational solution was needed. In the text below the
techniques used to solve the various systems of equations are discussed.
6.2 Transient premixed flame propagation

The governing equations for this problem have been described in chapter 4. In
this section the solution techniques used and the issues that were addressed will be
examined by using a generic one-dimensional equation with the same form. The

following equation will be examined,

92:33:

81 3x2 +f(¢), (6.1)

140

 

where 11) represents a vector of dependant variables, and the function f is a non-linear
function of the dependant variables.

There are several concerns that arise before numerical solutions begin. Primarily,
the numerical stability of the problem as time increases. Furthermore, this stability issue
is clouded with the nonlinear reaction term. It is clear that without the non-linear source
term, a numerical scheme could be used that would be unconditionally stable, such as the
Crank-Nicholson procedure. However, the need to linearize the equation for solution
also raises stability issues. The fact that 11) may contain more than one variable and be
coupled through the reaction term is addressed by solving all the equations in the system
simultaneously.

There are two ways to linearize equation (6.1): evaluate the non-linear term at the
previous time step, thus producing a source term with a known value, or expand the
function with a truncated Taylor series around the previous time step so that the function
is represented as a linear equivalent. The discretized equations for both of these methods
will be outlined below.

Along with the time lagging technique to linearize the system we apply a forward
differencing for the time derivative and a central differencing for the spatial derivative.
This scheme is first-order accurate in time and second in space; so in order to keep the
solution reasonably accurate a small time step will be needed. An additional restriction
on the size of the time step is introduced with this method. If the non-linear reaction term
is going to be accurately represented by the value at the previous time step, then the time

step cannot be too large. The discretized equations then become

n+1 n n+1 n+1 n+1
° ‘91 ¢1+1 “291 +¢i-1 n
‘ = +f 6.2

At Ax, 11> ) < )

141

where the i subscript represents the node number, while the superscript n represents the
time index. Furthermore, At and Ax are the time step and grid spacing respectfully. At
this point it should be pointed out that for problems where there are rapid changes in the
values of the variables both temporally and spatially such as combustion wave fronts or
shock waves, a non-uniform grid is usually desirable. In this way, there are more grid
points in the vicinity of the sharp changes so their characteristics are captured more
accurately. A non-uniform grid makes the discretization a bit more complex but does not
affect the solution procedure. In the schemes used in this text a uniform grid was used,
after the difference between a non-uniform and uniform mesh proved to be insignificant.

Equations (6.2) are then solved simultaneously in a block tri-diagonal format.
This system will produce a banded matrix of the form Ax=b, where the vector x
represents all the dependant variables of (1). Here A is a square matrix of order C*N,
where C is the number of dependant variables and N is the number of grid points.

The lagging method of linearization produced a stable scheme only when the time
step was very small. If the time step became too large, the non-linear term was unable to
be accurately modeled using previous time step data. As it turned out, for the flame
propagation problem, the time step using this method became so small that round off
error may have become significant. This is an unacceptable situation, both from an
accuracy standpoint as well as a practicality perspective. For this reason, the nonlinear
term was linearized to include information from both the previous and present time
indexes. This is discussed below.

The truncated Taylor series expansion of the reaction term is shown below,

142

111““) - f(¢” ) +A¢[-§%] (6.3)

where A¢=¢"+'-¢", and 3%¢ represents a n by n Jacobian matrix, with n being the

number of dependant variables in (1). The resulting set of algebraic equations is then
solved. By linearizing the system in this fashion and solving all n equations
simultaneously, a stable solution can be achieved by means of a much larger time step. It
is important to note that the time step must still be small to assure stability, but not as
small as the lagging technique would require.

The boundary conditions of the first kind were implemented by fixing the value at
the boundary node to equal that required by the boundary condition. The zero flux
boundary conditions, or boundary conditions of the second kind, i.e. a fixed gradient at
the boundary, were done in such a manner to ensure second order spatial accuracy. The
method of image points was used so that the following discretized equation would result.

«1131' = 2 * Ax * q"+¢1"_i' - (6.4)
In this case the i-l node represents the image point and is found as a function of the
interior node at 1+1 and q” is the heat flux at the boundary. Note that if this were to take
place at the boundary where i+l was the point outside the boundary, that would be the
image point and the exterior node would be found in terms of the interior node.

6.3 Two-dimensional steady state with co-flow reactants

The governing equations for this system have the form

34> 82¢ 82¢
c D +— =1 6.5

143

where C and D represent constants, and f is a non-linear reaction term. The reaction term
is handled identically to that of section 6.2.

Steady state schemes differ slightly from their transient counterparts in that there
is no temporal evolution. In steady schemes a fixed solution is desired and there is no
definite beginning and ending to the process. In this respect, the goal is to "find" the
solution via interative approaches, or in some instances a quasi-transient procedure. It
should be noted that although there is no temporal evolution, the iteration itself acts as a
pseudo-transient scheme so that numerical stability concerns can not totally be ignored.
Furthermore, with the second dimension, the set of algebraic equations will no longer
produce a block tri-diagonal form. If the set of equations is to be solved simultaneously,
then a block penta-diagonal matrix will result. Since solution methods for a penta-
diagonal scheme take significantly longer than those for a tri-diagonal matrix, it is
desirable to utilize a tri-daigonal method in solving. This is accomplished by using an
alternating direction implicit (ADI) scheme.

The ADI scheme is characterized by discretizing in one direction with an implicit
manner, and the other explicitly. This set of equations is solved and the solution
represents one half of a full iteration. The procedure is reversed and the formerly implicit
direction is then treated as implicit and visa-versa. The solution obtained then represents
one full iteration.

The discretization techniques are slightly different from those in section 6.1. The
second derivatives, i.e. the diffusion term, are still modeled by second order correct
central differences; but some care must be taken when examining the first derivative, i.e.

the convective term. If C is large enough, the convective term become dominant and a

144

 

wave type equation exists. In this case an upwind scheme is needed to accurately model
the problem. This is accomplished by using a second order accurate backward

differencing shown below

fl z 34%,) -4¢i-1.j “111—2,)
By 2Ay

 

(6.6)

Here, the i subscript refers to the nodes in the x-direction and the j subscripts to those in
the y-direction. It is also noted that the backwards direction refers to the direction
relative to the convective direction.

If C is less than some critical value, a second-order correct central differencing
can be used on the convective term with no numerically instability associated with wave
phenomena. This value of Ccritical can be found by examining the sign of the 411.1 and (1),-)
terms in the discretized equation. If the condition below is satisfied, then a central
differencing can be used with no loss of accuracy.

|C|Ay < 2|D| (6.7)
Since this condition is satisfied for the co-flowing reactant problem, a central
differencing was used for both the convective and diffusive terms.

The ADI method employed differed slightly from that described above. It was
found that the numerical scheme was much more stable when the central differencing
schemes involved evaluating the i,j term at the current iteration. This essentially resulted

in the following discritized equations.

 

 

k “2 k
C¢Cj+1'¢i‘.j-1_ ljfiz‘z‘l’tilflmfigz”RH-24)"; WU“ -f(¢"-)
2A 2 2 _ "J
y Ax A)’
(6.8)
k
k /2 af

+8115” -¢Cj{§b]

145

Solving with the y-direction derivatives evaluated at the (k+1) iteration step completes
the ADI. This scheme allows a block tri-diagonal scheme to be used for both directions,
and utilizes less CPU time than solving the block penta-diagonal system.

After each iteration each point is subjected to some relaxation to improve
convergence and help with numerical stability. The following equation represents how

this relaxation is done:

11",)“ =w¢“,}" +(1— out, (6.9)
where, the k+l* superscript indicates the undated next guess in the iteration procedure.
The value of 0) is chosen by examination of the residuals to equation (6.5). This will be
discussed in section 6.5. This iterative procedure will continue until the differences
between the updated solution and the previous solution is small.

6.4 Two-dimensional Transient with Co-flow reactants
The governing equations are similar that those seen in section 6.3, with the exception

that there is a art/8t is added to the left-hand side.

@ _a_1>_ 62¢ 32¢
at+Cay D[ax 2 +872}: f(¢) (6.10)

Furthermore, as previously stated, a block tri-diagonal scheme is desirable due to the
computational efficiency. Hence, an ADI scheme was implemented. The solution
procedure follows that described by Shih [l] in which approximate factorization was used
with a delta formulation of the difference equations. Equation (6.10) can be discretized

into the following formulation:

n+l_ n n+1
q)” _]_A__t¢1")_ :Kflj 4.6—4’] :I+O(At2) (6.11)

146

The right-hand side of equation( 6.11) is then re-written as

Ff Cad) 82¢ 82¢ \n -
+—-D—+D—+f(¢) +
RHs—l \ Cay ax 8y 1
—-§( 82 \M] (6.12)
-C-a—¢+ D 3¢+Da__2 ¢+f(¢)

 

 

 

 

The source term f(¢) is non-linear and hence cannot be handled directly at the n+1 time
step. The linearization is handled as in section 6.2. In the delta formulation, the M)“ =
¢“+‘-¢" which appears in the linearization is kept as the dependent variable and all
discretization is performed on this delta term. The variable B will represent the J acobian
matrix of Bf/Bq) evaluated at the nth time step. Substituting equation (6.12) into equation

(6.11) and introducing subscript notation to represent the derivatives yields

At
I+—(C¢ —D¢ -D¢ +B)]A¢" =
[ y x" W (6.13)
At(— Q), + D1,, + D1yy + 1(1))“
Here I represent the identity matrix. The approximate factorization is introduced to

produce a directionally dependent representation of equation (6.13).

At At
[N+-2—(c¢y— D¢W)]N [N+?(— D1),,.)]=RHS2 (6.14)

Where N = I +BAt/2, and ms; is the right-hand side of Eq(6.13). Equation (6.14) is
then solved by sweeping in the y-direction and then the x-direction via the following

procedure:

[N+At-(C¢y- -D¢yy)]A¢ =RHs2 (6.15)
At ,, ..
[N+—2—(— D¢xx)]A¢ =NA¢ (6.16)

147

 

$11.?! = ‘91".1 + “’31 (6.17)

Equations (615,616) are solved via second order correct finite differencing subject to the
constraints listed in section 6.3.
6.5 Error analysis

Code validation is an essential component to proper numerical analysis. This
validation is commonly done by comparing experimental data to the output from the
numerical method. This approach is useful when the governing equations used simulate
the actual physical phenomena rather than model it. A simplified model used to predict
only qualitative results will have a difficult time comparing well with experimental data.
However, because of the simplicity of such models, the governing equations are generally
well understood and several techniques can be used to ensure that the correct equations
are being solved.

Initially, the finite difference equations (FDE) are examined for consistency. That
is, a Taylor series expansion for each term is done about the i,j point. These are then
substituted into the FDE and check to ensure that the original equations are obtained as
the spatial or temporal terms approach zero. This will identify if the scheme is consistent
or conditionally consistent. See reference 2 for more details.

Once the FDE has been solved, it is necessary to check and see if the correct
equations have indeed been solved. The following describes how this was done for
steady problems. The first step was to identify the optimal convergence parameter for the
over-relaxation. Optimal in the sense that the residual sum of squares is minimized.
Secondly, the errors are estimated by assuming the equation solved has the order of

magnitude of error do the finite differencing as the original FDE.

148

It has long been known that when using a successive over-relaxation (SOR)
procedure that the choice of the relaxation parameter a) can greatly impact the numerical
schemes efficiency. However, it can also impact the accuracy of the numerical solution.
This is shown by looking at the residual sum of squares. The residual is evaluated at each
node and is effectively the imbalance between the right-hand side and the left-hand side
of the governing finite difference equations. Each residual is then squared and summed
over all interior nodes. If the numerical scheme was solving the FDE exactly, then the
residual sum of squares would be zero; however, this is not the case. So, a) is varied over
the range 0<m<2 to find the value that minimizes the residuals.

Examination of the error behavior is also an important check to determine if the
equations are being solved correctly. To do this it is necessary to know the accuracy of
the FDE used in the model. If the problems examined are second-order accurate in both
time and space, we can now assert that the true value can be given by the following
equation:

T=Rh +Eh2 (6.18)
where T is the true value, Rh is the numerical result at grid spacing h such that Ax=Ay=h,
and E is an order one constant. Taking results for three different grid sizes to be

h,h/2,and h/4 the error term can be estimated as follows.

2
E(%) = 5h1_2_3__R_h (6.19)

The true value can now be estimated for the grid spacing NZ and h/4 since the errors are

known. If the errors behave in a second-order correct fashion, then the difference

149

between the two estimated true values should be less than the error for the h/4 grid

spacing, and should be much less than order hz.

 

lTh/z‘Th/4l<iRh/4th/2I<<O(h2) (6.20)

If equation (6.20) holds true, then the numerical scheme behaves as expected and one can

have reasonable confidence that the correct equations are being solved.

6.6 References
1. Shih, T. and Chyu, W. AIAA Journal, 29, 1759-1760 (1991)

2. Anderson, D., Tannehill, J. and Pletcher, R. Computational Fluid Mechanics and Heat
Transfer, Hemisphere Publishing (1984)

150

spa——
. J

   

E

11111111111111111111111111