r—1
-3 F33 i

N
(x
g;

Ph.D.__

This is to certify that the
dissertation entitled

Diffusion Flame Stability

presented by

Amy B. Moore

has been accepted towards fulfillment
of the requirements for the

degree in i Applied Mathematics

 

ltd/LL

MMajor Professor's Signature
57/ l [ 2006

Date

MSU is an Afiinnative Action/Equal Opportunity Institution

__—.————,

LIBRARY
Michigan State
University

 

 

 

 

_._-—.—-—--..—4J

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 p:lC|RC/DaleDue.indd-p.1

 

DIFFUSION FLAME STABILITY
By

Amy B. Moore

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Department of Mathematics

2006

ABSTRACT
DIFFUSION FLAME STABILITY
By

Amy B. Moore

We analyze the solutions of a boundary value problem arising from a onedimensional
diffusion flame. We state properties satisfied by the linear operator of the system and
prove existence of the steady state solutions for certain parameter values. It has been
well established that the steady state solutions form an S-curve when no radiative
heat losses are included. We show that this S-curve transforms into an island shaped
curve and an ignition branch for large activation temperatures when radiation is in-
cluded, but islands never form for low activation temperatures. We also analyze the
stability of the steady state solutions by analyzing the eigenvalues of the linearized
system. The evolution of stable oscillations is seen for certain parameter values by
perturbing steady state solutions. Hopf bifurcation points are identified and classified
as subcritical or supercritical and the regions in which small perturbations lead to
stable oscillations are analyzed. A method for easily determining whether a Hopf
bifurcation point of a general system of differential equations is subcritical or super-
critical is developed. The method is used to more accurately identify the regions in

which stable periodic solutions exist in diffusion flames.

For Germ

iii

ACKNOWLEDGMENTS

I owe my deepest gratitude to my thesis adviser, Dr. Milan Miklavcié, for his knowledge,
guidance, and tremendous patience. I truly appreciate all he has done for me. I would
also like to thank the rest of my thesis committee: Dr. Indrek Wichman, Dr. Chichia
Chiu, Dr. Chang Yi Wang, and Dr. Keith Promislow for graciously dedicating their time.
I am incredibly grateful to my family and friends for all of their encouragement and for

never losing faith in me.

iv

TABLE OF CONTENTS

LIST OF TABLES ................................................................................ vi
LIST OF FIGURES ............................................................................... vii
SECTION 1

INTRODUCTION .................................................................................. 1
SECTION 2

MATHEMATICAL MODEL .................................................................... 4
SECTION 3

EXISTENCE OF THE STEADY STATE SOLUTIONS ..................................... 7
SECTION 4

PROPERTIES OF THE LINEAR OPERATOR .............................................. 18
SECTION 5

STEADY STATE SOLUTIONS ............................................................... 25
SECTION 6

STABILITY OF THE STEADY STATE SOLUTIONS ..................................... 29
SECTION 7

HOPF BIFURCATIONS ......................................................................... 35
SECTION 8

ANALYTICALLY IDENTIFYING SUPERCRITICAL HOPF BIFURCATION
POINTS OF A GENERAL SYSTEM ........................................................... 45
SECTION 9

IDENTIFYING SUPERCRITICAL HOPF BIFURCATION POINTS OF OUR
SYSTEM ............................................................................................ 71
SECTION 10

CONCLUSIONS .................................................................................. 78
REFERENCES .................................................................................... 8O

LIST OF TABLES

TABLE 1
Ta vs. values of R where the islands appear and disappear .................................... 28

TABLE 2

Ta vs. values of R where the interval of unstable steady solutions appears and
disappears and values of R at which the bifurcation points at the ends of the interval
become subcritical ................................................................................... 34

vi

LIST OF FIGURES

FIGURE 1

The one-dimensional diffusion flame between two porous walls .............................. 4
FIGURE 2

S-curve (R = O) of steady states when Le = 1 and Ta = 1.2 ................................... 25
FIGURE 3

Transformation of S-curve when Le = 1 and Ta = 1 ........................................... 26
FIGURE 4

Transformation of S-curve for small R when Le = 1 and T a = 1.2 ........................... 26
FIGURE 5

Formation and shrinking of islands when Le = 1 and Ta = 1.2 .............................. 27
FIGURE 6

Disappearance of islands when Le = 1 and T a = 1.2 ........................................... 28
FIGURE 7

Stability of the S-curve when Le = 1 < Lem, .................................................... 31
FIGURE 8

Stability of the S-curve when Le = 5 > Lem-t .................................................... 31
FIGURE 9

Real parts of the three leading eigenvalues at the beginning of the upper branch when

Le == 5 ................................................................................................. 32
FIGURE 10

Stability of the steady states when Le = 1 and Ta = 1 ........................................... 33
FIGURE 11

Stability of the steady states when Le = 1, Ta = 1.2, and R is low ............................. 34
FIGURE 12

Stability of the steady states when Le = 1, Ta = 1.2, and R is high ........................... 34

vii

FIGURE 13
Plot of the difference between the solution of the perturbed problem and steady

solution in which the perturbation leads to extinction .......................................... 36
FIGURE 14

Plot of the difference between the solution of the perturbed problem and steady

solution in which the perturbation approaches a stable periodic solution .................... 36
FIGURE 15

Plot of the difference between the solution of the perturbed problem and steady
solution in which the perturbation grows slowly, but does not approach a stable

periodic solution .................................................................................... 37
FIGURE 16

Stability of the islands when Le = 1, Ta =3, and R is low ...................................... 42
FIGURE 17 .

Stability of the islands when Le = 1, Ta = 3, and R is high ................................... 42
FIGURE 18

Plot of R vs. the perturbation interval size when Le = 1 and Ta = 3 .......................... 44
FIGURE 19

Plot of R vs. the stability coefficient when Le = 1 and Ta = 3 ................................. 77

viii

1 Introduction

Solutions of a boundary value problem arising from a one-dimensional, non-premixed,
film diffusion flame (in which the physics include only diffusion, chemical reaction,
and possibly convection or volumetric heat losses due to radiation) have been inves-
tigated, both theoretically and numerically. When no radiating heat loss is included,
Fendell [4] showed that plotting the maximum temperature of the steady state so-
lution versus the Damkohler number yields an S-curve. Vance et a1. [14] examined
the response of a diffusion flame to small perturbations by studying the eigenvalues
of the linearized system. In particular, Vance and colleagues discussed the impact of
the Lewis number on the stability of the upper branch of the S—curve. Sohn et a1. [11]
investigated when the Lewis number is larger than unity by integrating the conserva-
tion equations numerically. Kukuck and Matalon [8], Cheatham and Matalon [2], and
Kim et al. [6, 7] also explored the stability problem, but all used asymptotics instead
of analyzing eigenvalues. Kim and colleagues predicted pulsating instability when
the Lewis number is slightly greater than unity. Vance et al. [14] used eigenvalues
not only to determine stability, but also to predict when flame oscillations would oc-
cur. Oscillatory behavior in flames was observed in micro—gravity candle experiments
aboard the Mir space station, as well as during droplet burning experiments. Sohn
and colleagues [11] also observed both decaying oscillations and oscillations leading
to flame extinction.

In the later 19905, the effect of heat loss on diffusion flames was studied. Cheatham

and Matalon [1] included a linearized volumetric heat loss term, h(T — To), in the

energy equation. Kukuck and Matalon [8], using the same model, examined what
effect this heat loss had near the upper turn of the S—curve. Sohn and colleagues [12]
included heat loss due to radiation, which introduced a term of the form RD(T4 —T61)
to the temperature equation. This study showed that plotting the maximum tem-
perature of the steady state solution versus the Damkiihler number in this model
yields an island curve, or isola, instead of the usual S-curve. In addition, they found
a region of unstable steady flames which, when perturbed, evolved into stable oscil-
lations. This was the first time such behavior was reported. However, the behavior
seemed reasonable due to remarkably sustainable flame oscillations in a droplet flame
experiment [12]. Christiansen and colleagues [3] included heat loss due to radiation
in a way similar to Sohn et a1. [12], but also included variable properties, full species
multicomponent diffusion, and detailed multistep reaction chemistry. They, however,
were not able to find stable oscillations. They explain that this is due either to the
different geometry or to the fact that the stable oscillations appeared only in “an
extremely narrow regime.”

We will examine a system which includes a radiating heat loss term. We prove
the existence of the solutions of such equations for certain parameter values and
establish many important properties of the linear operator. We will also show how
an isola emerges from the usual S-curve at large activation temperatures, but never
develops for smaller activation temperatures. We demonstrate that when the ratio of
characteristic chemical and radiation time scales is large, only a lower ignition branch
remains. We then identify the parameter values affecting the stability of the steady
state solutions. In addition, we will Show when stable periodic solutions exist and

2

describe a method for easily determining the existence of such solutions.

2 Mathematical Model

We will consider a one-dimensional diffusion flame which lies between two porous
walls. Fuel issues from a large reservoir behind the wall at :1: = —1 and an oxidizer

diffuses from a free stream through the wall at a: =2 1, see Figure 1.

FUEL
FLAME
OXIDIZER

a
"eon-00.0000
N

x=—1

Figure 1: The one-dimensional diffusion flame between two porous walls

The equations governing the flame evolution over time t > 0 can be written as

 

 

8T 82T
E = m+w—Iw(14—751) (1)
2
av 62Y
“Ff = axzf T W’ (3)

Here T, the temperature, is a function of both the spatial and time coordinates, a:
and t, respectively. Y0 and Yf, the oxidizer and fuel mass fractions, respectively, are
also functions of x and t. Le is the Lewis number, taken to be the same for the fuel
and oxidizer, and R is the ratio of characteristic chemical and radiation time scales.

The reactivity term, W, which is a result of the chemistry, is given by

W = DYoYfe—Ta/T, (4)

4 .

where D is the Damkéihler number and T a is the activation temperature. Note that
all of these coordinates and functions have been nondimensionalized exactly as in

Sohn et al. [12]. We assume the following boundary conditions at the porous walls:

T=T0 Y0=0 Yf=1 at x=—1 (5)

T=T0 Y0=1 Yf=O at x=+l. (6)

This model is physically idealized. However, understanding this simple model can
aid in understanding more complicated situations. For this reason, this model has
been studied extensively before. In particular, the equations used by Vance et al. [14]
are identical to the above equations 1 - 6 when R = 0, however Vance’s equations
include a convection term. In our experience, the introduction of this convection term
does not change the resulting stability behavior substantially. The equations used by
Sohn et al. [12] are also equivalent to the above equations when Ta = 5, T0 = 0.1, and
Le = 1. A slightly difl'erent configuration is preferred by Kukuck and Matalon [8].
Christiansen et al. [3] added many complicated real-world influences. For a discussion
of radiative loss, see T’ien [13]. They all provide excellent physical descriptions of
diffusion flames and cite many references.

Throughout the analysis, we shall take T0 = 0.1, as Sohn et al. [12]. Except
for a brief discussion in Sections 6 and 7, we will take Le = 1, also as Sohn and
colleagues. In our experience, the introduction of other parameters, like the Peclet

number of convection, different starting fuel fractions Yf at :1: = —1, and different

wall temperatures To do not change the resulting stability behavior substantially.

3 Existence of the Steady State Solutions

We will begin by showing the existence of the steady state solutions for small values
of the Damkohler number, D. In Section 5 we will find the steady state solutions
numerically. We therefore know that these solutions exist, however, we will prove

existence for only small D values.
( u, I

For u = u2 cC’[—1, 1]3, define the function F(u) to be

\“3/

 

 

/ -Du2U3e-Ta/u1 + RD(u% — T61) \
F(u) : DuZu3e_Ta/u1
K Du2U3eTTa/u1 )

 

 

and the operator N u to be

 

 

 

 

 

 

 

 

(T0\ (Tod
Nu(x) = 0 1;:c_ 1 42—35
l 1 J l 0 1
._$ 17 — -—;1; 1
+1 2 /_1<-1-s)F<u(s))ds+ 1, [x (1—s)F<u<s))ds.

Claim 3.1 u is a fixed point of N if and only if u solves the boundary value problem

(T,\

0

 

 

l1)

U”($) = F W»

Proof 3.1 Let u be a fixed point of N. Hence,

So, u(—1)=

O

 

1—

+ 2

{T0\

0

 

 

\1/

U’(:v) = 0

 

\1
1
2

1

[To]

[1)

1’ L1H —s)F(u(s))ds+

 

 

 

and u(1) 2

MO] I

