Competition Effects in Combustion Chemistry and Th

and Their Role in Detonation Initiation∗

Shannon Browne†

California Institute of Technology

Pasadena,CA91125

May2004

Abstract

This study considers the issues of approximate modeling of chemical kinetics in detonation waves and the prediction of dynamic detonation parameters such as initiation energy.Several models

have been proposed for the initiation of detonation waves.Due to the complexity of the coupling

betweenfluid mechanics and chemistry,these models have relied on simple single-step reaction

mechanisms.Recently Short,Quirk,and Kapilla have proposed a three-step chemistry model that

includes initiation,branching,and termination components.They propose that initiation and failure

can be related to competition between branching and termination and the existence of a“crossover

temperature.”This chemistry model is a mathematical model and proposes to reduce a complicated

chemical mechanism to three generic reactions that do not depend on specific chemical properties.

A program of research is proposed to examine the validity of their model by comparing detonations

based on detailed chemistry with detonations based on their3-step model using realistic values of

the input parameters.

1Introduction

Direct initiation refers to any method of detonation initiation that uses a strong shock wave.Commonly, high explosives or electrical discharges provide the necessary energy to create a strong spherical blast wave in the test medium.If the exothermic chemical reactions and the wave remain coupled as the blast wave expands and decays,then a detonation has successfully been initiated.The coupling of the shock wave and exothermic chemical reactions involves both gas dynamics and chemistry.This study initially will focus on the role of the chemical reaction mechanism in determining the critical conditions for direct initiation success.

The combustion process in a detonation wave can be conceived of as a convected chain-branching/thermal explosion.The crucial feature of high temperature combustion chemistry that may play a role in this problem is the competition between branching chain reactions and chain termination reactions.The historical viewpoint[10]is that if termination reactions proceed more quickly than branching reactions, an explosion can be quenched.This effect may play a role in causing the detonation initiation process to fail.The critical parameter controlling this behavior can be expressed as a crossover temperature where the chain-branching and chain-termination reaction rates are equal.A simplified view of the situation is that above this temperature,branching dominates,while below,termination may quench the reaction. The situation with realistic combustion chemistry is more subtle and although chain termination may alter the explosion time,peroxide chemistry may ultimately still enable an explosion to occur[16].

∗Submitted in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy

†Candidacy Committee:Dr.Chris Brennen(chair),Dr.Joseph Shepherd(advisor),and Dr.David GoodwinThe proposed study plans to address three things.Determining whether this crossover temperature is a real chemical and physical phenomenon,establishing how well this three-step model reflects nature, and investigating if this model can be used tofind critical conditions for direct initiation in real systems.

In this report,first we present a discussion of previous work on detonation failure in Section2.1. Then we give two formulations of the ZND model of steady one-dimensional waves in Section2.2.As the focus of the study is the chemistry aspect of detonations,combustion chemistry is illustrated for the hydrogen-oxygen system in Section2.3.Once we have established a foundation,we describe more thoroughly the crossover effect in the context of specific fuel-oxidizer systems in Section3.Section4 explicitly gives the model that Short,Quirk,and Kapilla[17,18]propose as well as graphical results. Finally,the concluding Section5explains the work proposed to achieve the goals of this study.

2Background

2.1Previous Work

During the1950’s,Zeldovich et al.[21]and Manson and Ferrie[12]were thefirst to directly initiate a spherical detonation experimentally.In their experiments,they found that there is a minimum energy that must be provided in order to successfully detonate their test mixture.Because the source energy is an independent variable,knowing the value of this critical initiation energy would dictate under what circumstances there would be direct initiation success.

Critical initiation energy is termed by Lee as a dynamic detonation parameter[8]and is dependent on the structure of the detonation.Detonation structure is unsteady and multidimensional and for this reason,the most successful models for critical initiation energy are empirically based on extensive experimental data.These models include the Zeldovich et al.criterion[21],the Lee et al.model[9], and the Sichel Model[19].Litchfield et al.[11]first showed that the trend of the critical initiation energy as a function of equivalence ratio is a U-shaped curve.The minimum critical initiation occurs at approximately stoichiometric conditions while very rich or very lean mixtures require large amounts of energy to successfully initiate.

