S. Yungster and K. Radhakrishnan

Institute for Computational Mechanics in Propulsion, Cleveland, Ohio

Computational Study of Near-Limit

Propagation of Detonation in

Hydrogen-Air Mixtures

NASA/CR—2002-211889

October 2002

ICOMP–2002–06

AIAA–2002–3712

NI STITUTE FOR COMPUTATIONAL MECHAN CI S

RP NI OPULS OI N

Glenn Research Center

NASA

Ohio Aerospace

Institute
The NASA STI Program Office . . . in Profile

Since its founding, NASA has been dedicated to

the advancement of aeronautics and space

science. The NASA Scientific and Technical

Information (STI) Program Office plays a key part

in helping NASA maintain this important role.

The NASA STI Program Office is operated by

Langley Research Center, the Lead Center for

NASA’s scientific and technical information. The

NASA STI Program Office provides access to the

NASA STI Database, the largest collection of

aeronautical and space science STI in the world.

The Program Office is also NASA’s institutional

mechanism for disseminating the results of its

research and development activities. These results

are published by NASA in the NASA STI Report

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of

completed research or a major significant

phase of research that present the results of

NASA programs and include extensive data

or theoretical analysis. Includes compilations

of significant scientific and technical data and

information deemed to be of continuing

reference value. NASA’s counterpart of peerreviewed formal professional papers but

has less stringent limitations on manuscript

length and extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific

and technical findings that are preliminary or

of specialized interest, e.g., quick release

reports, working papers, and bibliographies

that contain minimal annotation. Does not

contain extensive analysis.

• CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

• CONFERENCE PUBLICATION. Collected

papers from scientific and technical

conferences, symposia, seminars, or other

meetings sponsored or cosponsored by

NASA.

• SPECIAL PUBLICATION. Scientific,

technical, or historical information from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific

and technical material pertinent to NASA’s

mission.

Specialized services that complement the STI

Program Office’s diverse offerings include

creating custom thesauri, building customized

databases, organizing and publishing research

results . . . even providing videos.

For more information about the NASA STI

Program Office, see the following:

• Access the NASA STI Program Home Page

at http://www.sti.nasa.gov

• E-mail your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA Access

Help Desk at 301–621–0134

• Telephone the NASA Access Help Desk at

301–621–0390

• Write to:

 NASA Access Help Desk

 NASA Center for AeroSpace Information

 7121 Standard Drive

 Hanover, MD 21076
S. Yungster and K. Radhakrishnan

Institute for Computational Mechanics in Propulsion, Cleveland, Ohio

Computational Study of Near-Limit

Propagation of Detonation in

Hydrogen-Air Mixtures

NASA/CR—2002-211889

October 2002

National Aeronautics and

Space Administration

Glenn Research Center

ICOMP–2002–06

AIAA–2002–3712

NI STITUTE FOR COMPUTATIONAL MECHAN CI S

RP NI OPULS OI N

Glenn Research Center

NASA

Ohio Aerospace

Institute

Prepared under Cooperative Agreement NCC3–542

Prepared for the

38th Annual Joint Propulsion Conference and Exhibit

cosponsored by the AIAA, ASME, SAE, and ASEE

Indianapolis, Indiana, July 7–10, 2002
Available from

NASA Center for Aerospace Information

7121 Standard Drive

Hanover, MD 21076

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22100

The Aerospace Propulsion and Power Program at

NASA Glenn Research Center sponsored this work.

Available electronically at http://gltrs.grc.nasa.gov
1

Abstract

A computational investigation of the near-limit propagation of detonation in lean and rich hydrogen-air mixtures is

presented. The calculations were carried out over an equivalence ratio range of 0.4-5.0, pressures ranging from 0.2 bar

to 1.0 bar, and ambient initial temperature. The computations

involved solution of the one-dimensional Euler equations

with detailed finite-rate chemistry. The numerical method is

based on a second-order spatially accurate total-variation-diminishing (TVD) scheme, and a point implicit,

first-order-accurate, time marching algorithm. The hydrogen-air combustion was modeled with a 9-species, 19-step

reaction mechanism. A multi-level, dynamically adaptive

grid was utilized in order to resolve the structure of the detonation. The results of the computations indicate that when

hydrogen concentrations are reduced below certain levels,

the detonation wave switches from a high-frequency, low

amplitude oscillation mode to a low frequency mode exhibiting large fluctuations in the detonation wave speed; that is, a

“galloping” propagation mode is established.

Introduction

Recent experimental and computational results indicate

that the specific impulse obtained in a hydrogen fueled,

air-breathing pulsed detonation engine (PDE) is highest

when it is operated under fuel lean conditions1-3. However,

when the fuel-air equivalence ratio, φ, reached values below

approximately 0.6, a sharp drop off in thrust and specific

impulse was observed experimentally, possibly due to having exceeded the detonability limits for hydrogen-air mixtures.

At fuel-rich conditions, the performance of an air-breathing PDE is significantly degraded. Nevertheless, a potential

application of a fuel-rich pulsed detonation rocket engine

(PDRE) has been proposed in combination with an ejector.

In such a PDRE-ejector engine, the excess fuel in the PDRE

tube is mixed and burned with entrained air and subsequently expanded through a nozzle in an unsteady version of

the ejector-ramjet concept4. Unsteady ejectors are potentially more efficient than their steady-state counterparts5.

There is therefore considerable practical interest in investigating the limits of detonability for both rich and lean mixtures. Published work on detonability limits is not extensive.

