DOE/ERy03077276
Courant Mathematics and
Computing Laboratory
U.S. Department of Energy
An Asymptotic Analysis of an
Expanding Detonation
J. Jones
CN (w
Research and Development Report
'^iioported by the Applied Mathematical Sciences
Drogram of the Office of Energy Research,
o c . Department of Energy under Contract
vo .'^ â€¢;J DEAC0276ER03077
t^ "(0 Â° :hematjcs and Computers
o o en
M H c
w e o a
o e a,
W W 0)
01 ca
uary 1987
il^s m NEW YORK UNIVERSITY
UNCLASSIFIED
DOE/ER/03077276
UC32
Mathematics and Computers
New York University
Courant Mathematics and Computing Laboratory
AN ASYMPTOTIC ANALYSIS OF AN EXPANDING DETONATION
J. Jones
January 1987
Supported by the Applied Mathematical Sciences
subprogram of the Office of Energy Research,
U. S. Department of Energy under Contract No.
DEAC0276ER03077
UNCLASSIFIED
1
DISCLAIMER
This report was prepared as an account of work .sponsored
by an agency of the United States Government. Melther
the United States Government nor any aprency thereof, nor
any of their employees, makes any warranty, express or
Implied, or assumes any lep;al liability or responsibility
for the accuracy, completeness, or usefulness of any
information, apparatus, product, or process disclosed, or
represents that Its use would not infrlnr^e privately
owned rights. Reference herein to any specific commercial
product, process, or service by trade name, trademark,
manufacturer, or otherwise, does not necessarily constitute
or imply its endorsement, recommendation, or favorlnf^ by
the United States Government or any ap;ency thereof. The
views and opinions of authors expressed herein do not
necessarily state or reflect those of the United States
Government or any agency thereof.
Printed In U.S.A.
Available from
National Technical Information Service
U.S. Department of Commerce
5285 Port Royal Road
Springfield, VA 22l6l
11
An Asymptotic Analysis of an Expanding Detonation
James Jones â– *â– â–
Courant Institute of Mathematical Sciences
New York University
New York, N.Y. 10012
ABSTRACT
An expanding cylindrically or spherically symmetric detonation
is analyzed in a regime in which the radius of the detonation is
much greater than the width of the reaction zone. Under this
assumption the fundamental equations may be approximated by
a system of autonomous ordinary differential equations for the
flow velocity and a reaction progress variable. The independent
variable of this system is a radial variable in the rest frame of
the detonation front. The radius of curvature and the detona
tion speed enter the system as parameters. At zero curvature
this system reduces to the plane wave equations of Zeldovich,
von Neumann and Doering. The plane wave equations possess a
degenerate bifurcation point with a nilpotent linear pjirt, which
bifurcates into a saddle node when the radius is finite. Any
smooth transonic solution must pass through the saddle node.
This fact determines the wave speed implicitly as a function of
radius. To leading order, the correction to the detonation speed
as a function of curvature is proportional to the curvature, on
the basis of formal and numerical considerations.
111
An Asymptotic Analysis of an Expanding Detonation
12 1
James Jones ' '
Qjurant Institute of Mathematical Sciences
New York University
New York, N.Y. 10012
Table of Contents
1. Introduction 1
2. Derivation of the Model 3
3. Reduction of Order for the System (2.9) 11
4. Statement of Main Results 18
5. The Proposed Normal Form 27
6. The Transonic Critical Point 33
7. Proof of Theorem 4.1 40
8. Conclusion 52
References 53
Appendix 55
1. Supported in pan by National Science Foundation Grant Nos. MC58207965 and MCS8301255.
2. Supponed in part by Army Research Office, Contraa No. DAAG2984K0130.
3. Supported in part by U.S. Dept. of Energy, Contraa No. DEAC0276ER03077.
 1 