The main difficulty in determining the critical initiation energy is that there is not a robust theoretical model for detonation structure.Currently the interaction between the shock and the reaction zone is not completely understood.Extensive effort has been invested in developing theoretical and numerical models that can accurately simulate detonation structure.

The complexity of the coupled chemical reaction andfluid mechanics requires that simplifications be made in order to analyze the initiation process.This is accomplished by neglecting certain aspects of the fluid mechanics and the chemistry.For example,manipulating the energy equation using a single-step reaction model,one can derive the temperature-reaction zone structure equation[3]for one-dimensional flows:

ηc p DT

Dt

=(1)−(1−γM2)

e k˙Ωk−

c2

γ

W

W k

˙Ω

k

(2)

+

j

R c−x

w2(D−w)(3) +w

dD

dt

−w

∂w

∂t

+

1

ρ

∂P

∂t

(4)

This relation shows that there are three major contributions to the temperature change(1);these are the energy released by chemical reaction(2),the wave curvature(3),and the unsteadiness(4).It has been shown by Eckett,Quirk and Shepherd[3]that the unsteadiness of the leading shock wave playsReaction A n E a

1.H+O2→OH+O5.13×1016−0.8216507

2.H+O2+M→HO2+M2.10×1018−1.000

3.H2+M→2H+M2.23×10120.592600

4.OH+H2→H2O+H1.17×10091.303626

5.O+H2→OH+H1.80×10101.008826

6.O+H2O→2OH6.00×10081.300

Table1:Partial Hydrogen Oxidation Mechanism and Rate Constants

a role in causing failure in direct initiation.They proposed a model,the“critical decay rate”theory, based on the balance between terms2and4.

Another theory of direct initiation is based on considering the structure of slightly curved quasi-steady detonation waves[6].In this model,the dominant balance is between chemical heat release and front curvature,terms2and4.They propose that excessive curvature prevents the appearance of a sonic point at the rear of the reaction zone.The absence of a rear sonic plane allows acoustic disturbances to weaken the leading shock wave and cause the initiation to fail.

Yet a third explanation,proposed by Short and Quirk[18],is that the chemical energy release,term 2,vanishes causing the shock wave and reaction zone to decouple.This model is discussed in more detail in Section4.

2.2ZND Model of Steady One Dimensional Waves

The primary governing equations for anyfluid mechanics problem are the conservation equations for mass,momentum,energy,and species.Zeldovich[22],von Neumann[20],and Doering[2]proposed that detonation structure can be modeled by considering a leading nonreactive shock followed by a steady flow with chemical reaction.This has come to be called the ZND model of detonation.The ZND model can be expressed in terms of a system of ordinary differential equations.A key part of this formulation is the thermicity,˙σ,which is the sum of two terms:the difference in the number of moles of products and reactants and the energy absorbed or released from chemical bonds[4,pg.78].The differential equations are less difficult to solve than the mixed set of algebraic expressions and ordinary differential equations that Zeldovich,von Neumann,and Doeringfirst proposed.Around the turn of the century,Chapman[1] and Jouguet[7]hypothesized that the steady detonation velocity is the minimum velocity consistent with the jump conditions.The speed is now termed the Chapman-Jouguet(CJ)detonation velocity, D CJ.The standard model of a propagating detonation used by most analysts is a ZND structure moving at the CJ speed in a uniform mixture.

2.3Detailed Reaction Mechanism

In a combustion system,collisions between different species cause reactions and the breaking and reform-ing of bonds,ultimately converting reactants to products.This is described by a reaction mechanism containing a set of reactions and reaction rates.Table11shows a selection of reactions and rates for a hydrogen-oxygen combustion mechanism[16]

An explosive chemical reaction is due to both thermal and chain reactions.There are four primary reaction types in a branching chain:chain initiation,chain branching,chain propagation,and chain termination.A chain initiation reaction is a collision resulting in a free radical.In the hydrogen mechanism(table1),reaction3is the initiation reaction.Here the hydrogen molecule reacts with another molecule to produce two hydrogen radicals.