Gordon et al.6carried out experiments in hydrogen-oxygen

mixtures with different diluents. Of particular interest to the

present study are their results for hydrogen-air mixtures

which exhibited large oscillations in detonation speed near

the detonability limits. Dupre et al.7 also showed that in

hydrogen-air mixtures large detonation velocity fluctuations

were obtained when hydrogen concentrations were lowered

below a certain level. Because of the large variation in propagation speed, these marginal detonations are sometimes

called “galloping” detonations.

Belles8 suggested that the detonability limits in hydrogen-air and hydrogen-oxygen systems are the mixture compositions for which the conditions at the von Neumann spike

lie outside the isothermal branched-chain explosion limits.

Belles’ analysis included only four chemical reactions. Dove

and Tribbeck9 extended the work of Belles by considering a

more complete combustion mechanism. They showed that

secondary reactions involving HO2 cannot be ignored, and

that the effects of exothermicity are also important and can

lead to a “chain-thermal” explosion. That is, explosions can

occur for compositions significantly outside the isothermal

limit. However, the numerical methods needed for solving

the stiff ODE’s describing the combustion process were not

sufficiently developed at the time, and they could not determine the boundaries outside the isothermal limits for which

explosive behavior would still exist.

In this paper, we investigate first the combustion of

hydrogen-air mixtures at the von Neumann spike conditions

using the LSENS code10,11. Subsequently, we examine in

more detail the structure and stability of the detonation wave

for near-limit mixtures by conducting one-dimensional computational fluid dynamic (CFD) calculations.

Reaction Kinetics at the von Neumann Spike

The reaction rate equations were integrated under the

assumption of a Zeldovich, von Neumann, Doring (ZND)

detonation. The actual calculation consisted of two parts:

(i) The determination of the state at the von Neumann

spike, which is used as the initial condition for the chemical

kinetic calculation. This state was found using the CEA

equilibrium code12.

(ii) The integration of the chemical kinetic rate equations

with the initial conditions found in step one. These calculaCOMPUTATIONAL STUDY OF NEAR-LIMIT PROPAGATION OF

DETONATION IN HYDROGEN-AIR MIXTURES

 S. Yungster and K. Radhakrishnan

Institute for Computational Mechanics in Propulsion

NASA/CR—2002-211889

Brook Park, Ohio 44142
2

tions were carried out using the LSENS chemical kinetics

code11.

The combustion mechanism used in the calculations in

part (ii) and in the CFD computations described later, is

based on Jachimowski’s model13, and is listed in Table 1. It

consists of 19 elementary reversible reactions among 8 reacting species, with N2 treated as an inert (i.e., nonreacting)

species.

Results are presented in figs. 1 and 2 in the form of ignition delay times as a function of mixture equivalence ratio

for two different initial fill pressures: p0 = 0.4 bar and

p0 = 1.0 bar. Here, the ignition delay is arbitrarily defined as

the time required to increase the temperature of the gas by

25 K. For p0 = 0.4 bar (fig. 1), ignition occurs within a few

microseconds for the equivalence ratio range 0.6 < φ < 3.0.

Outside this range, the ignition delay increases very rapidly.

For p0 = 1.0 bar (figure 2), a similar trend is observed, but

the curves are significantly steeper on the lean and rich sides

of stoichiometric than those at the lower pressure. Also, the

ignition delay times near stoichiometric conditions are

smaller at the higher pressure. At φ = 1, for example, the

ignition delay is 0.83 μs for p0 = 0.4 bar, and 0.36 μs for

p0 = 1.0 bar. Figure 2 also includes the value of the detonability limits computed by Belles8 as well as experimental

values obtained by Laffitte14, as reported in ref. 8.

It is clear from figs. 1 and 2 that the limits of detonability

are probably determined by the steep increase in the ignition

delay time for very lean and rich mixtures. That is, for mixtures significantly away from stoichiometric conditions, the

large ignition delay times may prevent the coupling between

the reaction front and the shock wave needed to form a stable

detonation wave. Based on this argument, it appears from

figs. 1 and 2 that the limits of detonability become further

apart as the pressure is lowered, a result also predicted by

Belles8. This somewhat counterintuitive behavior can be

explained from the variation of the ignition delay (or explosion limit) with pressure for hydrogen-air mixtures. Table 2

gives ignition delay times for a stoichiometric hydrogen-air

mixture at an initial temperature of 1250 K. It is observed

that the ignition delay time is not a monotonic function of

pressure. For pressures between 10 and 40 bars, the ignition

delay actually increases with pressure.

The precise value of the maximum ignition delay that

would allow the establishment of a stable detonation, or the

changes in the structure and stability of a detonation wave at

near-limit conditions cannot be determined from these type

of calculations. To answer these questions a more detailed

analysis, based on CFD computations that include the

dynamics of the flow and combustion process (and their

interaction) is needed. Such analysis is presented in the following sections.

CFD Simulations

As a first step in understanding the behavior of detonation waves at near-limit conditions, we carried out

one-dimensional CFD computations. It is well known that

detonation waves have a three-dimensional structure, and

therefore, a one-dimensional analysis will not be able to

investigate the effect of transverse waves. However, such a

study will be able to examine, in a relatively quick and inexpensive way (particularly when a detailed reaction model is

used to describe the chemistry) the longitudinal stability

characteristics of detonation waves. Moreover, a one-dimensional analysis can serve as a basis for subsequent