—1

_—

2

l K

 

 

ITO]

1

 

 

[0)

M0]

1

 

 

[0)

Tot

-1—x

 

ITO]

1

 

 

 

l0}

l—rc

. In addition,

1 ——

 

0}

1—2:

1 _
—-2-L (1 -—s)F(u(s))ds — 2

 

__ L1(—1 — s)F(u(s))ds + 1;—$(—1 — $)F(u(:1:))

(1 - $)F(U($))

1
lg; (l — s)F(u(s))ds.

 

 

0
— 1 —lx——s us s——1-—s uss
— a ,/_1(1 )F<())d [we )F(())d
l
[‘5/
and
Me) = —§<—1—z>F<u(z>)+§(1—x)F(ue»
= Fae»

Thus, u is a solution of the boundary value problem 7. To prove the converse, let u

be a solution of the boundary value problem 7. Then

 

 

 

 

 

 

(Tol (Tol
Nu(:c) : 0 1—2-2:_ 1 —12—a:
l1) l0)

l—a:

:1: _ _ 1
2 /I(-1—s)u”(s)ds+ 12 13/2: (l—s)u”(s)ds

 

 

 

 

 

 

 

 

 

 

(To \ (TO \
1 — 1: ‘1 - 33
= 0 2 - 1 2
t 1 I t 0 i
+1 2 [(—1 — s)u’(s)]__l + /_1 “[(Sldsl
+ 12 [(1— s)u’(s)]: +/a: “[(Sldsl

no] on

_ 1—1:_ —1—33
_ 0 2 1 2

[I] w

 

 

 

 

 

 

 

 

+1; [(—1 — x)u’(a:) — 0 + u(x) — u(—1)]
4; ““10 — <1 - qu’eI + no) - u(rr)l
= u(zr)

and u is a fixed point of N. Hence u is a fixed point of N if and only if u is a solution

of 7.

Fix M > max{T0, 1} and 0 < c < T0. Define

U = {MCI—1,1113% 3 ul 3 M, [122] g M, |u3| g M}.

{u\

For u = 112 cC[—1,1]3, define the norm

['13)

1

 

 

IIUII = 81tp$€[_1,1](IU1(-’B)I + IU2(~’B)I + |u3($)l)-

Claim 3.2 N has a unique fixed point in U if

10

 

D < min{

 

 

M—l M—m
8M2 “Ta/M’ RT4 ’
e _20_ + M26“Ta/M
c — T0 1 }
1274 ’ 2 —Ta/M '
—20— — Mze‘Ta/M —— RM4 48Me—Ta/M + 24M 8 2 T“ +32RM3
6

Proof 3.2 To prove that N has a unique fixed point, we will show that U is a com-

plete, nonempty metric space, N : U —> U, and N is a contraction mapping. First,

note that U is a closed, nonempty subset of the Banach space C [—1, 1]3. Hence, U is a

complete, nonempty metric space. If u =

 

“1

“2

\“3/

 

EU, then clearly NueC[-1, 1]3. To

show that N ucU , we must show that Nu satisfies the necessary bounds. Let xe[—1, 1],

then

 

 

|(Nu($))2|

~1—:I:

+ 2

 

17

|/\

|/\

|/\

|1+$+1_$/x
2 2 -1

(I:
1 + 2/ [Du2u3e—
1

(—1 — s)Du2u3e_Ta/u1ds

= 1+ 8DM2e_Ta/M.

11

1+4-2-DM28“Ta/M

l
f (1 — s)Du2u3e_Ta/u1ds]
T l
a/“1[d3 + 2/ |Du2u3e—Ta/u1|ds
:1:

l
1 + 4/ 1 |Du2u3e—Ta/u1|ds

 

 

 

. M — 1 . .
S'o, [(Nu(:t))2| S M ifD _<_ 8M23_Ta/M. Similarly,
I— :1: 1 — :1: :1: —T /u
[(Nu(x))3| = [ 2 + 2 1H — s)Du2u3e 0 1ds
‘1 " x 1 —T /u
+ 2 (1 - s)Du2u3e a 1ds
:1:

g 1 + 8DM2e_Ta/ M .

M—l

. < .
So, l(Nu(:r))3l s M 1I D —— 8M2e—Ta/M

 

Now, note that

1—1:

 

[31—1 — s)(—Du2u3e_Ta/u1

+RDu‘11 — RDTg)ds

(NU($))1 = T0+

_ _ 1
12 a: f (1 — 8)(-DU2’U3€_Ta/u1 + RDui — R076)“
:1:

x—l 1:2 1
2 (—2———+(a:+1))(—RDT3)

2
x T 4
/ 1(3 + 1)(—Du2u3e_ a/ul + RDu1)ds

+

 

 

35—1
2
1154-1 1 $2
+

4
2 (5 — ‘2— — (1 — 33))(—RDT0)

+

 

 

 

l
+ 1 —T
+16 2 f (s — 1)(—Du2u3e a/ul + RDul11)ds

 

_ x
(1 — 2:2) + ELIE—$2] 1(3 + 1)u2u3e_Ta/u1ds
1
+D—(x—2fl/ (1 - 8)U2U3e-Ta/u1d3
a:

+52(—:—ll/ 1(3 + 1)u111ds + W] (s — 1)u‘11ds.
_ :1:

12

Hence,

Ron} DU—x)

2
(1 — 2:)(1— x)M2e“Ta/M

 

(1 — 2:2) + (:1: +1)(:1: + 1)M2e_Ta/M

|/\

(Nu(:c))1 TO +
D(:I: +1)
+___—_

2

 

12014
= To + 2 0(1— $2) + D(1 — x2)M2e-Ta/M

l/\

RDT4
TO + -—29 + DM2e‘Ta/M.

M—TO

 

So, (Nu(:r))1 S M ifD g R . And,
20 +1‘42e—Ta/M

 

TO + RZTgu — 3:2) — D—(lzli)(z +1)(a: +1)M2e_Ta/M
—i$2+—1)-(1 — x)(1 — :I:)M2e—Ta/}M

RD —1

(N “(13))1

IV

RD(9: + 1)

2 (1 — 3:)(2: —1)M4

(:1: +1)(:1: + 1)M4 +

RT4
= To + D(1 — 2:2)(70 — M2e_Ta/M — RM4)

RT4
2 TO + ”To -— MZe—Ta/M — RM4).

—T
awNmoMZeng ‘ 0

 

4 . Therefore, N maps U into
3T0 2 —T /M 4
—— — e a — RM

2
U.

We will now show that N is a contraction mapping. Let u = ,z =

 

 

13

EU and :rr.[—1,1]. So, 6 S u1,zl g

 

 

M and hence, by the Mean Value Theo-

_ _ _ T
[e Ta/ul—e Ta/ZIISe Ta/Mé—llul—zl].

So,

|u2u3e_Ta/u1 — z2z3e_Ta/zl| =

|/\

|/\

 

|u2u3e—Ta/u1 — 22u3e_Ta/u1
+22u3e_Ta/u1 — z223e_Ta/u1
+2223e_Ta/u1 — 22z3eTTa/zl|
lu3e‘Ta/“1IIU2 - 22|
+|zge_Ta/u1[|u3 — 23]
+|22z3[|e_Ta/u1 - e_Ta/le
[We—Ta/Mlug - 22]
+Me—Ta/M|u3 — 23]

_ Ta
+M28 Ta/Mé—flul — 21]

S K1]]U — lea
. M2 —Ta/MT
where K1 = 2Me-Ta/M + e 2 a. In addition,
6
[ail—2411] 2 [(ul —z1)(u:13+u%zl +ulzf+zil)|

14

S 4M3I’UI - ZII

S KZUU' _ lev

where K2 = 4M3. Hence,

 

llNu—Nzll = “1‘2“” [_xII—I—s)<F(u(s))-F<z(s)>>ds

 

+ ‘12— z [:0 — s)(F(u(s)) - F(z(s)))d8[]

( —Du2u3e—Ta/u1 + RD(u‘% — T61) \

 

 

 

 

 

 

 

 

 

 

 

2 [[1 2$/_x1(_1 - S)( DU2U3e-Ta/"1
K Du2u3e—Ta/u1 )
{ -D2223e—Ta/zl + RD(2’11 — T61) \
— Dz223e-Ta/zl )ds
\ Dz223e—Ta/zl )
. I-Du.u.e-Ta/u1+RDe—Tgn
_12— 27/2: (1- s)( Du2u3e—Ta/ul
K Du2u3e—Ta/“1 )
I + — a) I
_ D22z3e—Ta/Zl )dsll
K Dz223e-Ta/Zl )
= sum-1,1]lll 23 [id—1 _ 3)(‘DU2U36_T“/”1

+RD(u‘11— 751) — (—D22z3e—Ta/ 21 + Roof] — T51)))ds

15

—1—:z:
2

 

l
+ / (1— s)(—Du2u3e_Ta/u1 + RD(-u£11— T6)
:1:

—(—02223e—Ta/ 31 + RD(z‘11— T61)))ds

 

 

_ :13
+]1 2 1:] 1(_1‘ S)(Du2u36—Ta/u1 — Dz223e-Ta/Zl)ds
_ _ 1
+ 12 33/ (1— s)(Du2u3e-Ta/u1 — D2223e_Ta/Zl )ds]
:1:

 

_ a:
+]1 2 13/ 1(_1- 8)(Du2u3e_Ta/u1 — D22z3e‘Ta/Zl)ds

—1—a:
2

1
4/ ] — DuQU3e—Ta/u1

 

I
+ f (1 — s)(Du2u3e_Ta/u1 — Dz223e_Ta/Z1)ds]]
a:

g
+RD(u‘11—- T61) —— (—-Dz2z3e—Ta/zl + RD(zl11 — T61))]ds
+4/11 [Du2u3e_Ta/u1 — DzQz3e_Ta/Z1]ds
+4 /_11 [Du2u36_Ta/u1 — D2223e—Ta/Zl [d3

S 8(3Dlu2u3e—Ta/u1 — z2z3e_Ta/Z1I+ RDlui1 — zill)

S 8(3DKIHU — z[1+ RDK2HU - 2H)

= Cllu - 2H,

where c = D(24K1 + 8RK2). Thus, N is a contraction mapping if

1

D
241g + 23ng

 

1

24MQe-Ta/MTQ
€2

 

48Me—Ta/M + + 32RM3

 

16

Hence, by the Contraction Mapping Theorem, N has a unique fixed point in U if

 

 

 

D < min{ M — 1 M — T0
8M2 -Ta/M ’ 1274 ’
8 -70— + M28_Ta/M
c - T0 1 }
R74 ’ 24 2 “Ta/M T '
—2Q — M2e‘Ta/M — RM4 48Me_Ta/M + M e 2 “ +32RM3
E

17

4 Properties of the Linear Operator

In this section, we will prove that the linear operator has important properties which

we will use in later sections.

We begin by transforming the system 1 - 6 into a system for a perturbation of a

steady state solution. Let the steady state solutions be denoted by T, Y0, and 16f.

Also, assume that T(a:) > 0 for all :rc[—1, 1]. By expressing T, Y0, and Yf as the sum

of the steady state solution and a perturbation, 1 - 3 can be written as

gee) + no, t))

gm + age. t))

~

gm) + u3($,t))

2 ..
gages) + u1 (x. t))

+D(Y0(2:) + u2(:r, t))

(Yfzx) + ”3(37 t))e—Ta/(T($) + "1(331 t))

—RD<(T(x) + u1($,t))4 - T3) (8)
1 62
2735?
—-2—6D(Y0~(:c) + u2(:1:, t))

(Yoix) + age, t))

~

(Yflx) + ”3(x1t))e—Ta/(T($) + u1(:r,t)) (9)

02 ~
L125$7(Yf($)+ u3(513,t))

—-L1;D(v0‘(z) + U2(:13, t))

(Yfix) + age, t))e—Ta/(Tlxl + “ll-”1 0). (10)

In order to linearize our nonlinear terms. we will define the following functions:

we) = D(Yol$) + cage, t))(Yflx) + cage. t))e-Ta/W) + “1WD

18

Z(c) = RD((T(2:)+cu1(2:,t))4— 4).

Then the nonlinear term of 8 is simply W(l) - Z (1) and the nonlinear term of 9 and
—1

10 is f—WU). In order to linearize these terms, we will first express W(l) and Z (1)
e