Once sufficient free radicals exist,chain branching and chain propagating reactions continue the chain.Chain propagating reactions produce one radical per input radical while chain branching reactions 1Units are moles,cubic centimeters,seconds,Kelvin,and calories/moleproduce more than one radical per input radical.Reaction4of the hydrogen mechanism is chain propagating while reactions1,5,and6are chain branching.In each of these branching reactions,two radicals are produced for each radical that reacts.

Chain termination reactions occur when one or more radicals react to form a stable species.Reaction 2of the hydrogen mechanism is a termination reaction.Unlike reaction1,the extra reactant molecule removes enough energy to H and O2to allow them to bond and form HO2.

The rate of reaction is proportional to the product of reactant concentrations.The proportionality constant,the reaction rate coefficient,has a modified Arrhenius form.The molecularity of the reaction is also a crucial factor in determining the reaction rate.Molecularity is the number of molecules that interact in the reaction step[5,pg.37].In the hydrogen mechanism,all of the reactions except reaction2 are bimolecular.Reaction2,the termination reaction,is trimolecular.The rate of termination reactions, which are trimolecular,increase more rapidly with pressure than bimolecular reactions.

3Crossover Effect

3.1Crossover Region

There exists a region of the temperature-pressure plane where the rates of reactions1and2(Table1) are equal.Above this region,reaction1,the branching reaction,proceeds more quickly than reaction2, the termination reaction.To determine this region in the specific case of hydrogen combustion,consider

the two reaction rates.

R1=A1T n1exp

−E a1

RT

[H][O2]R2=A2T n2exp

−E a2

RT

[H][O2][M]

[M]=

P

RT

R2=αR1⇒P=α

A1

A2

R T1+n1−n2exp

E a2−E a1

RT

(5)

Infigure1,the broken lines depict R1=αR2for several values ofα.The corresponding expression for the pressure as a function of temperature for these lines is

P=0.048αT1.18exp

−16507

2T

.

Figure1also shows a collection of gaseous hydrogen-air detonation post-shock states.Each point on the plot corresponds to a different overdrive value.The CJ detonation is noted.This plot indicates that the crossover effect may play a role in detonation initiation.For each case,the initial pressure was one atmosphere and the initial temperature was300K.

The post-shock states plotted infigure1were calculated with the Stanjan program[14].This program works with a thermodynamic datafile to calculate equilibrium states for species mixtures.It allows the user to specify the initial state including initial pressure,temperature,species mass fractions, and species phases.It then provides twenty-two different options for calculating an equilibrium state. The two important options for this study calculate the post-shock state and the Chapman-Jouguet detonation velocity.After choosing an option,Stanjan calculates the pressure,specific volume,specific energy,specific enthalpy,specific entropy,and species mass fractions of the equilibrium state.

3.2Induction Zone

Detonations have two major components:the shock wave and the reaction zone.The reaction zone usually consists of an induction zone that is almost thermally neutral followed by an exothermic recom-bination zone[15].

Figure1:Stoichiometric hydrogen-air detonation post-shock states with respect to the chemical crossover region

Figure2illustrates these elements of the detonation structure.In thisfigure,the induction length

.

coincides with the location of maximum thermicity,∆σ

Distance (cm)

Figure2:Induction length based on the maximum thermicity.It also shows the spatial temperature profile.Initial State was300K,1atm,and stoichiometric hydrogen-air.

The profiles infigure2were calculated with the ZND program[16].This program works with a chemical mechanism,the Chemkin library,an ordinary differential equation solver,and a rootfinder to spatially march from the post-shock state to the sonic plane where M=1.It allows the user to specifythe initial state including initial pressure,temperature,species mass fractions,and species phases as well as the shock speed and tolerances for the ordinary differential equations solver.It returns spatial profiles for the mach number,thermodynamic properties,velocity,thermicity,species mass fractions, and average molecular weight.Figure2was created with the GRI Mechanism version3.0[13].

3.3Crossover Temperature