multi-dimensional calculations.

Since detonation limits are not dependent on transport

phenomena6, we neglect all dissipation effects and therefore

based our study on the Euler equations. We use the conservation form of the unsteady Euler equations and replace the

global continuity equation with the 9 species conservation

equations. A detailed description of the set of differential

equations and additional state and constitutive equations

needed for system closure are given in ref. 15.

The numerical method used for solving the governing

equation set is based on Yee’s second-order spatially accurate total-variation-diminishing (TVD) scheme16, and a

point implicit, first-order time-accurate marching algorithm.

This method is a subset of a more general class of BDF

methods considered by Yungster and Radhakrishnan17. It is

obtained by setting the variables γ=0 and β=1 in equation 4

of ref. 17, and treating only the chemical source term implicitly. Further details about the numerical algorithm are given

in ref. 17.

In order to maintain adequate numerical resolution of the

detonation wave front, without the need to use hundreds of

thousands of grid points, a multi-level, dynamically adaptive

grid is utilized. Figure 3 shows a section of the grid at three

different times as the detonation wave moves from left to

right. The grid constantly adapts to keep the detonation front

within the finest grid level. An arbitrary number of levels can

be specified. Nine or ten grid levels were used in the present

study, and 100 points were included in the finest grid level.

CFD Results

There are two ways in which a detonation can be formed

in general: (1) by direct initiation, in which a strong shock

wave is generated in the tube (with a charge of solid explosive or by using a high pressure reservoir) and (2) by transition from a deflagration. In this paper, we consider only

NASA/CR—2002-211889
3

development of detonations with direct initiation. Figure 4

shows a schematic of the initiation process used in the

present study. A high pressure, high temperature driver gas,

consisting of nitrogen, was used in a small region next to the

head-end of the tube. As the computation is started, a shock

wave travels to the right and an expansion wave propagates

to the left, towards the head-end. The shock wave is strong

enough to initiate chemical reactions in the combustible mixture. Depending on the mixture properties, the shock wave

and the combustion front can subsequently merge and form a

detonation wave. The value of the driver pressure, pdriv, will

determine the degree of overdrive of the resulting detonation

wave. The degree of overdrive f is defined as

where D is the actual detonation propagation speed, and DCJ

is the theoretical Chapman-Jouguet detonation speed.

The first case considers a stoichiometric mixture of

hydrogen-air at p0 = 0.4 bar and T0 = 298 K and a driver

pressure ratio, rp = pdriv/p0, of 100. A grid refinement study

for this case is presented in figure 5. Numerical results are

presented for four successively refined grids having a minimum spacing, ∆xmin, indicated in each figure. The plots

show the variation in the detonation speed with time. Initially the shock wave travels at a constant speed (determined

by the initial conditions) slightly above 1500 m/s. Chemical

reactions behind the shock front start generating strong compression waves (i.e., a reaction shock) after approximately

10 μs . The reaction shock overtakes the shock front at

around 13 μs, strengthening it and sharply increasing its

propagation speed (figs. 5b-5d). Subsequently, the detonation wave speed decreases gradually.

It can be observed that the structure of the detonation

wave changes with grid spacing. On the coarsest grid

(fig. 5a), the detonation reaches a steady propagation speed

of D = 2020 m/s after the initial transient phase. For the next

finer grid (fig. 5b) the computed low frequency oscillation

appears to be a numerical artifact. With the finest two grids

used in this study (figs 5c and 5d), the detonation reaches a

high-frequency, low-amplitude propagation mode; however,

with the third grid (fig. 5c) a low frequency mode can still be

observed (in addition to the high-frequency one). The solution in Fig. 5d has an average propagation speed of D = 2025

m/s, nearly identical to the steady solution in fig. 5a. The frequency of the oscillations, ω, in fig. 5d is 0.90 ± 0.05 MHz.

Figure 6 shows pressure distributions at different times

for the grid spacing corresponding to fig. 5d. At t = 10.8 μs,

chemical reactions behind the shock front have produced a

sharp peak in pressure. No change in the shock front propagation speeds occurs at this time. At t = 13.7 μs, the reaction

shock has overtaken the shock front resulting in a spike in

pressure and propagation speed. After approximately

t = 38.3 μs, high frequency pressure oscillations become evident downstream of the detonation front.

The theoretical Chapman-Jouguet detonation speed for

this mixture, as computed from the chemical equilibrium

code CEA12 is DCJ = 1956.2 m/s. Therefore, the results presented in figs. 5 and 6 are actually those of an overdriven detonation, with a degree of overdrive f = 1.07. The frequency

of oscillation computed in fig. 5d is in close agreement with

the results of Sussman18, who obtained a frequency of

1 MHz on a similar piston-driven detonation computation for

a degree of overdrive of f = 1.1, slightly higher than that considered here.

In the ballistic range experiments of Lehr19, which consisted of spherical projectiles being fired into hydrogen-air

mixtures, a frequency of oscillation of 1.04 MHz was

obtained for a projectile speed of 2029 m/s, and for the same

mixture considered in Fig. 5, but at a slightly higher pressure

(0.42 bar).

The high-frequency mode of oscillation can be explained

by the McVey-Toong wave interaction mechanism20. Two

fundamental processes form the basis of this mechanism:

(i) when a new reaction front is created, compression waves

(reaction shocks) are generated that travel upstream and

downstream from the new reaction front; (ii) when an old

reaction front is extinguished, it must be accompanied by the

generation of upstream and downstream rarefaction waves