as the Taylor expansions of W(c) and Z (c), respectively, about c = 0 evaluated at

c = 1. Note that

W(l) = W(0) + W’(0) + éw”(0) + éw”’(0) + (11)
and
Z(1) = 2(0) + z’(0) + éz”(0) + éz”’(0) + (12)
where

W(0) = Drone—Tm
W’(0) = D(u )1? e-Ta/T + one )e-Ta/T + my Me—Ta/T
2 f 3 f T2
_ 7 ~ T "
w"(0) = 2D(u2)(u3)e Ta/T +2p(u2)yf_a:_f('2‘_1)e—Ta/T
~T ~
+2D(U3)Y0———aggl)e_Ta/T
T
~ ~ —2Ta(u1>2 Tim)2 4,)?
+DY0Yf(——:ra— + ‘77—)6
T ~
~ —2Ta(u1>2 Tim)? _T r
+3D(u2)Yf( T3 + T4 )e a/

19

- —2T 2 T2 2 ~
+3D(u3)yo(_a_('_‘1_)_ + fl)e—Ta/ T

T3 T4

 

 

~ ~ 6Ta('u1)3 673(111)?’ Td(u1)3 —Ta,/T
+DY0Yf( T4 — T5 + r6 )8
2(0) = RD(TA—T(31)
I _ ~3
Z(O) — 4RD(U1)T
z”(0) = 121’3’,D(u1)2fl~‘2
m» = 241mm

Using these expressions in 11 and 12 and using the fact that T, 170 and if are exact

steady state solutions, equations 8 - 10 can be written as

Bu

at

where

 

 

= Lu + f(u)
= Lu + N2(u, u) + N3(u, u, u) + 0(u4), (13)
( ul(:l:,t) \
“(117, t) = u2(:r, t) ,
Ku3(x,t)}
( W”(0)l— ZII(O) \
N2(U:U) = % —iW”(0) ,
l

 

 

K 173W”) l

20

{ “III/(0) __ ZIII(0) \
-Ll—ewm(0) ,

1
K _Ewlflm) )

N301, u, U)

curd

 

 

and the linear operator L is given by

:7 + DY},Y}% ‘Ta/ T —4RDT3 DYfe-Ta/ T DYoe—Ta/ T
D- Ta -Ta/T 1 {,2 D ~ —Ta/fi" D —Ta/T
L YoYf T28 Le 61:2 Lere 62126 Yoe
D ~ ~ Ta —Ta/T D 2' —T'a,/7~1 1 D ~ -Ta/T
Le YoYf Tze Le Yfe Le 61-2 Le Yoe

In addition, since the steady state solution satisfies the boundary conditions, the

perturbation must satisfy the zero boundary conditions,

u1 2 u2 = u3 = O at :c = :1:1. (14)

To establish the properties of L needed for out analysis, define the operators

 

 

 

 

[ 32— )
L0 _ 1_6 $02
M5221
K L6 61:2 j
and
T ~ ~ ~ ~ ~ ~
(DYOYf—T—Z— “ e‘Ta/T—4RDT3 DYfe"Ta/T DYoe—Ta/T )
DT ~ ~ T " D ~ ~ D ~ "
L = __ i —Ta/T __ —Ta/T __ —Ta/T
1 Le YoYf T2 6 Le Yfe Le Yoe
D ~ ~ Ta —-Ta/T D~ ——Ta/T D ~ —Ta/T
\ — Le YoYf T28 — Lere — Le Yoe )

21

Then L 2 L0 + L1 and L1 is a bounded operator. Let the domain of L be D(L) =
{u'eAC[—1,1] : u”cL2[—1,1],u(—1) = u(1) = 0}3, a subset of the Hilbert space
L2[—1,1]3.

Before examining L, we will first analyze another linear operator. Define S to be

II

with the domain of S, D(S) = {u’cAC[—1, 1] : u”cL2[—l,1],u(—1) = u(1)= 0}. It is
shown in [9, page 81] that —S is the operator associated with a sectorial sesquilinear
form. In addition, 5' is self-adjoint and has compact resolvent [9, page 83]. It is also

easy to see that for ucL2[—1,1],

 

_ :1: :r 1
5—1” = $21/_1(y+1)u(y)dv+—:l/x(y-1)u(y)dy

and the point spectrum of S is given by

—7r2k2

 

ap(3) = { :k=1,2,3,...}.

Note that

(so 0)

1
L = O—S 0
0 Le

1
\0 0 273/

 

 

22

and D(LO) = D(S)3. Clearly,

 

 

(8‘1 0 0 )
L61 = 0 LeS_1 0
( 0 0 Les—1)

The compactness of S “1 implies the compactness of L61 and hence L0 has compact
resolvent. Therefore, according to [9, page 38], the spectrum of L0 is equal to its
point spectrum and

42,62 _,,2k2
4 ’ 4Le

 

0(L0) = { :k = 1,2, 3, ....}

It is easy to see that L0 is symmetric. Thus, L0 is a symmetric, linear operator
in the Hilbert space L2[—1, 1]3 and zero, a real number, is in the resolvent set of L0.
Therefore, according to [9, page 71], L0 is sclf-adjoint.

Since L0 is self-adjoint, [9, page 73] states that [[(L0 — A)_1|] = dist(/\, 0(L0))”1
for all A in the resolvent set of L0. Hence, since L1 is bounded, we see that there
is a A0 in the resolvent set of L0 such that [[L1(L0 — A0)_1]| < 1. Therefore,
according to [9, page 25], (1 + L1(L0 — )\0)_1)T1 exists. It is also easy to see that
(L — A0)_1 = (L0 — AO)-1(1 + L1(L0 — A0)_1)_1. Since both (L0 — A0)‘1 and
(1 + L1(L0 — /\0)_1)_1 exist, A0 is in the resolvent set of L. Also, the compactness
of (L0 —2\0)—1 implies that (L-—/\())_1 is compact [9, page 31]. Thus, L has compact
resolvent.

Since —L0 is self-adjoint and hence [](—L0 — A)—1|| = dist()\, o(—L0))_1 for all

23

A in the resolvent set of —L0, it is easy to see that, by definition, —L0 is sectorial.
Therefore, since —L0 is sectorial, —L1 is bounded and linear, and D(—L0) is a subset
of D(—B), [9, page 194] states that -—L is sectorial. We therefore may apply what

we know about semilinear parabolic equations to equations of the form

= Lu + f(t,u).

SEIE’

24

5 Steady State Solutions

We find steady state solutions of 1 - 6 numerically by using Mathematica’s NDSolve.
We transform 1 - 6 into an initial value problem by specifying T’(—1), Yé(—1), and
Y;(—1) instead of T(1), Y0(1), and Yf(1). This is done in one of two ways. For a
fixed value of the Damkéhler number, D, we use Newton’s Method to find T’(—1),
Y0'(—1), and Y}(—1) such that T(1) = T0, Y0(1) = 1, and Yf(1) = 0. Alternatively,
for a fixed value of T’ (—1), we use Newton’s Method to find D, Yg(—1), and Yj’.(—1)
such that T (1) = To, Y0(1) = 1, and Yf(1) : 0. After finding a steady state solution,
we continue to use this shooting method to find others. This can be done fast if a
suitable continuation method is employed. We use either continuation in D or T’(—1)
and we must switch, often several times, to draw a curve of steady state solutions.

When radiation is neglected (R = 0), plotting the maximum temperature (Tmax)

versus D gives the classical S-curve, see Figure 2. Including a convection term, as in

Tmax

0.5 *
0.4 >
0.3 >
0.2 -

 

 

‘ ‘ * Log D
6 8 10 12
Figure 2: S-curve (R = 0) of steady states when Le = 1 and Ta = 1.2

Vance et al. [14], also yields this S—curve. When radiating heat loss is included, this

25

S-curve transforms. As R increases from 0, the back of the S-curve begins to push
down. For small activation temperatures (T a S 1), the back of the S-curve pushes

down and eventually flattens down into the ignition branch, see Figure 3.

Tmax
0.275 .

0.25 i
0.225 *
0.2 i
0.175 L
0.15 t
0.125 l

 

 

 

Figure 3: Transformation of S-curve when Le = 1 and Ta = 1

For large activation temperatures (Ta > 1), the transformation of the S-curve is
more complicated. When R is small, the transformation of the S-curve is similar to
the transformation seen for small activation temperatures. The curve looks like the S-

curve with the back of the curve pushed down, see Figure 4. However, as R continues

Tmax

0.5 -
0.4 >
0.3 .
0.2 *

  

 

 

L

‘ ‘ ‘ LogD
2.5 5 7.5 10 12.5 15

Figure 4: Transformation of S—curve for small R when Le = 1 and Ta = 1.2

26

to increase, the plot of the steady states looks quite different than those seen for small
activation temperatures. Sohn et al. [12] showed that with certain radiation values,
this plot gives an island curve, or isola. We see this behavior as well. As radiation
increases and the back of the S—curve pushes down, the curve eventually breaks into

an island and an ignition branch, see Figure 5. As R continues to increase, the island

0.225 i

0.175r

 

 

 

7 8 9 10 I I 12
Figure 5: Formation and shrinking of islands when Le = l and Ta = 1.2

shrinks and eventually disappears leaving only the ignition branch, see Figure 6. With
further increases in R, the ignition branch remains and approaches closer and closer
to To. This transformation occurs when Ta > 1. The R values at which an island
first appears and then disappears are summarized in Table 1 for various activation

temperatures.

27

0.24 .
0.22 .

0.2 '
0.18 l
0.16 ’

0.14 r R = 0.07
0.12 L R = 0.06

r

‘ ‘ * éLogD
7 8 9 10 11 12

R = 0.06

 

 

 

Figure 6: Disappearance of islands when Le = 1 and Ta = 1.2

 

Ta RAppear RDisappear
1.1 0.1111 0.117
1.2 0.05 0.07

2 3.9 x 10—5 1.9 x 10-3
3 3.1 x 10"”9 7 x 10'5
5 1.2 x 10—17 3.1 x 10‘7
6 6.7 x 10-22 2.9 x 10’8

Table 1: Ta vs. values of R where the islands appear and disappear

28

6 Stability of the Steady State Solutions

In this section we analyze the stability of the steady state solutions. A steady solution
is called stable if all solutions of 1 - 6 that start as small perturbations of the steady
state solution decay to the steady state. The steady state is unstable if there exists
a number 1‘ > 0 such that for every number 6 > 0 one can find a perturbation of the
steady solution which is initially closer to the steady solution than 6 yet eventually
differs from it by more than 7'.

We saw in Section 4 that every perturbation, u, of a steady state solution satisfies

the differential equation

5,3 = Lu + f (u),
where —L is sectorial. Hence, according to [9, page 265], the stability of a steady
state solution is determined by the spectrum of L. We also saw in Section 4 that L
has compact resolvent. Therefore, by [9, page 38], the eigenvalues of L have no finite
accumulation point and the spectrum of L is equal to its point spectrum. Hence,
the steady state solution is stable if all the eigenvalues of L have negative real parts.
On the other hand, if L has an eigenvalue with positive real part, the steady state
solution is unstable. Therefore, to determine the stability of the steady states, we

first solve the eigenvalue problem

29

u(ztl,t) = 0

 

 

To solve this eigenvalue problem, we discretized the equations using a second-order
central difference scheme and Mathematica solved the resulting matrix eigenvalue
problem. Increasing mesh size allowed for the discretization error to be kept small.
On the graphs we will use thick curves to denote stable steady states and thin
curves to denote unstable steady states. When the leading eigenvalue is complex,
the curve will be dashed. For example, in Figure 8, the leading eigenvalue is real
and positive on the middle branch. At point 1, the leading eigenvalue changes to
a pair of complex values with positive real part. At point 2, the real part of the
leading eigenvalue pair changes from positive to negative and the leading eigenvalue
becomes real and negative at point 3. The stability of the S—curve varies little with
varying parameter values. The lower branch seems to be always stable and the middle
branch is always unstable. However, Vance et a1. [14] showed that the stability of
the beginning of the upper branch depends on the Lewis number. There exists an
Lecrz't > 1 such that if Le < Lecrit then the entire upper branch consists of stable
steady state solutions, see Figure 7. When Le > Lecrit the beginning of the upper
branch is unstable, see Figure 8. We see in Figure 9 how the interactions of the three
leading eigenvalues determine the stability of the S-curve when Le > Lecrz't' At

point 1, the two positive leading eigenvalues join and form a complex conjugate pair

30

0.5 -
0.4 »
0.3 I
0.2 .

 

 

' * LogD
6 8 10 12

Figure 7: Stability of the S-curve when Le = 1 < Lecrit

 

 

 

Tmax
2r
0.4- I .
0.3 >
0.2 '
* ‘ LogD
6 8 10 12

Figure 8: Stability of the S-curve when Le = 5 > Lecrz't

of leading eigenvalues with positive real part. At point 2, the real part of this pair
becomes negative, implying the stability of the steady state solution. At point 3, the
real part of the complex conjugate pair of eigenvalues decreases below the negative,
constant eigenvalue.

To examine the effect of R on the stability, we set Le = 1. When R is nonzero, the
lower and middle branches seem to remain stable and unstable, respectively, however

the stability of the upper branch is affected by the value of R. Since 1 < Lecrit’

31

Real Parts of Leading Eigenvalues

 

 

 

1,1
0.5
. Lo D
0,, 6.2 6.4 W8 g
_.J[ .3\
—1

Figure 9: Real parts of the three leading eigenvalues at the beginning of the upper
branch when Le = 5

the entire upper branch is stable when R = 0. When R is nonzero, the beginning of
the upper branch is sometimes unstable, like the R = 0, Le > Lecrz't case. We will
discuss when the beginning of the upper branch is unstable only when the behavior
here proves to be interesting when examining the Hopf bifurcation points. The more
interesting result of radiation is the appearance and disappearance of an interval of
unstable steady solutions on the back of the S-curve, which we discuss in detail.
When Ta < .95, the upper branch remains stable for all radiation values. As
R increases and the S-curve flattens out, the unstable middle branch shrinks and
eventually disappears, creating a completely stable curve. But, when Ta 2 .95, we
see the effect of radiation on the stability of the upper branch. When R is small and
the back of the S-curve begins to push down, the stability is similar to the stability
of the S-curve. Then, at a certain R value, an interval of unstable steady solutions

appears on the back of the transformed S-curve. This interval appears for Ta 2 .95,

however the curves depend on the Ta values. Let us first consider when .95 S Ta S 1.

32

In this case, we also see that the beginning of the upper branch becomes unstable at
a certain R value. We saw in Section 5 that the curve flattens as R increases. As
the curve flattens, the regions of unstable steady solutions shrink and eventually the
entire curve becomes stable. An example of this stability behavior is seen in Figure

10 for Ta =1.

Tmax
0.275 ,

0.25
0.225 >
0.2
0.175
0.15 .
0.125 ’

 

 

 

Figure 10: Stability of the steady states when Le = 1 and Ta = 1

When Ta > 1, the interval of unstable steady solutions appears before the trans-
formed S-curve splits into the island and ignition branch, see Figure 11. For example,
when Ta = 1.2, the interval of unstable steady states first appears when R is ap—
proximately 0.0064 and the island does not form until approximately R = 0.05. The
interval persists on the island after the split, see Figure 12. As R increases further,
the island becomes a collection of completely unstable steady states and then even-
tually disappears. The stable lower branch persists. We see this transformation for
all Ta > 1. The first columns of Table 2 give the R values at which this interval of

unstable steady solutions appears and disappears for varying Ta values.

33

 

 

 

 

 

Log D

 

Figure 12: Stability of the steady states when Le = 1, Ta = 1.2, and R is high

Ta RAppear RDisappea'r RLeftSubcritical RRightSubcritical

 

.95 .28 .34 — —
.97 .19 .30 — —
l .12 .25 .202 ~—
1.1 .026 .12 .096 .11
1.2 .0064 .058 .048 .034

Table 2: T (1 vs. values of R where the interval of unstable steady solutions appears
and disappears and values of R at which the bifurcation points at the ends of the
interval become subcritical

34

7 Hopf Bifurcations

A Hopf bifurcation occurs when the leading eigenvalue is complex and its real part
changes sign. At a Hopf bifurcation point, according to E. Hopf [5], a branch of
time periodic solutions splits from the steady state solutions. These time periodic
solutions may be stable or unstable. If the time periodic solutions exist on the side
of the bifurcation point where the steady state solutions are stable, then the time
periodic solutions are unstable and the bifurcation is called subcritical. If the branch
of time periodic solutions exists where the steady state solutions are unstable, then
the time periodic solutions are stable and the bifurcation is called supercritical. In
this case, slightly perturbing an unstable steady solution near the bifurcation point
will lead to a solution which will approach the time periodic solution as t grows. Sohn
et al. [12] found stable periodic solutions on the island curves, implying the existence
of supercritical Hopf bifurcation points of the system 1 - 6 for certain parameter
values. We vary these parameter values to develop a better understanding of when
these stable periodic solutions exist.

One way of determining the type of Hopf bifurcation is to make a small initial
perturbation of an unstable steady state solution near the bifurcation point and solve
1 - 6 directly using a higher order finite difference scheme. Two different types
of flame behaviors may occur on the unstable side of a bifurcation point. Either
the perturbation diverges to extinction, see Figure 13, or the solution approaches a
periodic solution, see Figure 14. If the solution stabilizes and clearly approaches a

periodic solution, as in Figure 14, then the bifurcation is supercritical. However, we

35

dT
0. 04

0. 02

— 0.02
— 0.04
— 0. 06
— 0.08

 

 

t
2.5 5 7.5 10 12.5 15 17.5

Figure 13: Plot of the difference between the solution of the perturbed problem and
steady solution in which the perturbation leads to extinction

dT
0.015

0.01
0.005 ..1.‘ w I‘M

0
—0.005

 

L— t

2 4 6 8

Figure 14: Plot of the difference between the solution of the perturbed problem and
steady solution in which the perturbation approaches a stable periodic solution

may also see the solution diverge to extinction when the bifurcation is supercritical. If
the perturbation is too large or the unstable steady state solution is too far from the
bifurcation point, the solution may not approach an existing stable periodic solution.
For example, there is a supercritical Hopf bifurcation point at approximately D =
1.59 x 107 when Le = 1, Ta = 3, and R = 10—5. By making a small perturbation

near the bifurcation point, say at D = 1.62 x 107, we see the solution approach a

36

stable periodic orbit and hence this stable periodic solution exists. However, if we
make a larger perturbation at this point or move farther from the bifurcation point,
say to D = 1.622 x 107, the solution diverges to extinction. On the other hand, if the
bifurcation is subcritical, there is no stable periodic solution to approach and every
perturbation of every unstable steady state near the bifurcation point will eventually
lead to extinction. When looking for stable periodic solutions, it is important to
search near the bifurcation point and make small perturbations to prevent extinction
when a stable periodic solution exists. Unfortunately, this can be time consuming.
If one is too close to the bifurcation point or the perturbation is too small, then it
takes too long to determine whether a periodic solution evolves. It may appear as if

the solution is periodic when in fact it is just growing very slowly, see Figure 15. For

dT
0.0004 ~
0.0002 n {l

— 0. 0002 . U

 

 

 

 

 

 

 

 

 

 

 

 

 

- 0.0004 *

 

2 4 6 8 10 12 14

Figure 15: Plot of the difference between the solution of the perturbed problem and
steady solution in which the perturbation grows slowly, but does not approach a
stable periodic solution

example, there is a subcritical Hopf bifurcation point at approximately D = 703.64

when Le = 5, Ta = 1.2, and R = 0. An initial perturbation of size 5 x 10—5 leads to

a growth rate of only about 1% in amplitude after the first cycle. A first noticeable

37

increase in the growth rate in this case happens if we start with a 300 times larger
initial perturbation (or continue for about 5000 periods), and it leads to extinction
after about an 800 times larger initial perturbation. These small growth rates near the
bifurcation point make it difficult to determine the type of Hopf bifurcation. One way
of verifying the behavior of the solution is to examine a few oscillations and use the
frequency and growth rate to compute the apparent eigenvalue. If the solution seems
to be behaving exactly as the eigenvalues predict, then the solution will likely diverge
to extinction. However, if the solution seems to grow more slowly than expected,
then the nonlinear terms may be affecting the solution enough to stabilize it. Such
analyses aid in finding possible supercritical Hopf bifurcation points.

Our findings of oscillations persisting for hundreds or thousands of periods contrast
sharply with statements made by Sohn et al. in [12] and in the previous nonlinear
analysis (Sohn et al. in [11]), concluding that the oscillations should be terminated
after a few cycles. However, it is unlikely that an actual experiment would produce
these kinds of oscillations which would last virtually unchanged for hundreds of cycles.
The key reason for being able to produce graphs like Figure 15 is the ability to choose
D very close to the bifurcation value and make small perturbations of the unstable
steady state. For example, in the situation described earlier when Le = 5, Ta = 1.2,
and R = 0, if the Damkéihler number is dropped to D = 700 from D = 703.5, the
growth rate jumps from 0.0004 to 0.0100, i.e. a 0.5% drop in the Damkéihler number
causes a 25 fold increase in growth rate, which shortens the period of persistence of
oscillations by roughly a factor of 25.

When R = 0 we have the S—curve, as we saw in Section 5. Vance et a1. [14]

38

showed that when Le < Lecrit’ the middle branch is unstable and the upper branch
is stable. In addition, although they found flame oscillations when Le < Lecrz’tv they
found no Hopf bifurcation points. On the other hand, when Le > Leerz’t’ Vance and
colleagucs showed that the beginning of the upper branch is unstable and there is a
Hopf bifurcation point where the stability changes on the upper branch. However,
in this case, all perturbations of unstable steady solutions near Hopf bifurcation
points led to extinction. We made many calculations for many different values of
physical parameters, yet we always found the bifurcation to be subcritical when R = 0.
When Le = 1 and R is nonzero, we saw in Section 6 that there are sometimes Hepf
bifurcation points at the beginning of the upper branch. However, we found almost
all Hopf bifurcation points here to be subcritical. There is a very small range of
parameter values for which the Hopf bifurcation points here are supercritical. This
will be discussed below.

As we did in Section 6, we will set Le = 1 to discuss the effect of R on the
bifurcation points. For small R, the stability analysis is similar to the stability of
the S—curve and therefore no supercritical Hopf bifurcatiOn points were found. As R
increases, though, and the interval of unstable steady solutions on the back of the
S-curve appears, supercritical Hopf bifurcations appear at one or both ends of the
interval, see Figures 10, 11, and 12. Supercritical Hopf bifurcation points are circled
in the figures. Subcritical Hopf bifurcation points are not circled. We discuss the
Hopf bifurcations for various values of Ta below.

When Ta < .95, we saw no Hopf bifurcation points. When .95 S Ta S 1, there

are three Hopf bifurcation points which appear. There is a Hopf bifurcation point

39

at the beginning of the upper branch and one at each end of the interval of unstable
steady solutions on the back of the transformed S-curve. When these points first
appear, the Hopf bifurcation point at the beginning of the upper branch is subcritical
and the Hopf bifurcation points at the ends of the interval of unstable steady states
are supercritical. The Hopf bifurcation point at the right endpoint of the interval of
unstable steady states remains supercritical as long as it appears for all .95 S Ta S 1.
For .95 S T a < .98, the Hopf bifurcation point at the left endpoint of the interval of
unstable steady states remains supercritical as long as it appears. In this case, the
Hopf bifurcation point at the beginning of the upper branch becomes supercritical
right before it disappears. This is the only set of parameter values for which we see
a supercritical Hopf bifurcation point in this location and the bifurcation point is
supercritical for only a very small range of R. For example, when Ta = .97, this
bifurcation point is still subcritical at R = .255, however it is supercritical when it
disappears by R = .256. When .98 S T a S 1, the Hopf bifurcation at the beginning of
the upper branch is subcritical as long as it appears. The left endpoint of the interval
of unstable steady states remains supercritical until right before it disappears when
it changes to subcritical. This bifurcation point is subcritical for only a small range
of R. For example, when T a = 1, it is still supercritical at R = .201, however it is
subcritical when it disappears at around R = .204.

When T a > 1, the bifurcations are supercritical when they first appear. However,
as R increases, the Hopf bifurcations first become subcritical and then they disappear.
Also when T a > 1, at a certain R value, the S-curve splits into an island and an

ignition branch. Whether the bifurcations become subcritical before or after the split

40

depends on the Ta value. For lower Ta, values (Ta S 1.2), the bifurcation points
change from supercritical to subcritical before the curve splits into an island and
ignition branch. For example, when Ta, = 1.2 and R = 0.01, supercritical Hopf
bifurcations occur at D = 268973 and D = 843180 on the back of the S-curve.
However, when the island forms at approximately R = 0.05, these bifurcation points
are subcritical. The bifurcation points disappear altogether by R = 0.06. At higher
values of Ta (Ta > 1.2), the bifurcation is supercritical when the island forms. We
saw in Table 1 that when Le = 1 and Ta = 3, an island forms at approximately
R = 3.1 x 10"-9 and there are no longer islands by R = 7 x 10-5. In Figures 16 and
17, we see that there are stable periodic solutions on the islands when R = 10"7 and
R = 10-5. However, when R = 5 x 10-5, although there is still a Hopf bifurcation
point on the island, it is subcritical. When R = 6.4 x 10—5, the entire island is
unstable and hence there are no Hopf bifurcation points on the island. We will
further investigate the R value at which the bifurcation changes from supercritical to
subcritical later. The existence of stable oscillations on islands was first observed by
Sohn et a1. [12]. Sohn and colleagues always found stable oscillations on the islands
since they considered only Ta = 5, a relatively large value of Ta.

Slightly perturbing an unstable steady solution near a supercritical Hopf bifur-
cation point leads to stable oscillations. When far from the bifurcation point, per-
turbations may not lead to stable oscillations. How near the unstable steady state
must be to the bifurcation point in order for a small perturbation to lead to stable
oscillations depends on the bifurcation point. For smaller Ta values, perturbing any
of the unstable steady states in between two supercritical Hopf bifurcation points on

41

 

 

L

12 14 16 18 20 22 24

LogD

Figure 16: Stability of the islands when Le = 1, Ta = 3, and R is low

  

 

R=§-4*!0’5. - - . .
11 12 13 14 15 16 I7 18

 

Log D

Figure 17: Stability of the islands when Le = 1, Ta = 3, and R is high

the back of the curve leads to stable oscillations. For larger Ta values, perturbing
the unstable steady states near the bifurcation points leads to stable oscillations; but
perturbing an unstable steady state in the middle of this interval leads to extinction.
In this case, the size of the interval in which perturbations approach stable periodic
solutions varies. We showed that when Ta > 1, the Hopf bifurcation points on the
back of the S-curve (or island) change from supercritical to subcritical as R increases.

We see that as this R approaches the point at which the type of Hopf bifurcation

42

changes, the size of the interval in which small perturbations lead to stable oscillations
approaches zero. We can examine these perturbation interval sizes to approximate
the points at which the bifurcation points change from supercritical to subcritical.
For example, when Ta = 3, Le = 1, and R = 1 x 10—5, there is a supercritical
Hopf bifurcation point when D = 1.59 x 107. Perturbing steady solutions between
this point and D = 1.62 x 107 leads to stable oscillations, but perturbing steady
states at larger values of D leads to extinction. Therefore, the perturbation interval
size when R = 1 x 10—5 is approximately 3.0 x 105. However, when R increases
to 2.1 x 10"5 the size of this interval decreases to around 3.3 x 103. We estimated
the perturbation interval sizes at several R values between 1 x 10"5 and 2.1 x 10—5
and used those values to approximate a function, T(R), which gives the perturba-
tion interval size at R. We first expressed this function as its Taylor series centered
at R0, the point at which the bifurcation changes from supercritical to subcritical.
We then used the perturbation interval sizes computed at the largest four values of
R : 1.6 x 10—5, 1.8 x 10-5, 2 x 10—5, and 2.1 x10"5 to estimate the third order Taylor
polynomial. Using this approximation of T, we see that the perturbation interval size
becomes zero, and hence the bifurcation point becomes subcritical, at approximately
R = 2.3 x 10-5. This order three Taylor polynomial approximation of T is plotted
in Figure 18, along with the seven perturbation interval sizes computed. The last
two columns of Table 2 give the R values at which the two bifurcation points at the
ends of the interval of unstable steady states become subcritical for various activation
temperatures. However, these values were found by using the method developed in

Section 8, instead of by analyzing perturbation interval sizes.

43

Perturbation Interval Size

1 50000
1 00000
50000

 

 

R times 100000
1.5 2 2.5

Figure 18: Plot of R vs. the perturbation interval size when Le = 1 and Ta = 3

44

8 Analytically Identifying Supercritical Hopf Bi-
furcation Points of a General System

In the previous section we discussed classifying Hopf bifurcation points as either
subcritical or supercritical by perturbing unstable steady solutions near a bifurcation
point and examining the behavior of the solution of the perturbed problem. We
addressed some of the difficulties in using this method and how supercritical Hopf
bifurcation points may appear to be subcritical. In this section, we will develop
a method for determining whether a Hopf bifurcation point of a general system is
subcritical or supercritical. Instead of examining the apparent behavior of solutions,
this method will use the structure of the equations and the eigenvalues to explicitly
state whether stable periodic solutions exist. This method was partially developed by
Renardy in [10], however, his argument depended on the existence of inverse functions
which do not exist at the bifurcation point.

Let H be a Hilbert space of complex, vector-valued function whose inner product

has the property that the inner product of real functions is real. Then, for all U, 2161!,

(0,5) = (Re(u) — i1m(u), Re(v) — z'Im(v))

= (Re(u), Re(v)) + (Im(u), 1772(0)) + i((Re(u), Im(v)) — (Im(u), Re(v)))

 

= (Re(u) + i1m(u), Re(v) + i1m(v))

 

= (11,12).

45

Consider a general system

it = L(A)u + f(u; A), (15)

where, for each t, u(t)cH, A is a parameter, L is a linear operator, and f is nonlinear.
Assume that u = 0 is a steady solution of 15. We will examine real solutions of this
differential equation. Let the domain of f be denoted by D( f) Assume f (u; A) and

L(A)u are both real when u is real and f can be written as
f(u; A) = N2(u,u; A) + N3(u, u, u; A) + O(|u]4),

where N2 and N3 are the quadratic and cubic parts of f, respectively. Note that since
f(u; A) is real when u is real, both N2(u,u; A) and N3(u,u,u; A) are real when u is

real. Assume that f, N2, and N3 have the following properties for all “1 , 112, U3€D( f):

 

1- f(U1)= f(fi)

2- N2(u1,U2) = N2(’U2,U1)

3. N3(u1,u2,u3) = N3(U1,U3,U2) = N3(u2,u1,u3)
4. N2 and N3 are linear in each coordinate.

Note that these properties imply that for u, v, wcH,

1
N2(u,v) = Z(N2(u+v,u+v)—N2(u—-v,u—v)) (16)
l
N3(u,v,w) = 53(N3(u+v+w,u+v+w,u+v+w)

46

+N3(u—v—w,u—v—w,u—v—w)
+N3(—u+v—w,—u+v—w,—u+v—w)

+N3(—u—v+w,—u—v+w,—u—v+w)). (17)

Hence, we see that N2(u, v) and N3(u, v, w) are real if u, v, and w are real. This fact,

along with the linearity of N2 and N3, imply that for u, v, wcH,

N2(fi,5) = N2(Re(u) — i1m(u), Re(v) — i1m(v))
= N2(Re(u), Re(v)) — N2(1m(u),1m(v))

—z'(N2(Im(u), Re(v)) + N2(Re(u), Im(v)))

 

= N2(u, v)

and

N362, ‘0, 277) = N3(Re(u) — iIm(u), Re(v) — i1m(v), Re(w) — i1m(w))
= N3(Re<u>, Ree), Re<w)) — N3<Im<u>.1m<v),Re<w))
—N3<Im<u), Re(v).1m(w)) — N3<Re<u).1m<v).1m(w)>
+.-(N3(1m(u), 1mm. Im<w>> — N3<Im<u>. Ree), Re<w>>
—N3<Re<u).1m<v).Re<w)) — N3<Re(u), Ree), Im<w>>>

