Detonation_Reflection_ShockWaves

Analysis of Shock Wave Acceleration from Normal

Detonation Reflection

Donner T. Schoeffler

and Joseph E. Shepherd

California Institute of Technology, Pasadena, CA.

*Corresponding author(s). E-mail(s): dschoeff@caltech.edu;

Contributing authors: joseph.e.shepherd@caltech.edu;

Abstract

Normal detonation reflection generates a shock wave that exhibits complicated dynamics as it propa-

gates through the incident detonation and post-detonation flow. Ideal models have historically neglected

the influence of a finite detonation thickness on the reflected shock due to its small size relative

to laboratory scales. However, one-dimensional numerical simulations show that the reflected shock

accelerates to a large shock speed not predicted by ideal theory as it propagates through the inci-

dent detonation. Analysis with a derived shock-change equation identifies the principal role of the

highly nonuniform upstream flow on producing the large shock acceleration. Simulations of detonation

reflection show how a finite detonation thickness affects the entire trajectory of the reflected shock.

Keywords:

Detonation reflection, Shock acceleration, Shock-change equations, Heterogeneous fluids

1 Introduction

The propagation of gaseous detonation waves down

pipes is a common feature of laboratory detonation

physics experiments and possible industrial accident

scenarios. In nearly all cases, the detonation reaches

the end of a pipe and impacts a high-impedance

wall, where it reflects a shock wave that propa-

gates backwards through the incident detonation and

post-detonation flow. Normal detonation reflection is

a fundamental example of both detonation-structure

interactions and shock propagation through heteroge-

neous media. Despite the elementary nature of the

problem and its ubiquitous occurrence in most detona-

tion tube experiments, prior research has been sparse,

and basic questions relating to the reflected shock’s

evolution remain unanswered.

Nearly all analytical modeling of the problem

has focused on an idealization that the detonation

itself is small relative to the length of the pipe, and

so the reflected shock wave can be treated as if it

originated immediately from reflection of the equi-

librium, Chapman-Jouguet (CJ) state. From this for-

mulation several solutions have been obtained since

it considers shock reflection directly into the known

Taylor-Zel’dovich (TZ) wave. Despite this, all models

are inexact and make assumptions about the post-

reflected-shock flow. Stanyukovich [1] first solved the

problem by using the acoustic approximation for the

reflected shock and assuming that the post-shock ther-

modynamic state was uniform. Shepherd et al. [2]

performed numerical simulations of the fully non-

linear equations and also observed uniformity in the

post-reflected-shock thermodynamic state, corrobo-

rating Stanyukovich’s assumption. This observation

was used by Karnesky et al. [3] to develop a quasi-

empirical approach, where a pressure trace from the

detonation-reflecting pipe end wall was utilized with

Preprint, accepted for publication in Shock Waves, March 2023.

Published version available online at

https://doi.org/10.1007/s00193-023-01126-5

2 PROBLEM FORMULATION

the assumption that the post-reflected-shock pressure

was uniform. Strachan et al. [4] formulated a theory

for shock propagation through a generally heteroge-

neous upstream flow by extending methods developed

by Chisnell [5], which neglect influences behind the

shock. Strachan et al. evaluated their theory using

the idealized problem of detonation reflection and

compared their model to experimental results. Schoef-

fler and Shepherd [6] instead assumed that the post-

reflected-shock velocity gradient was uniform, which

enabled a more generally applicable model using inte-

gration of an otherwise-exact shock-change equation.

All of these models produce similar results since they

use identical assumptions for the upstream flow and

reflection from the CJ state.

Experiments have mostly been performed to inves-

tigate the structural dynamic loading induced by det-

onation reflection [2, 3, 7, 8]; they consequently

identified that some aspects of the observed reflected

shock conflict with what is predicted by ideal model-

ing. Inconsistencies were found in the post-reflected-

shock pressure by Karnesky et al. [3, 9]. Bifurcation

of the reflected shock was investigated as a possi-

ble explanation but was not observed in experiments

[10, 11]. Damazo and Shepherd [12, 13] performed

highly-resolved near-wall measurements of normal

detonation reflection and found that the reflected

shock speed was initially much larger than expected.

It was proposed that upon reflecting into the inci-

dent detonation’s induction zone, the twice-shocked

reactants explode, driving a blast wave into the flow

ahead, which results in the measurements of elevated

shock speed near the wall [13]. The authors used

a ”square-wave” detonation model to demonstrate

aspects of these dynamics, however further research

was warranted to clarify the evolution of the near-

wall reflected shock. The observation of anomalous

reflected shock speeds illustrated that the ideal detona-

tion assumed by modelers neglects important features