which have a strength comparable to the reaction shocks.

The interaction between these waves and the shock front

produces the high frequency oscillations, as described in

detail in refs. 20 and 17.

To verify that the high-frequency oscillations computed

in this case are described by this mechanism, the history of

pressure and density in the region immediately behind the

detonation front have been plotted in an x-t diagram in fig. 7.

This figure shows three pulses corresponding to three oscillation cycles. (The “staircase” appearance of this plot is due

to the periodic sliding of the grid to keep track of the detonation front). Careful examination of this figure reveals all of

the elements of the McVey-Toong mechanism.

Unless otherwise stated, all subsequent calculations presented below, were carried out on the same grid used in fig.

5d. Figure 8 shows plots of detonation speed as a function of

time for several lean mixtures. For φ = 0.8, the structure of

the detonation is very similar to that observed in the stoichiometric case, except that the frequency has decreased to

0.76 MHz. For φ = 0.6, the frequency decreases further to

0.52 MHz and an additional low frequency mode appears. At

φ = 0.5 the low and high frequency modes are present simultaneously, and at φ = 0.4 the detonation wave almost completely switches from a high-frequency, low-amplitude

f D2 DCJ

2 = ⁄

NASA/CR—2002-211889
4

oscillation mode to a low-frequency, high-amplitude one;

that is, a “galloping” propagation mode is established.

While there is no generally accepted criterion for defining limit behavior, galloping detonations represent the lowest

propagation mode of detonation in a tube, and it has been

suggested that the detonability limit should be based on the

onset of the galloping mode7.

Figure 9 presents the maximum pressure behind the

shock front as a function of time for the same hydrogen-air

mixtures as above. This figure shows the increase in the

amplitude of the oscillations and the switch to the low frequency mode as the fuel-air equivalence ratio is lowered.

Figure 10 shows the variation in detonation wave speed

with time for rich mixtures. Note that as φ increases, the frequency of oscillation decreases and the time it takes for the

detonation to develop also increases significantly. For

φ = 4.0, only a few oscillation cycles were captured before

the detonation wave reached the end of the computational

domain. The frequencies of oscillation are: 0.62 MHz for

φ = 2.0, 0.41 MHz for φ = 3.0, and 0.29 MHz for φ = 4.0.

Pressure effect

The effect of pressure on detonation wave structure was

examined and the results are presented in figs. 11 and 12, for

p0 = 0.2 and 1.0 bar respectively.

The theoretical Chapman-Jouguet detonation speed for a

stoichiometric hydrogen-air mixture at p0 = 0.2 bar and

T0 = 298 K is 1939.5 m/s. The average detonation speed

computed in fig. 11a after 150 μs is 1960.7 m/s, corresponding to a degree of overdrive of f = 1.02.(Note that this calculation used a driver pressure ratio of rp = 70).

The structure of the detonation wave is similar to that for

the p0 = 0.4 bar case, but the frequencies of oscillation are

much smaller. For φ = 1.0, the frequency is 0.34 ± 0.02 MHz,

and it decreases to 0.29 MHz for φ = 0.8. At φ = 0.6, the high

frequency oscillations are much less uniform and range from

approximately 0.18 MHz to 0.28 MHz. At this equivalence

ratio (φ = 0.6) the low frequency mode begins to become evident. At φ = 0.4, the only mode present is the low-frequency

one, except during maximum detonation speed when very

weak traces of the high frequency mode can be seen.

The results for p0 = 1.0 bar are shown in fig. 12. A finer

grid, having a minimum spacing of ∆xmin=7.324x10-5 cm,

was used in this case, due to the shorter ignition delay times

at this pressure near stoichiometric conditions. The theoretical Chapman-Jouguet detonation speed for a stoichiometric

hydrogen-air mixture at p0 = 1.0 bar and T0 = 298 K is

1976.7 m/s. The average detonation speed computed in fig.

12a after 15 μs is 2127.7 m/s, corresponding to a degree of

overdrive of f = 1.16. Note that this calculation used a driver

pressure ratio of rp = 150. In general, the higher the pressure,

the higher the degree of overdrive needed to obtain the

high-frequency, low amplitude mode of oscillation described

by the McVey-Toong mechanism.

The results for p0 = 1.0 bar show more variation in the

amplitude of oscillation than the p0 = 0.4 and 0.2 bar cases,

even though the degree of overdrive is higher in this case.

The frequencies of oscillation are also much higher at

p0 = 1.0 bar. For the φ = 1.0 case, the frequency of oscillation

is 3.2 ± 0.2 MHz. At lower equivalence ratios, the frequencies are 2.71 MHz for φ = 0.8, 2.08 MHz for φ = 0.6, and

1.28 MHz for φ = 0.4. Note that due to the high degree of

overdrive, the high-frequency mode persists at the smallest

value of the equivalence ratio (φ = 0.4).

Driver pressure ratio effect

The effect of driver pressure ratio, rp, was investigated

for hydrogen-air mixtures at p0 = 0.4 bar. Figure 13 shows

the detonation speed for four mixtures at a driver pressure

rp = 150. It can be seen that compared with the results of

fig. 8, (which were obtained with rp = 100) the oscillations

are much more uniform, and have a lower amplitude and a

higher frequency. The high-frequency mode is present even

for the φ = 0.4 case.

At φ=1.0, the average detonation speed after 50 μs is

2087.3 m/s, corresponding to a degree of overdrive f = 1.14.