The overall activation energy of a combustion system quantifies the relative ease with which the overall reaction will proceed.If the activation energy of the overall process decreases,the reaction proceeds more easily,but if the activation energy increases,the reaction will have more resistance and proceed more slowly.An Arrhenius plot graphically shows the activation energy.Induction length,∆,can be expressed as a function of temperature through the following simple model.

∆=t ind w2k=AT n exp

−E a

RT

dλ

dt

=k(1−λ)(6) In this model,λequals0in the unburned reactants and1in burned products.It is apparent in figure2that the induction zone occurs at approximately constant temperature equal to the post-shock temperature.Integrating equation6for afinite time with this in mind gives

ln(∆)=

E a1

2

+ln(A w2)

Now a linear relationship between the logarithm of the induction length and the inverse of the post-shock temperature is clear,and the slope of this line is the effective activation energy,E a/R.Different shock velocities will give different post-shock states.If there is a change in the slope of the Arrhenius plot, ln(∆)vs.1/T2,then the system’s effective activation energy has changed.This is an indication of a crossover effect from one chemical regime to another.

Figures3and4are two Arrhenius plots.Figure3shows a stoichiometric mixture of hydrogen and air while Figure4shows a stoichiometric mixture of methane and air.It is apparent that hydrogen exhibits a change of slope while methane does not.This change of slope for hydrogen-air is associated with a post-shock temperature of approximately1530K.

Thesefigures were created using∆σcalculated by the ZND program[16].D CJ and T2were given by the Stanjan program[14].Again the calculations relied on the GRI Mechanism version3.0[13]and a thermodynamics datafile that includes all of the pertinent species.

4Short,Quirk,&Kapilla Method

Short,Quirk,and Kapilla[17,18]propose that traditional one-step chemistry models are not adequate for detonation simulation.They propose an alternative reaction model that consists sequentially of a chain-initiation step,a chain-branching step,and chain-termination step.The Short,Quirk,and Kapilla model2is:

F→Y k∗I=exp

1

∗I

1

T∗I

−

1

T∗

(7)

F+Y→2Y k∗B=exp

1

∗B

1

T∗B

−

1

T∗

(8)

Y→P k∗C=1(9) 2nondimensional quantities are indicated with an*

Figure3:Arrhenius plot for stoichiometric hydrogen-air with varying overdrive values.

Figure4:Arrhenius plot for stoichiometric methane-air with varying overdrive values.

The Short,Quirk,and Kapilla[17,18]method of analyzing their three-step chemistry model uses dimensionless versions of the ZND equations with respect to the post-shock state.The three reaction rates they use are:

R∗I=fk∗I R∗B=ρ∗fyk∗B R∗C=y

With only three species,fuel,F,Y,and P,the product mass fraction can be expressed as a function of the other mass fractions,i.e.p=1−f−y.Consequently,there are only two species equations in thissystem.

d f dx∗=

−(R∗I+R∗B)

w∗

dy

dx∗

=

R∗I+R∗B−R∗C

w∗

After applying the algebraic constraints,Short,Quirk and Kapilla also solve the following non-dimensional algebraic equations.

P∗=−b±

b2−4ac

2a

1

ρ∗

=1−

(P∗−1)

γ

T∗=

P∗

ρ∗

w∗=

M2

ρ∗

where

a=−(γ+1)b=2(γ+1)c=−2γM22(γ−1)q+(γ−1−2γM22)

q=Q(1−f)−Qy

Besides the spatially varying thermodynamic properties,this system depends on certain constant values.γ,M2,T I,T B, I, B,Q,D CJ,and d,are also indirectly important to the system.Short,Quirk, and Kapilla specify Q rather than specifying M2,but M2can be determined from Q using the following system of equations.

M2=

(γ−1)D2+2

2γD2−(γ−1)

D=

√

dD CJ D CJ=

√

+

√

H H=

1

2

(γ2−1)

γ

∆h oρ1

P1

∆oρ1

P1

=

QP2ρ1

P1ρ2

P2

P1

=1+

2γ(D2−1)

γ+1

ρ1

ρ2

=1−

2

1−

1

D2

γ+1

The resulting detonation profiles for species mass fraction,temperature,pressure,and density are shown infigures5–8.The mixed system was solved by a Matlab script employing a forward Euler numerical method with the following values for the constants[17,18]:

