Hugh Douglas Perkins

Glenn Research Center, Cleveland, Ohio

Effects of Fuel Distribution on

Detonation Tube Performance

NASA/TM—2002-211712

August 2002
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

data bases, 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
Hugh Douglas Perkins

Glenn Research Center, Cleveland, Ohio

Effects of Fuel Distribution on

Detonation Tube Performance

NASA/TM—2002-211712

August 2002

National Aeronautics and

Space Administration

Glenn Research Center
Acknowledgments

Available from

NASA Center for Aerospace Information

7121 Standard Drive

Hanover, MD 21076

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22100

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

The Aerospace Propulsion and Power Program at

NASA Glenn Research Center sponsored this work.

I would like to acknowledge the patience and support of my thesis advisor Dr. C.J. Sung of Case Western Reserve

University in guiding this work. I would also like to acknowledge the NASA Glenn Research Center for its support

of my program of study and the development of this paper.
NASA/TM—2002-211712 1 

Effects of Fuel Distribution on Detonation Tube Performance 

Hugh Douglas Perkins 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Summary 

A pulse detonation engine (PDE) uses a series of high frequency intermittent detonation 

tubes to generate thrust. The process of filling the detonation tube with fuel and air for 

each cycle may yield non-uniform mixtures. Lack of mixture uniformity is commonly 

ignored when calculating detonation tube thrust performance. In this study, detonation 

cycles featuring idealized non-uniform H2/air mixtures were analyzed using the SPARK 

two-dimensional Navier-Stokes CFD code with a 7-step H2/air reaction mechanism. 

Mixture non-uniformities examined included axial equivalence ratio gradients, transverse 

equivalence ratio gradients, and partially fueled tubes. Three different average test 

section equivalence ratios (Φ ), stoichiometric (Φ = 1.00), fuel lean (Φ = 0.90), and fuel 

rich ( ) Φ = 1.10 , were studied. All mixtures were detonable throughout the detonation 

tube. It was found that various mixtures representing the same test section equivalence 

ratio had specific impulses within 1% of each other, indicating that good fuel/air mixing 

is not a prerequisite for optimal detonation tube performance. 

Nomenclature 

bi Fluid element body forces, ith dimension 

c Speed of sound 

di Dimensional characteristic boundary derivatives 

fi Mass fraction of the ith species 

hi Enthalpy of the ith species 

iref Centerpoint of X-grid spline 

k Thermal conductivity 

kb Backward reaction rate 

kf Forward reaction rate 

n Arrhenius expression temperature exponent 

ns Number of chemical species 

p Static pressure 

pL Local tube exit pressure 

po Ambient (atmospheric) pressure 

q Fluid element heat flux 

t Time 
NASA/TM—2002-211712 2 

tcycle Time to complete one detonation tube cycle 

u X – velocity 

ui Velocity in the ith direction 

uj Velocity in the jth direction 

uL Local tube exit velocity 

v Y – velocity 

wi Mass production of the ith species 

xˆ Unit vector in the X direction 

xi Spatial dimension in the ith direction 

xj Spatial dimension in the jth direction 

yˆ Unit vector in the Y direction 

Ai-Di Laminar viscosity coefficients 

A Arrhenius expression pre-exponential factor 

Cp Specific heat at constant pressure 

D Binary diffusion coefficient 

DTi Thermal diffusion coefficient of the ith species 

E Total energy of flow element 

Ea Arrhenius expression activation exponent 

F Vector of convective terms from governing equations 

Fˆ Vector of convective terms from governing equations, transformed coordinates 

G Vector of diffusive terms from governing equations 

Gˆ Vector of diffusive terms from governing equations, transformed coordinates 

H Detonation tube height 

H Vector of source terms from governing equations 

Hˆ Vector of source terms from governing equations, transformed coordinates 

I Total impulse 

Isp Specific impulse 

J Jacobian Matrix 

Keq Equilibrium species concentration 

L Detonation tube length 

Li Dimensional characteristic boundary wave amplitude terms 

Mi Molecular weight of the ith species 

Mj Molecular weight of the jth species 

Pcj Chapman-Jouget point pressure 

Po Detonation tube pressure prior to combustion 
NASA/TM—2002-211712 3 

Pp Plateau pressure 

Pr Prandtl Number 

Pvn Von Neumann spike pressure 

Ru Universal gas constant 

Sc Schmidt number 

T Static temperature 

U Vector of dependant variable from governing equations 

Uˆ Vector of dependant variable from governing equations, transformed coordinates 

V Velocity vector 

Vcj Chapman-Jouget equilibrium detonation velocity 

Vi Diffusion velocity of the ith species 

V j Diffusion velocity of the jth species 

Wfuel Detonation tube fuel weight per cycle 

X Length dimension 

Xi Mole fraction of the ith species or grid-point X-coordinate 

Xj Mole fraction of the jth species 

Y Height dimension 

ij δ Dirac delta function 

λi Dimensional characteristic boundary wave velocities 

μ Mixture viscosity 

μi Species viscosity of the ith species 

μ j Species viscosity of the jth species 

φ Local equivalence ratio 

φij Mixture viscosity weighting factor 

η 2nd transformed grid coordinate 

ρ Density 

ρi Species density of the ith species 

ρ j Species density of the jth species 

ξ 1st transformed grid coordinate 

ij τ Fluid element surface stress ( ) = τ

Φ Detonation tube test section average equivalence ratio (Alternately given as “Phi” 

in various numerical results plots) 
NASA/TM—2002-211712 4 

1.0 Introduction 

1.1 Motivation of Study 

The concept of a pulse detonation engine (PDE) for aerospace propulsion system 

applications is not new. The work at the University of Michigan in the 1950’s (Nichols 

et al., 1957) is a prime example. However, it was not until repetitive detonations with 

gaseous hydrocarbon fuels at frequencies at or above 25 detonations per second were 

demonstrated at the Naval Post-Graduate School in the mid-1980’s that it became 

apparent that practical devices might be possible, as described by Coleman (2001). 

These advances and others over the past decade have stimulated interest in the 

mainstream propulsion development community (Kailasanath, 1999; Kaemming, 2001). 

Both government and industry entities are currently considering the application of PDE’s 

and derivative configurations to missions as diverse as Single-Stage-To-Orbit (SSTO) 

launch vehicles, rocket upper stages, commercial aircraft, missiles, and tactical aircraft 

(Coleman, 2001). At the heart of each of these unique propulsion systems is a series of 

high frequency intermittent (pulse) detonation tubes. 

PDE’s are attractive for several reasons. First, Kailasanath and Patniak (1999), 

among others, have shown that the thermal efficiency of the detonation cycle can 

approach two times the thermal efficiency of the standard Brayton cycle found in gas 

turbine engines, 49% versus 27% in an example case where both combustion processes 

began at the same initial pressure. Second, since static pressure increases substantially 

during detonative combustion, instead of slightly decreasing as in deflagrative 

combustion, it is possible to decrease the amount of compression required before the 

combustion process, thus leading to simplified, lighter, lower cost engine architectures. 

For some missions, such as air-launched missiles, it is possible to eliminate all 

mechanical pre-compression and operate the PDE as a “supercharged” ramjet, even at 

relatively low subsonic Mach numbers. Lastly, a PDE detonation chamber does not have 

to be round, nor do the detonation chambers have to be grouped in circular arrays if 

turbomachinery components have been eliminated, thus leading to a more 

aerodynamically efficient airframe integration than with a gas turbine engine. These 

PDE advantages will become more apparent in the following sections as the detonation 

tube architecture and engine gas dynamics are explored. 

1.2 PDE Architecture Description 

 A variety of architectures have been proposed for pulse detonation engines, 

differing materially in tube geometry, air valve design, fuel injector placement, and other 

operational features. Figure 1 shows a generic configuration for an airbreathing pulse 

detonation tube. To begin a new cycle, the air valves on the left are opened to admit air 

from the inlet plenum. Fuel is injected into the incoming air stream as it enters the 
NASA/TM—2002-211712 5 

initiator tube and the main detonation tube. Typically, the start of fuel injection is 

delayed slightly after the air valves are opened in order to allow a thin buffer of unfueled 

air to form between the exhaust products from the previous cycle and the fresh charge, 

preventing pre-ignition. After the inlet fuel and air valves are closed, combustion is 

initiated using the spark plug in the initiator tube. The ensuing deflagration quickly 

accelerates into a detonation, aided by the small diameter of the initiator tube and 

typically by some amount of oxygen enrichment in the initiator tube. 

Figure 1 – Basic Pulse Detonation Combustor Tube Architecture 

 Once formed, the detonation wave travels down the initiator tube at a 

characteristic detonation wave speed in excess of 1500 m/s, depending on the fuel and 

oxidizer combination used and the initial conditions (pressure, temperature, 

stoichiometry). The detonation wave then propagates out of the initiator tube, detonating 

the main charge. Care must be taken in selecting appropriate geometric parameters (area 

ratio, rate of area change, cross section) that will allow the transition of the detonation 

wave into the larger tube. In some cases, the detonation may initially fail 

(shock/combustion decoupling) upon exiting the initiator tube, but will then re-initiate 

from hot spots in the main tube caused by intersecting shocks. 

 Two additional features of the detonation tube shown in Figure 1 that aid in the 

initiation and propagation of the detonation wave are the convergent-divergent (C-D) 

nozzle and the actively cooled combustor walls. The nozzle, when employed, is used to 

provide a higher initial (pre-detonation) pressure in the tube by restricting the outflow 

from the tube during the fill process. This higher initial pressure aids in the initiation of a 

detonation, since detonability is a function of both pressure and temperature (Kuo, 1986). 

A higher initial pressure also increases thrust potential due to the corresponding increase 

in peak and plateau pressures caused by the pressure multiplying effect of the detonation 

wave. In the common configuration, the detonation cycle is also enabled by the presence 

of the actively cooled combustor walls. The wall temperature must be maintained below 

the autoignition temperature of the fuel/air mixture being used to prevent deflagrative 

Actively Cooled Combustor 

Air Inlet Valves 

 Spark Ignitor 

Potential Fuel 

Injection 

Locations 

Fuel/Oxidizer 

Inlet Plenum 

Initiator 

Nozzle 

(Optional) 
NASA/TM—2002-211712 6 

mode ignition prior to the initiation of the detonation. However, any mixing or buffering 

scheme that keeps the flammable mixture away from the wall could also be used to 

prevent autoignition. 

1.3 Detonation Wave and PDE Wave Cycle Descriptions 

 Once fully initiated, the detonation wave assumes the characteristic Zeldovich, 

Von Neumann, Doring (ZND) structure, shown in the fixed, or laboratory, coordinates in 

Figure 2, taken in part from Kuo (1986). The leading shock compresses the fuel/oxidizer 

mixture (A) to a pressure and temperature sufficient for rapid combustion (B). This point 

is known as the Von Neumann spike, and is approximately equal to the pressure and 

temperature rise observed across a frozen shock moving at the same Mach number as the 

detonation wave. After a short kinetic delay, combustion occurs with a characteristic 

decrease in static pressure. Point (C) corresponds to the point at which combustion is 

complete, know as the Chapman-Jouget (CJ) point. The CJ point is uniquely defined by 

the pre-detonation conditions and represents the point at which the local Mach number in 

the combustion products equals 1 in a reference frame fixed to the detonation wave. 

While it is possible to have a detonation which does not end in the CJ point, such 

detonations are generally unstable, and so the leading shock will typically self-adjust its 

speed (and therefore strength) so as to create the CJ point at the end of the combustion 

zone. The resulting shock Mach number, referenced to the speed of sound in the 

unburned gas, is known as the CJ Mach number, and is again uniquely defined by the 

unburned conditions. 

Figure 2 – ZND Detonation Structure in a Tube (Laboratory Coordinates) 

Propagation 

Pressure 

X/L 1

Po

Pp

Pcj 

Pvn

D A

C

B

E
NASA/TM—2002-211712 7 

 

In order to maintain continuity in the closed tube, a series of expansion waves (D) 

propagate toward the closed end of the detonation tube from the detonation wave front to 

reduce the pressure from the CJ pressure to what is commonly referred to as the plateau 