The frequency of oscillation for this case is 1.19 ± 0.03

MHz. The frequencies for the other cases are: 1.09 MHz for

φ=0.8, 0.91 MHz for φ=0.6, and 0.63 MHz for φ=0.4.

Figure 14 shows the results for the same mixtures, but a

driver pressure ratio of rp = 40. Previous “steady” detonation

calculations3,4 (i.e. coarse-grid computations) using this

driver pressure ratio, have shown that the resulting detonation speeds agree closely with the Chapman-Jouguet theoretical speed. That is, for this driver pressure ratio, the degree of

overdrive should be very close to one (i.e., f ≈ 1.0). However,

for the high-resolution computations in the present study, the

detonation wave at this driver pressure ratio contains very

complex, irregular oscillations for the equivalence ratio

range 0.6 < φ < 1.0. (The φ = 0.4 case should be run for

longer times to determine the oscillation pattern). This

behavior for p0 = 0.4 bar, and degree of overdrive close to

one, agrees with the findings of Sussman18 for a stoichiometric hydrogen-air mixture at nearly identical initial pressure

and temperature. Bourlioux et al.

21 also showed that the

oscillations obtained with a one-step chemistry model

became nearly chaotic for low degrees of overdrive. It is not

clear wether the one-dimensional detonations considered

NASA/CR—2002-211889
5

here are inherently unstable at these conditions or wether

further grid refinement is necessary to obtain accurate solutions.

Conclusions

The one-dimensional Euler equations with finite rate

chemistry were solved to investigate the structure and stability of detonation waves for hydrogen-air mixtures at various

equivalence ratios and pressures. The results indicate that the

solutions to the governing equations do not include a steadily

propagating detonation wave. Instead, the solutions always

exhibit an unsteady oscillatory propagation mode. For mixtures close to stoichiometric conditions, the detonation wave

propagates in a high-frequency low amplitude mode that is

described by the McVey-Toong mechanism. For very lean

mixtures, the detonation wave switches to a low-frequency,

high-amplitude propagation mode similar to the “galloping”

detonations observed experimentally in near-limit mixtures.

The transition from the high-frequency mode to the

low-frequency mode depends on the degree of overdrive.

High degrees of overdrive tend to stabilize the detonation

wave and delay the transition to the low-frequency mode for

lean mixtures. The frequency and amplitude of the oscillations are very uniform for high degrees of overdrive. At low

pressures (p0 = 0.2 bar), solutions that behave according to

the McVey-Toong mechanism were possible for degrees of

overdrive close to one. For higher pressures (p0 > 0.4 bar),

these solutions were possible only for degrees of overdrive

above 1.07. Calculations with degrees of overdrive close to

one at the higher pressures resulted in complex, irregular

oscillations.

References

1. Schauer, F., Stutrud, J. and Bradley, R., “Detonation

Initiation Studies and Performance Results for Pulsed Detonation Engine Applications,” AIAA paper 2001-1129, Jan.

2001.

2. Povinelli, L.A. and Yungster, S., “Thermodynamic

Cycle and CFD Analyses for Hydrogen Fueled Air-Breathing Pulse Detonation Engines,” AIAA paper 2002-3629, July

2002.

3. Povinelli, L.A. and Yungster, S., “Airbreathing Pulse

Detonation Engine Performance,” NASA TM 2002-211575,

also JANNAF Conference, 26th Airbreathing Propulsion

Subcommittee, Destin, FL, April 2002.

4. Yungster, S. and Perkins, H.D., “Multiple- Cycle

Simulation of a Pulse Detonation Engine Ejector,” AIAA

paper 2002-3630, July 2002.

5. Porter, J.L. and Squyers, R.A., “A Summary/Overview of Ejector Augmentor Theory and Performance,” ATC

Report No. R-91100/9CR-47A, Sept. 1979.

6. Gordon, W.E., Mooradian, A.J. and Harper, S.A.,

“Limit and Spin Effects in Hydrogen-Oxygen Detonations,”

Seventh Symposium (Inter) on Combustion, P752-759, Butterworths, 1959.

7. Dupre, G., Knystautas, R. and Lee, J.H., “Near-Limit

Propagation of Detonation in Tubes,” 10th International

Colloquium on Dynamics of Explosions and Reactive Systems, Berkeley, CA, Aug. 1985.

8. Belles, F.E., “Detonability and Chemical Kinetics:

Prediction of Limits of Detonability of Hydrogen,” Seventh

Symposium (Inter) on Combustion, P745-751, Butterworths,

1959.

9. Dove, J.E. and Tribbeck, T.D., “Computational Study

of the Kinetics of the Hydrogen-Oxygen Reaction Behind

Steady State Shock Waves. Application to the Composition

Limits and Transverse Stability of Gaseous Detonations,”

Astronautica Acta, Vol. 15 pp. 387-397, 1970.

10. Radhakrishnan, K., “Combustion Kinetics and Sensitivity Analysis Computations,” in Numerical Approaches to

Combustion, E.S. Oran and J.P. Boris, Eds., AIAA, Washington, DC, 1991, pp. 83-128.

11. Radhakrishnan, K., “LSENS–A General Chemical

Kinetics and Sensitivity Analysis Code for Homogeneous

Gas-Phase Reactions. I. Theory and Numerical Solution Procedures,” NASA RP-1328, 1994.

12. McBride, B.J. and Gordon, S., “Computer Program

for Calculation of Complex Chemical Equilibrium Compositions and Applications. II. Users Manual and Program