= N3(u,v,w).

Now, let us examine 15 near a Hopf bifurcation point. Assume that the system

47

has a Hopf bifurcation point at A = 0. That is, assume L(A) has a complex conjugate

pair of leading eigenvalues 0(A) and 0( ),

Re(a(0)) = 0,

1m(0(0l) 3‘é 0,
d(Rea(A))

0 0.

d, ( ) 7e ,
and the rest of the spectrum of L(A) is in the left half plane. Let 022', where w > 0
denote 0(0). Therefore, iwi is the complex conjugate pair of leading eigenvalues of
L(0) at the bifurcation point. Let (L(A) and a*(A) denote the eigenvector and adjoint

eigenvector, respectively, corresponding to 0(A). So,

L(A)a(A) = a(A)a(A) and L*(A)a*(A) = a(A)a*(A).

Without loss, assume (a(A), a*(A)) = 1. Also, note that

 

mam) = (a*(A).%:—:]—:]a7ii)
= (001.5338)) 0))
= W “Tamra—9)) —_<a_<x>a*(A>:(/\))]
z EXT—T) powwow») — (o<A)a*(A),a(A)>]
= 0 (18)

48

Consider the projection P(A) given by

 

P(Alu = u-(“39*(A))G(A)-(u,a*(A))a(/\)

with the domain of P equal to H. The normalization of the eigenvector and adjoint

eigenvector and 18 imply that

P(A)a(A) = P(A)a(

V
II
CD
A
...;
to
v

Hence,

P(Aflu = P(A>(u—(u.a*(x>>a(x)—<u.a_*T5\7>Z(T))
= P(m - (u,a*(A))P(A>a(A) — (lemmas—(r;

= P(A)u

and P(A) is indeed a projection. We will now find P*(A).

Claim 8.1 P*(A)v = v — (v,a(A))a*(A) — (v,a(A))a*(A).

 

Proof 8.1 If u,veH, then

(P(Aluw) = (”u-(u,a*(A))a(A)-(uam)m,v)

= (nav) - (u,a*(r\))(a(A),v) - (med/Uav)

 

 

 

= (1M) - (u, (9(A),v)a*(A)) - (u, (a(/\),v)a*(/\))

 

: (uav _ (u(Alav)a*(A) _ (a(A),v)a*(A)).

49

Later we will need to know the null space of P*(A).

Claim 8.2 The null space of P*(A) is Span(a*(A),a*(A)).

Proof 8.2 We un'llfirst show that P*(A)a*(A) = P*(A)a*(A) = 0.

P*(/\)a*(r\) = 0*(A)-(9*(A),a(/\))a*(/\)-(9*(A),a( ))a*(A)

 

= a*(A) — 1 * a*(,\) — 0

= 0

and

 

P*(A)a*(A) = a*(/\) - (9*(A),a(/\))a*(>\) - (0*(A),a(A))a*(/\)

 

= m—0—1*a*(A)

=0.

So, Span(a*(A),a*(A)) is a subset of the null space of P*(A). Now, let u be in the

null space of P*(A). Then

P*(A)v = 0

v — (v,a(A))a*(A) - (v,a(A))a*(A) = 0

 

 

v = c1a*(A)+c2a*(A)

for constants c1 and c2. 30, the null space of P*(A) is a subset of Span(a*(A), a*(A))

50

and the null space of P*(A) is equal to Span(a*(A),a*(A)).

Throughout the rest of this section, assume that u is a real solution of 15. Now,

define
z = (u,a*(A)) and y = P(A)u.

Then 2 is a scalar and

 

za(A) + Ea(A) + y = (u, a*(A))a(A) + (u, a*(A))a(A) + u

 

-(u,a*(/\))a(/\) - (u, 9"‘(/\))a()\)

=71.

So, we may write u as the decomposition

u = za(A) + Ea(A) + y where z = (u, a*(A)) and y = P(A)u.

Note that since u is real and za(A) + Ea(A) = 2Re(za(A)), y is real as well. Now,

2 = (u,a;"(A))
: (u,a*(A))

= (L(A)u + f(za(A) + Za(A) + y; A), a*(A))

 

= (L(A)u, 0*(/\)) + (f(zaO‘) + 30W + y; A), 9*(A))

 

= (u, L*(A)a*(A)) + (f(zam + fun) + 0; mm»
= (wTMaVM) + (f(za(A) + 33(7) + y; A). a*(A))
= a<A>(u.a*(A)) + (mam + 2W) + y; arm»

51

= a(A)z + (f(za(A) + z@ + y; A), a*(A))

= 0(A)2 + 90.31;», (20)

where

gm,» = (f(za(A)+2‘Z(T)+y;A),a*(A)).
Similarly,

if = EUE+§(z,y A) (21)

where

9(z,y,/\) = (f(za(A)+3J/\_)+y;r\)am).
And,
y = P(A)u

= P(A)(L(/\)u + f(U; M)

 

= P(A)L(A)(za(A) + Ea(A) + y) + P(A)f(za(A) + 2a( )+ y; A)

= za(A)P(A)a(A) + Eo(A)P(A)a(A) + P(A)L(A)y + P(A)f(za(A) + Tia(A) + y; A)

 

= P(A)L(A)y + P(A)f(za(A) + 300‘) + y; A)

by 19. Hence,

9 = B(A)y+h(z,y;z\), (22)

where

and

h(z,y;/\) = P(A)f(za(A)+§a(/\)+y;z\)-

Claim 8.3 B(A) is singular.

Proof 8.3 Assume that B(A) has an inverse, B(A)—1. Then, noting that P(A)2 =

P(A) since P(A) is a projection,

B(A)B(A)"’1 = 1
P(A)L(A)B(A)_l = 1
P(A)2L(A)B(A)_1 = P(A)
P(A)L(A)B(A)_1 = P(A)

B(A)B<»\>‘1 = P(A),

a contradiction since P(A) is not the identity. Hence, B(A) is singular.

53

Claim 8.4 (B(A) — t))—1 exists for all a in the resolvent set of L(A), except zero.

Proof 8.4 We will denote the resolvent set of B(A) by p(B(A)) and the resolvent set
of L(A) minus zero by p(L(A)) — {0}. We will show that p(L(A)) — {0} is a subset of

p(B(A)). Let a be an element of p(L(A)) — {0}. Let ch and define

h : (L(A)_a)—1P(A)g_(9,aa(A))a(/\)_(9,9;(A))a(/\).
Then

(B(A) — 0012 = Penman) — a)‘1P(A)g — awe) — (Jo—1139);;

_WWWW) + (g. ammo)

 

-(—9—’3§34P<A)L<A>EW + (g,a*(,\))a(,\)
= Fame) — axm) — «no—lam + aP<A><L<A> — arlmm
—a<L(A> — arleg — Waowamm

meumm

 

+(9,a*(/\))0(A) -

 

 

+(9, a*(/\))a(/\)

= P<A>P(A)g + aP<A><L<A> — arlPam — cram - arleg

 

+(9,a*(A))a(/\) + (9, a*(/\))a(A)

= ng - «(L(A) — a)—1P<A>g,a*<A>)a(A)

 

 

—a<(L<A> — a)-1p(,\)g,a*(.\))a(,\)

 

+(9,a*(2\))a(A) + (g, a*(A))a( )
= ng — «may. «L(A) - a)-1)*a*(,\))a(,\)

54

 

_a(p(,\)g, ((L(A) — arlrfiifim‘)

+(g, a*(A))a(A) + (is—(rum (23)

According to [9, page 57/, ((L(A) — a)_1)* = (L*(A) — H)_1. In addition, since

(He) — mam) = W) — mam)

_ a*(A) = (L*(A)—a)‘1a*(A).

 

Similarly, (L*(A) — 6)—1a*(A) = a*(A). Hence, by 23,

 

0(A) — Er"

(B(A)—ooh = P(A>g-—-“—Q(P<A)g,a*a>>a(x)

0(A) -—

 

(P(A)9, a*(/\))a(/\)

O!

779‘) — a
+<g,a*(A))a(A) + (am—AW?)

 

 

= g — (ammo) — (g. a*(A)>a<A> — 75339 — (a. (fume)

0'
(I

-(9, 21_"‘—(/9)a—(l\—), 0*(A))0(A) - 30—) _ 0(9 - (9, Cl"‘(/\))a(z\)

-(9, a*(/\))a(/\), a"‘(A))a(x\) + (9, 9"‘(/\))a(/\) + (amfim

 

 

= g - mwaaflblldbl + mm,a*(A))(a(/\).a*(A))a(A)

 

0.(/\) _ 0(9aml(maa*(Alla(A) _ Em _ a(g,a*(A))m
+———“ (g.a*(A))(a(A),a*a>>Ea‘)

0(A) — a

__" ,a*A aA,a*AaA
1‘00)...” (DH) ())()

oz * or
0(A) _ a(g,a (0)90) + ——0(/\) _

+

 

 

 

 

 

 

 

: g—

55

'm

 

 

 

_m _ a(gva*(’\))a(’\) + ;Cf)‘ _ a

—_ a _—

(9,a*(r\))a(/\)

Thus, acp(B(A)) and p(L(A)) —— {0} is a subset ofp(B(A)). Therefore, (B(A) — 0:)"1

exists for all a in the resolvent set of L(A), except zero.

According to The Center Manifold Theorem, u approaches an invariant manifold.

On the manifold, there are constants {(110)};1: 1 such that

y = dl(A)z2 + dz<A>z2 + dun-22 + 03m; + d4(/\)z2§+ 311752232

 

+d3(/\)E3 + O(|z|4).
Note that some of the coefficients are conjugates due to the fact that y is real. Now,

3] = 2d1(A)Zi + d2(A)§2 + d2(A)Z§ + 2d1(z\)fi + 3d3(A)222 + 2d4(A)Z§2':

+d4(2\)22§ + (14092-22 + 2d4(A)zz_2 + 3d3(/\)‘22? + 0(Izl4).

We can solve for these coefficients {di(A)};1= 1 by setting this expression for y equal
to the expression given in 22. Note that g(z, y, A) = 0(|z|2) and 'g'(z,y, A) = 0(|z|2).

So, by 20 and 21,

B(A)d1(A):«:2 + B(A)d2(A)zE = 2d1(A)o(A)z2 + 2d1(A)zg(z, y; A)
wow—1(1):? + B(A)d3(A)z3 +d2(A)o(A)Ez + d2(A)Eg(z, y; A)

+B(A)d4(A)222 + B(A)d4(A)z"z’2 +d2(A)a(A)z'z' + d2(A)z§(z, y; A)

56

+B(A)E3TA)73

+ h(z, y; A) +2d1(A)o(A)E2 + 201(A)@(2, y; A)

+3d3(A)a(A)z3 + 204(A)0(A)222
+d4(A);(A)z2"z' + 04—(:\-)-0(A)z'z'2
+204(A)a(A)zz2 + 3d3(A)o(A)33

+O(|z|4). (24)

By analyzing f more closely and remembering that N2 and N3 are linear in each

coordinate,

f(za(A) + Ea(A) =
+d1(A)z2 + d2(A)zE
+m§2 + d3(A)23
+d4(A)222 + ERA—)2?2

+3317)? + 0044);»

N2(za(A) + :W + d1(A)z2 + d2(A)zE + (Tl—(T):2
+d3(A)z3 + 04(A)z22 + Ma? + d3(—A)§3 + 00244),
za(A) + Em + d1(A)z2 + d2(A)z‘z‘ + 07(7);2

+d3(A)z3 + d4(A)z2E + (1.1—(A7222 + my;

+0044); A) + N3(za(A) + 2m + d1(A)22 + d2(A)z'z‘
+m‘22 + d3(A)2:3 + d4(A)z2E + mag + @523
+0(|z|4), za(A) + ERA) + d1(A)z2 + d2(A)zE + We?
+d3(A)z3 + d4(A)22§ + ERA—)5? + @523 + 0(Izl4),
za(A) + 2m + d1(A)22 + 02mg: + me? + d3(A)z3
+d4(A)222 + mg? + 33(7):" + 00244); A) + O(|z|)4

N2(a(A), a(A); A)z2 + 2N2(a(A), a(A); A)z§

 

 

+N2(a(A>, am; 522 + (2N2<a(A>, c119); A)

57

 

+N3(a(A). a(A). a(A); A))z3 + (2N2(a(A). d2(A); A)

+2N2(a(A), d1(A); A) + 3N3(a(A), a(A), a(A); A) )2?

 

NI

+(2N2(a(A), d1(A); A) + 2N2(a(A), d2(A); A)

+3N3<a(A).E(_A)3‘JA‘); A))222 + <2N2<a(A).'oTlTA); A)

 

+N3<a(A). a(A). a(A); A»:3 + MW) (25)

By rewriting h and g with this representation of f and setting terms of equal
degree in 24 equal to each other, we can solve for the coefficients {dz-(A)};l___ 1. We
use Claims 8.3 and 8.4 to determine which operators are invertible and which are not,

noting that, except for a(A) and o( ), the spectrum of L(A) is in the left half plane.

The first coefficient d1(A) must satisfy

B(A)d1(A) = 20(A)d1(A)

+P<A)N2(a(A).a(A);A)

d1(A) = —(B<A)—20(A))‘1P<A)N2<a<A).a(A);A) (26)