of the underlying physics.

The purpose of the present work is to investi-

gate the shock dynamics from normal reflection of

finite thickness detonations by the analysis of numeri-

cal simulations. One-dimensional simulations of both

supported and unsupported detonations are performed,

where only in the latter case does the TZ form, and

so its effects can be isolated. A one-gamma, one-step

detonation model with reduced activation energy is

used to simplify exact analysis of simulation results

with methods that describe shock propagation through

nonuniform flows. The theory is an extension to the

shock-change equations, used previously to model

ideal detonation reflection [6]. The derived equation

describes shock acceleration caused by the equation’s

forcing, or source terms, which correspond to the net

effects of flow gradients ahead of and behind the

shock wave. The application of these methods to ana-

lyze simulation data enables clear conclusions to be

drawn and explain the observed reflected shock accel-

eration. Results show that the reflected shock exhibits

large near-wall acceleration, overshooting the speed

predicted by reflecting the CJ state. The large shock

acceleration is shown to be caused by shock prop-

agation through the highly nonuniform flow in the

incident detonation. The subsequent evolution of the

shock depends on the detonation thickness relative

to the pipe length. Simulations of unsupported det-

onations with varied pipe lengths illustrate how the

finite detonation thickness affects the reflected shock’s

entire trajectory compared to predictions by simpli-

fied modeling. It is shown that for sufficiently short

pipes, there is a transition in the type of reflected

shock behavior, where the shock velocity increases

monotonically until reaching the end of the TZ wave.

Factors underlying all of the above reflected shock

dynamics are discussed using the derived shock-

change methods.

The structure of this article is as follows. In

section 2, the problem of detonation reflection and

the assumptions inherent in the methods used here are

described. Numerical simulation methods and results

are given in section 3. Methods to describe shock

propagation through heterogeneous media are first

derived in section 4.1 and then applied to the simula-

tion results in section 4.2. Results from this study are

further discussed in section 5.

2 Problem Formulation

Normal detonation reflection produces a shock that

propagates through the incident detonation, hence

treatment of the problem requires first a formulation

of the incident detonation. Detonations are coupled

combustion and shock waves that have been well

established to be intrinsically multidimensional [14,

15] with significant influence from turbulence [16].

This is in contrast to normal shock waves that are

exactly treated by one-dimensional and inviscid mod-

els. Even in one dimension, the coupling between

chemical reaction and the shock front is well known

to be highly unstable [17], resulting in an oscillating

detonation speed, which leads to the cellular structure

2 PROBLEM FORMULATION

observed in multidimensional detonations [18]. High

levels of instability in the detonation front enables

compressible turbulence to accelerate reactant burn-

ing and energy release [19]. The detonation front

nonetheless exhibits periodicity such that a transverse

or span-wise average describes mean quantities of a

one-dimensional detonation like its speed and reaction

zone structure. These mean structures may be treated

using the highly-successful lower-dimensional and

inviscid theories, however they are not equivalent. The

models are only approximate but effective and useful

for describing a transverse-mean detonation structure.

Furthermore, although the one-dimensional detona-

tion can be unstable, only a quasi-steady structure is

considered important here.

The well-known classical theories describe the

detonation speed, equilibrium thermodynamic state,

reaction zone structure, and downstream flow.

Chapman-Jouguet (CJ) theory predicts the detona-

tion speed where chemical equilibrium coincides with

a sonic point, terminating the detonation. The sonic

point isolates the shock-driven combustion in the deto-

nation, described by the Zel’dovich-Neumann-D

oring

(ZND) equations, from the downstream fluid dynam-

ics. A zero-velocity downstream boundary condition,

as is typical in pipes, requires an unsteady expan-

sion, i.e., the TZ wave forms to bring the high-velocity

combustion products to a halt. A piston boundary

condition can instead match the combustion product

velocity, so that no expansion wave forms. These two

cases describe an unsupported and supported detona-

tion, respectively.

In addition to only using one-dimensional anal-

ysis, a highly-simplified one-gamma one-step deto-

nation model with reduced activation energy is used

throughout this article. The reduced activation energy

stabilizes the detonation [20], and the one-gamma

one-step model enables straightforward analysis with

theory, however as a consequence important charac-

teristics of real chemical reaction mechanisms are

lost. Because of the reduced activation energy, there

is no considerable induction zone, and coupling

between the post-reflected-shock thermodynamic state

and reaction rate is weak. Damazo and Shepherd

[13] discussed the possibility that prompt explosion

of twice-shocked reactants behind the initially-frozen

reflected shock may drive a substantially faster over-

driven detonation in the incident detonation’s induc-

