COMBUSTION A N D F L A M E 24, 319- 334 ( 1975) 319

A Computational Technique For The Evaluation

,
of Dynamic Effects of Exothermic Reactions

L. M. COHEN, J. M. SHORT and A. K. OPPENHEIM

University o f California, Berkeley, California 94 720

The paper presents a computational technique for the analysis of nonsteady flow fields generated by
exothermic reactions in a compressible medium. This is obtained by numerical integration of the set
of rate equations of chemical kinetics combined with the set of conservation equations of nonsteady
gasdynamics expressed in a Lagrangian form. Under the hypothesis that the power pulse of exothermic
energy is so short that the effects of diffusion, viscosity, and conductivity are, during its effective life
span, negligible, the place where it is generated can be restricted to a discrete Lagrangian cell, and,
consequently, the chemical kinetic set expressed in terms of ordinary differential equations. The com-
putational model devised in this manner is referred to as the exothermic center. Conceptually, it is a
simple flow field, which consists of a kernel, confined within a single Lagrangian cell around the center
where the exothermic reaction takes place, and of inert surroundings through which the pressure wave
generated by the expanding kernel propagates, the two separated from each other by an impermeable
interface.
The computations yield the characteristic features of the power pulse of work done by the kernel
on the surroundings. On the basis of this information, complex flow fields, where exothermic processes
occur, can be analyzed as a nonsteady (diabatic) flow subject to energy (heat) supply that is furnished
at a predetermined rate and given as a function of local thermodynamic state parameters, as demon-
strated in the companion paper presented at the 15th Symposium (International) on Combustion.

Introduction devoid of the dynamic effects that form their

Presented here is a computational technique for salient feature. Basic concepts of systems having
the analysis of the generation of nonsteady flow such properties have been first introduced as "ex-
fields under the influence of exothermic reactions plosive reaction centers" [1 ]. Chemical criteria
in compressible media. Under the hypothesis that governing their existence were studied by van
the power pulse of exothermic energy is so short Tiggelen [2] and, more recently, a stochastic
that the effects of diffusion, viscosity and conduc- theory of their formation has been presented by
tivity are, during its effective life span, negligible, Borisov [3]. Their utility in rationalizing physical
the place where it is generated is restricted to a phenomena occurring in combustion systems was
closed chemical system; a kernel where the exo- demonstrated by the coherence theory of the
thermic reaction takes place, separated by an inter- strong ignition limit [4], a characteristic feature
face, impermeable to mass and energy, from inert of auto-ignition in gaseous mixtures [5, 6]. Finally,
surroundings where the pressure wave generated direct experimental observations of their dynamic
by the expanding kernel propagates. Such a model effects were made in the course of our studies of
is referred to as the exothermic center. exothermic processes in shock ignited gases [7].
Exothermic centers are related to the classical Our technique permits the evaluation of the
concept of "hot spots." The latter were, however, dynamic properties of exothermic centers, mani-
fested by the work their kernels perform on the
*This work was supported by the United States Air
Force through the Air Force Office of Scientific Research surroundings, on the basis of a given chemical
under Grant AFOSR-72-2200, and the National Science kinetic reaction scheme. For the analysis of a com-
Foundation under Grant NSF GK--41614. plete flow field, the power pulse o f this work is

Copyright © 1975 by the Combustion Institute

Published by American Elsevier Publishing Company, Inc.
320 L.M. COHEN, J. M. SHORT and A. K. OPPENHEIM

evaluated for the whole range of initial thermo- quire a definite period of time (the so-called igni-
dynamic conditions that one may expect to en- tion delay or induction period) for the generation
counter, and then correlated with respect to an of chain carriers to usher, via chain branching
appropriate function of pressure and specific steps, the termolecular recombination processes
volume expressing the progress of reaction, while that constitute the major source of exothermic
representing, in effect, the fraction of exothermic energy.
energy it is capable of generating. On this basis, _ _ 1.I1
the whole field, where exothermic processes occur, ~.4' 4
can be analyzed as a nonsteady (diabatic) flow
subject to energy (heat) supply that is furnished
at a predetermined rate, and given as a function of
local pressure and density. Such an analysis, applied
to a particular case of strong ignition in an argon rn rn . ,
diluted hydrogen-oxygen mixture behind a re- 'm 2 . 4 . . . . n
flected shock wave is presented in a companion
paper [8].

The Model
The flow field of the exothermic center is described
on the time-space diagram in Fig. 11. It is pre- 1.4 ' 3+~r

sented there actually to scale for a specific case I i !