and d1(A) must satisfy

 

B(A)d1<A) + P(A)N2(a(A). a(A); A) 2a(A)d1(A)

211(7)" = —(B(A) — 2W)‘1P(A)N2(ZI(AAW; A).

58

 

 

 

 

The coefficient d2(A) must satisfy

309420) + 2P(/\)Nz(a(/\)» a(A); A) = 0006120) + 000420)-

This equation cannot be solved as the d1(A) and d1(A) equations since at the bifur-

cation point a(A) + o( ) = 0 and B (A) is singular. So, we must find d2(A) in another

way. Since B(A) = P(A)L(A), this equation becomes

H

P(A)L(A)d2(A) + 2P(A)N2(a(A),a—(A_); A) a(A)d2(A) + Maw)

P(A)(L(/\)d2(/\)+2N2(a(A),a(A);/\)) = (a(A)+0(/\))dz(/\)- (27)

In order for such an equation to have a solution, the righthand side must be in
the range of P(A). Thus, this equation will have a solution if and only if ((a(A) +
8(7))d2(A), 11,-) = 0 for all 1),- in the null space of P*(A). At the bifurcation point,
this will clearly held since 0(0) + T0) = 0, but if (d2(A), 12,-) = 0 for all u,- in the null
space of P*(A), then the equation will have a solution at each A, instead of just at