tion zone. This possibility cannot be investigated using

the present detonation model with reduced activa-

tion energy. Post-reflected-shock chemical reaction,

subsequent energy release, and its effect to acceler-

ate the reflected shock are considered here, however

the results are particular to the single-step model and

cannot be used to generalize for all cases of normal

detonation reflection. Analysis of the general problem

of detonation initiation in the incident detonation’s

induction zone and the ensuing dynamics remains for

future research.

The incident detonation shock reflects and prop-

agates first through the detonation’s reaction zone

structure, followed by the post-detonation flow. Prop-

agation of the reflected shock wave through the unsup-

ported detonation is diagrammed in Fig. 1, where

after traversing the detonation itself the shock trans-

mits into the TZ wave. The TZ wave is a self-similar

centered expansion that scales with the distance the

detonation has traveled. The one-dimensional solu-

tion for the TZ wave is well known [21]. Upon the

detonation’s collision with the wall, the TZ wave

extends the entire length of the pipe. Hence, there are

two essential length scales describing reflection of an

unsupported detonation: a length scale describing the

detonation thickness,

∆

, and the length of the pipe,

. If the detonation thickness is much smaller than

the pipe length such that

∆

→

, then all gra-

dients at the detonation sonic point go to zero, and

the detonation becomes supported. This shows that

the supported detonation is a particular case or the

limit of an unsupported detonation. Likewise, in the

same limit, the influence of the detonation structure

on the propagation of the reflected shock through the

TZ wave is negligible, and the reflected shock can

be approximated to have originated from reflection of

the CJ state. This suggests that, in the limit of large

pipe length, the reflected shock trajectory might be

formulated as a boundary layer problem, where the

outer problem describes the reflected shock propa-

gation through the TZ wave, and the inner problem

describes the shock’s propagation through the det-

onation reaction zone. The outer problem has been

modeled historically, as described above, where the

convenient idealized detonation was used. Here we

examine the inner problem, reflection of the sup-

ported detonation, and how the length scales in the

unsupported detonation reflection affect the reflected

shock’s evolution.

An important consequence of the one-step deto-

nation model is that the detonation is described only

by a single length scale. Typically, with more detailed

reaction mechanisms the ZND model exhibits two

principle length scales: one describing the induction

2 PROBLEM FORMULATION

0.2

0.4

0.6

0.8

1.0

0.6

0.8

1.0

1.2

1.4

TZ wave

zero-velocity

plateau state

reactants

characteristic

(

)

reflected

shock

detonation

wall

TZ wave

detonation

reactants

characteristic

shock front

reflected

shock

Fig. 1

Diagram of one-dimensional detonation reflection in a pipe

of length,

. The detonation, traveling at the Chapman-Jouguet (CJ)

speed,

, impacts the wall on the right, reflecting a shock wave

that travels to the left. The reflected shock traverses the detonation

itself followed by the Taylor-Zel’dovich (TZ) wave until it reaches

the uniform zero-velocity plateau state behind the TZ wave

length, and a second describing the exothermic energy

release length [15]. The single length scale in the

present analysis allows the detonation thickness to be

scaled unambiguously.

With a one-gamma detonation model, the perfect

gas Rankine-Hugoniot equations can be solved to find

useful preliminary results describing the relative shock

speeds reflected by the von Neumann (vN) and CJ

states. Upon impact with the wall, the initial speed

of the reflected shock is such that the incident gas

at the vN point is brought to rest by the reflected

shock to match the zero-velocity boundary condition.

∆

∼

, then after sufficient time for the shock to

propagate into the equilibrium asymptote and for non-

steady effects to vanish, the shock speed approaches a

limit where the CJ state is brought to rest. The magni-

tude of these two shock speeds can be obtained exactly

for the one-gamma detonation. The vN reflected shock

speed,

, is given by

1 + 2

−

+ 1