pressure (E), reconciling the high velocity flow following the detonation wave with the 

zero velocity requirement at the closed end wall. The tube then exhausts through the 

open end, with the blow down characteristics governed by the reflected expansion and 

compression waves propagating through the tube. Once the tube has reached a pressure 

sufficiently close to ambient pressure, the air inlet valves are opened to begin recharging 

the detonation tube leading to the initiation of a new cycle. The overall detonation cycle 

is shown graphically in Figure 3. 

2 Detonation Initiated

3 Detonation Transitions 

Into Main Tube

4 Detonation Propagates 

Through Length Of Main Tube 5 Tube Blows Down To 

Ambient Pressure 

1 New Fuel/Air Mixture Introduced,

Purging Out Remaining Hot Products

Figure 3 – Pulse Detonation Combustor Operational Cycle 

1.4 Scope of Numerical Study 

It is likely that in any application of detonative combustion for propulsion, the 

fuel and air will not be uniformly mixed, as is often assumed. Non-uniformity in fuel 

distribution may be created intentionally or may result from hardware limitations, and 

may be either axial or transverse, or both, depending on the source. Possible sources of 

non-uniform fuel distribution include, but are not limited to, 

• Combustion products/fresh charge buffering (purge) 

• Engine throttling 

• Emissions control 

• Wall heat transfer control 

• Detonation wave shape/strength control 

• Air and fuel valve transients 

• Non-uniform inlet air flow in space and time 

The study that follows examines the effects of various idealized fuel distribution 

non-uniformities on detonation tube thrust performance (specific impulse). Kailasanath, 

et al. (2000) studied the performance effects of unfueled purge fraction, alternately 
NASA/TM—2002-211712 8 

referred to as partial fill, showing an increase in Isp with increasing purge fraction 

(decreasing fill fraction). Their study was conducted using a transversely uniform 

stoichiometric mixture of ethylene and air in the fueled portion of the tube, leaving the 

exit end of the tube unfueled. No spatial transition (gradient) in equivalence ratio was 

used between the two regions. Purge fractions of up to 80% were studied. This current 

study builds upon and expands the consideration of fuel distribution effects by looking at 

the effects of equivalence ratio gradients in both the axial and transverse directions as 

well as comparable partial fill cases. The numerical results are intended to be 

comparative in nature, looking for trends and relative magnitude effects between the 

different cases. The actual values of performance are not validated by test data, and so 

may be subject to different levels of offset error. However, it is expected that the 

physics-based numerical modeling techniques employed will provide correct insight into 

the relative performance effects of the different fuel distributions. 

1.5 Paper Organization 

After briefly introducing the basic concept and operation of a PDE in the 

preceding sections, this thesis presents the specific geometry, operational conditions, and 

boundary conditions being addressed in the computational study, including all the unique 

formulations added to the existing SPARK code. The SPARK code itself is then briefly 

described, along with the basic implementation of the computational grids used in the 

study. Next, the test cases used to check the accuracy of SPARK for calculating H2/air 

detonations are described and the corresponding results presented. 

With the code thus validated, the results of the various cases used for the study are 

presented and compared. These results are then summarized and appropriate 

conclusions drawn. Lastly, an appendix is provided justifying the selection of the 

chemistry model chosen for this study. 

2. Problem Formulation

2.1 Geometry 

 The geometry of a detonation tube is generally simple outside of the transition 

section between the initiator tube and the main detonation chamber. While the geometry 

of the area transition is critical to operability of the device, this study will be confined to 

the constant-area main detonation chamber. This region represents the majority of the 

volumetric capacity and therefore energy release within the device, making it the 

dominant region in terms of thrust performance. To keep the geometry as simple as 

possible, the left hand side inflow plane (main air inlet valves – Figure 1) will be 

modeled as a solid wall, simulating a transitioned detonation with the inlet valves closed. 
NASA/TM—2002-211712 9 

Typical proof-of-concept PDE detonation tubes being used today are in the range 

of 5 to 12 cm in diameter and 120 to 180 cm in overall length (including the initiator 

tube). For this study a height of 3.0 cm and a length of 30.0 cm will be used in order to 

limit the computational time required. This basic two-dimensional geometry is shown in 

Figure 4. 

Figure 4 – Geometry for Numerical Study 

2.2 Simplifying Assumptions 

A number of simplifying assumptions are made in this study relative to the actual 

PDE cycle, as follows. 

2.2.1 DDT 

No attempt is made to model the deflagration to detonation transition (DDT) 

process, as the phenomena under study is primarily concerned with the steady state 

detonation propagation. As shown in Figure 4, the initiator tube itself is not modeled nor 

is the area change from the initiator to the main chamber modeled. A successfully 

transitioned detonation wave is numerically initiated through the imposition of a high 

temperature, high pressure ignition zone at the left-hand (inflow) side of the domain in a 

5.0 cm initiation/stabilization region. The high pressure, high temperature ignition zone 

is set uniformly to 150 atmospheres, 4000 K, and is the full height of the tube, but only 

0.6 mm long, with a composition of 25.5% H2O and 74.5% N2 by weight. The total 

energy contained in the ignition region is 335 J, which is approximately 1.3% of the heat 

release from the full 30 cm detonation tube when it is uniformly, stoichiometrically 

fueled (26.2 KJ). The remaining 4.94 cm stabilization region is maintained as a uniform, 

stoichiometric unburned mixture for each case to aid in establishing a successful 

detonation regardless of test case mixture properties. Although this detonation may be 

initially overdriven to insure propagation, it falls back to near the equilibrium ChapmanJouget condition by the end of the stabilization zone. H = 3 cm 

Stabilization Region ( 0 φ = 1. ) 

L = 30 cm 

Ignition Region (0.6 mm) 

Test Section – Tailored Mixtures 

 5 cm 
NASA/TM—2002-211712 10 

2.2.2 Partial Cycle 

Only the detonation propagation and blowdown portions of the full PDE cycle, 

portions 4 and 5 in Figure 3, are modeled. It is assumed that the purge and refresh 

cycles are identical for each case. While this might not be rigorously true in practice, 

such differences would be expected to be of significantly less importance. It is further 

assumed that the blowdown of the tube is complete when the tube exit velocity reaches 

zero anywhere along the exit plane. In a real cycle, some backflow is likely to occur due 

to the sub-ambient pressure in the tube from the last reflected expansion wave. This 

portion of the cycle is important when modeling the filling of the tube, both in terms of 

the fill time and the performance penalty, but this issue is beyond the scope of this study. 

2.3 Boundary Conditions 

 The left-hand, top and bottom walls of the two-dimensional detonation tube are 

all treated as impermeable, dictating that the pressure derivative, the species 

concentration derivatives, and the normal component of velocity be set to zero at each 

wall. No-slip wall and adiabatic conditions are set at each surface by setting the 

temperature derivative and parallel component of velocity to zero. The adiabatic wall 

condition was selected to simplify the analysis by not requiring resolution of the 

boundary layer region, however the no-slip wall was maintained to better capture any 

wave reflections that might occur. 

 The tube exit boundary conditions for this type of calculation have been a source 

of considerable study, with various approaches having been used by different individuals, 

such as Kailasanath and Patnaik (1999) and Wilson and Paxson (2001). Over the course 

of each cycle, there are six segments of time that need to be addressed when considering 

the outflow boundary conditions. 

1. During the wave propagation segment, the detonable mixture between the 

detonation wave and the tube exit is undisturbed by the approaching detonation 

wave. In all of the cases considered in this study, the tube fill pressure matches 

the ambient pressure, so there is no outflow during this time segment. In most, if 

not all, real cases, there will be some low velocity flow into or out from the tube. 

2. As the detonation wave passes out of the tube, the flow at the exit plane quickly 

accelerates to a supersonic condition, as the flow between the leading shock and 

the CJ point is characteristically supersonic. 

3. After this short period of supersonic flow, the high pressure behind the reaction 

zone chokes the flow at the tube exit (Mach number equal to 1). 

4. The flow remains choked until a region of supersonic flow generated by reflecting 

gas dynamic waves reaches the tube exit, causing the exit flow from the tube to be 

supersonic for a second time. 

5. The flow then returns to the choked flow condition until the pressure has dropped 

sufficiently for the flow to become subsonic. 

6. Subsonic outflow throughout the rest of the blowdown. 
NASA/TM—2002-211712 11 

Figure 5 shows this sequence of events at the detonation tube centerline from the baseline 

stoichiometric H2/air detonation calculation (Case 1) described in detail in Section 5. 

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

0 0.0002 0.0004 0.0006 0.0008 0.001

T im e (s)

Figure 5 – Time History of Detonation Tube Exit Flow Mach Number 

at the Tube Centerline 

This result is consistent with the results of Wilson and Paxson (2001) from their 

1-dimensional calculations. 

 The outflow boundary conditions required for each of these situations were 

formulated as follows. 

2.3.1 Subsonic Outflow 

For subsonic outflow, the methodology of Poinsot and Lele (1992) for subsonic 

reflecting outflows for the Navier-Stokes equations was adopted, from which the 

following notation is adapted. 

First, the basic characteristic velocity terms are calculated from the existing values at 

the boundary, with the downstream direction defined as the direction of flow. 

= u − c λ1 : The characteristic velocity of the upstream traveling wave (1) 

Choked Flow 