the bifurcation point. According to Claim 8.2, a basis for the null space of P*(A) is

{a*(A), a*(A)}. If (d2(A),a*(A)) = 0 and (020),?56) = 0, then

P(A)L(A)d2(A) = L(A)d2(A) — (L(A)d2<A).a*(A))a(A) — (L(A)d2(A)fa*—(A‘))E(T)

= L(A)d2(>\) — (6120‘), L*(A)a*(A))a(A) - (d209, L"‘(/\)a*(/\))a(/\)

 

 

= L(A)d2()\) - (€120), a(A)a*(/\))a(z\) - (d2(/\), a(A)a*(/\))a(/\

V

_ _—

= L(206120) - 0000120)), 0*(A))a(A) - 0090120), a*(/\))a(/\)

59

 

Hence, 27 becomes

L(A)d2(A)+ = (a(A)+0(/\))d2(z\)

2P(A)N2(a)(>\), 9(A); A)

<L<A)—(a(A)+a<A)))d2<A) = -2P(A)N2<a(A),a<A);A)
MA) = —2<L(A) — (a(A) + ETADH
P(A)N2(a(A).E(T);A). (28)

a quantity which exists both at and away from the bifurcation point.

3

By setting the coefficients on 2 equal and substituting in the values of d1(A) and

6120),

B(A)d3(A)+P<A)(2N2(a(A),d1(A);A) = 2d1<A><N2(a(A).a(A);A).a*(A))
+N3<a<A),a(A),a(A);A)) +d2(A)(N2(a(A).a(A);A),a—*(A—))
+30(A)d3(A)
MA) = (B(A)—3a<A))—1{—P(A)<2N2<a(A).
—(B(A) — 2or(A))‘1
P(A)N2(a(A),a(/\);x\);x\)

+N3((1(A), (1(0)) a(A); Al)

—2<B<A) — 2o<A))'1P(A)N2(a(A).

60

 

The 223 terms give

B(Ald4(/\) + P(A)(2N2(a(/\)A d2(r\); A)

—.

+2N2<a<A).d1(A);A)

+3N3(a(A), a(A), a(A); A))

61

a(A);A)(N2<a<A).a<A);A).a*(A))
—2<L<A) — (a(A) + WWW)

N2<a(A),a(A);A)

(N2(a(A), 9(A); A), a*(/\))l-

20(/\)al4(«\) + @6140)

+2€11 (”(21)/201(k), 3(7); A), a"‘(/\))

 