(

−

(1)

and the CJ reflected shock speed,

, is given by

−

(

+ 1)

−

γM

+ 1

(2)

where

and

are the CJ speed and Mach

number, respectively. These are shock speeds in the

laboratory-fixed frame, and

is the ratio of specific

heat capacities.

The ratio of these two reflected shock speeds is

plotted for

< γ <

and select CJ Mach num-

bers in Fig. 2. The ratio illustrates that for a gas the

shock speed ratio is always less than unity, i.e. the vN

reflected shock is always slower than the CJ reflected

shock. Since the reflected shock always originates

from the vN point, the reflected shock must acceler-

ate in order to reach the asymptotic CJ reflected shock

speed.

Fig. 2

Ratio of vN and CJ reflected shock speeds for a one-gamma

detonation

The ratio of vN and CJ reflected shock speeds is

also plotted for the strong detonation approximation in

Fig. 2, where the ratio is given by

≈

−

3 +

1 +

(17

+ 2)

(3)

and the ratio only equals unity when

9 +

√

≈

(4)

3 NUMERICAL SIMULATIONS

For the one-gamma detonation model used in this

article with

= 1

and

= 7

, the vN and CJ

reflected shock speeds are, respectively,

= 0

(5)

3 Numerical Simulations

3.1 Methods

Direct simulation of the one-dimensional reactive

Euler equations was performed using the finite-

volume CFD toolbox, OpenFOAM-9 [22], with the

additional library blastFoam-5 [23].

A one-gamma, one-step detonation model was

used for all simulations. An Arrhenius rate law was

used for the single first-order irreversible reaction.

The reaction equation, rate constant, and enthalpy are

given by

→

(

) =

exp(

−

)

−

Q ,

(6)

where

and

are reactant and product species,

respectively,

and

are their respective

enthalpies,

is the specific heat capacity at con-

stant pressure,

is the activation temperature,

is the pre-exponential factor, and

is the energy

release. The energy release was specified in the one-

gamma detonation through the CJ Mach number with

= 7

[24]. The gamma or ratio of specific heats

chosen was

= 1

. Although these values for the CJ

Mach number and gamma are typical, the parameters

are chosen for convenience and not to replicate the

properties for a particular mixture.

In order to simulate stable unsupported detona-

tions, a reduced activation energy,

= 3

, was

used, where

is the temperature at the von Neu-

mann point. The reduced activation energy substan-

tially decouples the reaction rate from the post-shock

thermodynamics, which eliminates instabilities [20] at

a cost of limiting the generality of results, as discussed

in section 2.

The rate law pre-exponential constant was used to

scale the ZND reaction zone size, which was taken

to be the half-reaction-zone thickness,

∆

, which

is the width of the ZND reaction zone where the

reaction progress variable equals one half. The ZND

reaction zone structure for this one-gamma, one-step

detonation model is shown in Fig. 3.

0 0

000

0 0

Fig. 3

ZND reaction zone structure for one-gamma, one-step deto-

nation model with

= 1

= 7

= 3

, and scaled by

half-reaction-zone width,

∆

/

. Quantities shown are (a) pressure,

(b) temperature, (c) thermicity, and (d) reaction progress variable,

with all appropriately nondimensionalized by the CJ state

The one-step reaction mechanism was constructed

and implemented first in Cantera [25] using the tech-

niques described by Kao [26]. Thermodynamics are

computed in both Cantera and OpenFOAM using

NASA-7 polynomials [27]. The single irreversible

reaction was translated to OpenFOAM’s native units

using the appropriate conversion factors. The ZND

equations were solved using the routines in the Shock

and Detonation Toolbox [28]. The temperature, pres-

sure, lab-frame velocity, and reaction progress vari-

able corresponding to the ZND solution were used to

prescribe the initial condition for all simulations in

OpenFOAM. The relevant initial condition and bound-

ary conditions for the supported and unsupported

detonation simulations are discussed in greater detail

below.

For all simulations, a rightward propagating deto-

nation is considered such that the detonation impacts

the right domain wall and reflects a shock wave to

the left. However, simulation results are flipped to

describe rightward reflected shock propagation with

positive lab-frame shock speeds.

For unsupported detonation simulations, results

for the reflected shock trajectory are compared with

3 NUMERICAL SIMULATIONS

those given by the idealized model with an infinites-

imal detonation thickness. The specific model used

is the one described by Schoeffler and Shepherd [6],

where a shock-change equation is integrated assum-

ing that the post-reflected-shock velocity gradient is

uniform.

Simulation data were post-processed to compute

the reflected shock position, velocity, and relevant

pre-shock and post-shock quantities. An algorithm

was implemented to consistently extract these quan-

tities by fitting the upstream and downstream flows

bounding the shock. The finite volume solver gener-

ates shocks that are approximately five cells thick. At

each saved time step, the maximum derivative of the

pressure data was used to initially locate the shock

wave. Using the pressure data, the shock wave, the

immediate upstream flow, and the immediate down-

stream flow were each fitted to lines. Twenty grid

points were used to fit the upstream and downstream

flows, respectively, except when the shock was near

the wall, where the known shock thickness was used

to approximate the upstream and downstream coor-

dinates. Using additional grid points in fitting the

upstream and downstream flows introduces nonunifor-

mity in fit residuals resulting from nonuniform flow

gradients, which would require higher-order poly-

nomials to properly fit. Since only the immediate

upstream and downstream states were of interest,

the linear fit to twenty grid points was found to be

adequate for evaluating these states and the veloc-

ity gradient. The two intersection points between the

shock-fitted line and the upstream and downstream

lines were used as coordinates for the pre-shock and

post-shock states, respectively. The position of the

shock was taken to be the average of these two points.

All other pre-shock and post-shock quantities were

then computed by interpolating simulation data at

these coordinates. An example of the results of this

post-processing algorithm is shown in Fig. 4.

The shock velocity was computed by numeri-

cally differentiating the position data. The post-shock

velocity gradient was extracted from the slope of the

line fitted to the downstream flow velocity. The time-

derivative of pre-shock quantities was also computed

by numerical differentiation.

All results presented below are nondimensional

using the following normalizations for pressure,

temperature,

, density,

, velocity,

, time,

, and

2.10

2.15

2.20

2.25

simulation data

fitted lines

shock position

pre-shock, state 1

post-shock, state 2

Fig. 4

Diagram of numerical shock features captured by post-

processing algorithm. Pressure is plotted, normalized by the CJ

pressure, against the spatial coordinate, normalized by the half-

reaction-zone width,

∆

/

distance,

, T

, ρ

, t

ref

, x

∆

(7)

where the tilde identifies dimensional variables, the

subscript

refers to the CJ state, and the reference

time scale is

ref

∆

(8)

Note that the velocity normalization uses the CJ speed,

, and not the gas velocity at the CJ state.

3.2 Simulation of Supported Detonation

Reflection

3.2.1 Simulation Details

Reflection of a supported detonation was simulated

by initializing the ZND reaction zone structure with

shock front coincident with end wall. The right

domain boundary is chosen to be the reflecting wall

with a zero-velocity boundary condition. The equi-

librium asymptote of the ZND model, i.e., the CJ

state, is extended to a sufficient domain length for the

reflected shock to propagate. The domain length cho-

sen was

= 80

. The prescribed initial condition for

the simulation is shown in Fig. 5.

Zero-gradient boundary conditions are used for all

other quantities, including the left domain wall veloc-

ity. The grid resolution was 300 cells per

∆

= 1

(i.e.,

300 cells per

∆

), and time steps are saved every

∆

= 0

077

3 NUMERICAL SIMULATIONS

0.0

0.5

1.0

1.5

2.0

Fig. 5

Initial conditions for simulation of supported detonation

reflection. The von Neumann point is located at the right wall

3.2.2 Simulation Results

Near-wall evolution of the pressure and lab-frame

velocity are plotted in Fig. 6. Fig. 6(a) shows a rapid

decrease in post-shock pressure as the left-running

shock wave accelerates away from the wall. Fig. 6(b)

shows that the velocity of the post-shock flow is nega-

tive with an approximately constant gradient between

the post-shock state to the zero-velocity boundary

condition.

77.0

77.5

78.0

78.5

79.0

79.5

80.0

(a)

77.0

77.5

78.0

78.5

79.0

79.5

80.0

0.0

0.5

(b)

Fig. 6

Near-wall evolution of pressure and lab-frame velocity from

reflected shock wave propagation. The wall is at

= 80

. The time

difference between each plot is

∆

= 0

The velocity of the reflected shock wave is plotted

over time in Fig. 7. The velocity is initially equal to

the shock speed reflected by the von Neumann point,

which is

= 0

. As the reflected shock prop-

agates away from the wall, it rapidly accelerates and

overshoots the CJ reflection limit,

= 0

. The

shock speed reaches a maximum velocity of

at time

= 3

and position

= 2

. The

shock speed subsequently decays over a substantially

longer time, first reaching the CJ reflection asymptote

at time

= 77

and position

= 37

. The shock

decays slightly below the CJ reflection asymptote to

a minimum of

= 0

at time

= 116

, before

accelerating again toward the asymptote.

100

150

0.2

0.4

0.6

0.8

simulation

vN reflection

CJ reflection

Fig. 7

Lab-frame reflected shock speed from supported detonation

reflection

Although the shock wave achieves a large peak in

lab-frame velocity, the shock Mach number, plotted in

Fig. 8, is at a maximum initially and then decays to

the CJ reflection limit. The initial decay of the shock

Mach number is very rapid and then subsequently

slower. The shock Mach number first decays below the

CJ reflection limit, before slowly increasing toward

the limit.

The effect of the shock acceleration is to displace

the space-time trajectory of the shock relative to the

straight-line trajectory of a shock traveling at the CJ

reflection speed, as shown in Fig. 9. The magnitude of

this displacement,

δx

, is given by integrating the dif-

ference between the reflected shock velocity and the

CJ reflection limiting value, i.e.,

δx

(

)

−

t ,

(9)

where

is the end of the simulation. Integrating the

simulation data yields

δx

= 5

. The significance of

3 NUMERICAL SIMULATIONS

100

150

1.5

2.0

2.5

3.0

simulation

vN reflection

CJ reflection

Fig. 8

Shock Mach number of reflected shock from supported det-

onation reflection

this displacement is that it manifests as a kinked space-

time diagram when time-of-arrival measurements are

used in experiments to examine the shock’s trajectory,

as illustrated by Damazo and Shepherd [13].

100

150

200

simulation

CJ reflection

Fig. 9

Space-time diagram for reflected shock from supported det-

onation simulation compared with constant shock speed from CJ

reflection

The pressure-transient at the end wall is plotted in

Fig. 10. The pressure is initially very large, although

the simulation does not quite resolve the vN reflec-

tion pressure. Similar to the shock speed and Mach

number, the pressure decays slightly below the CJ

reflection limit, before increasing toward it. The pres-

sure first intersects the CJ reflection limit at time

. Because of the finite detonation thickness there

is an excess impulse exerted on the pipe end wall,

which can be found from the integral,

wall

(

)

−

= 26

(10)

= 26

∆

(11)

where

is the excess impulse,

is the end wall area,

wall

is the pressure at the end wall, and

is the CJ

reflection pressure. The dimensional impulse is also

given in (11) to make the scaling explicit. The vast

majority of this excess impulse is exerted during the

initial decay from the peak pressure, i.e., over time

= 16

. This impulse increment is in addition to the

impulse exerted by CJ reflection, and, in general, the

impulse exerted by detonation reflection continues to

increase monotonically with time [29].

100

150

simulation

vN reflection

CJ reflection

Fig. 10

End-wall pressure over time following supported detona-

tion reflection

3.3 Simulation of Unsupported

Detonation Reflection

3.3.1 Simulation Details

An unsupported detonation was simulated by initial-

izing a width of the ZND model for a rightward

propagating detonation at the left domain wall, so that

upon simulating the wave’s evolution a zero-velocity

boundary condition at the left domain wall gener-

ates the TZ wave as the detonation propagates. The

width of the ZND model chosen was

5∆

. An ini-

tial kink in the propagating detonation between the

reaction zone and TZ wave smooths as the detona-

tion propagates. The evolution of this initial condition

is illustrated in Fig. 11. Zero-gradient boundary con-

ditions are used for all other quantities. The grid

resolution was 200 cells per

∆

= 1

(i.e., 200 cells per

∆

), and time steps are sampled every

∆

= 0

154

The unsupported detonation exhibits two length

scales: a length scale for the detonation thickness, here

∆

, and the domain length scale,

. Unsupported

detonation reflection was simulated for four domain

lengths,

50, 100, 200, and 500. All simulations

3 NUMERICAL SIMULATIONS

0.00

0.25

0.50

0.75

1.00

1.25

initial condition

Fig. 11

Velocity initial condition and its evolution for unsupported

detonation.

∆

= 7

between curves

begin with the same initial condition, such that by

extending the length of the domain the same detona-

tion is allowed to propagate further. Fig. 12 shows the

detonation pressure data just before the wave reaches

the domain end wall for the four simulation cases.

(a)

=50

100

(b)

=100

100

200

(c)

=200

250

500

(d)

=500

Fig. 12

Pressure data of unsupported detonation before impacting

right domain wall for four domain lengths

3.3.2 Simulation Results

Because there are two intrinsic length scales in this

problem, the simulation results are plotted using

both scales. Scaling by the reaction zone thickness

describes reflection of a given detonation in multiple

pipe lengths. Scaling by the domain length illustrates

how a varied detonation thickness produces different

reflected shock dynamics in a fixed length detona-

tion tube. When scaling by the domain length, the

coordinates in space and time are divided by the

nondimensional domain length, for example, time is

given by

∆

(12)

Results are only compared with the ideal model when

using the domain length scale, since the ideal model

only exists for an infinitesimal detonation thickness.

The reflected shock velocity is plotted over time

for the four domain lengths in Fig. 13. For all cases,

the shock velocity is initially equal to the velocity

of the shock reflected by the von Neumann point.

Fig. 13(a) shows that all reflected shocks exhibit a

period of rapid acceleration near the wall, which is

indistinguishable between the cases up until some

point in time. The subsequent dynamics vary distinctly

between the different simulated domain lengths. At

a later point in time, the shock trajectories for all

cases feature a local velocity maximum. This peak

corresponds to the reflected shock’s arrival in the zero-

velocity plateau state that exists behind the TZ wave.

In this region, all shocks begin to decay, hence the

maxima in shock velocity. The reflected shock behav-

ior between these two points, the early point where

trajectories first diverge and the later point where all

shocks reach the end of the TZ wave, is substantially

different. For the shorter domain length cases,

and 100, the reflected shock velocity increases mono-

tonically until reaching the end of the TZ wave. For

the longer domain length cases,

200 and 500,

the reflected shock velocity exhibits both a local max-

imum and minimum before the shock proceeds to

accelerate through the TZ wave.

The shock velocity trajectories are shown in Fig.

13(b) with time scaled by the domain length. As the

detonation thickness shrinks relative to a fixed domain

length, the shock velocity begins to approach the

ideal model with the exception of an early overshoot

resulting from the finite detonation thickness.

A consequence of these dynamics is that as

the domain length decreases, the maximum velocity

achieved by the reflected shock increases. The maxi-

mum shock velocity, which coincides with the shock’s

arrival at the end of the TZ wave, is plotted in Fig. 14

for each domain length.

The space-time trajectories for the four simulation

cases are plotted with the ideal model in Fig. 15. The

trajectories are only shown up until the shock reaches

4 SHOCK-CHANGE ANALYSIS

100

150

200

250

0.0

0.2

0.4

0.6

0.8

1.0

(a)

=50

=100

=200

=500

vN reflection

0.0

0.2

0.4

0.6

0.0

0.2

0.4

0.6

0.8

1.0

(b)

vN reflection

CJ reflection

ideal model

Fig. 13

Lab-frame reflected shock velocity from unsupported det-

onation for varied domain length plotted against time, scaled by

reaction zone thickness in (a) and scaled by domain length in (b)

100

200

500

0.80

0.85

0.90

0.95

s,max

Fig. 14

Maximum reflected shock velocity from unsupported deto-

nation reflection

the end of the TZ wave. Again, as the detonation thick-

ness shrinks the shock trajectory approaches the ideal

model.

Pressure at the detonation-reflecting end wall is

plotted over time for both scales in Fig. 16. Again,

0.0

0.1

0.2

0.3

0.4

0.0

0.1

0.2

0.3

0.4

=50

=100

=200

=500

ideal model

Fig. 15

Space-time trajectories of reflected shock from unsup-

ported detonation for four domain lengths. All trajectories are

terminated where the reflected shock arrives at the end of the TZ

wave

the simulation results do not quite capture the ini-

tial vN reflection pressure. In all cases, the pressure

decays monotonically. Fig. 16(a) shows that, as the

reflected shock propagates through the reaction zone,

the pressure decay is nearly indistinguishable between

the different cases, and the subsequent decay as the

shock propagates through the TZ wave scales with

the length of the domain. Fig. 16(b) illustrates how a

thicker detonation for a given detonation tube length

increases the duration of the pressure pulse at the

end wall. The excess impulse exerted against the end

wall resulting from the finite detonation thickness is

approximately equivalent to the result from simulation

of the supported detonation.

4 Shock-Change Analysis

Analysis of the numerical simulation results requires

methods that describe shock propagation through

nonuniform upstream media. Shock propagation the-

ories have been developed by numerous authors

over many decades. A common technique is to con-

sider the transport of shock-jump quantities, enabling

a combination of the Euler equations with the

Rankine-Hugoniot equations. The formalism, some-

times referred to as singular surface theory [30, 31],

enables derivation of equations referred to by sev-

eral names including shock-propagation and shock-

amplitude equations, however, in detonation physics,

they are most often known as the shock-change

equations, as presented by Fickett and Davis [14].

There have been many contributions to the methods,

some include those by Chen and Gurtin [32], Bowen

4 SHOCK-CHANGE ANALYSIS

100

150

200

250

(a)

=50

=100

=200

=500

vN reflection

CJ reflection

0.0

0.2

0.4

0.6

(b)

=50

=100

=200

=500

vN reflection

CJ reflection

Fig. 16

End-wall pressure plotted over time following unsupported

detonation reflection

and Chen [33], Nunziato and Kennedy [34, 35], and

most recently by Radulescu [36, 37].

Shock propagation into nonuniform media has

also been examined by several authors, however

shock-change methods have typically only been used

for uniform upstream flows. A notable exception is the

work by Nunziato and Walsh [38–40]. Other methods

were developed by Chisnell [5], and later generalized

by Strachan et al. [4], to describe shock accelera-

tion from upstream gradients by neglecting nonsteady

influences behind the shock wave. The propagation

of two-dimensional wave fronts through regions of

nonuniformity has been treated with extensions to

Whitham’s theory [41], for example by Collins and

Chen [42] and Catherasoo and Sturtevant [43, 44].

Methods using shock-change equations are dis-

tinct because they provide an exact description of how

shock propagation is affected by upstream and down-

stream flow gradients. However, the gradients must

be known

a priori

in order to make predictions of

shock trajectories. To use the equations to compute

shock evolution additional models are required and,

although inexact, can yield highly accurate results

[6, 37]. Methods have been developed to make higher-

order approximations, as described by Best [45] and

Sharma and Radha [30].

A powerful application of the shock-change

equations is instead to use them to analyze simula-

tion or experiment data, and this is the strategy we are

taking here. The appropriate shock-change equation

describing shock propagation through a nonuniform

flow is first derived generally and then specified for a

one-gamma perfect gas. The result is an ordinary dif-

ferential equation for the shock velocity, which makes

explicit how shock acceleration originates from forc-

ing terms in the equation that describe the influence

of flow gradients both ahead of and behind the shock

wave.

4.1 Derivation of Shock-Change Equation

The equations of motion for a chemically reacting

flow in one dimension, neglecting diffusion of mass,

momentum, and heat, are

−

∂u

∂x

(13)

−

∂P

∂x

(14)

(15)

(

= 1

,...,N

)

(16)

where species

has mass fraction

, molecular weight

, and net production rate

. An equation of state

of the form

(

P,ρ,

Y

)

can be expanded as

∂h

∂ρ

∂h

∂P

ρ,

∂h

∂Y

ρ,P,Y

̸

(17)

where

Y

is a vector containing the

species mass

fractions.

Expressing the total differentials in (17) as mate-

rial derivatives and combining with (15) gives

−

∂h

∂P

ρ,

∂h

∂ρ

∂h

∂Y

ρ,P,Y

̸

(18)

4 SHOCK-CHANGE ANALYSIS

After further manipulation, the following result is

obtained

ρa

(19)

where

is the frozen sound speed [46],

∂h

∂ρ

−

∂h

∂P

ρ,

(20)

and

is the thermicity,

−

∂h

∂Y

ρ,P,Y

̸

(21)

is the coefficient of thermal expansion, and

is the

specific heat capacity at constant pressure.

The result (19) is well known and sometimes

referred to as the adiabatic-change equation [14]. It

can be used to eliminate the material derivative of

density from the continuity equation (13),

−

ρa

∂u

∂x

ρa

σ,

(22)

resulting in a system of two equations with the

momentum equation (14).

Fluid elements travel with velocity

relative

to a consistent reference frame, often termed the

laboratory-fixed frame. The time variation of flow

quantities, such as pressure, density, and velocity,

along particle paths is described by the material

derivative. Similarly, the time variation of flow quan-

tities along the path of a shock wave can be computed

with the corresponding total derivative. An arbitrary

field quantity,

(

x,t

)

, evaluated along the path

of a shock wave is given by

(

)

, where

(

)

is the position of the shock wave. Then, the time

variation of

along the path of the shock is

∂q

∂t

∂q

∂x

∂q

∂t

∂q

∂x

(23)

where

is the lab-frame shock velocity. The time

derivative of a quantity along the path of the shock

will be subsequently referred to as the quantity’s

shock

derivative

The material derivative can be written in terms of

the shock derivative with

+ (

−

)

∂q

∂x

(24)

which can be used to transform (14) and (22). This

change of variables gives

+ (

−

)

∂u

∂x

−

∂P

∂x

(25)

+ (

−

)

∂P

∂x

−

ρa

∂u

∂x

ρa

σ .

(26)

Equations (25) and (26) can be combined to either

eliminate the pressure-gradient term or the velocity-

gradient term to arrive at a shock-change equation.

Here we consider only the velocity gradient and elim-

inate the pressure gradient, and the corresponding

shock-change equation is

−

∂u

∂x

(27)

where the

is the sonic parameter,

= 1

−

(28)

and

−

is the flow velocity in the shock-

fixed frame. The subscript

was also introduced to

indicate that the variables are post-shock quantities,

whereas a subscript

refers to a pre-shock or upstream

quantity. The subscripts make explicit the relationship

between the shock-change equations and the shock-

jump equations, which describe how state

quantities

are related to state

quantities. The relevant Rankine-

Hugoniot or shock-jump relations can be expressed as

(

,ρ

Y

)

(29)

(

,ρ

Y

)

(30)

i.e., the post-shock pressure and velocity are functions

of the shock speed, the upstream flow velocity, and the

upstream thermodynamic state. The thermodynamic

state is here defined by the set of independent variables

{

P,ρ,

Y

}

, however any equivalent set can be used.