1. Introduction. The diverging detonation is interesting because the non
linear interaction between the chemistry and the hydrodynamics is displayed
in a particularly sharp form. The theoretical study of detonations in a diverg
ing geometry began in the late 1940's with the work of Eyring, Powell, Duf
fey, and Parlin [9] and of H. Jones [12]. These papers were largely con
cerned with the relationship between the detonation velocity and the finite
diameter of the charge which has come to be called the "diameter effect".
Wood and Kirkwood (1954) [19] were the first to elucidate the relationship
between the wave velocity and the radius of curvature of the detonation
wave. This fundamental work has been extended by a number of authors [3,
7, 8, 15]. In all of these works the geometry considered was a cylindrically
symmetric charge with the detonation propagating parallel to the axis of sym
metry. An excellent review of this problem can be found in Fickett and
Davis [10], pp. 199  229. The axial geometry is useful to both experimental
ists and theorists because steady state solutions exist. For the experimentalist
this means that stationary diverging detonations can be produced in the
laboratory. For the theorist this reduces the problem to an analysis of a sys
tem of ordinary differential equations.
Although a substantial body of work on steady state, diverging detona
tions has accumulated, there has been little progress in developing an under
standing of the dynamics of diverging detonations. One reason for this is the
presence of nonuniformities that appear in straightforward attempts to per
form an asymptotic analysis of this problem. In the present work, we
present an analysis of a spherically symmetric detonation, or a radially pro
pagating cylindrically symmetric detonation (independent of the axial
coordinate). Although this is not a steady state problem, we will show that
to first order in the ratio of reaction zone width to radius of curvature, the
detonation may be modeled by a system of autonomous ordinary differential
equations in space, with time entering as a parameter. The nonuniformities
are exhibited as unboundedness of the vector field along the sonic locus, at
which the flow velocity relative to the shock equals the sound speed. A
bounded transonic solution, i.e. a solution crossing the sonic locus, must pass
through a critical point of the vector field. We will employ the methods of
bifurcation theory to show that an appropriate critical point exists. The evo
lution in time of the detonation is then determined by a shooting problem
connecting the state at the shock interface to the critical point. Our main
result, derived in conjunction with the numerical results of Bukiet [5, 6], is a
demonstration on the level of numerical computation that the leading order
correction to the detonation wave velocity is proportional to the curvature.
A second main result is the formal derivation of ordinary differential equa
tions for the leading order large radius asymptotics for the internal structure
of the expanding detonation wave. A third main result is the analysis of this
system of ordinary differential equations. We show that the plane wave
sonic critical point bifurcates into a unique critical point which is a saddle
point, and that this critical point is a smooth function of the shock curvature.
A normal form is proposed for the bifurcation point, and a one parameter
unfolding is presented illustrating the effect of curvature on the phase por
trait.
In a recent work, Bdzil and Stewart [4] have independently presented a
related asymptotic theory for the dynamics of curvilinear detonations. They
3
employ the strong shock approximation and a specialized rate law for which
the reaction terminates at a finite distance behind the shock, and coincides
with the sonic transition. These idealizations permit them to study two
dimensional detonation dynamics and to deal with the characteristic nonuni
formities in a more traditional fashion. Although Bdzil and Stewart
emphasize that their analysis cannot be expected to apply to more general
reaction kinetics, it is interesting to note that they aJso find that the leading
order correction to the speed of the detonation wave is linear in the shock
curvature.
2. Derivation of the Model. The equations for a cylindrically symmetric,
transport free, reactive, polytropic gas are
(2.1)
m
III f\
m,
E,+
m
+ P
mJE  p)
+
m~
pr