+d2(A)(N2(a(/\), a(A); A), a"‘(/\))

+d2(A)(2N2(a(A), a(A); A), a*(*))

 

+2d1(A)(N2(a(A).a(A);am)

(B(A) - (20(A) +a—a)))—1[—P(A)

(2N2<a(A) — 2(L(A) — (a(A) + a( )))—1

P(A)/V2010), a(A); A); A) + 2N2(a(A),

 

—<B<A) — 2a<A))—1P(A)N2(a(A),

a(A); A); A) + 3N3(a(«\), a(A), a(A); A))

a

—2(B<A) — 2a(A))—1P(A)N2<a<A),

a(A);A)(2N2(a(A),W;A).a*(A))

—2(L(A) — (a(A) + o(A)))—1P(A)

N2(a(/\), a(A); A)(N2(a(/\), a(A); A),

2

and the 2? terms give

3mm

+P(A)(2N2(a(/\)A 4190; A)

 

+2N2(a(A), 02

a*(A)) - 2(L(A) —— (a(A) + ETA—n)“1
P(A)N2(a(A).E(_AS;A)(2N2(a(A).
300; MW)

—2(B(A) - arm-1m)

N2( (A) (—a/\ ) W)(N2(a(z\) a(A )A),

 

a"‘(>A))l

a(AfiZA’) + 2W
+2d1(A)(N2(Zm?J‘);A).a*(A))
+d2(A)(2N2(a(A).‘cXA_);A),a*(A))
+d2(/\)(N2(m, 3(7); ALE-(70)
+2m(2)v2(a(i),m;,\),m)

(B(A) — (a(A) + 2m))-1[_p(,\)(21v2(a(,\),
—(B(A) - 23m)‘1P(A)N2(5(_A),E(T);A);A)

+a( )))‘1

P(AlN2(a(/\),a—(x\_); A); A) + 3N3(a(/\),9—(/\_),3(7); A))

+2N2(&(_A)1—2(L(A) — (a(A)

—2(B(A) — 20(A))“1P(A)N2(a(A). a(A); A)(N2(Z("A‘),

a(A);A).a*(A))—2(L(A)— (0 (A A)+o(A)))‘1P(A)

N2(a(/\). WA); A)(2N2(a(/\), a(—/\); A), 9"‘(M)

62

 

N.

—2(L(A) - (a(A) + 397))‘1P(A)N2(a(A).

 

a(A);A)(N2(a(A),a( );A),a*(A))

 

—2(B(A) — 2m)‘1P(A)N2(a_(A_),m; A)

(2N2(a(A), W; A), ETA»)

Lastly, the 33 terms give

 

B(A)d3(A) = 30(A)d3(A)

 

 

+P(A)(2N2(a( )» (11 (A); A) +d2(A)(N2(9_(T), VA); A), 9’"(A))

 

 

 

+N3(a(A),a(A),a(A); A)) +2W(N2(E(A_),m; A).a*(A))

em = (B(A)—3a(A))—1(—P(A)(2N2(a( ).

—(B(A) - 20(A))‘1P(A)N2(a—(A). W); A);A)

 

 

+N3(a(A), a(A), a(A); A))

—2(B(A)—2o_(A)) 1P(A)N2(fl.&(—) A) ‘
(N2(a—(A_), a(_A); A), 75W). E“
Now, by equation 20,
= a(A)Z+(f(m(A)+?(A)+d1(A)Z2+~+(1—3(A)33020 l4 ) A) (A ))

= o(A)z + A1(A)z2 + A2(A)zz + A3(A)22 + B1(A)z3 + 320)}: + B3(A)zE2

63

+B4(A)E3 + 0(|z|4), (29)

where, by 25,
41(A) = (N2(a(A),a(A);A),a*(A)) (30)
42(A) = (2N2(a(A)am;A),a*(A)) (31)

A3(A) = (N2(a(A),a(A);A),a*(A))

B1(A) = (2N2(a(A),d1(A);A)+N3(a(A),a(A),a(A);A),a*(A))

32(A) = (2N2(a(A)»d2(A);A)+2N2(a(A),d1(A);A)

+3N3(a(A), a(A), 9(A); A), 0"‘(A)) (32)

 

B4(A) = (2N2(a( AWN+N3(a(A).a(A),a(/\);A),a*(A)).

To understand the behavior of the solution 2 of 29, we will find a related equation
which is in Poincare' normal form. Thus, we set, for some values {ai(A)}z-7: 1 which

we will choose later,

11) = z + a1(A)z2 + a2(A)zE + a3(A)E2 + cz4[(A)2:3 + a5(A)223 + a6(A):-:E2

+a7(A)E3 + 00:44). (33)

64

 

We would like to write

2 = w + 01(A)w2 + b2(A)ww + 03mm? + b4(A)w3 + b5(A)w2w + 06(A)ww2

+b7(A)m3 + 0(]w|4). (34)
Note that

22 = 1122 + 201(A)w3 + 202(A)w2‘w‘ + 2130mm? + 0(IwI4)

23 = w'w + 010020052 + W102i; + b3(A)w3 + 01(A)w2w
+02(A)ww2 + 03(A)w3 + 0(le4)
22 = a? + 25173??? + 2027520102 + 2337510210 + 0(|w|4)
23 = w3 + 0(le4)
:22- : 10% + 0(lwl4)
232 = wfi2 + O(]w|4)

:3 = m3 + o(|w|4).
By substituting these expressions into equation 33, we see that

w = 11) + w2(b1(A) + a1(A)) + ww(b2(A) + a2(A)) + w2(b3(A) + a3(A))

+w3(b4(A) + 231 (A)))1 (A) + 920%?) + a4(A)) + 102127050)

+2(1L1(A)’>2(A) + a2(A)b2(A) + a2(A)b1(A) + 203(A)b3(A)

Ms (A)) + @2969) + 2a1(A)b3(A) + a2(A)b1(A)

65

+a2(A)b2(A) + 2a3(A)E(‘/fi + a6(A)) + @3079)

+a2(A)b3(A) + 203(A)01(A) + a7(A)) + 0020(4).

In order for this equation to hold, we must have that the coefficients on the terms

of degree higher than one are zero. We can write our {bi(A)}7

z = 1 in terms of the

{ai(A)};-7: 1 to ensure that the necessary terms are zero, although we will not write

those values of {bl-(AH: 1 here. So, we see that we can write 2 as 34 by having our

{bl-(A))z-f: 1 defined in this way. Now,

11) = Z + 2a1(A)zz' + (12(A)25 + a2(A)Ez' + 2a3(A)E? + 3a4(A)z2z' + 205(A)::Ei'

+a5(A)32§ + a6(A)E22 + 2a6(A)z?§ + 3a7(A):~23 + 002(4)

= o(A)z + c1(A)z2 + 02(A)zE + C3(A)z2 + C4(A)Z3 + C5(A)222 + (.~6(A).~.:2
+c7(A)r3 + 00204. (35)
where, by 29,

CM) = A1(A)+2a(A)a1(A)

C29) = A2(A) +a(A)a2(A)+a(A)a2(A)

c3(A) = A3(A)+2a(A)a3(A)

C4(A) = 81(A) + 2al(A)A1(/\) + 02(A)A3(A) + 3a(A)a4(A)

05(A) = 32(A)+291(A)42(A)+02(A)A2(A)+

02(A)A1 (A) + 2a3(A)A3(A) + 20(A)<15(A) + a(A)05(A)

66

06(A) = B3(A) + 2611(A)/13(A) + a2(A)Al(A) + 02(A)42(A)

+2a3(A)A2(A) + a(A)a6(A) + 20(A)a6(/\)

C7(A) = B4(A) + a2(A)A3(A) + 203(A)A1(A) + 30(A)a7(A).
By using 34 to express each 2 in terms of w, 35 becomes

2) = a(A)w + (a(A)b1(A) + c1(A))w2 + (a(A)b2(A) + c2(A))ww + (a(A)b3(A)
+C3(A))m2 + (a(A)b4(A) + 201(A)01(A) + image) + C4(A))’w3
+(a(A)b5(A) + 2b2(A)c1(A) + Wow) + b1(A)c2(A) + 2b:T(A_)C3(A)
+cs(A))w2w + (a(A)b5(A) + 2b3(A)c1(A) + maze) + b2(A)c2(A)

+2b2(A)c-3(A) + c6(A))ww2 + (a(A)b7(A) + b3(A)c2(A) + 2Wc39)

+C7(A))w3 + 0(le4). (36)

Since both the {bi(A)}27= 1 and {Ci(’\)}Z= 1 are expressions of {ai(A)}Z= l’
we can choose our {ai(A)}Z= 1 to make all of the above coefficients on terms of

2

degree higher than one, except the coefficient on the w 117 term, to be zero. At the

bifurcation point, 0(0) + 3(0) 2 0, which causes the a5(0) terms in the coefficient of
w2fi to vanish. Therefore, we cannot eliminate this term in the differential equation

at the bifurcation point. We will not include the values of all the coefficients here,

but instead will give only those values which prove useful. By substituting in the

67

chosen values of {ai(A)};-7: 1, the above equation 36 becomes

11') = o(A)w + (a(A)b5(A) + 202(A)c1(A) + mew) + b1(A)c2(A)

+2b3_(—A_)C3(A) + C5(A))w2-u7 + 0010(4).

Now, the expressions for these b,-(A) and ci(A) values are

 

 

 

 

 

 

 

 

_ 41(A)
b1(A) — 0(A)
A2(A)
b A = -=
2” a(A)
_ A3(A)
b3“) ‘ —o(A)+2a—(A)
341(A)42(A) A2(A)A2(A)
b A = —a A __ —:
5H 5( H a(A)a(A) + a(A)o(A)
@0430) _
(20(A) - a(A))(-0(A) + 2<7(A))
61(A) = -41(A)
42(A)0(A)
C A = ———_—__—
2( ) 0(A) _
_ _ 2A3(A)o(A)
C _ _241(A)42(A) a a _Ai(A)42(A)_A2(A)A2(A)
a 0—_ 2143(A)A3(A)
+ 5(A) (A) _aqua—(E.

Hence, by substituting in these values and using the fact that 11121? = |w|2w,

1b = o(A)w + 0|w|2w + 0(|w|4), (37)

68

where

_ 2 2
—A1(A)A2(A)(20(A)+0(A))_lA2(’\)l +32(,\)+ 2'A3(’\)'_

Ia(A)I2 3(7)‘ a(A) - 20(A)

+a5(A)(o(A) + a(A)).

 

 

fl:

This is in Poincaré normal form. Thus, its stability is fully analyzed in [5]. We know
that the periodic solutions which bifurcate from the steady solution w = 0 are stable

if and only if Re(,0) < 0 at the bifurcation point. Now, at the bifurcation point since

 

 

0(0) = m,
0: A1(0);42(0)z' _ @(Fmfi + 32(0) _ _migigfi
So,
Rem) = Re(A1(O):2(O)i + 32(0))

 

A 0 A 0 '
and the bifurcating periodic solutions are stable if and only if Re( 1( )w 2( )2 +

82(0)) < 0. Let us now find the values needed to compute this value. By 26, 28, 30,

31, and 32,

A1(0) = (N2(a(0),a(0);0),a*(0))

A2(0) = (2N2(a(0),a(0);0),a*(0))

32(0) = (2N2(a(0),d2(0);0)+2N2(m,d1(0);0)+3N3(a(0),0(0),9—(0_);0),a*(0))

69

«(0) = —(B(0)—2iw)‘1P(0)N2(a(0).a(0);0)

d2(0) = —2L(0)—1P(0)N2(a(0).a_(6')‘;0).

Now, remember that

N
I

— w + b1(A)w2 + b2(A)wfi + b3(A)'1'u'2 + b4(A)w3 + b5(A)w2E + b6(A)wfi2
+b7(A)w3 + 0(le4),
_ 2 ~- _—2 3 2— _ —2

y — d1(A)z + d2(A)2.z + d1(A)z + d3(A)z + d4(A)z z + d4(A)zz

+139)?3 + 0(lzl4),

and

 

u = za(A) + 'z'a(A) + y.

Hence, the bifurcating periodic solutions of the system 37 are stable if and only the
bifurcating periodic solutions of the original system 15 are stable. So, the periodic
solutions of 15 which bifurcate from u = 0 are stable if only if Rem) < 0. For this

reason, we call Rem) the stability coefficient of 15.

70

 

9 Identifying Supercritical Hopf Bifurcation Points

of Our System

In this section we will use the method established in Section 8 to develop another
method for classifying the Hopf bifurcation points of the system 1 - 6. We saw in