Q=3 I=1/20 B=1/8T I=3T B=0.82γ=1.2d=1.2 Once the publishedfigures could be reproduced,the thermicity equations were non-dimensionalized. The profiles were constructed again with the dimensionless thermicity equations.

dP∗dx∗=−

γρ∗w∗˙σ∗

η

dρ∗

dx∗

=−

ρ∗˙σ∗

ηw∗

dw∗

dx∗

=−

˙σ∗

η

We have looked at both the initial state3and the post-shock state as the reference state.

3Q=8.33 I=1/37.5 B=1/10T I=3P2ρ1

P1ρ2T B=0.82P2ρ1

P1ρ2γ=1.2d=1

Figure 5:Mass fraction profiles based on the three step model.

0.911.11.21.31.41.5t e m p e r a t u r e position

Figure 6:Temperature profile based on the three step model.The reference state is the post-shock state.Results of solving the mixed algebraic-differential formulation of the ZND model.

5Proposed Work

The goal is a comprehensive model of direct initiation of detonation that can quantitatively predict the critical initiation energy.To do this we must construct simple but robust models for the chemical kinetics and include these in models of detonation structure.

This study proposes to examine how the chemistry,both realistic and modeled,affects the shock wave/reaction zone interaction.The chemical regime of the post-shock state is one factor in determining if the reaction will proceed or decouple from the shock wave.The post-shock temperature and pressure will allow us to quantify which regime the combustion system is in.We propose to study the existence and usefulness of this crossover temperature in real chemical and physical situations.

We plan to carry out simulations using the three-step model with parameters appropriate to realistic chemical mixtures.We then plan to compare these results with experimental data to determine how well the model emulates nature.This will allow us to comment on the possible applications of this reduced

1position

p r e s s u r e

Figure 7:Pressure profile based on the three step model.The reference state is the post-shock state.Results of solving the mixed algebraic-differential formulation of the ZND model.

1d e n s i t y position

Figure 8:Density profile based on the three step model.The reference state is the post-shock state.Results of solving the mixed algebraic-differential formulation of the ZND model.

chemistry mechanism.

Initially,the three-step method needs to be related back to dimensional physical parameters so dimensional profiles can be created.In addition,Arrhenius plots using this chemistry model need to be created to see if the change in slope is similar to those seen in Figures 3and 4.One of the most challenging aspects of this model is the generic nature of the species.The numerical results from the ZND and Stanjan programs depend on specific properties of the involved species –properties not included in the model.These properties include molecular weights,a variable gas constant,and a variable ratio of specific heats.

In addition,to accommodate this generic three-step chemistry model,we must analyze simulation results to see if there exist common threads between different combustion systems.An extensive study must then include different fuels,oxidizers,stoichiometric conditions,initial shock speeds,and initial

In the long term,we propose to compare unsteady simulations of initiation using the three step model with experimental results and simulations using other approximate reaction mechanisms.

Nomenclature

a,b,c constants

A pre-exponential factor

A modified A

[A]concentration of species A

c frozen soun

d speed

c p specific heat capacity

d overdriv

e factor

D detonation wave velocity

D CJ Chapman-Jouguet detonation

wave velocity

e specific internal energy

E a activation energy

f fuel mass fraction

F fuel

h specific enthalpy

∆h o heat of reaction

H dimensionless heat of reaction

based on initial state

j geometry constant

k reaction rate coefficient

k B branching reaction rate coefficient k C termination reaction rate coefficient k I initiation reaction rate coefficient

M mach number

M generic reactant molecule

n empirically determined power

p product mass fraction

P pressure

P product

q dimensionless chemical energy

Q dimensionless heat of reaction

based on post-shock state

R gas constant

R c position of the shock

in thefixed reference frame R reaction rate

R B branching reaction rate

R C termination reaction rate

R I initiation reaction rate

t time

t ind induction time

T temperature

T B branching crossover temperature

T I initiation crossover temperature

v specific volume

wfluid velocity in shock-still frame

W k molecular weight of species k

W average molecular weight of mixture x distance

y radical mass fraction