Description,” NASA RP-1311, 1996.

13. Jachimowski, C.J., “An Analytical Study of the

Hydrogen-Air Reaction Mechanism with Application to

Scramjet Combustion,” NASA TP-2791, Feb. 1988.

14. Laffitte, P. F., “Flames of high-speed detonation,” Science of Petroleum, Vol. IV, pp. 2995-3003, 1938.

15. Yungster, S., “Numerical Study of Shock-Wave

Boundary Layer Interactions in Premixed Combustible

Gases,” AIAA Journal, Vol. 30, No. 10, 1992, pp. 2379-2387.

16. Yee, H.C., “Construction of Explicit and Implicit

Symmetric TVD Schemes and Their Applications,” Journal

of Computational Physics, Vol. 68, 1987, pp. 151-179.

17. Yungster, S. and Radhakrishnan, K., “A Fully Implicit

Time Accurate Method for Hypersonic Combustion: Application to Shock-Induced Combustion Instability,” Shock

Waves, Vol. 5, 1996, pp. 293-303.

NASA/CR—2002-211889
6

18. Sussman, M.A., “A Computational Study of

Unsteady Shock Induced Combustion of Hydrogen-Air Mixtures,” AIAA paper 94-3101, June 1994.

19. Lehr, H.F., “Experiments on Shock-Induced Combustion,” Acta Astronautica, Vol. 17, 1972, pp. 589-597.

20. McVey, J.B. and Toong, T.Y., “Mechanism of Instabilities of Exothermic Hypersonic Blunt Body Flows,” Combustion Science and Technology, Vol. 3, pp. 63-76, 1971.

21. Bourlioux, A., Majda, A.J. and Roytburd, V., “Theoretical and Numerical Structure of Unstable One-Dimensional Detonations,” SIAM Journal of Applied Mathematics,

Vol. 51, pp. 303-343, April 1991.

TABLE 1.—H2-Air Reaction Mechanisma

No. Reaction A b Θ

1 H2 + O2 = HO2 + H 1.0×1014 0.0 28197.38

2 H + O2 = OH + O 2.6×1014 0.0 8459.21

3 H2 + O = OH + H 1.8×1010 1.0 4481.37

4 H2 + OH = H + H2O 2.2×1013 0.0 2593.15

5 OH + OH = O + H2O 6.3×1012 0.0 548.84

6b H + OH + M = H2O + M 2.2×1022 -2.0 0.0

7b H + H + M = H2 + M 6.4×1017 -1.0 0.0

8b H + O + M = OH + M 6.0×1016 -0.6 0.0

9b H + O2 + M = HO2 + M 2.1×1015 0.0 -503.52

10 O + O + M = O2 + M 6.0×1013 0.0 906.34

11 HO2 + H = OH + OH 1.4×1014 0.0 543.81

12 HO2 + H = H2O + O 1.0×1013 0.0 543.81

13 HO2 + O = O2 + OH 1.5×1013 0.0 478.35

14 HO2 + OH = H2O + O2 8.0×1012 0.0 0.0

15 HO2 + HO2 = H2O2 + O2 2.0×1012 0.0 0.0

16 H + H2O2 = H2 + HO2 1.4×1012 0.0 1812.69

17 O + H2O2 = OH + HO2 1.4×1013 0.0 3222.56

18 OH + H2O2 = H2O + HO2 6.1×1012 0.0 720.04

19b H2O2 + M = OH + OH + M 1.2×1017 0.0 22910.37

a

Forward rate coefficient kf = ATbexp(-Θ/T); units are moles,

 seconds, centimeters, and Kelvins.

bThird-body efficiencies:

 Reaction 6: H2O = 6.0; Reaction 7: H2O = 6.0;

 Reaction 8: H2O = 5.0; Reaction 9: H2O = 16.0, H2 = 2.0;

 Reaction 19: H2O = 15.0.

TABLE 2.—Ignition Delay Time

(tig) for φ=1.0, T1=1250.0 K

p1 (bar) tig (μs)

0.5 50.3

1.0 25.1

10.0 4.8

15.0 11.1

20.0 41.0

25.0 54.7

30.0 60.3

40.0 62.3

50.0 60.1

75.0 52.0

100.0 44.9

NASA/CR—2002-211889
7

0.3 0.5 0.7 0.9 1.1

Equivalence ratio, φ

0

5

10

15

20

25

30

Ignition Delay (μsec)

1

Aφ + B -----------------

A= 1.98519

B= -0.74646

0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5

Equivalence ratio, φ

0

5

10

15

20

25

30

35

40

Ignition Delay (μsec)

Fig. 1. Ignition delay for Chapman-Jouguet detonations in H2-Air; p0=0.4 bar, T0=298 K.

0.3 0.5 0.7 0.9 1.1

Equivalence ratio, φ

0

5

10

15

20

25

30

Ignition Delay (μsec)

1

Aφ + B -----------------

A= 5.65195

B= -2.80115

Detonability limits

Theory (Belles 1959)

Experiment (Laffitte 1938)

0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5

Equivalence ratio, φ

0

5

10

15

20

25

30

35

40

Ignition Delay (μsec)

Detonability limits

Theory (Belles 1959)

Experiment (Laffitte 1938)

Fig. 2. Ignition delay for Chapman-Jouguet detonations in H2-Air; p0=1.0 bar, T0=298 K.

NASA/CR—2002-211889
8

H2+Air

p=p0; T0=298 K