Section 4 that a perturbation of a steady solution must satisfy

0a

 

 

— = Lu + u
a, f( )
= Lu + N20), u) + N3(u, u, u) + 0(04), (33)
where
( ”1(3)” \
“(37:15) = u2(x, t)
\ “3(3) t) j
and the linear operator L is given by
2 ~ ~ .. ~ ~ - ~ -
62—2- + DYoYf-ge—Ta/T — 4RDT3 DYfe—Ta/T DYoe—Ta/ T
:1:
_2~ ~33. 42.)? is: _ 2~ 42.)? _2- 4m
LeYOY-f T26 Le 332 Le Yf e 62Leyoe
_2 - ~ T9. —Ta/i' -2 ~ —Ta/7" .1._ _ 2 ~ —Ta/T
Leyoy.f ~26 Ler e Le 63:2 LeYOe

71

 

 

Also,

and

where

W”(0)

WIII(0)

{ WII(0) _ ZII(O) \

Nth—t

~20)» = -géwvm)

1
A 721”“) /

 

 

 

 

K WIII(0) _ 2111(0) \
N3(u,u,u) = $- —-L};Wm(0) ,
1 m
A '22W (0) l

— ~ ~ Tofu ) —T T

2D(u2)(u3)e Ta/T+2D(u2)Yf—7—2—1—e 9/
~ T _ ”

+2D(u3)Y0—“—l'2‘—1—)e Ta/T

T
+0“; (4700102 +T§(91)2)e—Ta/tf

0 1 f3 T4

T _ ..

60(u2)(u3)—“T5—‘§1—)-e Ta/T

 

—2:ra(ul)2 + T3(u1)2)e_Ta/r

T3 T4
—2rr..(u1)2 + T3902) —Ta/T

- ~ 6

T3 T4

~~6T 36T23T3u3_~
+DYoYf( a???) _ #17511) + a§61))e Ta/T
12RD(u1)2T2

+3D(u2)Y~f(

 

24RD(u1)3T.

72

And, the perturbation must satisfy the zero boundary conditions,
“1 = u2 = u3 = 0 at x 2 21:1. (39)

Define H = {ueAC[—1, 1] : u'eL2[—1,1],u(—1) = u(l) = 0}3 as a subset of
L2[—1,1]3. It follows from [9, page 81] that H is a Hilbert space of complex, vector-
valued functions. Let [I . || denote the norm in H. Note that u = 0 is a steady solution

of 38 in H. By examining the forms of W”(0) and Z”(0), we see that we can write

)

K 21 S id S 3ai,j"iuj

 

 

N2(u,u) = 21 S iaj S 3bi,juz-uj
A 53131.1 .<_ 301,1"in /
and
21 S M, k S 3di,j, kuiujuk \
N301, 7‘, u) = 21 S i,j, k S 363,j, kuz'UjUk A

 

 

\ Z1 £23231: 3 3f1,j.k"z'"j"k )

where aid" bi,jv Ci,j) 1 S i,j S 3 and di,j,k) ei,j,k) fi,j,lc’ 1 S i,j,k S 3 are

functions of x. Also note that these functions may be chosen so that a,- j. = aj i)

bi.) = biz" 62'.) = Cm) di,j,k = dais.) : disk) 81,231: = 623m” = exile and

fi,j, k = fi,k,j = fjfl, k for all 1 S i,j, k S 3. Therefore, it follows from 16 and 17

73

that,

f 21 :1.) s 301,1“in ]

 

 

”20"”) Z ’31 s M s 3%"in
K 21 s 2'.) s 3%.)“in /
and
f 21 3 my I: 3 3(1),), kUz'ijk ]
N3(“’”’“’) = 21 Si,j,ks3ei,j,k“i"jwk

 

 

K ’31 s 23).): s 3fi.j,k“i"jw1c )

It is easy to see in these forms that N2 and N3 are linear in each coordinate and that
N2(u,1)) = N2(1),u) and N3(u,1),w) = N3(u,w,1)) = N3(v,u,w) for all u,v,wcH.
Also, note that f(17) = T(1),). Hence f, N2, and N3 have properties 1 - 4 in Section 8.

Suppose that 38 has a Hopf bifurcation point at the Damkéihler number D0. Define
our bifurcation parameter, A, to be A = D — D0. Thus, 38 has a Hopf bifurcation
point at A = 0. Note that the linearized system L depends on D and hence L = L(A).
In addition, let us use the same notation for the leading eigenvalue, eigenvector,

projection, and so on, as we did in Section 8. The following theorem is a result of

Section 8.

Theorem 9.1 Assume that 38 has a Hopf bifurcation point at D0 and Rem) < 0

74

where

 

Re([i) = Re(A1(0)32(0)i + 32(0))

and

A1(0) = (N2(a(0),a(0);0),a*(0))

A2(0) = (2N2(a(0),a(0);0).a*(0))
32(0) = (2N2(a(0).d2(0); 0) + 2N2(9_(0_), d1(0);0) + 3N3(a(0). mm 0), a*(0))
d1(0) = —(B(0)-2w)‘1P(0)N2(a(0).a(0);0)

(12(0) = —2L(0)‘1P(0)N2(a(0),a(0);0).
Then there is a time periodic solution 0(x,t) of I - 6, a critical Damko'hler number

Dcritr and an e > 0 such that if Yo is a steady state solution of 1 - 6 at D,

 

 

(r)

where D is between D0 and Dcrit and p(x) is an initial perturbation of 370

K91

 

 

(r)

with ||p(x)|| < c, then Y0 + u(x,t), the solution of 1 - 6 with u(x,0) = p(x),

K5)

approaches 0(x, t) as t —+ 00.

 

 

75

Note that all the values needed for computing the stability coeflicient, such as the
eigenvalue and eigenvector, can be found (or very closely approximated) by Math-
ematica. As in the method for determining the stability of the steady states (see
Section 6), this is done by discretizing and therefore changing the operator L(A) into
a real matrix. Hence, finding the adjoint eigenvector, a*(A), can even be done easily
by finding an eigenvector of the transpose of L(A).

Using Mathematica to compute the stability coefficient, 6, is much more straight-
forward than searching for stable periodic solutions as we did in Section 7. Determin-
ing where a Hopf bifurcation point changes from supercritical to subcritical is now
done with case. For example, in Section 7, we examined when T a = 3 and Le = 1.
We saw that for smaller R values, the size of the interval in which small perturbations
lead to stable oscillations is relatively large. However, as R increases, the size of this
interval decreases and approaches zero as R approaches the point at which the bifur-
cation point becomes subcritical. By analyzing the sizes of these intervals for various
R values, we found that the bifurcation changed from supercritical to subcritical at
approximately R = 2.3 X 10-5. Finding this value of R required much time and
many computations. However, this result is easily verified in Figure 19. We see that
the stability coefficient changes signs, and hence the bifurcation point changes from
supercritical to subcritical, at approximately R = 2.3 x 10—5. It is in this way that
we computed the R values at which the bifurcation points changed from supercritical

to subcritical in Table 2.

76

Stability Coeflicient
100

 

/ R times 100000
1.9 2.1 2.3 2.5

      

— 100L
-200
- 300
- 400

 

Figure 19: Plot of R vs. the stability coefficient when Le = 1 and Ta = 3

77

10 Conclusions

We organize our conclusions by examining different ranges of Ta and explaining the
transformation of the steady state solutions, the stability, and the Hopf bifurcations
as R increases.

When Ta S 1 and R increases from zero, the back of the S-curve pushes down
to T 0 and no islands form. When Ta < .95, the upper branch remains stable for all
values of R. When .95 S T a S 1, the beginning of the upper branch becomes unstable
and an interval of unstable steady solutions appears on the back of the transformed
S-curve as R increases. The Hopf bifurcation point at the right endpoint of the
interval of unstable steady solutions on the back of the curve is always supercritical.
In addition, the left endpoint of the interval of unstable steady states on the back
of the curve is a supercritical Hopf bifurcation point when it first appears and the
Hopf bifurcation point at the beginning of the upper branch is subcritical when it
first appears. When .95 S Ta < .98, the Hopf bifurcation point at the left endpoint
of the interval of unstable steady states remains supercritical as long as it appears.
The Hopf bifurcation at the beginning of the upper branch becomes supercritical
right before it disappears. When .98 S Ta S 1, the Hopf bifurcation point at the
beginning of the upper branch is always subcritical. The Hopf bifurcation point at
the left endpoint of the interval of unstable steady solutions becomes subcritical right
before it disappears.

When Ta > 1 and R is small, the transformation of the S-curve is similar to the

Ta S 1 case. As the Ta S 1 case, an interval of unstable steady solutions appears on

78

the back of the transformed S—curve and the Hopf bifurcation points at the ends of the
interval are supercritical when they first appear. After the appearance of this interval,
the transformed S-curve breaks, creating an island and a lower branch. The lower
branch is completely stable. The island contains a Hopf bifurcation point. When
1 < Ta S 1.2, the bifurcation points change from supercritical to subcritical before
the island forms. When Ta > 1.2, the bifurcation point on the island is supercritical
when the island forms, but changes to subcritical as R increases further. As R
increases, the island shrinks and the lower branch pushes down closer to T0. At a
certain value of R, as the island shrinks, the island becomes a collection of completely
unstable steady solutions and then disappears. The stable lower branch persists.
The method established in Section 8 made determining the type of Hopf bifur-
cation a straightforward process. Although we can often correctly identify whether
a Hopf bifurcation point is subcritical or supercritical by perturbing unstable steady
solutions near the Hopf bifurcation point, it is time consuming and at times it may
appear as though a supercritical Hopf bifurcation point is subcritical. By finding the

stability coefficient, we clearly establish when stable periodic solutions exist.

79

References

[1] CHEATHAM S AND MATALON M, Heat Loss and Lewis Number Effects on
the Onset of Oscillations in Diffusion Flames, Proceedings of the Combustion
Institute, 26 (1996), 1063-1070

[2] CHEATHAM S AND MATALON M, A General Asymptotic Theory of Diffusion
Flames with Application to Cellular Instability, J. Fluid Mechanics, 414 (2000),
105-144

[3] CHRISTIANSEN E W, TSE S D, LAW C K, A Computational Study of Os-
cillatory Extinction of Spherical Diffusion Flames, Combustion and Flame, 134
(2003), 327-337

[4] FENDELL F E, Ignition and Extinction in Combustion of Initially Unmixed
Reactants, J. Fluid Mechanics, 21 (1965), 281-303

[5] HASSARD B D, KAZARINOFF N D, WAN Y H, Theory and Applications of Hopf
Bifurcation, London Mathematical Society Lecture Note Series, 1981 (Cambridge:
Cambridge University Press)

[6] KIM J S, Linear Analysis of Diffusional-Thermal Instability in Diffusion Flames
with Lewis Numbers Close to Unity, Combust. Theory Modelling, 1 (1997), 13-40

[7] KIM J S, WILLIAMS F A, RONNEY P D, Diffusional-Thermal Instability of
Diffusion Flames, J. Fluid Mechanics, 327 (1996), 273-302

[8] KUKUCK S, MATALON M, The Onset of Oscillations in Diffusion Flames, Com-
bust. Theory Modelling (2), 5 (2001), 217-240

[9] MIKLAVCIC, M, Applied Functional Analysis and Partial Differential Equations,
1998 (Singapore: World Scientific)

[10] RENARDY, M, Hopf Bifurcation,

http : //ww . math . vt . edu/people/renardym/ class_home/nova/bif s/node20 . html ,

1998

[11] SOHN C H, CHUNG S H, KIM J S, Instability-Induced Extinction of Diffusion
Flames Established in the Stagnant Mixing Layer, Combustion and Flame, 117
(1999), 404-412

[12] SOHN C H, KIM J S, CHUNG S H, MARUTA K, Nonlinear Evolution of Dif-
fusion Flame Oscillations Triggered by Radiative Heat Loss, Combustion and
Flame, 123 (2000), 95-106

[13] T’IEN J S, Diffusion Flame Extinction at Small Stretch Rates: The Mechanism
of Radiative Loss, Combustion and Flame, 65 (1986), 31-34

80

[14] VANCE R, MIKLAVCIC M, WICHMAN I S, On the Stability of One-Dimensional
Diffusion Flames Established Between Plane, Parallel, Porous Walls, Com-
bust. Theory Modelling (2), 5 (2001), 147-161

81

 

 

 

u“[(1)[((j[][j)(3)21)