computed using our technique. The center consists 0 I 2
of a kernel where the chemical reaction takes place r(mm)
in bulk (i.e., only as a function of time throughout
Fig. 1. Time-space diagram o f an e x o t h e r m i c cen-
its volume, as it has to for a single Lagrangian cell) ter. Presented to scale f o r a 2H 2 + O 2 m i x t u r e in a
and of the surroundings, which behave effectively kernel o f an i n i t i a l radius o f 1.0 m m , w h i c h , at t = 0,
as an inert gas, set in gasdynamic motion by the was suddenly brought to a pressure of 3.236 atm and
expansion of the kernel. The diagram shows the a temperature of 1336°K, while the surroundings
behaved as a perfect gas with specific heats ratio
path of the interface between the kernel and the
")"= 1.14 in a plane-symmetrical flow field (/" =0).
surroundings, and of the boundaries of the first
three Lagrangian cells of the surroundings. In an
actual system, chemically the surroundings are, of Because of the extremely short effective life
course, not inert, but undergo an induction process. span of the center [ 1], as pointed out in the intro-
The outward energetic properties of this process duction, the interface between the kernel and the
are, however, negligible during the time interval of surroundings can be considered as impermeable to
the major expansion of the kernel, and for this both mass and heat. Consequently, the exothermic
reason it can be disregarded in the analysis of the energy is equal to the change of internal energy in
dynamic behavior of the center. The existence of the kernel and is manifested in the form of work
the induction process is certainly essential for this performed by it in displacing the surroundings or,
purpose. Its occurrence in combustion systems what amounts to the same thing, in creating a pres-
represents, on the other hand, one of their most sure wave that propagates through the surroundings
characteristic features, as a consequence of the and sets them in motion.
fact that combustion processes are associated, as
a rule, with chain reaction mechanisms that re- Fundamental Equations
For the present analysis, the substance in the kernel
is considered as a mixture of perfect gases, the rate
l l t should be noted that we use the convention, em-
constants for elementary chemical kinetic steps are
ployed throughout the paper, of indicating instants of
time for Lagrangian coordinates by superscripts and posi- expressed in terms of Arrhenius equations, while
tions in space by subscripts. the gasdynamic flow field is computed in a Lagran-
DYNAMIC EFFECTS OF EXOTHERMIC REACTIONS 321

gian frame of reference. Presented here, respectively, tions for an inviscid gas. In Lagrangian form these
are the equations used for the computation of the are, for mass, momentum, and energy, respectively,
rates of production of chemical species, and for the as follows:
evaluation of gasdynamic variables.
_a_~ : ,, [_!_ a (6)
1. Kinetic Rate Equations Ot rl ~r (rJu)],
For a set o f r ( r = 1,2 . . . . . R) elementary steps
DU_ =_v ~P-, (7)
L kr L Ot Or
air A l ~ ~ blr A 1 , (1)
l=1 l=1
a_e . . . . p _O_v, (8)
the net rate of production of species A t (l = 1 , 2 , Ot Ot
. . . . i . . . . . L) can be written as follows [9] :
where t and r are the time and space coordinates
of the flow field, v is the specific volume, u, the
0 [Ai]- = ~
Ot r=l l L f
(bir - air)k r l~ 1 [Al]alr , (2) particle velocity, p, pressure, and e, internal
energy, while/" = 0, 1, or 2 for planar, cylindrical,
or spherical symmetry, respectively.
where the square bracket denotes concentration
expressed in terms of number density of moles, Computational Scheme
and k r is the reaction rate constant for step r, as According to the postulates of the model, chemical
indicated in Eq. (1). reactions are confined within a single Lagrangian
Since mass, rather than volume, is conserved cell. Consequently, the kinetic rate equations, Eqs.
in the course of the processes under study, in the (4) become de facto reduced to a set of ordinary
analysis, species concentrations have to be ex- differential equations. In the computational scheme
pressed in terms of the number of moles per unit they can be, therefore, de-coupled from the gas-
mass: dynamic equations. For this purpose, one must
P
take k r = const over an elementary step, h. Thus,
Yl = [ A l ] V " (3) in effect, time is allowed to progress only at a con-
stant temperature and specific volume. The funda-
Hence, the appropriate form of the kinetic mental principle of our computational technique
rate equations for our purpose is: is to adopt this feature as the basis for a stepwise
numerical procedure. Thus, at the end of each
ay i - ~ time step, the gasdynamic parameters are read-
0t r=l
l (bir-air)k
'
r/=1
H yl alr
1
, (4)
justed, while the medium is considered to be at
frozen composition, so that appropriate conditions
where of constraint, imposed by the energy balance for
the kernel and the gasdynamic equations for the
!
k r = k r v la, (5) flow field, are satisfied. The process paths of the
elementary computational steps involved in this
while procedure are shown on the Le Chatelier diagram,
Fig. 2, representing an actual sequence of such in-
L termediate steps, plotted to scale, for the example
Ia= 1 - ~, a i r . of Fig. 1.
1=1
As a consequence of the fact that numerical
coefficients in the kinetic rate expressions, Eqs.
2. Gasdynamic Equations (4), differ from each other by orders of magni-
According to the basic assumption of our analysis, tude, variables that exhibit very small changes in
the flow field is governed by conservation equa- comparison to others, may have a commanding
322 L. M. COHEN, J. M. SHORT and A. K. OPPENHEIM

ment that its pressure be equal to that in the first

Lagrangian cell of the surroundings, i.e.,
3.19

p~+l = p nl+V2
+l (10)

"1" 3.17
//j,J / The enthalpy in Eq. (9) is evaluated by the
use of JANAF Tables [13] according to the fol-
lowing expression

L
v J / _ 1 ~ Xl hi ' (11)
Hk PovoM(xl) l-1

where x I is the mole fraction and h l is the molar

3.15
enthalpy o f A l at a temperature
I I I
1.340 1.3411 1.356
Po Vo
® rk - 6~ Pk VkM(Xl) (12)
Fig. 2. The Le Chatelier diagram representing the actual
computational steps taken in the vicinity of instant n
When the pressure matching condition, Eq.
shown in Fig. 1.
(10), is satisfied, Eq. (9) yields the specific volume,
influence on the results of numerical integration.
V~ ÷ 1, and the finite difference relation based on
Equations having such properties are known as
Eq. (6) establishes the motion of the interface,
stiff [9, 10], and in computations they have to
be treated with particular care. Rkn+l -Rn'k At the same time the equation of
We adapted for this purpose the linear multi- state, Eq. (12), gives Tkn+l , and thus all the param-
step scheme developed by Gear [10] and recently eters at the subsequent point become specified. In
formalized by Hindmarsch [ 11] ; its salient features order to make sure that the step-wise integration
are described here in Appendix A. path is sufficiently accurate, we use also, as a step-
The gasdynamie flow field is evaluated by limiting restriction, the condition
means of a Lagrangian finite difference technique
associated with the use of the yon Neumann-
</st,
Richtmyer artificial viscosity according to the
formalism of Wilkins [ 12] ; significant aspects of
this method are summarized here in Appendix B. where specifically, in the numerical example re-
The procedure for matching the chemical kinet- ported here, we have taken/5 T = 10°C.
ic computations with those for the gasdynamic Flow diagrams of the computational scheme
flow field are as follows. are presented in Figs. 3 and 4. The first describes
In terms of symbols, specified in Figs. 1 and 2 the executive program, and the second the com-
and defined in the Nomenclature, the energy bal- putational process for the exothermic center, re-
ance for the adiabatic change of state in the kernel ferred to in the former as EXOCENTER. In each
can be expressed in the following finite difference block of the diagram of Fig. 4, the results are listed
form. in the order in which they are actually evaluated.
The computer program for the exothermic
HIn =_21(Vff+l+ V;)(pn+I
H n +I_.,. pn),k
k (9) center handles also its extreme cases:
(1) the process at constant volume;
while the restriction imposed by the flow field is, (2) the process at constant pressure;
for the kernel, reduced, in effect, to the require- (3) the process along a Rayleigh line.
DYNAMIC EFFECTS OF EXOTHERMIC REACTIONS 323

©
IREAD OR iNITIALIZE I
P ~ INDICES I rich < h c hch = h c
~ ytrs

r
IlilTIgl.
IREAI) INPUT PARAMETERS
[MID DEFINE CASE
I ESTNT
READRESTARTDATA laRCH
j I CHEMICAL
.

Xk
n+l
KINETICS
/hn+=~
~ ~"ch #
L
I-
I
sTART OF I"
GASDY NAMICS
1
IRUNTERMiNATIONIYES
CRITERION MET? I
PRINTAPPROPRIATE
TERMINATIONMESSAGE
l Un+i Rn+l .a+l _n+l
I ' i I V i + ~ t I"i+11r
NO

CEll ADDITION ~
CRiTERiONMET?I
ADD ONE CELL TO I
I FLOWFIELDMESH ADJUSTMENT
L
1"
CefUT~ PR(X:ESSOF
I ESTART ]
kRll'ERESTARTI~,TA
Rn+l
I ,U~
+lit • n+l _n+l
Vi +½, Pi+ i , Hk ,
_n+l

BOTI'B~IC CF.NTER I FOR NEXT RUN n+l n+l _n+l En+l

Pk ' V k , T k i k
I
COMPUTENEXTTIMESTEPBASEl)
ONCOIJR/~ITSTABILITYCRITERION
I _n+l _n - n+l
Ilk - h , I < 8T I " ' ~ ' l "÷l=hcll hcl~/2
i
6'
Fig. 4. Flow diagram of the EXOCENTER routine.
Fig. 3. Flow diagram of the computational program.
Pkn+l-1 _ Pkn - 1
While, of course, none of these !nvolves gas-
l-V; +' 1-v
dynamic computations, all are subject to the energy
balance condition expressed by Eq. (9). Thus the
only difference between them is in the principle Example
for ADJUSTMENT, whose position in the pro- The example used here as an illustration of the
gram is shown in Fig. 4, in that, instead of the computational technique is associated with the
pressure matching condition of Eq. (10), one has acquisition of fundamental background for the
for (1): interpretation of experimental observations of the
dynamic effects of exothermic reactions described
v ; + l = V; ; in a paper presented at the Fifteenth Symposium
(International) on Combustion [8]. There the
medium was 2H2 + 02 + 27A, and the initial state
for (2): corresponded to a pressure of 3.326 atm and tem-
perature of 1336°K established behind a reflected
pn+l shock. Due to a significant dilution by argon, the
k =P~c ; results pertaining to exothermic center were very
close to those for a reaction at constant pressure.
and for (3): For comparison, as well as to emphasize the effects
324 L . M . COHEN, J. M. S H O R T and A. K. O P P E N H E I M