N2

p=pdriv

L

T=2500.0 K

Fig. 4. Schematic of detonation wave initiation.

0 20 40 60 80

Time (μs)

1500

2000

2500

3000

3500

Detonation speed (m/s)

0 20 40 60 80

Time (μs)

1500

2000

2500

3000

3500

Detonation speed (m/s)

0 20 40 60 80

Time (μs)

1500

2000

2500

3000

3500

Detonation speed (m/s)

0 20 40 60 80

Time (μs)

1500

2000

2500

3000

3500

Detonation speed (m/s)

∆xmin= 1.953x10-3 cm

∆xmin= 9.766x10-4 cm

∆xmin= 4.883x10-4 cm

∆xmin= 2.441x10-4 cm

Fig. 5. Detonation speed as a function of time; H2-Air;

p0 = 0.4 bar, φ=1.0; rp = pdriv/p0 = 100; DCJ = 1956.2 m/s;

f = 1.07.

(a)

(b)

(c)

(d)

Fig. 3. Computational grid at three different times.

D 2025 m/s=

D 2020 m/s=

ω = 0.90 MHz

NASA/CR—2002-211889
9

Fig. 8. Detonation speed as a function of time; H2-Air;

p0 = 0.4 bar, rp = 100, ∆xmin=2.441x10-4 cm.

0 20 40 60 80

Time (μs)

1200

1700

2200

2700

3200

Detonation speed (m/s)

0 20 40 60 80

Time (μs)

1200

1700

2200

2700

3200

Detonation speed (m/s)

0 20 40 60 80

Time (μs)

1200

1700

2200

2700

3200

Detonation speed (m/s)

0 20 40 60 80

Time (μs)

1200

1700

2200

2700

3200

Detonation speed (m/s)

φ = 1.0

φ = 0.8

φ = 0.6

0 20 40 60 80

Time (μs)

1200

1700

2200

2700

3200

Detonation speed (m/s)

φ = 0.4

φ = 0.5

x

t

x

t

Pressure

Density

Fig. 7. History of pressure and density (x-t diagram)

showing the Mcvey-Toong mechanism; rp=150.

Shock front

New reaction front

Induction zone

Reaction shock

Fig. 6. Pressure ratio as a function of distance at different

times; H2-Air; p0 = 0.4 bar, φ=1.0; rp = pdriv/p0 = 100.

0 3 6 9 12 15

Distance (cm)

0

10

20

30

40

50

p/p0

t=10.8 μs 13.7

38.3 62.8

NASA/CR—2002-211889
10

40.0 45.0 50.0 55.0 60.0 65.0 70.0

Time (μs)

26.0

28.0

30.0

32.0

34.0

36.0

pmax/p0 

40.0 45.0 50.0 55.0 60.0 65.0 70.0

Time (μs)

26.0

28.0

30.0

32.0

34.0

36.0

pmax/p0 

40.0 45.0 50.0 55.0 60.0 65.0 70.0

Time (μs)

24.0

26.0

28.0

30.0

32.0

34.0

pmax/p0 

40.0 45.0 50.0 55.0 60.0 65.0 70.0

Time (μs)

22.0

24.0

26.0

28.0

30.0

32.0

34.0

pmax/p0 

40.0 45.0 50.0 55.0 60.0 65.0 70.0

Time (μs)

20.0

25.0

30.0

35.0

40.0

45.0

pmax/p0 

φ=1.0

φ=0.8

φ=0.6

φ=0.5

φ=0.4

Fig. 9. Maximum pressure as a function of time; H2-Air;

p0 = 0.4 bar, rp = 100, ∆xmin=2.441x10-4 cm.

0 20 40 60 80 100

Time (μs)

1500

2000

2500

3000

3500

Detonation speed (m/s)

0 20 40 60 80 100

Time (μs)

1500

2000

2500

3000

3500

Detonation speed (m/s)

0 20 40 60 80 100

Time (μs)

1500

2000

2500

3000

3500

Detonation speed (m/s)

Fig. 10. Detonation speed as a function of time; H2-Air;

p0 = 0.4 bar, rp = 100, ∆xmin=2.441x10-4 cm.

φ=2.0

φ=3.0

φ=4.0

ω = 0.62 MHz

ω = 0.29 MHz

ω = 0.41 MHz

ω = 0.90 MHz

ω = 0.76 MHz

ω = 0.52 MHz

NASA/CR—2002-211889
11

0 20 40 60 80 100 120 140 160 180 200

Time (μs)

1000

1500

2000

2500

3000

Detonation speed (m/s)

0 20 40 60 80 100 120 140 160 180 200

Time (μs)

1000

1500

2000

2500

3000

Detonation speed (m/s)

0 20 40 60 80 100 120 140 160 180 200

Time (μs)

1000

1500

2000

2500

3000

Detonation speed (m/s)

0 20 40 60 80 100 120 140 160 180 200

Time (μs)

1000

1500

2000

2500

3000

Detonation speed (m/s)

φ=1.0

φ=0.8

φ=0.6

φ=0.4

Fig. 11. Detonation speed as a function of time; H2-Air;

p0 = 0.2 bar, rp = 70, ∆xmin=2.441x10-4 cm.

0 5 10 15 20

Time (μs)

1500

2000

2500

3000

3500

Detonation speed (m/s)

0 5 10 15 20

Time (μs)

1500

2000

2500

3000

3500

Detonation speed (m/s)

0 5 10 15 20

Time (μs)