^(^ ^ P) ^
(P^)<
P
(p^)m _
P''
pR
= 0.
Here m = pM is the mass flax, p is the density, u is the radial velocity and p
is the pressure, i is the time and r is the radial coordinate. \ is the reaaion
progress parameter, which varies from (all reactant) to Mail product). The
total energy density is E = pe ~ m / 2p, where the specific internal energy e
4
is given by
(2.2)
e â€”
(1  M^ ,
P(71)
7 is the polytropic gas constant and q is the heat released per unit mass by
the complete reaction. R is the reaction rate function, which we assume has
the Arrhenius form
k{\
X) exp
A
r
,r, Â»o  J â€”j^ â€” ail.
where/ ^ (0. 1). The relevant observation here is that w scales as â€” . Let e
be the ratio of the reaction zone width n^^ to a typical radius. Since the wave
speed is constant in the plane wave limit, : and t are nearly proportional for
sufficiently large i. If we assume that uq and the plane wave speed i,) are
0(1), then : and / are O
I
â‚¬
V /
. We now transform the differential equations
(2.3) to an appropriate form for asymptotic analysis by translating the radial
variable r to the rest frame of the shock and rescaling t and : to be 0(1) for
large times. Define
(2.4) T s 6f
Ut)  â‚¬r(r)
X = :(t)  r = ^^  r
e
\=i\x = z\x
where t â€” j â€¢ Note that the transformation r  r involves a reversal of
ax
orientation, as shown in Fig. 2.1. This convention was chosen so that x will
be positive behind the shock in the reaction zone. The flow velocity relative
to the shock, oriented parallel to jc, is v. Now change variables from r, r, r,
and u to T, X, ^, and v respectively to obtain
(2.5) ep,  vp, + v,p = ~y^j^/^
Px
ev + vVj + â€” = e^
P
epx + v/p, + 7pv, = ^(7  l)p/?(X,r)  '^^^j ^/^
7
â‚¬X, + v\,  /?(x,r).
Figure 2.1
There are two alternate derivations of (2.5) which are worth mentioning.
First, one could transform to dimensionless variables and set all 0(1) parame
ters to 1. Alternately, rather than taking w and x to be (7(1) and : and t
large, one could assume that r and t are 0{l) and that u and x are small by
rescaling R and doing an "inner expansion" in v behind the shock.
The system (2.5) has the form
{2.5a) â‚¬W + 0(w) â€¢ w^  h(w, x, e).
8 
It will be seen in Section 3 that the desired solution of (2.5) must possess
both subsonic and supersonic regions with a smooth transition between. At
such a transition point the quasilinear operator O becomes singular. This
singularity may be interpreted as a turning point of the system. If we
attempt to solve for w^ in (2.5a) near the turning point we find that the x
derivatives blow up unless certain solvability conditions are satisfied. Specifi
cally h  ew must lie in the range of O. We will not attempt here to prove
asymptotic convergence of solutions of (2.5) as â‚¬  to the planar (e = 0)
solution, but will proceed in the spirit of formal perturbation theory by
assuming on physical grounds that a smooth transonic solution of (2.5) exists
which is differentiable with respect to e. The turning point structure will be
apparent in the final version of the model which we will derive in Section 3.
We thus assume that a representation
p(.x, t, â‚¬) = po(.r, t)  ep(.r, t)  PrcJx, t. e)
exists in e for p(x, f, e), where
,. PremC^' ^ ^
lim ; ;^ â€” ^
Â£0 6pi(x,0
uniformly in x and t. We also assume analogous representations for \\ p . \
and i. Now expand (2.5) formally in powers of â‚¬. grouping together all
terms of each power. The zeroth order equations are
(2.6) VqPo, ^ Vo,p =
Po
These are the equations for a one dimensional steady state detonation studied
by Zeldovich, von Neumann and Doering (ZND). The hydrodynamic equa
tions may be integrated to obtain p, v, and p as functions of \. The rate
equation then constitutes an ordinary differential equation for \ as a function
oi X,
(2.7) PqVq  m { = constant)
Po ~ Ps _ â– ,
dx Vq '
"v â€” 1 1
where a = ^ r. and V = â€” is the specific volume. The s subscript indi
^7+1 P
cate variables evaluated immediately behind the shock. The first of equations
(2.7) states that the mass flux is constant in the shock frame. The second
equation defines a line of slope m in the p, V plane, and is referred to as
the Rayleigh line. The third of equations (2.7) defines a family of
Hugoniot curves in the p , V plane. The intersections of the Rayleigh line
with the \ = Hugoniot curve determine the possible shock transitions. The
solution terminates at point where the Rayleigh line intersects the \  1
Hugoniot curve. In general, there is a one parameter family of solutions,
parameterized by the mass flux m, or for a given ambient state ahead of the
10
shock, by the wave speed ^o
The first order equations are
Po(to  Vo)
(2.8) po^ + VoPi,  v,po, ^ PoV;, ^ PiVo,
PoVo^ + PoVqViv ^ Po^iVo,, + PiVqVo., ^ Pu = Po4o
7Po(to  vo)
^(7  l)(po^i ^ pi^o) 
4o
Since the zeroth order equations are steady state, the time derivatives drop
out of the first order equations, and ^ = 0. Thus the first order solution is
"quasisteady" in the sense that the solution is given by steady state equations
with time entering as a parameter (through ^o ^d ^q) Note also that the tr
term which appears in the denominator of the geometric source terms in (2.5)
does not appear through first order. Thus, in the reaction zone, the deviation
of the flow divergence from the value at the shock is a second order effect.
If the time derivatives and e.r terms are dropped from (2.5) one obtains a
system of autonomous ordinary differential equations
^P(C  V)
(2.9) vp, + v,p =
w, + â€” =
C
11 
v\,  R{\J) .
It is easily verified that (2.9) possesses the same formal zeroth and first order
equations as does the original system (2.5). On the basis of the assumed
regularity of the asymptotic expansion for (2.5) we now take (2.9) as our
first order model for the expanding detonation and turn our attention to an
analysis of this system. Of course one could also pursue the more usual
course of studying the linear nonautonomous perturbation equations (2.8),
but the autonomous system (2.9) permits a much cleaner treatment by the
technique of phase space analysis.
3. Reduction of Order for the System (2.9). Note in (2.9) that e and I
occur only in the ratio â€” so that â‚¬ is really a redundant parameter. This
redundancy may be removed by setting â‚¬ = 1. This inverts the scale
transformation, so that t  t and i  z. These identifications will be made
henceforth.
The steady state energy equation is
The velocity equation in (2.9) can be written as
\iy\ ' vp.  .
These two equations may then be added to obtain
 12
1 .
Vp
= 0,
which integrates to yield Bernoulli's Law:
(3.1)
^w ^ e + Vp ^ fit)
Denote values upstream of the shock by the subscript a (for "ambient"). It
will be assumed that the ambient state is constant and unreacted (X^  0),
and that the ambient flow velocity is zero (u^  0, or v^ = z). For our
polytropic equation of state, (3.1) becomes
(3.2)
1
_v + â€”
2 1
1
1
1
c" â€” X.^
1 .1
V
1
Here c  (yp/p) is the sound speed. Note that we are able to connect
across the shock to the ambient state since Bernoulli's law is one of the
RankineHugoniot jump conditions. This fact determines the function /(/) in
(3.1).
Now eliminate/?, between the velocity and pressure equations to obtain
(3.3fl)
q{y  [)R
v. =
> >
c~ â€” V
^(7  1)^(1  X)exp
A7
iz  ^â– )c'
c  V
By (3.2) c~ is a known function of v and X. Therefore the right hand side of
the equation above is also a known function of v and X. This form of the
13
velocity equation may be combined with the rate equation
(3.3Z;) X.  
to obtain a self contained system of two equations for u and \.
For a steady undriven plane detonation, the reaction zone terminates at
the Chapman Jouguet. or CJ point v ^ c in the Hugoniot diagram. At this
point the Rayleigh line is tangent to the \ = 1 Hugoniot curve. This cri
terium determines a unique solution of (2.7). The diverging detonation is
weakened by rarefactions produced behind the shock front and terminates
below the sonic point on the weak detonation branch of the X ^ 1 Hugoniot
curve. This means that a sonic transition must occur in the reaction zone
from the subsonic flow behind the shock to the supersonic flow at termina
tion. However, the denominator in the velocity equation (3.3a) vanishes at a
sonic transition, so for a smooth sonic transition to occur, the numerator
must vanish simultaneously. The sonic transition is just the turning point
mentioned in Section 2, and the condition that the numerator vanish is
equivalent through first order to the solvability condition for w, in (2.5a).
We may use substitute v for c~ in (3.2) to obtain the sonic locus
i3Aa) v2 = ^  ^j^f + 2y.'qK.
7+1 ^
A solution of (3.3) which crosses the sonic locus will be called transonic. We
thus seek a transonic soluuon of (3.3) which satisfies
(3.4/7) q(y  1)/?  ic(r  v) =
at the sonic transition. A solution (v^., XJ of the system (3.4) will be
14
referred to as a sonic critical point. Note that when z â€” x, equation (3.4b)
yields /?  at the sonic critical point. If the detonation has not failed this
implies X â€” 1 at the sonic critical point, consistent with the aforementioned
result that a steady undriven plane detonation terminates at the CJ point.
The wave speed Jq and finaJ (X  1) flow velocity v^y for the plane wave
may be determined from equations (2.7) and (3.4a) and the condition
Vcj  ccj^ where ccj is the sound speed at the CJ point. The result is
(3.5)
^0
vcy =
[iT  l)qc; + 1 + (((7  l)qc;' + 1)"  1)^'=
2d
.L'2
_L "'â€¢2,12
y + 1
The positive solution in (3.5) must be chosen in order to satisfy the jump
conditions at the shock.
The system will (3.3) be easier to analyze if transformed into a more
conventional form. Define the singular change of variable
exp
(3.6)
, = ;.
_ Ay
cix'Y}
dx'
{cix'y  v(x')2)v(.r') â€¢
Now change variables from x to y to obtain
(3.7)
v,  qiy  \)k(l  X)v
X^  ^(l  X)(c  V).
(z  v)v ,
^ ^â€” c exp
2
Ax
Observe that the structure of the phase curves in the (v, X) plane is unal
tered by this change of independent variable since the transformed vector
15
field is proportional to the initial vector field. The integral curves have sim
ply been reparameterized to eliminate the singularity in the denominator of
the velocity equation. The right hand side of (3.7) is now bounded in the
region of interest, and the sonic critical point defined by (3.4) is recognized
as a stationary point of the system.
Several observations about (3.7) can be made immediately. Since the
first equation is proportional to v, the \ axis is a phase curve of the system.
Likewise, the \  1 line is a phase curve, since the second equation has a fac
tor of 1  \. These two phase curves intersect in a fixed critical point
(0, 1). When z â€” ^, the vector field is proportional to 1  X, so that the
entire X = 1 line is stationary. The z  ^ sonic critical point is thus embed
ded in a manifold of critical points. We will see that the z = ^c sonic critical
point is a bifurcation point for the system, i.e. a point in the phase plane
where the topology of the phase curves is unstable to small perturbations of
the vector field. These notions will be made more precise in the following
section.
For a fixed ambient state ahead of the shock, the RankineHugoniot
jump conditions determine a one parameter family of solutions behind the
shock, parameterized by the wave speed :. The critical point equations (3.4)
define v^, and X^ as functions of z and i. The desired solution of (3.7) will
connect the ambient state (via the jump conditions) to a sonic critical point;
this shooting problem defines the functional relationship between z and z.
The method of determining the full dynamic solution of (3.7) is now
clear. We first analyze the system allowing z and z to vary as independent
parameters, thus determining the phase plane structure throughout some
16
domain in parameter space. We then may recover the leading order dynamics
by solving the shooting problem. This paper is primarily concerned with the
first half of this procedure, i.e. determining the phase plane structure of
(3.7). The shooting problem has been solved numerically by Bukiet [5, 6].
After eliminating c^ between (3.2) and the second of equations (3.5), we
obtain
(3.8) c^  v2 = qi^  1)(1  \)  ^^(v2  vl),
where
V,  (v^, + ^ii'  zi)r~
is the flow velocity at the z = ^c sonic bifurcation point. With z and z
independent, v^ is now a function of z and zq.
We may obtain a single differential equation for v(\) by dividing the
first equation in (3.7) by the second to obtain
(3.9) (1  \){c'  v)^ = ^(7  1)(1  X)v  ^~^ c' exp ^
This equation clearly shows the (nonlinear) turning point character of the
equation. When z = ^c^ we may use (3.5) and (3.8) to write equation (3.9)
in the form
1
v^  2q\ir{\  X) ^
dv r
Â±mÂ£
â– dX =
The left hand side of this equation is the differential of the function
v:  2flix(l  X)
/(v, X) = V + ^ ^ ^ + C .
 17
Denoting by (\v, Xj any fixed reference point, we have
(3.10) V ^ ^^^^^'(^^^) ^ , . '^ ^ ~^^'^'  ^)
The choice v^ = v^, \ ^ 1 yields the separatrix solution
(3.11) (v  v,)=  2q^iK  I) =
for the plane wave. This result could also have been obtained by solving
(2.7) forv(X).
Denote by v  ^(X, r, i) the solution of the shooting problem for (3.9)
connecting the ambient state to a sonic ciitical point. Evaluating g at the crit