(Mach # =1)

Time (s) 

Tube Exit Local Flow Mach Number 
NASA/TM—2002-211712 12 

λ2 = u : The characteristic convection velocity (2) 

λ3 = u : The characteristic velocity at which the y velocity component is 

advected downstream (3) 

= u + c λ5 : The characteristic velocity of the downstream traveling wave (4) 

Next, terms associated with the variation of the wave amplitude (

t

u

c

t

p

∂

∂ − ∂∂

ρ ) for 

each velocity term are calculated using the 1-dimensional Euler equations. The inviscid 

Euler equations are the best available approximation for the viscous (Navier-Stokes) case, 

and insure compatibility between the physical boundary conditions being calculated and 

the amplitudes of waves crossing the boundary. 

dt

dp L1 = −L5 − 2 : The amplitude term for the characteristic upstream traveling wave (5) 





 





∂

∂ − ∂∂ = Xp

X

L c ρ λ 2

2 2 : The amplitude term for the wave associated with the 

characteristic convection velocity (6) 





 





∂

∂ = X

v L3 λ3 : The amplitude term for the wave associated with the 

characteristic velocity at which the y velocity component is advected downstream (7) 





 





∂

∂

+

∂

∂ = X

u

c

X

p L5 λ5 ρ : 

The amplitude term for the characteristic downstream traveling wave (8) 

The λ4 and L4 terms are omitted as they are used in reference to the 3rd spatial 

dimension, not required for this analysis. 

Next, the following collection of derivatives at the boundary are cast in terms of 

the wave amplitude associated terms for use in simplifying the governing equations. 

( ) 2 2 5 1 1 2

( ) 1 1 L L L d

X c

u = 



 



 = + +

∂

∂ ρ (9) 

( ) 5 1 3 2

1 1 L L d

X c

p

X

u

u = − = ∂

∂

+

∂

∂

ρ ρ

 (10) 
NASA/TM—2002-211712 13 

L3 d4 X

v

u = = ∂

∂ (11) 

The derivatives d2 , associated with the energy equation, and d5 , associated with the 

third spatial dimension, are not required for the subsonic reflecting outflow case and are 

therefore omitted here. 

Lastly, these derivatives are used in the reformulated version of the 2-dimensional 

Navier-Stokes equations below to advance the desired quantities to the next time step at 

the tube exit, assuming constant tangential stresses at the boundary. 

( )

1 = 0

∂

∂

+ +

∂

∂

Y

v d

t

ρ ρ (continuity equation) (12) 

( ) ( )

Y X

uv ud d

t

u

∂

∂ = ∂∂

+ + +

∂

∂ 11

1 3

ρ τ

ρ ρ (X-momentum equation) (13) 

( ) ( )

Y Y

p

Y

v vd d

t

v

∂

∂ = ∂∂

+

∂

∂

+ + +

∂

∂ 22

2

1 4

ρ τ

ρ ρ (Y-momentum equation) (14) 

where 





 



 



 





∂

∂

+

∂

∂ − ∂∂ = Yv

X

u

X

u

3

2

τ 11 μ 2 (15) 

and 





 



 



 





∂

∂

+

∂

∂ − ∂∂ = Yv

X

u

Y

v

3

2

τ 22 μ 2 (16) 

By first solving for the density using equation (12), the u and v velocity terms can 

then be found using equations (13) and (14), respectively. For subsonic outflow, the 

static pressure is known, and is taken to be constant at 1 atmosphere for all the cases run 

during this study. Wilson and Paxson (2001) have shown that for tubes with a 

length/diameter (L/D) ratio greater than approximately 10, this is a valid assumption. 

This is because the residual effects of the blast wave on the static pressure in the vicinity 

of the tube exit have dissipated by the time the detonation tube unchokes during the last 

part of the blowdown, provided the tube pre-detonation pressure is approximately equal 

to the ambient static pressure. Lastly, the temperature at the exit boundary is determined 

by extrapolation from interior points assuming a constant local heat flux normal to the 

tube exit boundary. 
NASA/TM—2002-211712 14 

2.3.2 Supersonic Outflow 

 The values of p, T, u, and v are linearly extrapolated from the 2 immediate 

interior points, as described by Anderson (1995). The density, ρ , is then calculated 

using the equation of state. 

2.3.3 Choked Outflow 

 The values of p, T, ρ and v are handled as in the supersonic outflow case, but 

with u limited to the local speed of sound. 

2.4 Detonation Tube Performance 

 There are two fundamental methods for calculating the impulse generated by the 

detonation tube, as taken from Oates (1984). The first is to integrate the pressure along 

the left hand wall over the partial cycle modeled, subtracting the frictional pressure loss 

along the two side walls, 

dX dt

Y

u dX

Y

u I pdY

t

bottom

L

wall

top

L

wall

left

H

∫ ∫ ∫ ∫ 































∂

∂ − 















∂

∂ − 













 =

0 0 0 0

μ μ (17) 

The second method is to integrate the momentum flux two-dimensionally at the tube exit 

over the partial cycle, correcting for non-ideal expansion, 

I u dY p p dY dt

t H

L

H

∫ ∫ ∫ L 



 



 = − −

0 0

0

0

2

ρ ( ) (18) 

Both methods, shown per unit width of detonation chamber, will be implemented in the 

computer code and the relative results will be compared. The calculated impulse can 

then be used to calculate the specific impulse of the detonation tube. 

fuel

sp W

I I = (19) 

where Wfuel is the weight of fuel in the tube for one cycle. 
NASA/TM—2002-211712 15 

2.5 Fuel-Air Combinations 

 A stoichiometrically fueled 5.0 cm detonation initiation section is used to provide 

a stabilized detonation wave to the test section, as described in Section 2.2. The test 

section equivalence ratio is varied in either the X-direction or the Y-direction from case 

to case. Three test section average equivalence ratios, alternately denoted as Φ or “phi”, 

are used, one stoichiometric, one fuel-lean, and one fuel rich. Non-stoichiometric 

equivalence ratios of 0.9 and 1.1 were selected, representative of moderate levels of nonuniformity, but still well within expected stable detonation limits throughout the 

detonation tube. Different buffers or gradients of fuel/air mixture are used to achieve the 

desired average test section equivalence ratio. Transverse gradient cases are 

representative of a tube centerline fuel injection configuration, while axial gradient cases 

are representative of a non-constant tube filling process brought on by time varying tube 

inlet and exit conditions. Partial fill cases are representative of fuel lead (H2 buffer) or 

fuel lag (air buffer) injection timing and provide a point of comparison to the previously 

referenced work of Kailasanath, et al. (2000). Baseline uniformly fueled cases are used 

to normalize the results of the subsequent cases. The test cases are grouped as follows, 

with Φ referencing the test section average equivalence ratio and φ referencing the local 

equivalence ratio within the test section. 

Stoichiometric Combustion (Φ = 1.00 )

 Case 1 – Baseline – uniformly fueled throughout test section 

Case 2 – Linear transverse gradient - φ = 1.10 at the centerline; φ = 0.90 at each wall 

Case 3 – Linear axial gradient - φ = 1.10 at the closed end (after the initial 

stabilization zone); φ = 0.90 at the open end 

Fuel-Lean Combustion (Φ = 0.90 )

Case 4 – Baseline – uniformly fueled throughout test section 

Case 5 – Linear transverse gradient - φ = 1.00 at the centerline; φ = 0.80 at each wall 

Case 6 – Linear axial gradient - φ = 1.00 at the closed end (after the initial 

stabilization zone); φ = 0.80 at the open end 

Case 7 – 2.50 cm air buffer at the end of the test section; φ = 1.00 up to the buffer air 

in the test section 

Fuel-Rich Combustion (Φ =1.10 ) 

Case 8 – Baseline – uniformly fueled throughout test section 

Case 9 – Linear transverse gradient - φ = 1.20 at the centerline; φ = 1.00 at each wall 

Case 10 – Linear axial gradient - φ = 1.20 at the closed end (after the initial 

stabilization zone); φ = 1.00 at the open end 

Case 11 – 0.07 cm H2 buffer at the end of the test section; φ = 1.00 up to the buffer 

H2 in the test section 
NASA/TM—2002-211712 16 

Figures 6 and 7 show graphically the different axial and transverse fuel distribution 

schemes to be used in the “test section” of the detonation tube, excluding the baseline 

uniformly fueled cases (Cases 1, 4, and 8) and the partial-fill cases (Cases 7 and 11). 

Test cases performed using the opposite orientation of fuel distribution nonuniformity (transverse distributions that were fuel lean at the centerline and fuel rich at 

the walls, and axial distributions that increased in equivalence ratio along the length of 

the tube instead of decreasing) yielded nearly identical detonation tube performance 

results as the cases shown above, and so were not included in this study. 

X (cm)

0 5 10 15 20 25 30

E

quivale

nce Rati

o

0.7

0.8

0.9

1.0

1.1

1.2

1.3

Case 3 - Stoichiometric - Phi=1.0

Case 6 - Fuel Lean - Phi=0.9

Case 10 - Fuel Rich - Phi=1.1

Figure 6 – Test Section Axial Fuel Distributions 
NASA/TM—2002-211712 17 

Y (cm)

0123

Equivalence Ratio

0.7

0.8

0.9

1.0

1.1

1.2

1.3

Case 2 - Stoichiometric - Phi=1.0

Case 5 - Fuel Lean - Phi=0.9

Case 9 - Fuel Rich - Phi=1.1

Figure 7 – Test Section Transverse Fuel Distributions 

 

3. Numerical Tool Description 

 The numerical tool chosen for this study was the SPARK 2D Navier-Stokes Code 

developed at the NASA Langley Research Center. The SPARK code was developed 

primarily for the study of high speed reacting flows, particularly those found in scramjet 

engines, under the National AeroSpace Plane (NASP) program. SPARK is a versatile 

time accurate code that has been used in a number of supersonic combustion and 

detonation studies (Drummond and Mukunda, 1988; Drummond, 1988; Carpenter, 1989). 

The code has been previously described in detail by Drummond (1990), so this 

description will be limited to a general overview of the principal features salient to this 

study. 

 

3.1 Governing Equations 

SPARK uses the strong conservation form of the Navier-Stokes, energy, and species 

continuity equations as derived by Williams (1985). 
NASA/TM—2002-211712 18 

Continuity: 

+ ∇ ⋅( ) V = 0

∂

∂

ρ ρ

t

 (20) 

Momentum: 

( ) ( ) ∑=

+ ∇ ⋅ = ∇ ⋅ +

∂

∂ ns

i

i i VV f b

t

V

1

ρ τ ρ ρ (21) 

where τ is here defined as 

















∂

∂

+

∂

∂ = = − +

i

j

j

i

ij ij x

u

x

u

τ τ δ p μ (22) 

Energy: 

( ) ( ) ( ) ∑ ( ) =

+ ∇ ⋅ = ∇ ⋅ ⋅ − ∇ ⋅ + +

∂

∂ ns

i

i i VE V q fi b V V

t

E

1

ρ τ ρ ρ (23) 

where 

( ) i j

ns

i

ns

j i

T i

u

ns

i

i i i V V

M D

XjD

q k T h f V R T G G G − 













 = − ∇ + ∑ + ∑∑ =1 = = 1 1

ρ (24) 

2

2 2

1

p u v E h f

ns

i

i i

+ = ∑ − + = ρ

 (25) 

and 





 





∂

∂

+

∂

∂ = − y Y

f

x

X

f

f

D V i i

i

i ˆ ˆ G (26) 

Species Continuity: 

( ) ( ) ( )i i i i

i V f w f V

t

f ρ ρ ρ + ∇ ⋅ = − ∇ ⋅ ∂∂ (27) 
NASA/TM—2002-211712 19 

The expression for diffusion velocity, Vi

G

, assumes that N2 is the major species present, 

with all other species present in “trace” amounts. It can then be assumed that only the 

diffusion of each species with N2 be considered, and further that all binary diffusion 

coefficients are equal, as described by Drummond (1990). 

 Anderson (1995) has pointed out that when the governing equations are cast in the 

strong conservation form as above, there is no issue with capturing shocks in the flow 

field. Since the equations are set up in terms of fluxes instead of the primitive variables 

( ) ρ,u,v, p,etc. , there are no discontinuities across shock waves within the governing 

equations that require separate consideration for closure. 

 

 

3.2 Detailed Chemistry 

 Following the work of Cambier and Adelman (1988) and Lynch and Edelman 

(1994), a 7-step, 7-species hydrogen-air detailed chemical mechanism was chosen for this 

study instead of the more standard 18-step, 9-species mechanism of Jachimowski (1988). 

For computational efficiency, it is desirable to utilize the simplest possible mechanism 

that can with reasonable accuracy provide the correct time evolution and chemical 

composition of the detonation wave and the resultant combustion products. The 7-step 

mechanism is a simple reduction of the 18-step mechanism through the elimination of 

HO2 and H2O2 from the chemical system, leaving H2, O2, H2O, OH, H, O, and N2 (inert). 

This simplification can be justified when it is considered that in the basic ZND detonation 

model, ignition occurs after the leading shock, at which point the flow temperature is 

elevated to a point where the chemistry of HO2 and H2O2 is no longer dominant. 

The Appendix provides further analysis to justify the use of the simplified chemical 

mechanism. 

 Each chemical reaction is modeled using the standard Arrhenius expression for 

the forward reaction rate, 

R T

Ea

n

f u k A T e

−

= ∗ ∗ (28) 

For the 7-step mechanism used in this study the values of A, n, and Ea are as follows. 
NASA/TM—2002-211712 20 

Table 1 – 7-Step H2-Air Mechanism 

Reaction 

Ref. # 

Reaction A n Ea 

(cal/g-mole) 

1 H2 + O2 ⇔ 2OH 0.1700e14 0.00 48150 

2 H + O ⇔ OH + O 2 0.1420e15 0.00 16400 

3 H2 + OH ⇔ H2O + H 0.3160e08 1.80 3030 

4 O + H2 ⇔ OH + H 0.2070e15 0.00 13750 

5 OH + OH ⇔ H2O + O 0.5500e14 0.00 7000 

6 OH + H + M ⇔ H2O + M 0.2210e23 –2.00 0 

7 H + H + M ⇔ H2 + M 0.6530e18 –1.00 0 

The 3rd body efficiencies of all species for reactions 6 and 7 are assumed to be one. The 

backward rate constant is then found from the simple expression, 

eq

f

b K

k

k = (29) 

where Keq is found using a Gibbs free energy minimization scheme similar to those found 

in chemical equilibrium codes such as the NASA Glenn Research Center’s CEA code 

developed by Gordon and McBride (1994). The chemical system involving the 6 

reactive species is solved at each time step, thus providing the source term for the species 

equation (27). The mass of nitrogen present is determined by summing the mass 

fractions of each of the other species and then subtracting the result from the necessary 

total of 1. Mass conservation is maintained throughout the computational domain at each 

time step. 

3.3 Thermodynamic and Transport Properties 

 SPARK utilizes previously developed analytical expressions for the required 

thermodynamic and transport properties used throughout the code, which, as with the 

previously discussed aspects of the code, have been related in detail in other publications, 

such as Drummond (1990). Species specific heat at constant pressure and species Gibbs 

energy are calculated using 4th and 5th order polynomial curve fits of temperature, 

respectively. The species laminar viscosity, μi , is calculated using an expression known 

as Sutherland’s law, 
NASA/TM—2002-211712 21 

i

i i

i

i i T D

C D

B

T A

+

+















 =

1.5

μ , (30) 

where Ai, Bi, Ci, and Di are constants which vary with species. 

The mixture value of viscosity is then calculated using Wilke’s law. 

∑

∑ 













 +

=

= ≠

ns

ij

ns

j j i

j

i

i

X

X

1

1,

1 1 φ

μ μ (31) 

where 

0.5

2 0.5 0.25

1

2

4

1















 +

















































+

=

j

i

j

i

j

i

j

i

ij

M

M

M

M

ρ

ρ

μ

μ

φ (32) 

 Lastly, the binary diffusion coefficient, D, and the thermal conductivity, k, are 

calculated from the Schmidt number (0.22) and Prandtl number (0.72), respectively, 

which are provided as inputs to the code. 

Sc

D

ρ

μ = , 

Pr

Cpμ k = 

All properties were calculated assuming laminar flow throughout. 

3.4 Computational Grid 

Looking back at the detonation structure shown in Figure 2, it is easily seen that 

the region in the immediate vicinity of the detonation wave front will determine the 

required physical grid spacing for the computational model. If the goal of this study were 

to examine the detailed structure of the shock wave coupled to the reaction zone, it would 

be necessary to provide a computational grid spacing on the order of the width of the 

leading shock wave (shown by Shapiro (1953) to be on the order of 1.4×10–7 m), and 

probably finer. Fortunately, it is not necessary to resolve this structure in order to obtain 

accurate flow properties required to calculate the performance values that are the goal of 
NASA/TM—2002-211712 22 

this study. It is sufficient to provide enough grid resolution to capture the shortest 

reaction times of the chemical system as the flow behind the shock is convected 

downstream. For the 7-step chemical system and conditions investigated in this study, 

this worked out to be on the order of 1×10–5 m. Grid spacing significantly larger than 

this, on the order of 1×10–4 m , resulted in oscillations in the combustion zone. When 

even courser grids were tried, the calculation went completely unstable. From a practical 

standpoint, it is desirable to use the largest acceptable grid spacing in order to reduce the 

number of required grid points. Also, since the system of equations being solved are 

numerically “stiff”, Radhakrishnan (1991) has shown that going to a finer grid spacing 

results in a proportionally smaller required time step, leading to very long computational 

times. 

 Since an X-grid spacing of 1×10–5 m would yield 30,000 grid points in the xdirection for each Y-grid location for a 30 cm tube, it becomes immediately apparent that 

this grid spacing needs to be confined to the region where it is required, namely the 

leading section of the detonation wave where steep gradients in all the principal flow and 

composition variables are present. Since the detonation wave is moving in space, the 

fine grid must move with it. This was accomplished for this study by starting with a 

uniform 300-point grid in X at the fine grid spacing (1×10–5 m). As the detonation wave 

approaches the right-hand edge of the computational domain, the grid is stretched by 

approximately 30 times the fine grid spacing at the left hand side where the reaction has 

already occurred, lengthening the overall domain in X. The fine grid section and the 

coarser grid section are connected through a spline fit, 

( )

1 0.05555( ) 1 ref i i

coarse fine

i i coarse

e

X X

X X X − − − +

∆ − ∆ = + ∆ − , (35) 

where ref i is the center grid point of the spline, which is set to 60 for all the simulations in 

this study. This value of ref i gives approximately 100 grid points at the fine grid spacing 

at the right hand side of the domain. The wave structure is thus always maintained within 

the fine grid section, which is effectively shifted to the right every time the detonation 

wave travels 30 times the fine grid spacing. The current result is then interpolated onto 

the new grid and the calculation is continued. Once the desired length of the full 

detonation tube is reached, the fine grid is maintained at the right hand (exit) end of the 

tube. Figure 8 shows the physical X-grid spacing at detonation wave positions (domain 

lengths) of 3 cm, 10 cm, and 30 cm. 
NASA/TM—2002-211712 23 

0.000001

0.00001

0.0001

0.001

0.01

0 50 100 150 200 250 300

X Grid P oint

Grid S

p

a

cin

g (m)

Do main Len gth = 30 cm

Do main Len gth = 10 cm

Do main Len gth = 3 cm

Figure 8 – Variable X-Grid Spacing 

 The Y-grid spacing remains the same throughout each calculation. Grid points 

are clustered near the wall in order to capture with reasonable fidelity the severe gradient 

of velocity generated by the detonation wave in that region. As with the X-grid, 

significantly finer grid spacing yields proportionally smaller time steps, leading to 

excessively long computational times. Grid spacings as small as 5×10–6 m did not 

require time steps smaller than those required for the finite-rate chemistry. The general 

Y-grid distribution used throughout the study is shown in Figure 9. Minor adjustments to 

this grid were used case-to-case to minimize run times while maintaining detonation 

wave stability. 

0.000001

0.00001

0.0001

0.001

0.01

0 20 40 60 80 100

Y Grid Point

Grid S

pacin

g (

m)

Figure 9 – Y-Grid Spacing 
NASA/TM—2002-211712 24 

Once the grid has been generated, the method of Smith and Weigel (1980), 

referred to as the “two-boundary technique,” is used to transform the mesh into a uniform 

2-dimensional boundary-fitted rectangular grid ( ) ξ, for computation. The disconnected η

domain boundaries are defined by analytic functions, which are then connected together 

using parametric polynomials to generate the grid. This transformation allows for more 

efficient computation of the time-accurate solution, which will be further elaborated on in 

the next section. The solution is then transformed back into the original physical space 

for output. 

3.5 Calculation Implementation and Numerical Solution 

Once the geometric domain has been discretized into the computational grid, the 

governing equations are vectorized for solution in the form, 

H

Y

G

X

F

t

U = ∂∂

+

∂

∂

+

∂

∂ (36) 

where the U vector represents the dependant variables, the F and G vectors the 

convective and diffusive terms in X and Y, respectively, and the H vector the source 

terms. In the transformed coordinate system, equation (36) becomes, 

H GF

t

U ˆ ˆˆ ˆ

= ∂

∂

+

∂

∂

+

∂

∂

ξ η

 (37) 

where 

UJU

G G

ˆ = , (38) 

G X F Y F G G G

η ∂η

∂ − ∂∂ ˆ = , (39) 

F Y G X G G G G

ξ ∂ξ

∂ − ∂∂ =ˆ , (40) 

HJH

G G

ˆ = (41) 

and where 
NASA/TM—2002-211712 25 

ξ η ξ ∂η

∂

∂

∂ − ∂∂

∂

∂ = X Y Y X J (42) 

 This vectorized equation is then solved using a MacCormack explicit predictorcorrector scheme as described by Drummond et al. (1991), that is 4th order accurate in 

space and 2nd order accurate in time. While the code allows for both partially and fully 

implicit schemes to be used, these options were not selected for this study in order to 

maintain the greatest accuracy possible. The time step was limited so as to allow a 

maximum change in any species mass fraction of 0. 8% per iteration at any grid point. 

With these options selected the typical simulation contained in this study required 60 

hours to complete on a Silicon Graphics Octane2 workstation with a 400 MHz R12000A 

processor. 

4. Code Validation 

4.1 Introduction 

Before beginning to utilize the SPARK code for the configuration and conditions 

of interest, it is necessary to demonstrate that the code can accurately compute one or 

more known solutions in similar flow regimes. The SPARK code has been previously 

validated for such cases as 

• Non-reacting 

o Laminar spatially developing shear layer 

o Turbulent spatially developing mixing layer 

o Mach 14 laminar 15 degree compression corner 

o Transverse jet behind a rearward facing step 

• Reacting 

o H2-air turbulent mixing and reaction 

o Burrows-Kurkov Mach 2.44 H2 fueled vitiated scramjet combustor 

All of these cases are described in Drummond (1991). 

 Additionally, prior to initiating the detonation computations included in this 

study, several standard non-reacting test cases were performed, including a supersonic 

expansion corner, a supersonic compression corner, and a H2-N2 shock tube problem. In 

all cases the results were consistent with the corresponding classical analytical solution. 

In addition to the above validations, it was necessary to perform validation calculations 

specifically for detonation. Two cases were selected for which known solutions with 

detailed kinetic mechanisms exist or could be easily calculated using other validated 

codes, as follows. 
NASA/TM—2002-211712 26 

4.2 H2-Air Oblique Detonation Wave 

Thaker and Chelliah (1997) used the SPARK code to calculate the detonation 

structure of an oblique detonation wave created from a Mach 9.3 stoichiometric H2-air 

mixture flowing over a 27 degree wedge. The geometry and computational domain are 

shown in Figure 10. 

Detonation wave

Non-reactive shock Deflagration waves

Slip line

Computational domain

Wedge

Flow

X

Y

Figure 10 – Oblique Detonation Geometry and Computational Domain 

As shown, the computational domain was tilted 27 degrees to coincide with the 

wedge surface. A 200 by 150 uniform grid was used to model the 5 mm by 2.5 mm 

region. 0.5 mm of open boundary was included in the lower boundary before the leading 

edge of the wedge. A slip-wall, adiabatic boundary condition was used along the wedge 

wall. 

The flow entering the domain was at the following conditions, 

Pressure = 1 atmosphere Temperature = 300 K 

X-direction velocity = 3430 m/s Y-direction velocity = -1748 m/s 

H2 mass fraction = 0.0283 O2 mass fraction = 0.2265 

N2 mass fraction = 0.7452 

 The baseline calculation of Thaker and Chelliah (1997) was first repeated using 

the 12-step H2-air mechanism used in their study. This mechanism was taken from the 

full 18-step mechanism given in the Appendix by removing the H2O2 reactions. The 

resulting temperature and H2O mass fraction results are shown in Figure 11 and 

Figure 12. These results are consistent with those of Thaker and Chelliah (1997). 
NASA/TM—2002-211712 27 

X mm

01234

Y mm

0

1

2

Figure 11 – 12-Step Mechanism Oblique Detonation Temperature Contours 

X mm

01234

Y mm

0

1

2

Figure 12 – 12-Step Mechanism Oblique Detonation H2O Mass Fraction Contours 

 Next, the calculation was repeated using the 7-step mechanism chosen for this 

study. As shown in Figures 13 and 14, the initial shock/reaction zone is very similar to 

that obtained with the 12-step mechanism. However, it is apparent that the final stages of 

heat release are somewhat slower with the 7-step mechanism. This is consistent with the 

1-D results to be shown in the next section. It is not expected that this slight delay in heat 

release will influence the comparative results of this study. 

0.199 

0.161 

0.180 

Y (mm) 

X (mm) 

3310 

2710 

3010 

X (mm) 

Y (mm) 
NASA/TM—2002-211712 28 

X mm

01234

Y mm

0

1

2

Figure 13 – 7-Step Mechanism Oblique Detonation Temperature Contours 

X mm

01234

Y mm

0

1

2

Figure 14 – 7-Step Mechanism Oblique Detonation H2O Mass Fraction Contours 

3310 3010 2710

0.161 

0.180 

0.199 

X (mm) 

X (mm) 

Y (mm) 

Y (mm) 
NASA/TM—2002-211712 29 

4.3 Uniform H2-Air Mixture ZND Detonation 

As a final code validation case, the detonation structure computed by SPARK in a 

uniform stoichiometric H2-air mixture initially at 1 atmosphere pressure and 300 K (Case 

1 of this study, detonation tube centerline results) was compared to that obtained from the 

1-dimensional ZND code developed by Shepherd (1986). To run the ZND code, the 

desired kinetic mechanism is first input to the CHEMKIN code, described by Kee et al. 

(1989), that provides the kinetic and thermodynamic properties required. Next, the 

Gordon and McBride (1994) CEA code, or any other equilibrium detonation code, is used 

to determine the detonation (CJ) velocity. With these inputs, ZND calculates the 1-D 

detonation structure by first assuming a normal shock based on the detonation velocity 

and then reacting the mixture via the kinetic mechanism and gas properties from 

CHEMKIN using the post-shock conditions as the initial conditions. ZND performs 

these calculations in the reference frame fixed to the detonation wave, shown in Figure 

15, so the velocity must be transformed back to the stationary reference frame by adding 

the detonation wave velocity. All of the static properties are unaffected by the change in 

reference frames. 

Figure 15 – Stationary and Detonation Wave Reference Frames 

Figures 16 through 22 show the values of pressure, temperature, velocity, and 

H2O, OH, O, and H mass fractions calculated by SPARK against those calculated by 

ZND using both the 7-step kinetic mechanism and the full 18-step mechanism. The 18-

step ZND results are shown to additionally quantify the accuracy of the 7-step 

approximation beyond what is shown in the Appendix. The results from ZND begin after 

the leading shock. The results shown for SPARK are taken from the baseline case, Case 

1, which will be discussed further in Section 5. As shown, the ZND and SPARK results 

agree very well. The slight delay in completing the heat release that was discussed in the 

previous section is again evident in the temperature and H2O mass fraction plots. This is 

a very minor concern. Based on these results and those preceding, it appears that SPARK 

with the 7-step H2-air mechanism is adequate for the desired detonation calculations. 

Unburned Gas 

 V1=0 1 

Stationary 

Detonation 

Front 

1 2 

Burned Gas 

 V2

Vcj

Detonation 

Front 

Stationary Reference Frame 

(SPARK) 

Burned Gas 

 V2 

Reference Frame Fixed to Detonation Wave 

(ZND) 

Unburned Gas 

V1=Vcj

2 
NASA/TM—2002-211712 30 

Distance (cm)

0 1

Pressure (atm)

5

10

15

20

25

30

ZND 7 step 

ZND 18 step 

SPARK 7 step 

Figure 16 – Comparison of Pressure Profiles Between SPARK and ZND Codes 

Distance (cm)

0 1

Tem

perature (

K)

500

1000

1500

2000

2500

3000

3500

ZND 7 step 

ZND 18 step 

SPARK 7 step 

Figure 17 – Comparison of Temperature Profiles Between SPARK and ZND Codes 
NASA/TM—2002-211712 31 

Distance (cm)

0 1

Velocity (m/s)

200

400

600

800

1000

1200

1400

1600

1800

ZND 7 step 

ZND 18 step 

SPARK 7 step 

Figure 18 – Comparison of Velocity Profiles Between SPARK and ZND Codes 

Distance (cm)

0 1

H O Mass Fraction 2

0.00

0.05

0.10

0.15

0.20

0.25

0.30

ZND 7 step 

ZND 18 step 

SPARK 7 step 

Figure 19 – Comparison of H2O Mass Fraction Profiles Between SPARK 

and ZND Codes 

H O Mass Fraction 2
NASA/TM—2002-211712 32 

Distance (cm)

0 1

OH Mass Fraction

0.000

0.005

0.010

0.015

0.020

0.025

0.030

ZND 7 step 

ZND 18 step 

SPARK 7 step 

Figure 20 – Comparison of Detonation OH Mass Fraction Profiles Between SPARK 

and ZND Codes 

Distance (cm)

0 1

O Mass Fraction

0.000

0.005

0.010

0.015

0.020

0.025

ZND 7 step 

ZND 18 step 

SPARK 7 step 

Figure 21 – Comparison of Detonation O Mass Fraction Profiles Between SPARK 

and ZND Codes 
NASA/TM—2002-211712 33 

Distance (cm)

0 1

H Mass Fraction

0.000

0.001

0.002

0.003

0.004

0.005

ZND 7 step 

ZND 18 step 

SPARK 7 step 

Figure 22 – Comparison of Detonation H Mass Fraction Profiles Between SPARK 

and ZND Codes 

5. Computational Results 

Having shown the general applicability of the SPARK code with the 7-step H2/air 

mechanism for detonation calculations, we will now look at the results of each of the 11 

cases outlined in Section 2.5. Figure 4 is repeated below as Figure 23 as a reminder of 

the geometry and detonation initiation/stabilization approach used throughout the study. 

Figure 23 – Geometry for Numerical Study 

 H = 3 cm 

Stabilization Region ( 0 φ = 1. ) 

L = 30 cm 

Ignition Region (0.6 mm) 

Test Section – Tailored Mixtures 

 5 cm 
NASA/TM—2002-211712 34 

5.1 Introduction 

Before moving on to compare the results of the subject test cases, there are a 

number of introductory questions that need to be addressed. These questions deal with 

the stability of both the numerical solution and the detonation itself, the appropriate 

impulse calculation to use in the Isp calculation, and the effect of equivalence ratio on Isp

for uniform mixtures. It will also be useful to look at simplifications that can be made in 

the analysis of the 2-dimensional numerical results. These issues can all be addressed by 

looking at elements of the solutions that will be more fully presented in later sections. 

5.1.1 Dimensionality of Results 

It is appropriate to focus on the last issue first, namely, any simplifications that 

can be made in analyzing the 2-dimensional results obtained from SPARK. Detonation 

tubes are often simulated using 1-dimensional codes, as demonstrated by Kailasanath and 

Patniak (1999), Wilson and Paxson (2001), and Sterling et al. (1995). If the assumption 

of primarily 1-dimensional behavior for uniform and axially varying mixtures can be 

verified, then those cases can be visualized and compared by looking at 1-dimensional 

slices along the length of the tube, typically taken along the centerline. It is obvious that 

the transversely varying mixtures will still need to be examined to some extent using 2-

dimensional methods, but even these may be looked at 1-dimensionally for some flow 

variables. 

Figures 24 to 28 show 2-dimensional plots of pressure, transverse velocity and 

temperature from Cases 1, 2, 3, 4, and 7 as examples of a completely uniform mixture, a 

transversely varying mixture, an axially varying mixture, a uniform mixture with a step 

change in equivalence ratio along the tube, and a partially filled tube, respectively. All 

the plots in Figures 24 to 28 represent the flowfield just prior to the detonation wave 

exiting the tube. While all the cases show small traces of transverse velocity, only Case 

2 shows any signs of significant multi-dimensionality, seen in the temperature plot for 

that case. However, the pressure plot for Case 2 is still very 1-dimensional, which is to 

be expected, since there can be no net flow in the transverse direction. Thus for all the 

cases, one-dimensional comparisons of pressure should be valid. 

It is of interest that there is no observed boundary layer in any of the plots. A 

boundary layer is present in all cases, but it is so thin as to be unobservable in the scale of 

the plots, and it has essentially no effect on the flow features of the cases under study. It 

should also be noted that in Case 7, a partial fill case, the detonation wave has propagated 

into the air only region at the end of the tube, leading to the region of lower temperature 

at the end of the tube seen in that plot. 
NASA/TM—2002-211712 35 

Figure 24 – Case 1 Pressure, Transverse Velocity, and Temperature Contours 

Figure 25 – Case 2 Pressure, Transverse Velocity, and Temperature Contours 

Temperature (K)

Pressure (Pa)

Transverse Velocity (m/s) 

Temperature (K)

Pressure (Pa)

Transverse Velocity (m/s) 
NASA/TM—2002-211712 36 

Figure 26 – Case 3 Pressure, Transverse Velocity, and Temperature Contours 

Figure 27 – Case 4 Pressure, Transverse Velocity, and Temperature Contours 

Pressure (Pa)

Transverse Velocity (m/s) 

Temperature (K)

Pressure (Pa)

Transverse Velocity (m/s) 

Temperature (K)
NASA/TM—2002-211712 37 

Figure 28 – Case 7 Pressure, Transverse Velocity, and Temperature Contours 

5.1.2 Numerical and Detonation Stability 

 It is of great importance to establish that each simulation is both numerically and 

physically stable. This is most easily accomplished by looking at the detonation wave 

speed. The wave speed was determined by tracking the wave position, defined as the 

tube centerline position with the steepest X-velocity gradient, as a function of time during 

the wave propagation period of the calculation. The wave speed was found to be very 

sensitive to both concerns during the validation stages of this study. The equilibrium 

value of the Chapman-Jouget wave speed is easily calculable using the CEA code of 

Gordon and Mcbride (1994). In Figure 29, the three baseline case wave speeds are 

compared to each other, from which it can be seen that each simulation is stable, and that 

wave speed is proportional to equivalence ratio. Figures 30 through 33 show the wave 

speed for each baseline case compared to several reference detonation wave speeds 

calculated using the CEA code. After being initially overdriven to insure detonation 

initiation, the wave relaxes to just under the CJ wave speed in each case, and then slowly 

converges on the equilibrium CJ wave speed throughout the “steady state” propagation. 

The steady state wave speed is generally within 2% to 3% of the equilibrium value, also 

indicative of a stable detonation. 

Pressure (Pa)

Transverse Velocity (m/s) 

Temperature (K)
NASA/TM—2002-211712 38 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

2000

2200

2400 Wave Speed - Case 1 - Test Section Phi=1.0

Wave Speed - Case 4 - Test Section Phi=0.9

Wave Speed - Case 8 - Test Section Phi=1.1

Figure 29 – Detonation Wave Speeds for Cases 1, 4, and 8 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 1 - Uniform Mixture

CJ Wave Speed - Phi=1 

CJ Wave Speed - Phi=0.9 

CJ Wave Speed - Phi=1.1 

Figure 30 – Case 1 (Φ = 1.0 ) Wave Speed Compared to Equilibrium CJ Wave Speeds 
NASA/TM—2002-211712 39 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 4 - Uniform Mixture

CJ Wave Speed - Phi=1.0 

CJ Wave Speed - Phi=0.9 

CJ Wave Speed - Phi=0.8 

Figure 31 – Case 4 (Φ = 0.9 ) Wave Speed Compared to Equilibrium CJ Wave Speeds 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 8 - Uniform Mixture

CJ Wave Speed - Phi=1.2 

CJ Wave Speed - Phi=1.1 

CJ Wave Speed - Phi=1.0 

Figure 32 – Case 8 (Φ = 1.1) Wave Speed Compared to Equilibrium CJ Wave Speeds 
NASA/TM—2002-211712 40 

5.1.3 Total Impulse Calculation 

 As discussed in Section 2.4, there are two methods for calculating the impulse 

generated by the detonation and subsequent blowdown. The first is to integrate the 

pressure acting on the head-end of the tube over the cycle time, subtracting off the 

viscous drag of the walls, as given in equation 17. It is also acceptable to integrate the 

momentum flux exiting the tube over the cycle time, correcting for the incomplete 

expansion to ambient pressure, as given in equation 18. Figure 33 shows the results of 

both impulse calculations versus time for Case 1. It is not possible to differentiate the 

impulse generated by the initiation and stabilization zones, so calculated total impulse 

includes the effects of these regions. The weight of H2 in the stabilization zone is 

included in the total fuel weight in order to calculate the overall specific impulse. 

0

1

2

3

4

5

6

7

8

9

10

0 0.0002 0.0004 0.0006 0.0008 0.001

Time (s)

Total Impulse (N-s)

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

Shear Impulse (N-s)

Total Impulse (Head End Pressure)

Total Impulse (Exit Momentum Flux)

Total Shear Impulse

Figure 33 – Impulse Calculations for Case 1 – Uniform Stoichiometric Mixture 

 From Figure 33 it is seen that the head end pressure calculation has come close to 

its maximum value, while the exit momentum flux calculation is still rising. The two 

impulse values are about 6% apart at the termination of the calculation due to flow 

reversal at the tube exit boundary. This result is understandable in that flow is still 

exiting the tube at a considerable velocity at this point (see Figure 5). However, the 

pressure at the head-end has returned close to the initial pressure. The head-end pressure 

will in fact go below the ambient pressure if the calculation is allowed to continue due to 

reflected expansion waves in the tube. It is expected that if the calculation were run out 

far enough so that the tube returned to a quiescent state, these two impulse values would Total Shear Impulse (N-s) 
NASA/TM—2002-211712 41 

equalize. However, the point at which the head-end pressure begins to go subatmospheric is precisely the point at which it would be expected that the inlet air valves 

would open, drawing in fresh fuel and air and preventing the head pressure from going 

significantly negative. Since the thrust wall pressure integration has approached a 

maximum, it is evident that this is the value to use in the Isp calculation, being the closest 

to the expected maximum value. All further references to Isp values are based on the 

head-end pressure integration method of equation 17. 

 Before leaving Figure 33, it is interesting to note that the wall shear force impulse 

penalty also included with the figure shows that the wall drag only represents a penalty of 

approximately 0.75% of the total impulse. While the viscous drag is included in the total 

impulse calculation, it could be readily ignored without significantly altering the results 

of any of the study cases. 

5.1.4 Comparison of Baseline Uniformly Fueled Cases 

 It will be useful in understanding the effect of equivalence ratio gradients if the 

general effect of equivalence ratio in uniformly fueled mixtures in the range of interest is 

known. Figure 34 shows that Isp is inversely proportional to equivalence ratio, with the 

lean case achieving the highest Isp (4706 s), the stoichiometric case in the middle (4450 

s), and the fuel rich case yielding the lowest Isp (4177 s). If these maximum Isp values 

are plotted against overall equivalence ratio, as in Figure 35, it is seen that the Isp

dependence on equivalence ratio in this range is approximately linear. Inclusion of the 

stabilization zone raises the overall equivalence ratio of the fuel lean cases to 0.92, and 

lowers the overall equivalence ratio of the fuel rich cases to 1.08. Since the effect of the 

stabilization zone cannot be isolated within the calculation, this equivalence ratio bias is 

included in all Isp calculations. 

 Since Isp is calculated from the force generated by the fluid pressure pushing 

against the closed end of the tube, the head-end pressure can be plotted against time to 

get more insight into the source of the variation in Isp with equivalence ratio, as in 

Figure 36. It is evident that decreasing the equivalence ratio below an equivalence ratio 

of one decreases the plateau pressure achieved, but not in proportion to the decrease in 

the amount of fuel used, thus causing an increase in Isp for Case 4. However, no 

significant increase in plateau pressure was achieved by over-fueling the tube as in Case 

8, causing a decrease in Isp due to the increase in fuel used. It is of interest to note that 

the duration of the pressure plateau and the slope of the blowdown did not change 

significantly as a function of equivalence ratio. 
NASA/TM—2002-211712 42 

Time (s)

0.0000 0.0002 0.0004 0.0006 0.0008 0.0010

I

sp (s)

0

1000

2000

3000

4000

5000

I

sp

 - Case 1 - Phi=1.0

I

sp

 - Case 4 - Phi=0.9

I

sp

 - Case 8 - Phi=1.1

Figure 34 – Isp for Baseline Uniformly Fueled Cases 1, 4, and 8 

4100

4200

4300

4400

4500

4600

4700

4800

0.9 0.95 1 1.05 1.1

Overall Equivalence Ratio

Isp (s)

Figure 35 – Isp as a Function of Overall Equivalence Ratio, including Stabilization Zone, 

for Cases 1, 4, and 8 
NASA/TM—2002-211712 43 

Time (s)

0.0000 0.0002 0.0004 0.0006 0.0008 0.0010

Head End Pressure (atm)

0

2

4

6

8

Head End Pressure - Case 1 - Phi=1.0

Head End Pressure - Case 4 - Phi=0.9

Head End Pressure - Case 8 - Phi=1.1

Figure 36 – Detonation Tube Head-End Pressure for Cases 1, 4, and 8 

 Having dealt with these introductory issues, we can now move on to the analysis 

of the performance results for each group of mixtures. 

5.2 Stoichiometric Mixtures 

 Before looking at the performance results, it is of value to verify the numerical 

and physical stability of the detonations. Figures 37 through 40 show the wave speed 

results for the first 3 cases, first compared to each other and then to several reference 

equilibrium values calculated using the CEA code. Each simulation is clearly stable. It 

is of interest to note that the transverse fuel distribution, Case 2, reaches the same wave 

speed as the uniformly fueled mixture, Case 1, after initial regions of overshoot and 

undershoot, and that this wave speed is stable. This indicates that the detonation wave is 

underdriven in the rich mixture and overdriven in the lean mixture. We will come back 

to this point after looking at the performance results. The wave speed of the axially 

varying fuel distribution case follows the initial equivalence ratio distribution fairly well, 

deviating from the equilibrium value slightly more at the step change in equivalence ratio 

located at the 5 cm point (~2.2×10–5 s) as there is not enough time for the detonation 

wave to react to the change in mixture at that point. 
NASA/TM—2002-211712 44 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 1 - Uniform Mixture

Wave Speed - Case 2 - Transverse Gradient

Wave Speed - Case 3 - Axial Gradient

Figure 37 – Detonation Wave Speeds for Stoichiometric Cases 1, 2, and 3 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 1 - Uniform Mixture

CJ Wave Speed - Phi=1 

CJ Wave Speed - Phi=0.9 

CJ Wave Speed - Phi=1.1 

Figure 38 – Case 1 (Uniform Mixture) Wave Speed Compared to 

Equilibrium CJ Wave Speeds 
NASA/TM—2002-211712 45 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 2 - Transverse Gradient

CJ Wave Speed - Phi=1 

CJ Wave Speed - Phi=0.9 

CJ Wave Speed - Phi=1.1 

Figure 39 – Case 2 (Transverse Gradient) Wave Speed Compared to 

Equilibrium CJ Wave Speeds 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 3 - Axial Gradient 

CJ Wave Speed - Phi=1 

CJ Wave Speed - Phi=0.9 

CJ Wave Speed - Phi=1.1 

Figure 40 – Case 3 (Axial Gradient) Wave Speed Compared to 

Equilibrium CJ Wave Speeds 
NASA/TM—2002-211712 46 

Figure 41 shows the specific impulse traces versus time for each of the 

stoichiometrically fueled cases. In the scale of Figure 41, there is essentially no 

difference in Isp across the three simulations. Figure 42 provides a closer view of the 

final section of the Isp curves from Figure 41. We can see that there is some difference in 

the final values of Isp, but that the results are all within 1% of each other, with the axial 

distribution slightly outperforming the other 2 distributions. 

Figure 43 shows that the head-end pressure is essentially unchanged for the three 

stoichiometric cases. Case 3 shows a slight increase in the head end pressure during the 

plateau region, consistent with the slightly improved performance shown in Figure 42. 

The plateau duration and blowdown contour are virtually identical for all three cases. 

Time (s)

0 2e-4 4e-4 6e-4 8e-4

Isp (s)

0

1000

2000

3000

4000

5000

I

sp

 - Case 1 - Uniform Mixture

I

sp

 - Case 2 - Transverse Gradient

I

sp

 - Case 3 - Axial Gradient

Figure 41 – Specific Impulse of Stoichiometric Cases 1, 2, and 3 
NASA/TM—2002-211712 47 

Time (s)

7e-4 8e-4 8e-4 9e-4 9e-4

Isp (s)

4000

4100

4200

4300

4400

4500

I

sp

 - Case 1 - Uniform Mixture

I

sp

 - Case 2 - Transverse Gradient

I

sp

 - Case 3 - Axial Gradient

Figure 42 – Enlarged View of Final Specific Impulse of Stoichiometric 

Cases 1, 2, and 3 

Time (s)

0 2e-4 4e-4 6e-4 8e-4 1e-3

He

a

d-En

d Pre

s

sure (atm) 0

2

4

6

8

Head-End Pressure - Case 1 

Uniform Mixture

Head-End Pressure - Case 2 

Transverse Gradient

Head-End Pressure - Case 3 

Axial Gradient

Figure 43 – Head-End Pressure as a Function of Time for Stoichiometric 

Cases 1, 2, and 3 
NASA/TM—2002-211712 48 

 Lastly, we can look at plots of H2 or H2O mass fraction to determine if any 

significant mixing or post-detonation combustion is occurring in the detonation tube. 

Figure 44 shows the hydrogen mass fraction in the tube at five different times during the 

simulation. The first three plots are during the detonation propagation, and the last two 

plots are during the blowdown. The post-detonation H2 distribution remains essentially 

stratified throughout the entire propagation and blowdown, with the composition 

distribution at the left of the tube stretching to the right to fill the tube as the blowdown 

progresses. The products of the uniformly fueled stoichiometric stabilization zone 

eventually push out most of the stratified combustion products by the end of the 

simulation. The relatively high levels of residual H2 at the left-hand (closed) end of the 

tube seen in this time sequence are a result of the thermal breakdown of H2O from in the 

high pressure, high temperature region used to initiate the detonation and are not a 

product of the detonation process itself. It can be assumed that negligible mixing occurs 

in the axial direction if little mixing is observed in the transverse direction, since the axial 

composition gradients are approximately 1/16th as steep as the transverse composition 

gradients. This is born out by the one-dimensional time history of H2O mass fraction for 

Case 3 shown in Figure 45. 

Figure 44 – Two-Dimensional H2 Mass Fraction Distributions for Case 2 – 

Transverse Gradient 

t=0.88 ms 

t=0.47 ms 

t=0.15 ms 

t=0.08 ms 

t=0.03 ms 
NASA/TM—2002-211712 49 

X (cm)

0 5 10 15 20 25 30

H O Mass Fraction 2

0.00

0.05

0.10

0.15

0.20

0.25

0.30

time = 0.03 ms

time = 0.08 ms

time = 0.15 ms

time = 0.51 ms

time = 0.89 ms

Figure 45 – One-Dimensional H2O Mass Fraction Distributions for Case 3 – 

Axial Gradient 

5.3 Fuel Lean Mixtures 

 The wave speed results for the fuel lean cases, Figures 46 through 50, are much the 

same as the stoichiometric cases. The uniformly fueled mixture, Case 4, and the 

transversely varying mixture, Case 5, follow a similar profile of wave speed. The axially 

varying mixture, Case 6, follows the initial mixture profile at near the equilibrium 

detonation speed. The added case, a uniform stoichiometric mixture with an air buffer at 

the end of the tube, shows the expected drop off in wave speed as the detonation fails in 

the unfueled air. The wave speed falls off precipitously as the detonation wave enters the 

unfueled region, limiting the transfer of energy to the buffer air region. 
NASA/TM—2002-211712 50 

 Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1200

1400

1600

1800

2000

2200

2400

Wave Speed - Case 4 - Uniform Mixture

Wave Speed - Case 5 - Transverse Gradient

Wave Speed - Case 6 - Axial Gradient

Wave Speed - Case 7 - Partial Fill (Air)

Figure 46 – Detonation Wave Speeds of Fuel Lean Cases 4, 5, 6, and 7 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 4 - Uniform Mixture

CJ Wave Speed - Phi=1.0 

CJ Wave Speed - Phi=0.9 

CJ Wave Speed - Phi=0.8 

Figure 47 – Case 4 (Uniform Mixture) Wave Speed Compared to Equilibrium 

CJ Wave Speeds 
NASA/TM—2002-211712 51 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (

m/s) 1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 5 - Transverse Gradient

CJ Wave Speed - Phi=1.0 

CJ Wave Speed - Phi=0.9 

CJ Wave Speed - Phi=0.8 

Figure 48 – Case 5 (Transverse Gradient) Wave Speed Compared to 

Equilibrium CJ Wave Speeds 

 Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (

m/s) 1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 6 - Axial Gradient

CJ Wave Speed - Phi=1.0 

CJ Wave Speed - Phi=0.9 

CJ Wave Speed - Phi=0.8 

Figure 49 – Case 6 (Axial Gradient) Wave Speed Compared to 

Equilibrium CJ Wave Speeds 
NASA/TM—2002-211712 52 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1200

1400

1600

1800

2000

2200

2400

Wave Speed - Case7 - Partial Fill (Air)

CJ Wave Speed - Phi=1.0 

CJ Wave Speed - Phi=0.9 

CJ Wave Speed - Phi=0.8 

Figure 50 – Case 7 (Partial Fill (Air)) Wave Speed Compared to Equilibrium 

CJ Wave Speeds 

 The specific impulse results, Figure 51, are also similar to the stoichiometric results, 

allowing for the difference in magnitude due to the change in overall equivalence ratio as 

documented in Section 5.1. The total variation in specific impulse is on the order of 1%. 

We again see the axial distribution slightly outperforming the uniform and transverse 

distributions. The new “partial fill” case slightly under-performs the others. This 

indicates that it would be slightly better to throttle a PDE by decreasing the equivalence 

ratio through reduced fuel flow rather than by partial filling, if that is a practical approach 

when other factors are considered. 

 The head end pressure time histories in Figure 52 show more case-to-case 

variation than did the stoichiometric mixtures. The axially varying mixture has a 

continually decreasing plateau pressure, and the partial fill case has a shortened pressure 

plateau. While there are differences, it is evident from the specific impulse results that 

the positive and negative aspects of these differences work to cancel each other out. 
NASA/TM—2002-211712 53 

Time (s)

0.00075 0.00080 0.00085 0.00090 0.00095

Isp (s)

4400

4500

4600

4700

I

sp

 - Case 4 - Uniform Mixture

I

sp

 - Case 5 - Transverse Gradient

I

sp

 - Case 6 - Axial Gradient

I

sp - Case 7 - Partial Fill (Air)

Figure 51 – Enlarged View of Final Specific Impulse of Fuel Lean 

Cases 4, 5, 6, and 7 

 

Time (s)

0.0000 0.0002 0.0004 0.0006 0.0008 0.0010

Head End Pressure (atm)

0

2

4

6

8

10

Head End Pressure - Case 4 - Uniform Mixture

Head End Pressure - Case 5 - Transverse Gradient

Head End Pressure - Case 6 - Axial Gradient

Head End Pressure - Case 7 - Partial Fill (Air)

Figure 52 – Head-End Pressure as a Function of Time for Fuel Lean Cases 4, 5, 6, and 7 
NASA/TM—2002-211712 54 

 Again checking for mixing in the detonation products for the transversely varying 

mixture, Case 5, Figure 53 shows a time series of two-dimensional plots of O2 mass 

fraction. Since there are no fuel rich zones in this distribution, it is expected that there 

will be no more than trace amounts of H2 present in the combustion products, so O2

serves as a better indicator of mixing and/or combustion within the tube. The same 

characteristic behavior is seen here as was seen for the stoichiometric case. The stratified 

post-detonation composition distribution remains stratified as it is pushed out of the 

detonation tube during the blowdown process. No significant mixing or diffusion leading 

to post-detonation combustion is observed. 

Figure 53 – Two-Dimensional O2 Mass Fraction Distributions for Case 5 – 

Transverse Gradient 

5.4 Fuel Rich Mixtures 

The wave speed results, Figures 54 through 58, for the fuel rich cases also show 

similar behavior to those for the previous two overall equivalence ratios. The uniformly 

fueled mixture wave speed, Case 8, and the transversely varying mixture wave speed, 

t=0.08 ms 

t=0.16 ms 

t=0.48 ms 

t=0.90 ms 

t=0.03 ms 
NASA/TM—2002-211712 55 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 8 - Uniform Mixture

Wave Speed - Case 9 - Transverse Gradient

Wave Speed - Case 10 - Axial Gradient

Wave Speed - Case 11 - Partial Fill (H2)

Figure 54 – Detonation Wave Speeds of Fuel Rich Cases 8, 9, 10, and 11 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 8 - Uniform Mixture

CJ Wave Speed - Phi=1.2 

CJ Wave Speed - Phi=1.1 

CJ Wave Speed - Phi=1.0 

Figure 55 – Case 8 (Uniform Mixture) Wave Speed Compared to Equilibrium 

CJ Wave Speeds 
NASA/TM—2002-211712 56 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 9 - Transverse Gradient

CJ Wave Speed - Phi=1.2 

CJ Wave Speed - Phi=1.1 

CJ Wave Speed - Phi=1.0 

Figure 56 – Case 9 (Transverse Gradient) Wave Speed Compared to 

Equilibrium CJ Wave Speeds 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed-Case 10 - Axial Gradient

CJ Wave Speed - Phi=1.2 

CJ Wave Speed - Phi=1.1 

CJ Wave Speed - Phi=1.0 

Figure 57 – Case 10 (Axial Gradient) Wave Speed Compared to 

Equilibrium CJ Wave Speeds 
NASA/TM—2002-211712 57 

Time (s)

0.0 2.0e-5 4.0e-5 6.0e-5 8.0e-5 1.0e-4 1.2e-4 1.4e-4 1.6e-4

Wave Speed (m/s)

1800

1900

2000

2100

2200

2300

2400 Wave Speed - Case 11 - Partial Fill (H2)

CJ Wave Speed - Phi=1.2 

CJ Wave Speed - Phi=1.1 

CJ Wave Speed - Phi=1.0 

Figure 58 – Case 11 (Partial Fill (H2)) Wave Speed Compared to Equilibrium 

CJ Wave Speeds 

Case 9, do not agree quite as well as before, but are still relatively close. The axially 

varying mixture wave speed, Case 10, still follows the initial mixture distribution. The 

partial fill distribution, Case 11, this time with H2 at the end of the tube instead of air, 

shows a marked increase in wave speed as the detonation wave enters the end buffer 

region. This is due to the change in the local speed of sound that is significantly higher in 

hydrogen than in air. 

The specific impulse results (Figure 59) are, if anything, more tightly clustered 

than in either of the previous groups of mixtures. The uniformly fueled case and the 

axially varying case are almost identical, with the transversely varying case slightly 

lower. The partial fill case has gone from being the worst performer with an air buffer 

for the fuel lean cases, to being the best performer for the fuel rich cases. However, again 

the total variation in results is on the order of only 1%. 

Looking at the head end pressure time histories in Figure 60, we see almost no 

variation at all, despite the slightly increased variability in wave speed observed in 

Figure 54. 

Lastly, we see in Figure 61 the same lack of mixing in the post-detonation 

composition distribution of the 2-dimensional plots for the transverse mixture in Case 9. 

H2 mass fraction has again been used for these plots since there are fuel rich zones in the 

pre-combustion mixture.
NASA/TM—2002-211712 58 

Time (s)

0.00070 0.00075 0.00080 0.00085 0.00090

Isp (s)

3800

3900

4000

4100

4200

4300

I

sp

 - Case 8 - Uniform Mixture

I

sp

 - Case 9 - Transverse Gradient

I

sp

 - Case 10 - Axial Gradient

I

sp - Case 11 - Partial Fill (H2)

Figure 59 – Enlarged View of Final Specific Impulse of Fuel Rich 

Cases 8, 9, 10, and 11 

Time (s)

0.0000 0.0002 0.0004 0.0006 0.0008 0.0010

Head End Pressure (atm)

0

2

4

6

8

10

Head End Pressure - Case 8 - Uniform Mixture

Head End Pressure - Case 9 - Transverse Gradient

Head End Pressure - Case 10 - Axial Gradient

Head End Pressure - Case 11 - Partial Fill (H2)

Figure 60 – Head-End Pressure as a Function of Time for Fuel 

Rich Cases 8, 9, 10, and 11
NASA/TM—2002-211712 59 

Figure 61 – H2 Mass Fraction Distributions for Case 9 – Transverse Gradient 

6. Summary and Conclusions 

6.1 Summary of Results 

 A total of eleven 2-dimensional detonation tube simulations were conducted. Of 

these, three were run with stoichiometric H2/air mixtures, four were run with fuel lean 

H2/air mixtures, and four were run with fuel rich H2/air mixtures. An approximately 

linear relationship between Isp and equivalence ratio was observed between the uniformly 

fueled cases (1, 4 and 8). Very little variation (<1%) in Isp was observed between cases 

of the same overall stoichiometry. These results are summarized in Table 2. 

t=0.87 ms 

t=0.49 ms 

t=0.14 ms 

t=0.08 ms 

t=0.03 ms 
NASA/TM—2002-211712 60 

Table 2 – Summary of Test Case Performance Results 

Case 

Number 

Overall 

Stoichiometry 

Type of 

Distribution 

Isp

(s) 

1 Stoichiometric Uniform 4450 

2 Stoichiometric Transverse 4429 

3 Stoichiometric Axial 4459 

4 Fuel Lean Uniform 4706 

5 Fuel Lean Transverse 4699 

6 Fuel Lean Axial 4722 

7 Fuel Lean Partial Fill (Air) 4676 

8 Fuel Rich Uniform 4177 

9 Fuel Rich Transverse 4168 

10 Fuel Rich Axial 4175 

11 Fuel Rich Partial Fill (H2) 4204 

 Transversely and axially varying fuel/air mixtures remained basically stratified 

throughout the detonation propagation and blowdown process. 

6.2 Conclusions 

The primary conclusion to be drawn from this study is that a lack of fuel-air 

mixing up to a moderate level in a hydrogen fueled airbreathing detonation tube has 

almost no effect on the thrust performance of the system. Put another way, thrust 

performance is almost completely independent of combustion efficiency in the vicinity of 

an overall equivalence ratio of one. This is an encouraging result for the design of such 

systems, as it appears to be unnecessary to go to extremes to achieve good mixing in the 

fuel injection process. These results also imply that it is probably unnecessary to match 

air and fuel valve opening and closing profiles to maintain constant equivalence ratio. 

These results are limited to those cases where the entire fuel/air charge, while not fully 

mixed, is nonetheless fully detonable. 
NASA/TM—2002-211712 61 

There are a number of secondary conclusions that may also be drawn for H2/air 

systems. 

1. Decreased equivalence ratio and partial-filling of the detonation tube are 

essentially equivalent. If anything, it appears slightly more advantageous to lower 

the equivalence ratio as long as the mixture remains detonable. Obviously, a 

point is reached where it is impossible to lower the equivalence ratio and maintain 

a detonable mixture. At this point, further throttling can be achieved by partially 

filling the detonation tube with a lean mixture. 

2. Non-transversely varying mixtures may be treated as one-dimensional, ignoring 

wall shear losses. This means that future studies with such mixtures may be 

safely conducted using one-dimensional codes without fear of significant loss of 

accuracy. This will save considerable time in computational analysis. 

3. Transversely and axially varying mixtures remain essentially stratified throughout 

the detonation propagation and blowdown process. One benefit of this 

observation is that it will be possible to qualitatively evaluate mixing within a 

detonation tube by looking at the distribution of products coming out of the tube, 

thus reducing the need for costly optical test sections in detonation tubes. 

4. H2/air detonation tube systems can be adequately computed using simplified 

chemistry, thus reducing required computational time. Early test runs using the 

SPARK code with the full chemical mechanism found in Table A1 showed a 

greater than one order of magnitude increase in computational time required for a 

detonation tube calculation. This was due not only to the increased number of 

species, but also the finer axial grid spacing required to capture with precision the 

trace concentrations of HO2 and H2O2 in the reactions subsequently eliminated. 

6.3 Recommendations for Future Work 

There are a number of areas worthy of further inquiry that remain at the 

conclusion of this study. These include, 

Increasing the range of equivalence ratios studied. It would be of great value to study 

levels of fuel/air mixture variation that included marginally-detonable and non-detonable 

mixtures. Such mixtures would result from detonation tube film cooling schemes or poor 

mixing configurations. 

1. Performing calculations using a representative hydrocarbon fuel. Hydrocarbon 

fuels are likely candidates for most non-access-to-space applications due to their 

high density and ease of use. It is likely that such systems will behave 

significantly differently. Hydrogen is an excellent monopropellant because of its 

low molecular weight. Therefore, unburned hydrogen raised to a high 

temperature by the leading shock wave of the detonation front still provides 

significant impulse. A hydrocarbon fuel will not behave in the same manner and 

will likely penalize the system to a greater extent. 
NASA/TM—2002-211712 62 

2. Further development of PDEs is likely to bring forward a number of injection 

patterns not addressed by this study. As realistic patterns become available, 

further analysis should be performed. 

3. Simulations including turbulence should be performed to assess the impact of 

turbulence, both before and after the passage of the detonation wave, on the 

overall system performance. 
NASA/TM—2002-211712 63 

Appendix 

Justification of Truncated Kinetic Mechanism 

 The standard H2-air kinetic mechanism of Jachimowski (1988) is given in 

Table A1 below. For this study, reactions 8 through 18 have been eliminated. These 

reactions are those that contain HO2 and H2O2, species that generally govern the low 

temperature ignition process. However, ignition in a detonation wave occurs 

immediately after the moving shock wave, so the combustible mixture at ignition is not at 

a low temperature. For an equilibrium stoichiometric H2-air detonation wave with initial 

conditions of 1 atmosphere and 300 K, the Gordon and McBride (1994) CEA code 

predicts a detonation wave Mach number of 4.8. This corresponds to a temperature and 

pressure behind the shock of approximately 1650 K and 27 atmospheres, respectively, as 

calculated using the CEA code shock option. It is therefore in the neighborhood of this 

temperature and pressure that the necessity of incorporating the HO2 and H2O2 reactions 

should be considered. 

Table A1 – 18-Step H2-Air Mechanism 

Reaction 

Ref. # 

Reaction A n Ea 

(cal/g-mole) 

1 H2 + O2 ⇔ 2OH 0.1700e14 0.00 48150 

2 H + O ⇔ OH + O 2 0.1420e15 0.00 16400 

3 H2 + OH ⇔ H2O + H 0.3160e08 1.80 3030 

4 O + H2 ⇔ OH + H 0.2070e15 0.00 13750 

5 OH + OH ⇔ H2O + O 0.5500e14 0.00 7000 

6 OH + H + M ⇔ H2O + M 0.2210e23 –2.00 0 

7 H + H + M ⇔ H2 + M 0.6530e18 –1.00 0 

8 H + O2 + M ⇔ HO2 + M 0.3200e19 –1.00 0 

9 OH + HO2 ⇔ H2O + O2 0.5000e14 0.00 1000 

10 H + HO2 ⇔ H2 + O2 0.2530e14 0.00 700 

11 H + HO2 ⇔ OH + OH 0.1990e13 0.00 1800 

12 O + HO2 ⇔ O2 + OH 0.5000e12 0.00 1000 

13 HO2 + HO2 ⇔ O2 + H2O2 0.1990e13 0.00 0 

14 H2 + HO2 ⇔ H + H2O2 0.3010e12 0.00 18700 

15 OH + H2O2 ⇔ HO2 + H2O 0.1020e14 0.00 1900 

16 H + H2O2 ⇔ OH + H2O 0.5000e15 0.00 10000 

17 O + H2O2 ⇔ OH + HO2 0.1990e14 0.00 5900 

18 H2O2 + M ⇔ OH + OH + M 0.1210e18 0.00 45500 
NASA/TM—2002-211712 64 

 The chain-branching reactions in the H2/air system are reactions 2, 3, and 4 in 

Table A1. Kuo (1986) indicates that reaction 2 is the rate-controlling step in the chainbranching explosion, as it is the most endothermic. Looking at reactions 2-4 as a standalone system wherein H2 and O2 are consumed and H2O and H radicals are produced, we 

see that for each H radical consumed, 3 H radicals are produced, for a net of 2 H radicals. 

Thus the net production rate of H radicals in the chain-branching system is 2 times the 

forward rate constant of reaction 2, or 2kf2, multiplied by the concentrations of H and O2. 

The balancing reaction that prevents the mixture of H2 and O2 from running away into an 

explosion is reaction 8, wherein HO2 is produced. HO2 is relatively non-reactive, thus 

terminating the branching reaction chain. Since reaction 8 contains a 3rd body, it is 

dependent on the number density of molecules, and thus on pressure as well as 

temperature. Therefore, for a given temperature, the importance of the HO2 producing 

chain terminating reaction is dependant upon the pressure. If we eliminate HO2 from the 

chemical mechanism, we must determine whether or not we are significantly influencing 

the ignition process. In order to determine this, we need to calculate the temperature at 

which the H radical consumption reaction just balances the H radical production system 

of reactions. If this temperature is significantly below the post-shock temperature, then 

we can safely ignore HO2, and consequently H2O2 that is produced from the HO2. 

The equality that must be solved for the critical temperature is 

[ ][ ] [ ][ ][ ] MOHkOHk ff = 82 222 (43) 

where [ ] TR

p M

u

= , as defined by Kuo (1986), the concentration of molecules available 

for collisions, i.e. molar density. 

As mentioned previously, each reaction rate constant is of the form 

TR

Ea

n

f u eTAk

−

∗∗= (28) 

So we can substitute the appropriate values from Table A1 into equation (43) and 

simplify. 





















 



 ×=×

−

TR

p

T

e

u

uTR 1 )10320.0(*)10284.0( 19

16400

15 (44) 

Solving for temperature at a pressure of 27 atmospheres, the pressure immediately behind 

the shock wave, we find that the creation and consumption of the H radical equilibrate at 

about 1300 K, 350 degrees K below the expected post-shock temperature. Thus we can 

conclude that the HO2 reactions should not be a significant factor for this study. 
NASA/TM—2002-211712 65 

References 

Anderson, Jr., J.D. Computational Fluid Dynamics, McGraw-Hill, Inc., New York, 1995. 

Cambier, J.-L.; Adelman, H.G. “Preliminary Numerical Simulations of a Pulsed 

Detonation Wave Engine,” AIAA Paper 88–2960, July 1988. 

Carpenter, M.H. “Three-Dimensional Computations of Cross-Flow Injection and 

Combustion in a Supersonic Flow,” AIAA 89–1870, June 1989. 

Coleman, M.L. “Overview of Pulse Detonation Propulsion Technology,” CPIA Report 

No. CPTR 70, April, 2001. 

Drummond J.P.; Mukunda, H.S. “A Numerical Study of Mixing Enhancement in 

Supersonic Reacting Flow Fields,” AIAA Paper 88–3260, July 1988. 

Drummond, J.P. “A Two-Dimensional Numerical Simulation of a Supersonic, 

Chemically Reacting Mixing Layer,” NASA TM–4055, 1988. 

Drummond, J.P. “Supersonic Reacting Internal Flow Fields,” NASA TM–103480, 1990. 

Drummond, J.P.; Carpenter, M.H.; White, J.A. “Development and Application of the 

SPARK Combustion Codes,” 10th National Aero-Space Plane Technology Symposium, 

April 1991. 

Drummond, J.P. “Supersonic Reacting Internal Flowfields,” in Numerical Approaches to 

Combustion Modeling, E.S. Oran and J.P. Boris, Eds., AIAA, Washington, D.C., 1991. 

Gordon S.; McBride, B.J. “Computer Program for Calculation of Complex Chemical 

Equilibrium Compositions and Applications,” NASA Reference Publication 1311, 1994. 

Jachimowski, C.J. “An Analytical Study of the Hydrogen-Air Reaction Mechanism with 

Application to Scramjet Combustion,” NASA Technical Paper 2791, February 1988. 

Kaemming, T.A. “Integrated Vehicle Comparison of Turbo-Ramjet Engine and Pulsed 

Detonation Engine (PDE),” ASME Paper 2001–GT–0451, June 2001. 

Kailasanath, K. “Applications of Detonations to Propulsion: A Review,” AIAA Paper 

99–1067, January, 1999. 

Kailasanath, K.; Patnaik, G. “Pulsed Detonation Engines – What is its Performance,” 

Proceedings, JANNAF Combustion Meeting, November 1999. 
NASA/TM—2002-211712 66 

Kailasanath, K.; Patnaik, G.; Li, C. “On Factors Controlling the Performance of Pulsed 

Detonation Engines”, in Control of Detonation Processes, G. Roy, et al, editors. ElexKM, Moscow, 2000. 

Kee, R.J.; Rupley, F.M.; Miller, J.A. “Chemkin-II: A FORTRAN Chemical Kinetics 

Package for the Analysis of Gas-Phase Chemical Kinetics,” Sandia National Laboratories 

Report SAND89–8009, 1989. 

Kuo, K.K. Principles of Combustion, John Wiley & Sons, New York, 1986. 

Lynch, E.D.; Edelman, R. “Computational Fluid Dynamic Analysis of the Pulse 

Detonation Engine Concept,” AIAA Paper 94–0264, January 1994. 

Nichols, J.; Wilkinson, H.; Morrison, R. “Intermittent Detonation as a Thrust Producing 

Mechanism,” Jet Propulsion Journal of the American Rocket Society 27 [5] (1957), 

534–541. 

Oates, G.C. Aerothermodynamics of Gas Turbine and Rocket Propulsion, AIAA, 

Washington D.C., 1984. 

Poinsot, T.J.; Lele, S.K. “Boundary Conditions for Direct Simulations of Compressible 

Viscous Flows,” Journal of Computational Physics, 101, (1992), 104–129. 

Radhakrishnan, K. “Combustion Kinetics and Sensitivity Analysis,” in Numerical 

Approaches to Combustion Modeling, E.S. Oran and J.P. Boris, Eds., AIAA, Washington 

D.C., 1991. 

Shapiro, A.H. The Dynamics and Thermodynamics of Compressible Fluid Flow, Ronald 

Press Company, New York, 1953. 

Shepherd, J.E. “Chemical Kinetics of Hydrogen-Air-Diluent Detonation,” in Progress in 

Astronautics and Aeronautics, Vol. 106, p. 263–293, AIAA, New York, 1986. 

Smith R.E.; Weigel, B.L. “Analytic and Approximate Boundary-Fitted Coordinate 

Systems for Fluid Flow Simulation,” AIAA Paper 80–0192, January 1980. 

Sterling, J.; Ghorbanian, K.; Humphrey, J.; Sobota, T. “Numerical Investigations of Pulse 

Detonation Wave Engines,” AIAA Paper 95–2479, July 1995. 

Thaker A.A.; Chelliah, H.K. “Numerical Prediction of Oblique Detonation Wave 

Structures Using Detailed and Reduced Reaction Mechanisms,” Combustion Theory 

Modelling, [1], (1997) 347–376. 
NASA/TM—2002-211712 67 

Williams, F.A. Combustion Theory, Addison-Wesley Publishing Company, Inc., 

Menlo Park, 1985. 

Wilson, J.; Paxson, D.E. “On the Exit Boundary Condition for One-dimensional 

Calculations of Pulsed Detonation Engine Performance,” Proceedings, International 

Colloquium on Detonation, Explosion, and Reactive Systems, July 2001. 
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

Technical Memorandum

Unclassified

National Aeronautics and Space Administration

John H. Glenn Research Center at Lewis Field

Cleveland, Ohio 44135–3191

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

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

August 2002

NASA TM—2002-211712

E–13463

WU–708–48–13–00

74

Effects of Fuel Distribution on Detonation Tube Performance

Hugh Douglas Perkins

Pulse detonation engines; Detonation; Combustion; Performance; Mixing; Computational

fluid dynamics

Unclassified -Unlimited

Subject Category: 07 Distribution: Nonstandard

Responsible person, Hugh Douglas Perkins, organization code 0142, 216–977–7414.

A pulse detonation engine (PDE) uses a series of high frequency intermittent detonation tubes to generate thrust. The

process of filling the detonation tube with fuel and air for each cycle may yield non-uniform mixtures. Lack of mixture

uniformity is commonly ignored when calculating detonation tube thrust performance. In this study, detonation cycles

featuring idealized non-uniform H2/air mixtures were analyzed using the SPARK two-dimensional Navier-Stokes CFD

code with 7-step H2/air reaction mechanism. Mixture non-uniformities examined included axial equivalence ratio gradients, transverse equivalence ratio gradients, and partially fueled tubes. Three different average test section equivalence

ratios (Φ), stoichiometric (Φ = 1.00), fuel lean (Φ = 0.90), and fuel rich (Φ = 1.10), were studied. All mixtures were

detonable throughout the detonation tube. It was found that various mixtures representing the same test section equivalence ratio had specific impulses within 1 percent of each other, indicating that good fuel/air mixing is not a prerequisite

for optimal detonation tube performance.