1500

2000

2500

3000

3500

Detonation speed (m/s)

0 5 10 15 20

Time (μs)

1500

2000

2500

3000

3500

Detonation speed (m/s)

φ=1.0

φ=0.8

φ=0.6

φ=0.4

Fig. 12. Detonation speed as a function of time; H2-Air;

p0 = 1.0 bar, rp = 150, ∆xmin=7.324x10-5 cm.

(a)

(b)

(c)

(d)

(a)

(b)

(c)

(d)

ω = 3.20 MHz

ω = 2.71 MHz

ω = 2.08 MHz

ω = 1.28 MHz

D 2128 m/s=

Dcj = 1977 m/s}f=1.16

ω = 0.34 MHz

ω = 0.29 MHz

0.18 < ω < 0.28 MHz

D 1961 m/s=

Dcj = 1940 m/s}f=1.02

NASA/CR—2002-211889
12

0 20 40 60 80

Time (μs)

1200

1700

2200

2700

3200

Detonation speed (m/s)

0 20 40 60 80

Time (μs)

1200

1700

2200

2700

3200

Detonation speed (m/s)

0 20 40 60 80

Time (μs)

1200

1700

2200

2700

3200

Detonation speed (m/s)

0 20 40 60 80

Time (μs)

1200

1700

2200

2700

3200

Detonation speed (m/s)

Fig. 13. Detonation speed as a function of time; H2-Air;

p0 = 0.4 bar, rp = 150, ∆xmin=2.441x10-4 cm.

0 20 40 60 80 100

Time (μs)

1000

1500

2000

2500

3000

Detonation speed (m/s)

0 20 40 60 80 100

Time (μs)

1000

1500

2000

2500

3000

Detonation speed (m/s)

0 20 40 60 80 100

Time (μs)

1000

1500

2000

2500

3000

Detonation speed (m/s)

0 20 40 60 80 100

Time (μs)

1000

1500

2000

2500

3000

Detonation speed (m/s)

Fig. 14. Detonation speed as a function of time; H2-Air;

p0 = 0.4 bar, rp = 40, ∆xmin=2.441x10-4 cm.

φ=1.0

φ=0.8

φ=0.6

φ=0.4

φ=1.0

φ=0.8

φ=0.6

φ=0.4

D 2087 m/s=

f = 1.14

ω = 0.63 MHz

ω = 1.19 MHz

ω = 1.09 MHz

ω = 0.91 MHz

NASA/CR—2002-211889
This publication is available from the NASA Center for AeroSpace Information, 301–621–0390.

REPORT DOCUMENTATION PAGE

2. REPORT DATE

19. SECURITY CLASSIFICATION

 OF ABSTRACT

18. SECURITY CLASSIFICATION

 OF THIS PAGE

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson

Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

Form Approved

OMB No. 0704-0188

12b. DISTRIBUTION CODE

8. PERFORMING ORGANIZATION

 REPORT NUMBER

5. FUNDING NUMBERS

3. REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE

6. AUTHOR(S)

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

11. SUPPLEMENTARY NOTES

12a. DISTRIBUTION/AVAILABILITY STATEMENT

13. ABSTRACT (Maximum 200 words)

14. SUBJECT TERMS

17. SECURITY CLASSIFICATION

 OF REPORT

16. PRICE CODE

15. NUMBER OF PAGES

20. LIMITATION OF ABSTRACT

Unclassified Unclassified

Contractor Report

Unclassified

1. AGENCY USE ONLY (Leave blank)

10. SPONSORING/MONITORING

 AGENCY REPORT NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546–0001

Institute for Computational Mechanics in Propulsion

22800 Cedar Point Road

Brook Park, Ohio 44142

Available electronically at http://gltrs.grc.nasa.gov

October 2002

NASA CR—2002-211889

ICOMP–2002–06

AIAA–2002–3712

E–13571

WU–708–48–13–00

NCC3–542

18

Computational Study of Near-Limit Propagation of Detonation

in Hydrogen-Air Mixtures

S. Yungster and K. Radhakrishnan

Unstable detonations; Detonability limits; Computational fluid dynamics

Unclassified -Unlimited

Subject Category: 34 Distribution: Nonstandard

Prepared for the 38th Propulsion Conference and Exhibit, cosponsored by the AIAA, ASME, SAE, and ASEE,

Indianapolis, Indiana, July 7–10, 2002. Project Manager, Charles Trefney, Turbomachinery and Propulsion Systems

Division, NASA Glenn Research Center, organization code 5880, 216–433–2162.

A computational investigation of the near-limit propagation of detonation in lean and rich hydrogen-air mixtures is

presented. The calculations were carried out over an equivalence ratio range of 0.4 to 5.0, pressures ranging from

0.2 bar to 1.0 bar and ambient initial temperature. The computations involved solution of the one-dimensional Euler

equations with detailed finite-rate chemistry. The numerical method is based on a second-order spatially accurate totalvariation-diminishing (TVD) scheme, and a point implicit, first-order-accurate, time marching algorithm. The hydrogenair combustion was modeled with a 9-species, 19-step reaction mechanism. A multi-level, dynamically adaptive grid was

utilized in order to resolve the structure of the detonation. The results of the computations indicate that when hydrogen

concentrations are reduced below certain levels, the detonation wave switches from a high-frequency, low amplitude

oscillation mode to a low frequency mode exhibiting large fluctuations in the detonation wave speed; that is, a

“galloping” propagation mode is established.