of g e o m e t r y , presented here are results c o m p u t e d t h e r m i c energy o u t p u t manifests itself in the f o r m

for an u n d i l u t e d m i x t u r e o f 2H2 + 02 correspond- of w o r k p e r f o r m e d in displacing the surroundings
ing to the same initial conditions as before. The and is equal to the change in internal energy.
nonreactive m e d i u m in the surroundings is assumed Hence, for the finite difference analysis, the exo-
in this case to behave as a perfect gas w i t h 3' = 1.14. thermic energy i n c r e m e n t is, in accordance w i t h
This value was obtained by fitting a rectangular Eq. (9)
h y p e r b o l a to an equilibrium H u g o n i o t curve for an
2H2 + 02 m i x t u r e , corresponding to the same ini-
AW~c+l/z = ~(Pk
1 n+l n n + l - Vkn )1
+Pk)(Vk
tial c o n d i t i o n s as those q u o t e d above. The hyper-
bola based on this value o f 7 a p p r o x i m a t e s quite
closely the actual C h a p m a n - J o u g u e t states, as well
I
(13)
as the end p o i n t s o f processes at c o n s t a n t pressure
and at constant volume.
The kinetic rate constants are given in Tables 1
and 2, the first specifying the Arrhenius factors and
H n_Hn+l_O
k k
n +On+l
k x
t
!

the second, the third b o d y efficiencies.

Since, b y definition, the kernel o f the exother- and the corresponding e x o t h e r m i c p o w e r is t h e n
mic center undergoes an adiabatic process, its exo- simply

Table 1
Kinetic Scheme for the Hydrogen-Oxygen System
and the Parameters of the Reaction Rate Constants k = A X T n exp(-E ~-/RT).ab

Elementary step A n E Refs.

01. I42 + M ~ H+H+M + (2.08E15) (0.07) (103.83) [14]

- 1.00El5 0.0 0.0
02. O2+M ~ O+O+M + 1.85Ell 0.50 95.56 [15]
- (4.51 E08) (0.85) (- 23.60)
03. H2 + 0 2 ~ OH+OH + 1.70E13 0.0 47.78 [16]
- (2.80E10) (0.36) (28.58)
1. H + O2 ~ OH + O + 1.22E17 -0.907 16.62 [17]
- (3.52El 4) (-0.52) (-0.64)
2. O+H 2 ~ OH+H + 4.76E17 -0.907 16.62 [18]
- (2.71E17) (-0.94) (14.69)
3. OH + H 2 ~ H + H20 + 5.20E13 0.0 6.50 [19]
- (1.41E15) (-0.22) (22.06)
4. H + O2 + M ~- HO 2 + M + 2.00E15 0.0 -0.87 [20]
- (1.78E16) (-0.23) (46.18)
5. H+OH+M ~ H20+M + 7.50E23 -2.60 0.0 [21]
- (4.24E25) (-2.75) (119.40)
6. H+HO 2 ~- OH+OH + 6.00El 3 0.0 0.0 [22]
- (2.30E10) (0.66) (37.58)
7. O+HO 2 ~ OH+O 2 + 1.00El3 0.0 0.0 [22]
- (1.33E12) (0.27) (54.85)
8. OH + OH ~ O + H20 + 1.70E06 2.03 -l.19 [23]
- (8.11E07) (1.84) (16.30)
9. H2 + HO2 ~ OH + H20 + 2.00El 1 0.0 24.0 [24]
- (2.08E09) (0.44) (77.15)

ak is in (cm3/mole)/sec or (cm3/mole)2/sec, T in °K, and E ~ in kcal/mole.

bNumbers in parentheses have been computed to satisfy the relation, k+/k_ = Kc, assuring proper approach to final
equilibrium state.
DYNAMIC EFFECTS OF EXOTHERMIC REACTIONS 325

Table 2
Third Body Efficiencies

Step H2 02 H 0 OH H20 HO 2 A Refs.

01. 2.5 1 20 1 1 6 1 1 [25-27]

02. 1 10 1 30 1 1 1 1 [15]
4. 1 1 1 1 1 25 1 1 [20]
5. 1 1 1 1 1 20 1 1 [21]

/ 2.0
wn+~
rr k
=A
~ ' [4/n+l/z
k I
I hn+l (14) V

1.6
Time profiles of the exothermic energy and
power are presented, for our example, in Figs. 5
and 6, respectively. Included there are results ob- 1.2
tained for the plane symmetrical (/=0), line sym-
metrical (j=l), and point symmetrical (]'=2) exo-

i.o
ta 0 . 8
thermic center 2 , and for processes at constant
volume (II) and constant pressure (P), as well as
the Rayleigh line process (R) corresponding to a o,
Chapman-Jouguet deflagration.
For comparison, corresponding results obtained
for the diluted mixture are shown in Fig. 7. The
data (represented by points) for plane symmetrical
exothermic center are in this case so close to those -0.4
for the constant pressure process, that there is no
I 1 I 1 I
need to consider intermediate cases of other geome- 0 2 4 6 6 le
tries. t(p.s,c)
The fundamental feature of our method is assoc- Fig. 5. Time profiles of internal energy in the kernel.lni.
iated with careful introduction of the artifice of a tial conditions: 2H2 + 02 atPo = 3.236 atm and To =
o 0
perfect gas with constant specific heats whose dy- 1336 K, r k = 1.0 mm. Surroundings: perfect gas with "7
namic behavior, under the influence of exothermic = 1.14. Legend." Vrefers to constant volume;P, to con-
energy would be essentially the same as that of the stant pressure;R, to the Rayleigh line for deflagration;
j = 0, 1 and 2 for plane-, line-, or point symmetrical flow
reacting mixture. For this purpose we use the con-
field, respectively.
cept of heat input into a perfect gas equivalent to
the exothermic energy generated in the reacting whence, taking into account Eqs. (13) and (14) the
medium, and we refer to it as the equivalent exo- equivalent specific power is evaluated from the fol-
thermic energy. Its increment corresponding to lowing finite difference equation
Awkn+ V2 of Eq. (13) is
On+% p~+l+p~ (V~ +1 V~)
k ~ . . . . . . . .
At3n+lA=En p n + l _ 1_ o n O n+l~ (15) 2h n+1
~k k-~k T-1 ( k- k /'
on+l n )
2It should be noted that, as in Fig. 1, all the results + 1_ k -Ok (16)
presented here were computed for the initial kernel radius 7-I hn+l- "

fly= 1 mm. The effect of the size o f the kernel on the

namic behavior o f the exothermic center is outside the
scope o f this paper and will be discussed in a forthcoming On the other hand, the cumulative equivalent
publication. exothermic energy
326 L.M. COHEN, J. M. SHORT and A. K. OPPENHEIM

R (18)

j:z
where
j=l

.= /3=(3'-1)/(3'+ 1). (19)

@.=
Consequently, to express the "reactedness" of

"'t
the exothermic process we introduce as the prog-
ress parameter

o
i f - ~O0 Q
I I I I I I 'I' = - (20)
0 2 4 6 8 10 ~f-~O QT
t(Flec)

Fig. 6. Power pulses of work done by the kernel in where

displacing the surroundings. (All the conditions and
legend are the same as in Fig. 5.)
= (P +/7) (V-/3), (21)

I.I while f d e n o t e s the end equilibrium state and QT

is the total amount of exothermic energy, a unique
function of the end state for given initial condi-
1,0
q tions. The value of qJf, obtained by the hyperbolic
fit to the equilibriurri Hugoniot curve described
0.9 30 previously, for which it was established that 3' =
UJ
1.14 (corresponding to/3 = 0.0654), is 2.22 (cor-
responding to QT = 9.36).
0.8 zo ~"
l, The variation of the progress parameter, qs,
in
E with time is, for our example, shown in Fig. 8.
0,7 I0 .~

1.0
0.6 / ~ o
I I i I i
o 2o 40 60 80 0.8
t (Fsec)

Fig. 7. Time profiles of internal energy and power pulses 0.6

for the case of 2H 2 + 02 + 27A. (All the other conditions
are the same as in Fig. 5. Points represent data for a plane-
symmetrical exothermic center.) 0.4
\p
pn+l + Pflc
Q;+I _ k (Vkn+l _ V 0 ) 0.2
2 k

O
+ _ 3'-1
II (~¢+1_(~¢) (17) ; I I I I l
0 2 4 6 8 10
t(Fsec)
can be considered as a specification of the Hugoniot Fig. 8. Time profiles of the progress parameter. (All the
constant conditions and legend are the same as in Fig. 5.)
DYNAMIC EFFECTS OF EXOTHERMIC REACTIONS 327

The equivalent specific power, {~, evaluated by square and cube roots, denoted by P1 and P2, cor-
the use of Eq. (16), and the parameter ~ , defined responding, respectively, to/---0, 1, and 2.
by Eq. (20), are of central importance to the
analysis of an exothermic center, in that they de-
fine the coordinates of its phase plane, that is a
plane on which a single curve describes all the dy-
namic properties of the system. It is this feature 2.4
EQUILIBRIUM
G ~ H U G O N I O T
that serves, in fact, as the key to the gasdynamic
analysis of nonsteady flow fields associated with 2.0
exothermic processes presented in our companion
paper [8]. a. 1.6
For the present case of an undiluted, stoichio-
metric hydrogen-oxygen mixture, the pulses of 1.2
{2 = {~(q~) are depicted in Fig. 9. They cover the
full range from the Rayleigh line to the constant 0.8
volume process, as shown on the P- V diagram of
Fig. 10, and, in particular, they demonstrate the 0.4
effects of geometry, from spherical (j=2) and cy- l I I I I I I I
1.0 1.4 1.8 2.Z 2.6 3.0 3.4 3.8
lindrical Q'=I) to planar (/=0) symmetry. These v
effects are also quite prominent in the trajectory
of the interface on the time-space diagram, Fig. Fig. I0. Pressure-specific volume diagram of the processes
under study. (All the conditions and legend are the same
11, on the corresponding pressure pulses, Fig. 13,
as in Fig. 5.)
and on the time profiles of specific volume, Fig.
12. Asymptotic values for interface trajectories
in Fig. 11 are given by the solution obtained for
the process at constant pressure. This is depicted
by broken lines representing relative variation of
I 0 r-
i=Z ,:o/,
/j
the specific volume, denoted by Po, as well as its
8-
50

V
---6
0
40 j:O Q
U
:L , /
"~ 3O
N,
j:l

J:Z
P
4 // ~/
;'
v
•o 20 2

I0
0
I I I I
1.0 1.4 1.8 2.2
R
L I I I I I
0 0.2 0.4 0.6 0.8 1.0
Fig. 11. Interface trajectories. (All the condi-
tions and legend are the same as in Fig. 5,
Fig. 9. Phase diagram of the equivalent power pulse in while subscripts 0, 1 and 2 denote, respec-
terms of the progress parameter. (All the conditions and tively, the specific volume, its square root,
legend are the same as in Fig. 5.) and its cube root.)
328 L.M. COHEN, J. M. SHORT and A. K. OPPENHEIM

velop the concept of partial equilibrium as a com-

plement to exothermic center.

j=2 Appendix A
> j=l
Numerical Integration of Kinetic Rate Equations
j=O
For the purpose of numerical integration, the
kinetic rate system of Eqs. (4) may be expressed
l in vector form as
I I I I I I
o 2 4 6 6 IO ~=f(y), (A1)
t(p.sec)
Fig. 12. T i m e profiles o f specific volume. (All the con- where, as a consequence of confining the chemis-
ditions a n d legend are the same as in Fig. 5.) try to a single Lagrangian cell,

/._ _d (y).
- dt

2.0
The eigenvalues of Eqs. (A 1) differ from each
J=O other by orders in magnitude. In any conventional
O. 1.5 J- I integration scheme this causes an uneconomically
large number of steps to be taken in order to take
proper account of the slowly decaying constituents.
This property has become known in the literature
as characteristic of the so-called stiff equations
0.5 [10, 11]. Such equations are treated best by the
method of Gear whose stability conditions are so
I I I I I I
0 2 4 6 8 I0 specified as to permit rapidly decaying constituents
t[Fsac] to be properly accounted for by time steps of rea-
Fig. 13. Pressure pulses. (All the conditions and legend sonable size, assuring thus that the total relaxation
are the same as in Fig. 5.) time is spanned by a manageable number of inte-
grations.
The above results give complete information on The Gear method is illustrated on the mnemonic
all the dynamic features of an exothermic center diagram in Fig. A 1, where it is presented in com-
in a given chemical system under specified initial parison to the predictor-corrector techniques of
conditions and size covering all the possible forms Adams. In particular, the Gear corrector, in con-
that they can exhibit. One should note that the trast to that of Adams-Moulton, is numerically
only parameters of which they are comprised are stable for the relaxation process with step sizes
pressure and specific volume-the dynamic variables in the range of 0 + < h ~< oo. Hence step sizes are
related to force and displacement. These were used constrained only by the accuracy requirements,
also as primary variables in our computations, which are usually much less restrictive than those
while temperature was treated as their function, of stability.
subject to the variable composition of the reacting The details of the Gear method are described
mixture. The latter had to be, of course, evaluated in the report of Hindmarsch [ i 1]. Presented here
as well in the course of our computations, but the are its salient features.
corresponding results are omitted from here for The method is set up in such a way that the
lack of space. We are planning, however, to present order of accuracy is equal to the number of steps,
them in a future publication, on which we are now p. This can be described as follows. The solution
working, where, on their basis, we intend to de- up to and including t n is expressed in terms of a
DYNAMIC EFFECTS OF EXOTHERMIC REACTIONS 329

GEAR PREDICTOR

GEAR CORRECTOR

t ADAMS--BASHFORTH PREDICTOR

ADAMS--MOULTON CORRECTOR

n-p+l n-p+l n-p+S n n+!

r'l Y from interpolation polynomial

0 ~" from Interpolation polynomial
• Y unknown at next time instant
• f unknown at next time instant

Fig. A 1. Mnemonic diagram of the Gear method in comparison to the Adams predictor-
corrector techniques.

set o f p vectors,_Y n-p +i(i = 1 , 2 , . . . ,p), and optimized with respect to the usual criterion for
the vector of derivatives yn. In order to evaluate the truncation error. If the step size is changed, a
_Yn +1 one expands all these vectors into a set of polynomial in time is fitted to the previous solu-
Taylor series with respect to the earliest element, tion and re-evaluated at equal time intervals to
__yn-p +1, and demands that all the terms in this satisfy the fundamental requirement that for a
set containing derivatives of an order o f p and given integration, previous values of _Yat equally
smaller vanish. This yields a set of numerical co- spaced points in time be used.
efficients t~i , for the set of vectors _Y and a value The solution at t n + 1 is obtained by predicting
of the coefficient/3n for the vector of derivatives first an initial estimate, Q_yn+l)o, according to
I;'n . One thus gets the following expression for an Eqs. (A2), and then iterating Eqs. (A3) until a
explicit predictor desired criterion is met. In principle, the mth
iterate, (__yn~1)m' is determined from the (m-1)th
iterate, (yn 1)m -1' by the use of an iterative
yn+l= ~ ai Yn-p+i +hn+l B yn (A2)
--
i=1 -- ' n - - ° form of Eqs. (A3):

The corrector is based on the same principle, P ot;yn-p+i n+l * n+l

except for the use of the derivative at t n+l instead
(-yn+l)m =i~1 +h /3n+l(f )m-1
of that at t n , yielding the value of the coefficient (A4)
/3n*+1 for the vector _fn + 1. Thus the implicit cor-
rector is given by the expression Specifically this is accomplished by the use of the
vector of the difference error functions
P ,
_yn+l = ~ oti _yn-p+ i +hn+l~*-n+l_J'n+l _~_n+l)m --- ( y n + l ) m _h n+l/3:+l(fn+l)m
,'=1 (A3)

for which the exact values o f f n+l are specified by -

P
Z
* yn-p+i ,
~,_ (A5)
the given set of Eqs. (A 1). i=1
The multistep scheme is started by holding the
order at p = 1 until a sufficient backlog of informa- which should vanish if the differences between
tion is stored. Then, the order and step size are the mth and the (m-1)th iterates are annihilated.
 330 L. M. COHEN, J. M. SHORT and A. K. OPPENHEIM

Thus one obtains Appendix B

_(gn+l)m_l = (yn+l)m - (yn+l)m_l Numerical Integration Scheme

for Gasdynamic Equations.
The conservation system described by Eqs. (6,
+hn+lfl:+l [(fn+l)m_ _(fn+l-~v-'m-l]" (A6) 7, 8) may be cast into the following nondimen-
sional form:

The mean value theorem expressed in terms of the or= v a (#v)],

Jacobian matrix, J(jplj), where (B1)

Jolt = ~:+1 (3fp/3Yp)m-1 (A7) _Off_= _ V a P , (B2)

3T 3R
permits Eqs. (A6) to be transformed into the
Newton-Raphson form: _~E__=_p a V , (B3)
Or Or
cyn+l)m = (__yn+l)m_l - IS] ~_1 (gn+l)
m_ _ -1,(A8)
As a consequence of differentiation along the
particle path, the mass velocity is given by the
where relation

U = O~R. (B4)
[S]m_l - [i+hn+lj], Or
while I is an L × L identity matrix; L expressing Rapid changes of state produced by gasdynamic
the number of constituents. discontinuities are accommodated within a rela-
The iteration is terminated when, using the L 2 tively narrow region by replacing P with P + II in
norm, Eqs. (B2) and (B3) where

I OU dR + c°~ ( -3u
R d R ) 2]
n = 1_
v [ q 4-o
(-zn+l)m 2 ~< e, (A9)
(B5)
is the von Neumann-Richtmyer artificial viscosity
as modified by Wilkins [12]. The coefficents CL
while the components of (Z n+ 1)m are
and Co are constants of an order of unity, selected
so as to minimize numerical instabilities.
(Zl)m =-hn+1 [(f~+l)m - (f~+l)m-l]/(Yl)max The computational grid for the finite difference
(AIO) scheme is shown in Fig. B 1.
In order to maintain proper relationship be-
where tween the derivatives, as demanded by the conser-
vation equations, velocity is evaluated, as a rule,
0
(Vl)max ------max Lv7, Yt'
1 n
Yl' Y~ . . . . . Yl ] (A 11) at full steps in radius, that is at cell boundaries,
and half steps in time, i.e., one may have Uin+- 21,
and, in our case, 10"s -<-e ~< 10-4 . but not un_+1_. On the other hand, all the thermo-
One should note the inclusion ofy~, the equi- 2
librium value of component l, in Eq. (A 11). This dynamic parameters are evaluated at full steps in
is done to avoid complications that may arise if time, that is at given instants, and half steps in
y~------0(/=0,1,2 . . . . . n). radius corresponding to conditions in the middle
DYNAMIC EFFECTS OF EXOTHERMIC REACTIONS 331

(U~ <P~ tUJ tp) UJ

feel
P~-}
n+l tP~

uZ/ + u.".'~
4 u~":t
C~ ql

I"I

u;:~, . + u;-* + u~ tU~

IF-JI
o;:i ¢P~

I
,;.+~ I+ Z R

Fig. B 1. Computational grid for the Lagrangian finite difference scheme

of gasdynamic analysis.

of a cell, i.e., one may have pn+_½' V ni+ ½ and A~ - [pn + Ii n-½1Ji+½ _ [en + II n-½ ] i_½,(B 7)
En+ l/~,but not P7 +-½, vn +-l/2,or En +-l/~.How-
ever, since II is a parameter involving both the velo- FR n R n R n _R n -
city and pressure, it is evaluated at both half steps ~b~_l 1~ i+1-
~.-- i + i i-1 , (BS)
in space and in time. One can have, therefore, only 2 L i+½ vni-½
II n -+'/2 i.e., at points denoted by crosses in Fig.
i+_lA ,
B1.
The computational scheme is explicit if one IIn-½= 1 ( 1 + ) (CIB+C2oC),
treats first the equation of motion, Eq. (B 2), fol- i+½ 2 Vn V n-1 "
lowed by the kinematic equation, Eq. (B4), con- i+-'½ i-½ (B9)
tinuity equation, Eq. (B 1), and energy equation,
Eq. (B 3). Their finite difference forms are pre-
sented then here in this sequence, including the
stability criterion for the evaluation of the time
B= ¢( )( n-1½
Pi+ n ½
Vl~+_ +v.-,,
i+½} AU ,(B10)

step.

Equation of Motion c = (Auy, (U 11)

In finite difference form, Eq. (B2) is expressed as
while
follows:
+ F,,,,-½ '/~] (B12)
~u-=_ L~,_+l - u ; -
uin+½ = uin-½ _ A m A ~ / ~ , (B6)
It should be noted that the artificial viscosity
is required only to smooth down the effects of
where Ar n is determined below, Eq. (B 28), and excessive pressure gradients that appear in the flow
332 L.M. COHEN, J. M. SHORT and A. K. OPPENHEIM

field as a consequence of shock waves, that is rapid

R n+lA=! [Rn+l+Rn l, (B19)
changes of state associated with compression and z 2
velocity decrease. Thus in the computations, one
introduces the condition
x = ]q-~ [a~" +vq [(uT**lV2y - (u~"+Y2)~1,
24
I n-½
i+½
= O, (B20)

while, as before,] = 0, 1, or 2 depending upon

if whether the flow is plane-, line-, or point-sym-
metrical.
Vn ~ Vn- 1 (B 13)
i+-½ i+y2, Energy Equation
In finite difference form, Eq. (B 3) is expressed as
follows:
or
E n -W
AU>O En+l _ i+V2 (B21)
i+ V2
1 +(7-1)AV/V n+l
/ i+V2
Kinematic Equation
In finite difference form, Eq. (B 4) is simply where

R n+l = R n + U F + % A r n + ½ , (B14) 1 pn (B22)

w= [-2 + ~]i+w A v,
where
A V = Vn+1 _ Vn (B23)
i+lA i+*A,
A r n + Y z = r n + l - r n. (B15)

Continuity Equation
Following Wilkins [12], in order to maintain ~ = 1 [Iin + l/2 + l-ln- g2 (B24)
2 t i+V2 i+% ]'
proper mass balance, this equation is treated with
particular care. Thus Eq. (B 1) is expressed in finite
difference form as follows:
En+l=on+l / (B25)
i + '/2 i + V2 (7-1),

A T n + 1/2
vn+l = V n + - - [A(URJ)+X],(BI6)
i + lA i + lA Mi + V2
and

®n+l=pn+l vn+l (B26)

where i + V2 i + V2 i +1/2

[R°+ll ]+1 - [ROl 1+1

= , (B17) Stability Criterion
Mi+V2 (]+ 1) V°+V2
The size of the time step is adjusted so as to com-
ply with the Courant condition, which, according
to Wilkins [12] is expressed as follows:
A(UR]) = Vi+lHn+V2 t[Rn+g2]
j-i+1 uin+lAJR7+IA]]'
(B 18) A r n+l/2 = min [Ar~ +'/2 , 1.4Arn-V2], (B27)
DYNAMIC EFFECTS OF EXOTHERMIC REACTIONS 333

5. Meyer, J. W., and Oppenheim, A. K., Thirteenth

whence
Symposium (International) on Combustion, The Com-
bustion Institute, 1971, pp. 1153-1164.
Arn=! [Arn+½ + Ar n-½] (B28) 6. Vermeer, D. J., Meyer, J. W., and Oppenheim, A. K.
2 Combust. Flame 18, 327-336 (1972).
7. Meyer, J. W., Cohen, L. M., and Oppenheim, A. K.,
Combustion Science and Technology 8, (4) 185-
In the above 198 (1973).
8. Oppenheim, A. K., Cohen, L. M., Short, L M., Cheng,
An+½ =2 [ ARn+V2 ] R. K., and Horn, K., Dynamics of the Exothermic
Process in Combustion, Proceedings of the Fifteenth
min over all i's
Symposium (International} on Combustion, Tokyo,
(B29) Japan, August 1974 (in press).
9. Curtiss, C. F., and Hirschfelder, J. O., Proc. of the
Nat. Acad. of Sci., 38 (3) 235-243 (1952).
where 10. Gear, C. W., Numerical Initial Value Problems in
Ordinary Differential Equations, Prentice-Hall, Inc.,
1971, XVII + 253 pp.
A = X/-7 o n (B30) 11. Hindmarsch, A. C., Ordinary Differential Equation
i+½ '
System Solver, Lawrence Livermore Laboratory Re-
port # UC1D-30001, Rev. 2, Livermore, California,
III+ 40 pages, 1972.
B = 8Co v~ v, (B31) 12. Wilkins, M. L., Calculation of Elastic-Plastic Flow,
Lawrence Livermore Laboratory Report #UCRL-
7322, Rev. 1, Livermore, California, IV + 101 pages,
1969.
13. Stull, D. R., Prophet, H., et al., JANAF Thermochemi-
ARn+ ½ = RT+ 1 - R 7 , (B32)
cal Tables, U.S. Department of Commerce, NSRDS-
NBS 37, June 1971, 2nd Edition.
14. Mallard, W. G., and Owen, J. H.,Intern. J. Chem.
V = ( Vnt+½ -vn-1)/AT"n-V2'i+½ (B33)
Kinetics 6,753 (1974).
15. Watt, W. S., and Myerson, A. L., J. Chem. Phys. 51,
1639 (1969).
16. Jachimowski, C. J., and Houghton, W. M., Combust.
Flame 17, 25 (1971).
~_ - - 1
2
[ V i+½
n + v n-lq
i+l/2j (B34) 17. Schott, G. L., Combust. Flame 21,357 (1973).
18. Schott, G. L., Getzinger, R. W., and Seitz, W. A.,
164 Amer. Chem. Soc. Meeting, New York, 1972.
19. Gardiner, W. C. Jr., Mallard, W. G., and Owen, J. H.,
One should note that, as a consequence of the J. Chem. Phys. 60, 2290 (1974).
condition of Eq. (B 13) that annihilates the effects 20. Gutman, D., Hardwidge, E. A., Dougherty, F. A.,
of artificial viscosity for rarefactions, in the above and Lutz, R. W.,J. Chem. Phys. 47, 4400 (1967).
21. Homer, J. B., and Hurle, I. R., Proc. Roy. Soe. (Lon-
don) A314, 585 (1970).
B=0 whenever V~>0. (B35)
22. Browne, W. G., White, D. R., and Smookler, G. R.,
Twelfth Sym. (Int'l.) on Combustion, The Combus-
tion Institute, 1969, p. 557.
References
23. Rawlins, W. T., and Gardiner, W. C. Jr., J. Chem.
1. Zajac, L. J., and Oppenheim, A. K., AIAA Journal Phys. 60 (12) 4676--4681 (1974).
9 (4) 545-553 (1971). 24. Voevodsky, V. V., Seventh Symp. (lnt'l.) on Com-
2. Van Tiggelen, P. J., Combustion Science and Tech- bustion, Buttersworth, 1958, p. 34.
nology 1,225-232 (1969). 25. Jacobs, T. A., Giedt, R. R., and Cohen, N., J. Chem.
3. Borisov, A. A., On the Origin of Exothermic Cen- Phys. 47, 54 (1967).
ters in Gaseous Mixtures. Paper presented at the 26. Getzinger, R. W., and Blair, L. S., Combust. Flame
Fourth International Colloquium on Gasdynamics 13,271 (1969).
of Explosions and Reactive Systems, Acta Astro- 27. Gay, A., and Pratt, N. H., Procs. Eight lnt'l Shock
nautica 1(7-8 909-920 (1974), Tube Sym., Chapman and Hall, 1971, 13 pp.
4. Meyer, J. W., and Oppenheim, A. K., Combust.
Flame 17, 65-68 (1971). Received 8 August 19 74
 334 L.M. COHEN, J. M. SHORT and A. K. OPPENHEIM

Nomenclature T temperature
u flow velocity
Symbols u/v%o vo
a stoichiometric coefficient of a constituent of v specific volume
reactants V V/Vo
A chemical species x mole fraction
b stoichiometric coefficient of a constituent of
(7-1)/(7+1)
products
7 ratio of specific heats
e specific internal energy
® PV
E e/povo l-I artificial viscosity
F end equilibrium state of a constant pressure
r t/to
process (,°+~) (V- ~)
G end equilibrium state of a constant volume
(~ -¢o)/(~i-¢o)
process
h computational time step; molar enthalpy
H h/po Vo Subscripts
/ geometrical parameter (0 for planar, 1 for c Courant condition
cylindrical, and 2 for spherical symmetry) ch chemical kinetic
k reaction rate constant f end equilibrium state
K end equilibrium state for the Chapman- i inner boundary of the ith cell of the sur-
Jouguet deflagration roundings
L total number of species i + 1_ thermodynamic parameter of the ith cell
M molecular weight 2 of the surroundings
P pressure k kernel
P p/po 0 initial state
r radius r elementary reaction step number
R r/r°k; number of elementary reaction steps
universal gas constant
Superscripts
t time
0 initial time
to r°k/X/po vo n time instant