y k mass fraction of species k

Y Radical

αproportionality constant

∆induction length

∆σinduction length based

on maximum thermicity

I initiation reaction activation energy B branching reaction activation energy γratio of specific heats

λrate of progress variable

ηsonic parameter

˙ωk molar rate of

production of species k

˙Ω

k

mass rate of

production of species k

ρdensity

˙σthermicity

()1initial state

()2post-shock state

()k species k

References

[1]D.L.Chapman.On the rate of explosion in gases.Philos.Mag.,14:1091–1094,19.3

[2]W.Doering.On detonation processes in gases.Annals of Physics,43:421–436,1943.3[3]C.A.Eckett,J.J.Quirk,and J.E.Shepherd.The role of unsteadiness in direct initiation of gaseous

detonations.Journal of Fluid Mechanics,421:147–183,2000.2

[4]W.Fickett and W.Davis.Detonation Theory and Experiment.Dover Publications,INC.,1979.3

[5]I.Glassman.Combustion.Academic Press,3edition,1996.4

[6]L.He and P.Clavin.On the direct intiation of gaseous detonations by and energy source.Journal

of Fluid Mechanics,277:227–248,1994.3

[7]E.Jouguet.On the propogation of chemical reactions in gases.J.de Mathematiques Pures et

Appliquees,1:347–425,1905.continued in2:5-85(1906).3

[8]J.H.S.Lee.Dynamic parameters of gaseous detonations.Annual Review of Fluid Mechanics,

16:311–336,1984.2

[9]J.H.S.Lee,R.Knystautas,and C.Guirao.Proceedings of the First Sepcialist Meeting on Fuel-Air

Explosions,chapter SM Study No.16,pages157–187.University of Waterloo Press,Waterloo, Canada,1982.2

[10]B.Lewis and G.von Elbe.Combustion,Flames and Explosions of Gases.Academic Press,second

edition,1961.1

[11]E.L.Litchfield,M.H.Hay,and D.R.Forshey.Direct electrical initiation of freely expanding

gaseous detonation waves.Proceedings of the Ninth Symposium(International)on Combustion, pages282–286,1963.2

[12]N.Manson and F.Ferrie.Contribution to the study of spherical detonation waves.Proceedings of

the Fourth Symposium(International)on Combustion,pages486–494,1953.2

[13]University of California Berkeley.http://www.me.berkeley.edu/gri mech.GRI-Mechanism3.0.6

[14]W.Reynolds.The element potential method for chemical equilibrium analysis:Implementation in

the interactive program stanjan.Technical report,Mechanical Engineering Department,Stanford University,1986.4,6

[15]J.E.Shepherd.http://www.galcit.cal-tech.edu/EDL/projects/JetA/Glossary.html.Explosion Dy-

namics Laboratory Glossary.4

[16]J.E.Shepherd.Chemical kinetics of hydrogen-air-diluent detonations.Progress in Astronautics

and Aeronautics series,106:263–293,1986.1,3,5,6

[17]M.Short, A.K.Kapila,and J.J.Quirk.The chemical-gas dynamic mechanism of pulsating

detonation wave instability.Philosophical Transactions of the Royal Society of London,357:3621–3637,1999.2,6,7,8

[18]M.Short and J.J.Quirk.On the nonlinear stability and detonability limit of a detonation wave

for a model three-step chain-branching reaction.Journal of Fluid Mechanics,339:–119,1997.2, 3,6,7,8

[19]M.Sichel.A simple analysis of the blast initiation of detonation.Acta Astronaut,4:409–424,1977.

2

[20]J.von Neumann.John von Neumann,collected works,volume6,chapter Theory of detonation

waves.Macmillian,New York,1942.3

[21]Y.B.Zeldovich,S.M.Kogarko,and M.N.Simonov.An experimental investigation of spherical

detonation gases.Sov.Phys.Tech.Phys.,1(8):16–1713,1956.2

[22]Ya.B.Zeldovich.On the theory of the propagation of detonations in gaseous systems.Zh.Eksp.

Teor.Fiz.,10:542–568,1940.(English Translation:NACA TM1261,1960).3