Numerical Computation of Detonation Stability - KAUST Repository

0 downloads 0 Views 5MB Size Report
library of the scientific Python stack), hence, determined our choice of the ...... use solver ode45 with relative and absolute tolerances set to 10−10 to solve the.
Numerical Computation of Detonation Stability

Dissertation by Dmitry I. Kabanov

In Partial Fulfillment of the Requirements For the Degree of Doctor of Philosophy

King Abdullah University of Science and Technology Thuwal, Kingdom of Saudi Arabia

April, 2018

2

EXAMINATION COMMITTEE PAGE

The dissertation of Dmitry Kabanov is approved by the examination committee:

Committee Chairperson: Aslan Kasimov Committee Co-Chair: Athanasios Tzavaras Committee Members: Ravi Samtaney, David Ketcheson, Rodolfo Ruben Rosales

3

© April, 2018

Dmitry I. Kabanov All Rights Reserved

4

ABSTRACT

Numerical computation of detonation stability Dmitry I. Kabanov Detonation is a supersonic mode of combustion that is modeled by a system of conservation laws of compressible fluid mechanics coupled with the equations describing thermodynamic and chemical properties of the fluid. Mathematically, these governing equations admit steady-state travelling-wave solutions consisting of a leading shock wave followed by a reaction zone. However, such solutions are often unstable to perturbations and rarely observed in laboratory experiments. The goal of this work is to study the stability of travelling-wave solutions of detonation models by the following novel approach. We linearize the governing equations about a base travelling-wave solution and solve the resultant linearized problem using high-order numerical methods. The results of these computations are postprocessed using dynamic mode decomposition to extract growth rates and frequencies of the perturbations and predict stability of travelling-wave solutions to infinitesimal perturbations. We apply this approach to two models based on the reactive Euler equations for perfect gases. For the first model with a one-step reaction mechanism, we find agreement of our results with the results of normal-mode analysis. For the second model with a two-step mechanism, we find that both types of admissible travelling-wave solutions exhibit the same stability spectra. Then we investigate the Fickett’s detonation analogue coupled with a particular reaction-rate expression. In addition to the linear stability analysis of this model, we demonstrate that it exhibits rich nonlinear dynamics with multiple bifurcations and chaotic behavior.

5

ACKNOWLEDGEMENTS

First and foremost, I would like to thank my supervisor, Prof. Aslan Kasimov, for advising me and suggesting ideas for projects and numerical experiments throughout my research. He created a friendly atmosphere governed by freedom and scientific curiosity in his research group which helped me to complete my PhD studies, and I am truly grateful for that. I would also like to thank my co-supervisor, Prof. Athanasios Tzavaras, who supported and encouraged me during the late stages of my PhD work. Additionally, I would like to thank the Committee members, Profs. Ravi Samtaney, David Ketcheson, and Ruben Rodolfo Rosales, for helpful suggestions and comments. I gratefully acknowledge everybody in the Department of Applied Mathematics and Computational Sciences, for all the enjoyment that I had while studying here and taking a lot of excellent classes on numerical methods, fluid dynamics, and high-performance computing. The support and care of all my friends in KAUST, especially, not in any particular order, of Vera Solovyeva, Slava Korneev, Naeemullah Khan, Sultan Albarakati, Dasha Bolbot, Yana Ukhova, and Alena Denisova, made my life easier and more enjoyable during my PhD journey. I am much obliged for that. For everybody else who I know in KAUST and have not mentioned here explicitly, I thank you too. Last, but not least, I am deeply grateful to King Abdullah University of Science and Technology for an opportunity to study here, in this amazing place on a shore of the Red Sea, combining a research institute with a resort.

6

To my parents

7

TABLE OF CONTENTS

Examination Committee Page

2

Copyright

3

Abstract

4

Acknowledgements

5

List of Abbreviations

10

List of Figures

11

List of Tables

13

1 Introduction 1.1 CJ and ZND theories . . . . . . 1.1.1 CJ theory . . . . . . . . . 1.1.2 ZND theory . . . . . . . 1.2 Stability computations . . . . . 1.3 Overview of the present work . 1.4 Outline . . . . . . . . . . . . . .

. . . . . .

14 16 16 19 24 29 32

2 Reactive Euler equations with one-step chemistry 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 ZND solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Linear unsteady problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Numerical algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Computational domain and integration of the ZND solution 2.5.2 Integration of the linearized system . . . . . . . . . . . . . . . .

34 34 36 38 40 41 41 43

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

8 2.5.3

Postprocessing via dynamic mode decomposition . . . . . . 2.5.3.1 Description of the DMD algorithm . . . . . . . . . . 2.5.3.2 Application of DMD to the time series of detonation velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.3.3 Testing the method with synthetic data . . . . . . . 2.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6.1 ZND and perturbation solutions . . . . . . . . . . . . . . . . . 2.6.2 Comparison with normal-mode results . . . . . . . . . . . . . 2.6.3 Neutral stability curves . . . . . . . . . . . . . . . . . . . . . . . 2.6.4 Code verification by convergence study . . . . . . . . . . . . . 2.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Reactive Euler equations with two-step chemistry 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 ZND solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Linear unsteady problem . . . . . . . . . . . . . . . . . . . . . . 3.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 ZND and perturbation solutions . . . . . . . . . . . . . 3.5.2 Migration of spectra with increasing endothermicity 3.5.3 Neutral stability curve . . . . . . . . . . . . . . . . . . . . 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Fickett’s detonation analog 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Equations in laboratory and shock-attached frames 4.2.2 Ambient and Rankine–Hugoniot conditions . . . . . 4.2.3 Shock-evolution equation . . . . . . . . . . . . . . . . 4.3 ZND solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Linear unsteady problem . . . . . . . . . . . . . . . . . . . . . 4.5 Nonlinear simulations by Godunov method . . . . . . . . . 4.5.1 Description of the method . . . . . . . . . . . . . . . . 4.5.2 Riemann problem for the Fickett’s model . . . . . . . 4.5.3 Computation of detonation velocity . . . . . . . . . . 4.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6.1 ZND solutions . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . .

. . . . . . . . . . . . .

. 46 . 49 . . . . . . . .

52 54 57 58 63 64 67 68

. . . . . . . . .

. . . . . . . . .

70 70 72 73 78 79 79 85 87 90

. . . . . . . . . . . . .

91 . 91 . 93 . 93 . 94 . 96 . 97 . 98 . 100 . 100 . 102 . 105 . 109 . 109

9 4.6.2 Demonstration of the stable and unstable solutions . . . . . 4.6.3 Comparison of linear and nonlinear solutions . . . . . . . . . 4.6.4 Migration of linear spectrum as activation energy is varied . 4.6.5 Neutral stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6.6 Comparison with normal-mode analysis . . . . . . . . . . . . 4.6.7 Nonlinear dynamics of the model . . . . . . . . . . . . . . . . 4.6.8 Analysis of stable limit cycles . . . . . . . . . . . . . . . . . . . 4.6.9 Code verification . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 109 . 112 . 112 . 113 . 115 . 117 . 122 . 122 . 129

5 Conclusions 132 5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 5.2 Future Research Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 Appendices

136

A Details of the linearization of the reactive Euler equations

137

B Conversion between different scales used in the detonation stability research 140 C Details of the normal-mode analysis for the Fickett’s detonation analog C.1 Linearization with normal-modes . . . . . . . . . . . . . . . . . . . . . C.2 Formulation of stability problem in reaction variable . . . . . . . . . C.3 Derivation of the boundedness condition . . . . . . . . . . . . . . . . C.4 Formulation of the stability problem . . . . . . . . . . . . . . . . . . .

143 . 143 . 144 . 145 . 147

D Software and computational environment

149

References

151

10

LIST OF ABBREVIATIONS BDF

Backward-Differentiation Formulae

CJ

Chapman–Jouguet

DMD

Dynamic Mode Decomposition

ENO

Essentially NonOscillatory

FFT

Fast Fourier Transform

MUSCL Monotonic Upwind Scheme for Conservation Laws ODE

Ordinary Differential Equation

PDE

Partial Differential Equation

SSA

Singular Spectral Analysis

SVD

Singular Value Decomposition

VODE

Variable-coefficient ODE

WENO

Weighted Essentially NonOscillatory

ZND

Zel’dovich–von Neumann-Döring

11

LIST OF FIGURES

1.1 Schema of the flow in CJ theory. . . . . . . . . . . . . . . . . . . . . . . . 17 1.2 Diagram of possible intersections of the Rayleigh–Michelson line and the equilibrium Hugoniot curve. . . . . . . . . . . . . . . . . . . . . 23 1.3 Schema of the flow in a self-sustained detonation wave according to the ZND theory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.1 Schematics of the numerical grid. . . . . . . . . . . . . . . . . . . . . . . 42 2.2 Comparison of FFT and DMD algorithms for the analysis of time series. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.3 Steady solutions for the one-step Euler model . . . . . . . . . . . . . . 59 2.4 Computed time series of the perturbation of detonation velocity . . 61 2.5 Perturbation profiles for the one-step Euler model . . . . . . . . . . . 62 2.6 Comparison of neutral stability boundary with normal-mode results 64 2.7 Neutral stability for the fundamental mode . . . . . . . . . . . . . . . . 66 2.8 Example time series, in which mode 0 is neutrally stable while mode 1 is unstable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 2.9 Neutral stability curves for for modes 0 and 1 . . . . . . . . . . . . . . 67 3.1 Total heat release versus steady reaction progress variable . . . . . . 77 3.2 ZND subsonic-supersonic profiles . . . . . . . . . . . . . . . . . . . . . 81 3.3 ZND subsonic-subsonic profiles . . . . . . . . . . . . . . . . . . . . . . . 82 3.4 Perturbation profiles for the subsonic-supersonic case . . . . . . . . . 83 3.5 Perturbation profiles for the subsonic-subsonic case . . . . . . . . . . 84

12 3.6 Growth rate and frequency versus endothermicity for subsonicsupersonic case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.7 Growth rate versus endothermicity for subsonic-subsonic case . . . 86 3.8 Neutral stability curves for two-step model . . . . . . . . . . . . . . . . 89 3.9 Thermicity profiles for two-step model . . . . . . . . . . . . . . . . . . 89 4.1 Solution of the Riemann problem, shock wave case . . . . . . . . . . . 106 4.2 Solution of the Riemann problem, rarefaction wave case . . . . . . . 107 4.3 ZND profiles for the Fickett’s model . . . . . . . . . . . . . . . . . . . . 110 4.4 Time series of the perturbation of detonation velocity for Fickett’s model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 4.5 Perturbation profiles for the Fickett’s model . . . . . . . . . . . . . . . 111 4.6 Comparison of solutions of linearized and nonlinear problems . . . 113 4.7 Migration of the linear spectrum for the Fickett’s model . . . . . . . . 114 4.8 Neutral stability for the Fickett’s model . . . . . . . . . . . . . . . . . . 115 4.9 Carpet plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 4.10 ZND solution and perturbations obtained via normal-mode analysis118 4.11 Bifurcation diagrams for the Fickett’s model . . . . . . . . . . . . . . . 121 4.12 Analysis of nonlinear dynamics for θ = 0.95 . . . . . . . . . . . . . . . . 123 4.13 Analysis of nonlinear dynamics for θ = 1 . . . . . . . . . . . . . . . . . . 124 4.14 Analysis of nonlinear dynamics for θ = 1.004 . . . . . . . . . . . . . . . 125 4.15 Analysis of nonlinear dynamics for θ = 1.055 . . . . . . . . . . . . . . . 126 4.16 Analysis of nonlinear dynamics for θ = 1.065 . . . . . . . . . . . . . . . 127 4.17 Analysis of nonlinear dynamics for θ = 1.089 . . . . . . . . . . . . . . . 128

13

LIST OF TABLES

2.1 Results of applying the postprocessing algorithm to synthetic data . 56 2.2 Unstable spectra for γ = 1.2, Q = 50, and E = 35 or E = 40 . . . . . . . . 60 2.3 Comparison of the unstable spectra between the present work and the normal-mode results for θ = 1.2 and Q = 50 . . . . . . . . . . . . . . 65 2.4 Verification of the unsteady linear solver for the one-step Euler model 68 4.1 Critical values of activation energy and the frequency of oscillation on the neutral stability boundary for the Fickett’s model . . . . . . . . 115 4.2 Comparison of the growth rates and frequencies of perturbations obtained via unsteady linear computations and normal-mode computations for the Fickett’s model . . . . . . . . . . . . . . . . . . . . . . . 118 4.3 Bifurcation points and corresponding approximations of the first Feigenbaum constant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 4.4 Convergence study for the linear solver with the Fickett’s model . . . 130 B.1 Scales used in the detonation stability research . . . . . . . . . . . . . 141 B.2 Conversion factors for conversion between different scales . . . . . . 142

14

Chapter 1

Introduction

The goal of the present work is to study the stability of the ZND solutions of detonation models using an approach that combines ideas from linear stability analysis with direct numerical simulations. Detonation models usually consist of a system of conservation laws of fluid mechanics complemented by the relations connecting hydrodynamics to thermodynamics and chemistry. Mathematically, such a system admits a ZND solution—a solution which is a one-dimensional planar wave propagating with a constant velocity. However, such solutions are unstable to small perturbations and almost never observed in physical experiments. Why is it important to study detonation stability? During detonation, the chemical energy of an explosive is released after compression of the explosive by a shock wave, so that the whole process occurs with a significant increase of pressure—typically, 20–60 bars for gaseous mixtures. Deviation of detonation process from ZND regime can drastically change the maximum pressure attained, therefore, it is important to know the stability of the ZND regime in engineering applications, such as construction of chemical factories. Detonation phenomenon gained attention at the end of the nineteenth century after the series of explosions happened in coal mines in France. Mallard & Chatelier (1881) and, independently, Berthelot & Vieille (1881) studied detonations in glass tubes and found that detonation can be distinguished from combustion in that the speed of propagation is supersonic with respect to the initial gaseous mixture.

15 For much time it was believed that the detonation propagates as a plane wave with a constant velocity. Chapman (1899) and Jouguet (1905) developed the first successful theory based on the idea of an instantaneous chemical reaction occurring inside the detonation front. Agreement of the detonation velocity obtained from this theory with the results of the physical experiments of that time determined the longtime success of the Chapman–Jouguet (CJ) theory. Later, in 1940s, Zel’dovich (1940), von Neumann (1942), and Döring (1943), independently of each other, elaborated on the Chapman–Jouguet theory, and, with an assumption of chemical reactions with finite rates, created a theory (which is known now as ZND theory, after its creators) of detonation, in which a detonation wave is represented as a complex consisting of a shock wave and a following reaction zone, such that the length of the reaction zone is constant. The shock wave propagates into the initial gas mixture, compresses and ignites the gas, so that the gas undergoes chemical transformations in the reaction zone. The downstream end of the reaction zone moves at the same speed as the detonation front, so that the thickness of the reaction zone is constant. ZND theory is discussed in detail in Courant & Friedrichs (1948) and Zel’dovich & Kompaneets (1955). However, new physical experiments conducted to confirm the ZND theory found discrepancies between measurements and the theory (White, 1961; Voytsekhovskiy et al., 1963). It was found that a detonation front has a complex three-dimensional structure with new transverse shock waves appearing inside the reaction zone. So, although ZND theory predicts one-dimensional steady wave solutions which are mathematically feasible, these solutions are unstable and generally not observable in practice. The outline of this chapter is the following: in Section 1.1 we give an overview of the Chapman–Jouguet and Zel’dovich–von Neumann–Döring theories; in Section 1.2 we give an overview of previous research on detonation instability; in

16 Section 1.3 a high-level overview of our approach to study detonation instability is given; in Section 1.4 we give an outline of the chapters that follow.

1.1 CJ and ZND theories 1.1.1 CJ theory CJ theory developed by Chapman (1899) and Jouguet (1905) was the first successful theory describing detonation phenomenon. Its success was determined by the good agreement of the theoretical detonation velocity and the velocity observed in the physical experiments. The main properties of this theory are: • a detonation occurs in the motionless ambient mixture at a uniform ther-

modynamic state with no chemical reactions; • a detonation wave is a shock wave propagating with a constant velocity into

the mixture; • detonation velocity is very fast, and the main driver of energy transfer is

convection, so the effects of heat conduction and viscous dissipation are neglected; • the detonation front is a jump discontinuity; • chemical reaction is caused by compression in the detonation front; • the chemical reaction is instantaneous, so that chemical equilibrium is

attained right behind the shock. Figure 1.1 shows the schema of the flow according to the CJ theory—the top part of the figure shows the flow in the laboratory reference frame, and the bottom part shows the flow in the reference frame attached to the detonation front (we

17 Burned mixture

Detonation front Ambient mixture

Expansion of products

Laboratory reference frame

All reactions occur inside the front Burned mixture Ambient mixture Expansion of products

Reference frame attached to the detonation front

Detonation front

Figure 1.1: Schema of the flow in CJ theory. will call it shock-attached frame further). In the laboratory reference frame an ambient mixture—the mixture ahead of the shock—is kept at rest, and no chemical reactions occur inside. Right after the detonation front the burned mixture moves with speed u b . In the shock-attached frame the detonation front remains at the same position for all times, while ambient mixture flows into the front, and burned mixture flows away from the front (see the bottom part of Figure 1.1). The hydrodynamic state right behind the front is defined by Rankine–Hugoniot conditions, which express conservation of fluxes of mass, momentum and total specific enthalpy across the detonation front:

ρ b (u b − D ) = −ρ a D,

(1.1a)

p b + ρ b (u b − D )2 = p a + ρ a D 2 ,

(1.1b)

p b (u b − D )2 pa D 2 + eb + = ea + + , ρb 2 ρa 2

(1.1c)

18 where ρ , u , p , and e denote density, velocity, pressure, and internal energy, respectively; D is the detonation velocity; subscript “a” denotes quantities in the ambient mixture ahead of the detonation wave; subscript “b” denotes quantities in the burned mixture behind the detonation wave. The ambient mixture is defined by known quantities ρ a , p a , and e a . Therefore, we have four unknowns (ρ b , u b , p b , and D ) and only three equations (e b is found using the thermodynamic relation e = e (p, ρ )), thus, the closure condition is required. It appears in this theory as the

hypothesis based on the following observation. The detonation wave is stable, thus, the burned state is isolated from the acoustic waves that can appear behind the wave due to the expansion of reaction products. Therefore, at the burned state local flow velocity relative to the detonation wave must be sonic:

u b − D = −c b ,

where c b is the sound speed in the burned mixture. This condition, coupled with conditions (1.1) and with an explicit equation of state of the calorically perfect gas

eb =

p −Q, (γ − 1) ρ

where γ is the ratio of specific heats and Q the heat release of the reaction, allows us to compute Chapman–Jouguet detonation velocity of steady reaction D = D CJ : s D CJ =

c a2 +

γ2 − 1

and hence, compute the burned state.

2

s Q+

γ2 − 1

2

Q,

19

1.1.2 ZND theory CJ theory was successful in the first decades of the 20th century, because of the agreement between detonation velocity observed in the physical experiments and theoretical minimal possible velocity. However, assumption of the “instantaneous” chemical reaction inside the detonation front was weak under the following consideration: the thickness of the front is approximately several mean free paths, so that the molecules of a gas do not have much possibility to collide; at the same time, chemical reactions need many thousands of collisions of molecules to change the chemical composition significantly, so that chemical equilibrium can be attained only on distances many times longer than the front thickness. This consideration led researchers to extend the CJ theory, and in the 1940s, independently of each other, Zel’dovich (1940), von Neumann (1942), and Döring (1943) developed the new theory, now known as the ZND theory. Main assumptions of this theory are: • detonation wave propagates with a constant velocity into the motionless

ambient mixture at a uniform thermodynamic state with no chemical reactions; • the main driver of energy transfer is convection, while the effects of heat

conduction and viscous dissipation are neglected; • the detonation front is a jump discontinuity; • chemical reaction is caused by compression in the detonation front and the

rate of chemical reaction is finite, such that chemical equilibrium is attained at some distance behind the shock; • the whole process is steady in the reference frame attached to the front.

20 So, in the ZND theory, a detonation wave is considered as a complex consisting of a shock wave, which compresses the gas to the state in which chemical reaction can occur, and a following reaction zone. To continue with a concrete example, we, as previously, consider calorically perfect gas, in which one irreversible unimolecular reaction “A → B” occurs with the ratio of specific heats of A and B being the same; the degree of completeness of the reaction is denoted by the reaction progress variable λ: λ = 0 in the ambient mixture and λ = 1 in the end of the reaction zone where chemical equilibrium is attained. The state right behind the shock (in the point known as the von Neumann point in the ZND theory) is still unreacted and connected to the ambient state by the Rankine–Hugoniot conditions

ρ s (u s − D ) = −ρ a D, p s + ρ s (u s − D )2 = p a + ρ a D 2 , es +

pa D 2 p s (u s − D )2 + = ea + + , ρs 2 ρa 2 λs = 0,

where subscript “s” denotes von Neumann point, and internal energy e s does not contain any chemical energy released. The flow in the reaction zone is described by steady reactive Euler equations written in the shock-attached frame: d(ρ (u − D )) = 0, dx ¡ ¢ d ρ (u − D )2 + p = 0, ³ ³ dx ´´ 2 p d ρ (u − D ) e + ρ + u2 = 0, dx d(ρ (u − D ) λ) = ρω (p, ρ, λ) , dx

(1.2a) (1.2b) (1.2c) (1.2d)

where ω (p, ρ, λ) is the production rate of the product. The first three equations

21 can be easily integrated, while the last one can be simplified by using product rule for differentiation and continuity equation (1.2a). This yields

ρ (u − D ) = const1 ,

(1.3a)

ρ (u − D )2 + p = const2 ,

(1.3b)

e+

u2

p + ρ 2 dλ dx

= const3 , =

ω (p, ρ, λ) , u −D

(1.3c) (1.3d)

which can be integrated, and the flow state can be found once we know the detonation velocity, as it enters the equations as a parameter. Equations (1.3a–1.3b) can be manipulated , which yields the equation of the Rayleigh–Michelson line p − p a = ρ 2a D 2 (v a − v ) ,

(1.4)

where v = 1/ρ is a specific volume. The Rayleigh–Michelson line is a straight line with a slope −ρ 2a D 2 in the (v, p ) plane and depends on the parameter D . Consequently, equations (1.3a–1.3c) can be manipulated to exclude velocity u and obtain the equation for a family of Hugoniot curves

e (p, v, λ) − e a (p, v, 0) =

1 (p + p a )(v a − v ) , 2

(1.5)

which depend on parameter λ. As both equations (1.4) and (1.5) follow from the same conservation laws, any given Rayleigh–Michelson line must intersect with the whole family of Hugoniot curves for every 0 ≤ λ ≤ 1, otherwise the ZND detonation wave is not possible. From the whole family of Hugoniot curves the most important one for intersection is the equilibrium Hugoniot curve, namely, the curve corresponding to λ = 1. Looking at the Figure 1.2, we note that the Rayleigh–Michelson line and equilibrium

22 Hugoniot curve has two possibilities to intersect: 1. D = D CJ . In this case the gas state “jumps” abruptly from point A to point N and then continuously changes along the Rayleigh–Michelson line towards point C , where the Rayleigh–Michelson line and the equilibrium Hugoniot curve are tangential to each other. 2. D > D CJ . In this case the gas state “jumps” from point A to point N and then continuously changes along the Rayleigh–Michelson line towards point S . The significance of the case 1, where D = D CJ , is that the velocity D CJ is the minimal velocity that makes steady self-sustained detonation wave possible, at least theoretically. It can be shown that in the end of the reaction zone products move with a local sound speed c CJ relatively to the detonation front:

u CJ − D CJ = −c CJ ,

(1.6)

which, in terms of the characteristics of the system, means that the velocity u + c of the characteristic emanating from the end of the reaction zone and travelling towards the shock is the same as the lead shock velocity D . The end of the reaction zone is called the Chapman–Jouguet, or sonic, plane (point CJ on Figure 1.2, and the condition (1.6) is called the Chapman–Jouguet condition. The flow inside the reaction zone is subsonic relative to the shock, therefore acoustic waves generated by released chemical energy can reach the detonation front and support its strength. The presence of the sonic plane exactly at the end of reaction zone isolates the whole detonation wave from any rarefaction generated behind. Therefore, these two factors—support of the front by chemical energy and isolation from rarefaction waves—make self-sustained detonation waves possible. In contrast, when D > D CJ , the flow at the end of the reaction zone (point S on Figure (1.2)) is subsonic, therefore, the detonation wave can be destroyed by a

23

Rayleigh-Michelson lines

Shock Hugoniot

Equilibrium Hugoniot

Figure 1.2: Diagram of possible intersections of the Rayleigh–Michelson line and the equilibrium Hugoniot curve. rarefaction wave coming from behind. Figure 1.3 shows the schema of the flow in a self-sustained detonation wave according to the ZND theory. The von Neumann spike is denoted on the figure right behind the shock, where the ambient mixture is in shocked but unreacted state. Subsequently, the chemical reaction progresses in the direction away from the shock (arrow with the beginning λ = 0 and the ending λ = 1). The length of the reaction zone is constant, and the whole flow inside reaction zone is steady, that is, in the reference frame attached to the detonation front all flow variables are functions of only a spatial coordinate, not time. The following flow can be unsteady and depends on the boundary conditions on the left. The situation shown on Figure 1.3 is when a rarefaction wave in the reaction products travels between the left boundary and the reaction zone to eliminate the differences between high pressure in the products and low pressure on the left boundary.

24 Detonation front moves Rarefaction wave in reaction products

Reaction zone

Laboratory reference frame

von Neumann spike right behind the front

Sonic (CJ) plane

Rarefaction wave in reaction products

Ambient mixture

Reaction zone

Ambient mixture

Reference frame attached to the detonation front

Detonation front rests

Figure 1.3: Schema of the flow in a self-sustained detonation wave according to the ZND theory.

1.2 Stability computations As soon as it was discovered experimentally (White, 1961; Voytsekhovskiy et al., 1963) that ZND detonation in gases was not stable, the question of the stability analysis arose. Shchelkin (1959) formulated criterion of the stability of ZND detonations based on a simplifying assumption that the detonation front is followed by a region of constant state coupled with an infinitesimally thin reaction zone (square-wave detonation model). Consequent work on the detonation stability focused mainly on three branches: linear stability analysis, asymptotic analysis, and direct numerical simulations. Below we review some of the most important works of the first and third branches for the one- and two-dimensional detonations. Linear stability analysis. Erpenbeck (1962) formulated the problem of detona-

tion stability as a hydrodynamic stability problem for three-dimensional flow with several reactions and proposed an approach for the solution applying his

25 approach to the flow with the one-step irreversible Arrhenius-type reaction (Erpenbeck, 1964). Erpenbeck’s technique consisted of linearization of the equations about a steady-state solution and application of the Laplace transform in time to reduce the problem to a system of linear inhomogeneous ODEs and then analysis of the poles of the transform function in a complex plane. The extensive overview of the Erpenbeck’s work, as well as of other early studies, is provided in Fickett & Davis (2000). Later, Lee & Stewart (1990) studied the linear stability of one-dimensional detonation with the normal-mode approach by posing the stability problem as a two-point eigenvalue problem with the linearized Rankine–Hugoniot conditions at the shock and with the boundedness condition in the end of the reaction zone. The latter condition (also known as the radiation condition) serves as a dispersion relation and mathematically means that eigenfunctions remain bounded or, from a physical point of view, that no information propagates from the rear end of the reaction zone to the lead shock. To solve the boundary-value problem they used the shooting method coupled with the “carpet search” to find the domain of attraction of the boundedness-condition function. The authors obtained neutral stability curves in different parameter spaces. In the space “heat release vs. overdrive factor” they found qualitative agreement with the results of Erpenbeck (1964). Their main conclusion was that the increase of the degree of overdrive stabilizes the ZND solutions, while the increase of activation energy destabilizes the ZND solution. Later, Sharpe (1997) modified the approach of Lee & Stewart by deriving asymptotic value of flow state at the end of the reaction zone and consequently (using this state as an initial value) integrating towards the shock. As integration in the direction of the shock is better conditioned, it allowed consideration of not only the range of the parameters of Lee & Stewart (where his results agrees with theirs), but also extension of the analysis to the larger activation energy range,

26 where he found pure exponential growth of instability in contrast with the results of Buckmaster & Neves (1988). Short & Stewart (1998) extended normal-mode analysis of (Lee & Stewart, 1990) to study two-dimensional instability. The authors considered the dependence of the stability spectrum on the degree of overdrive and found that the number of unstable modes decreases when the degree of overdrive increases. They also investigated stability of detonation with large heat release and varying activation energy. Finally, they investigated how detonation stability depends on heat release and found that the lower the heat release, the more stable detonation becomes. They also demonstrated that even in the cases when detonation is stable to longitudinal perturbations, it can still be unstable to transverse perturbations. Several works also applied normal-mode analysis to models with several reactions. Short & Dold (1996) considered detonation in a calorically perfect gas with three-step chain-branching kinetics, while Gorchkov et al. (2007) formulated the stability problem for detonations with arbitrary thermodynamic and kinetic models. Sharpe (1999) considered the case of pathological detonation, in which two reactions occur consecutively and the second reaction is endothermic, and found that, due to the embedded sonic point, the stability properties did not depend on the part of the reaction zone behind this point. Kasimov & Stewart (2002); Kasimov (2004) applied normal-mode analysis to detonation in a circular tube of finite diameter with the one-step reaction. Stewart & Kasimov (2005) considered an extension of the theory of sonic locus of ZND detonations to the case of multidimensional unsteady detonations and provided a reformulation of the boundedness condition of the linear stability theory. Tumin (Tumin, 2007a,b,c) employed multidomain spectral collocation methods to revisit Erpenbeck’s method and found an agreement with the normal-mode approach of Lee & Stewart (1990).

27 Numerical simulations. Fickett & Wood (1966) employed the model of Erpen-

beck and studied the long-time behavior of shock pressure using the first-order explicit method of characteristics for several cases of overdriven detonation which were analyzed previously by Erpenbeck. Numerical solutions confirm the result of linear analysis by Erpenbeck. Also, the authors studied the growth of instability by specifying approximately steady piston motion which generated a ZND wave and then allowed the truncation error of the numerical method to perturb the wave. They found that when they halved the grid step, the amplitude of oscillations in the exponential growth stage is smaller, however, when the limit cycle is obtained, the system no longer “remembered” the amplitude of these oscillations, that is, long-time behavior was the same. They also varied the rear boundary condition using different piston motions and found that, again, although initial transient behavior depended on the piston movement, long-time behavior was the same. Fickett et al. (1972) extended their code to simulate simplified two-step branching kinetics. Colella et al. (1986) applied a fractional-step method to two models of detonation: simplified Majda’s model and reactive compressible Navier–Stokes equations and found that, for heat-release values less than the critical value or for a fast reaction rate, solution bifurcated and separated into two distinct waves: a selfsustained (weak) detonation wave in which reaction completed to equilibrium followed by an inert shock wave, and this solution was stable in unsteady computations with viscous effects included such that the distance between two waves was increasing. Furthermore, the authors conducted unsteady simulations of the reactive Euler equations and found that for coarse meshes a similar pattern can be observed in which the precursor moves with a speed of one computational cell per time step. The authors gave a theoretical explanation to this purely numerical phenomenon as it disappeared with increased grid resolution. The

28 authors also emphasized the importance of the simplified models for theoretical and numerical purposes. Bourlioux et al. (1991) developed a numerical method that combines the thirdorder piecewise-parabolic method (Colella & Woodward, 1984) with adaptive mesh refinement and shock tracking and applied it to the Erpenbeck–Fickett–Wood problem (Fickett & Wood, 1966) pointing out that shock-capturing methods smear physical picture of detonation dynamics if used on low-resolution grids. Cai et al. (1996) used spectral methods with shock tracking and by employing Fourier transform plotted a neutral stability curve in the plane “heat release–degree of overdrive”. Xu et al. (1997) developed code based on the fifth-order WENO interpolation (Jiang & Shu, 1996) and applied this to a test problem of detonation initiation. Several researchers investigated the required grid resolution for obtaining convergent solutions for the model with one-step chemistry: Sharpe & Falle (2000); Sharpe (2001); Sharpe & Falle (1999) used a second-order Godunov method with grid adaptation and concluded that they had to use 2.5 times finer resolution than was reported previously (Bourlioux et al., 1991) to obtain a convergent solution for stable and weakly unstable ZND solutions, while Hwang et al. (2000) used the third-order ENO scheme (Harten et al., 1987; Fedkiw et al., 1997) with local Lax–Friedrichs flux and concluded with similar needs for grid resolution. Short et al. (1999) studied the detonation dynamics with the model of three-step chemistry (Short & Dold, 1996; Short, 1997) using the adaptive Godunov method and gave an explanation of the interaction of the lead shock with expansion and internal detonation waves; additionally, their model allowed them to demonstrate detonation quenching, in which the shock decouples from the reaction zone completely. All of the previously described numerical simulations used the methods based on either shock capturing or front tracking to analyze detonation time series.

29 Kasimov & Stewart (2004) and Kasimov (2004) conducted an analysis of detonation dynamics in the reference frame attached to the lead shock, such that detonation velocity was explicitly computed on each time step by the first-order method of characteristics. The authors also demonstrated, by employing a model with finite length of reaction zone, that the sonic surface at the end of the reaction zone acts as an interface between those acoustic waves that can reach the detonation shock and those that cannot. Ng et al. (2005a,b) used the tools of the theory of dynamic systems to construct a bifurcation diagram and point out the possibility of chaos in reactive Euler equations. Henrick et al. (2006) and Henrick (2008) improved the implementation of the shock-attached simulations (Kasimov & Stewart, 2004) by simulating one-dimensional pulsating detonation with fifth-order method and plotted detailed bifurcation diagram for a weakly unstable case. Taylor et al. (2009) and Taylor (2010) conducted shock-attached simulations of two-dimensional weakly unstable detonations. Stewart & Kasimov (2006), Lee (2008), Shepherd (2009), and Ng & Zhang (2012) give extended overviews of the works on detonation instability based on linear stability analysis and direct numerical simulations.

1.3 Overview of the present work In the present work we develop a novel method for linear stability analysis of detonations, namely, the stability of steady travelling-wave solutions (which we refer to as ZND solutions in this work) of different detonation models varying from the models based on the reactive Euler equations with one- and two-step reaction mechanisms to the Fickett’s detonation analog. Here we give an overview of the method. We consider one-dimensional flow of a gas which is ignited by a detonation shock with consequent chemical reaction or reactions. Gas is assumed to be

30 inviscid such that the flow is described by a system of hyperbolic balance laws:

z t + f ( z ) x = s ( z ),

where z ∈ Rn is a system state, f ∈ Rn the flux function, and s ∈ Rn the source term. The Rankine–Hugoniot conditions connect the state ahead of the shock with the postshock state. The method is as follows: 1. Switch to the system moving with the detonation wave and write equations in a nonconservative form:

z t + ( A ( z ) − D I ) z x = s ( z ),

(1.7)

where D is the detonation-wave velocity, A ∈ Rn×n is the Jacobian of f , and I ∈ Rn the identity matrix. Now the Rankine–Hugoniot conditions become

boundary conditions. 2. Compute the ZND solution for the given model by omitting time dependence: ¡

¯I A (z¯) − D

¢ dz¯

dx

= s (z¯),

(1.8)

which involves finding speed D¯ , which is an eigenvalue of (1.8). 3. Linearize the system (1.7) about the ZND solution. Linearize the Rankine–Hugoniot conditions. This results in an linear unsteady system in which unknowns are the perturbations of the ZND solution:

z 0t + A (z¯)z x0 + B (z¯)z 0 −

dz¯ 0 ψ = 0, dx

where “prime” denotes perturbations with ψ0 being the perturbation of detonation velocity D¯ .

31 4. Pose the linear unsteady problem consisting of the linear system along with the linearized Rankine–Hugoniot conditions as a boundary condition on the upstream boundary and outflow boundary condition on the downstream boundary. Set initial condition using functions of magnitude much smaller than the magnitude of ZND quantities. 5. Solve the obtained linear unsteady problem numerically recording the time series of the perturbation of detonation velocity ψ0 . 6. Finally, postprocess the resultant time series to extract the growth rates and frequencies (stability spectrum), assuming that the perturbations can be decomposed into a linear decomposition of complex exponential functions:

ψ0 ( t ) =

N X i =0

¡ ¡ ¢¢ c i exp(αi t ) + exp α∗i t ,

where c i is the amplitude of the i th mode, αi ∈ C, and α∗i is the complexconjugate of αi for i = 0, . . . , N . Correspondingly the real and imaginary parts of αi are the growth rates and the frequencies of the linear perturbations. It is crucial that the last step (postprocessing) is carried out as efficiently as possible. In the early stages of this work we planned to achieve this goal using the nonlinear least-squares (NLS) method. However, in the case when the time series of the perturbation of detonation velocity contains multiple modes growing and decaying at different rates, it becomes difficult to choose an initial guess required by the NLS method: optimization solver can easily stop at a local minimum far from the global minimum for some particular initial guess. Hence, NLS method was discarded as not robust enough for the present work. Further research of the methods applicable for extracting the growth rates and frequencies led us to the Prony method (Hildebrand, 1987), which is a special method (actually predating the Fourier method) for decomposition of the time series into a linear combina-

32 tion of complex-valued exponential functions. However, its application did not yield much success as we were not able to extract multiple exponential functions due to numerical instabilities. Finally, this thesis was saved by discovering that in the fluid dynamic community a new method had been developed to extract stability spectra from the time series recorded in numerical simulation or physical experiments. The method is called dynamic mode decomposition (Schmid, 2010). It is at the core of our postprocessing algorithm and is described in details along with the solution of a linear unsteady problem in Chapter 2. The main novelty of this work is that we combine together linearization of the governing equations, their numerical solution, and postprocessing via dynamic mode decomposition for analysis of detonation stability. The idea of studying stability of hyperbolic systems via numerical solution of the linearized equations, to the best of our knowledge, was first described by Godlewski & Raviart (1999) and was applied in other areas of fluid dynamics, for example, in magnetohydrodynamics (Samtaney, 2009).

1.4 Outline The following chapters consider analyses of three different detonation models. In Chapter 2 we develop our approach and analyze the model based on the reactive Euler equations with one-step chemistry, which is the model that was extensively analyzed by the normal-mode approach by previous researchers. This is the central chapter of this thesis that contains a full detailed description of the numerical algorithms used to solve the linearized problem and postprocess the resultant time series of detonation velocity using the dynamic mode decomposition algorithm. In Chapter 3 we apply our approach to the model based on the reactive Euler equations with two-step chemistry with the second reaction being

33 endothermic. In Chapter 4 we study Fickett’s detonation analog by our method and also study the nonlinear dynamics of the model. Finally, in Chapter 5 we summarize the most significant results achieved in this work. A set of appendices contains details augmenting the main body of the thesis. Appendix A contains the details of the linearization of the reactive Euler equations. Appendix B contains the details of conversion between different scales used by researchers for analysis of detonation stability. Appendix C contains details of the algorithm used for the normal-mode analysis of the Fickett’s model. Appendix D describes the software used to implement the present work.

34

Chapter 2

Reactive Euler equations with one-step chemistry

2.1 Introduction In this chapter we consider linear stability of a detonation model which is based on the one-dimensional reactive Euler equations for a perfect gas that has a constant ratio of specific heats and undergoes single first-order irreversible monomolecular reaction reactant → product.

(2.1)

such that the molar masses of both species are the same. This model is ubiquitously used in research on stability of detonations; it was first used by Erpenbeck (1964) due to its relative simplicity and was studied extensively using normalmode approach (Lee & Stewart (1990); Short & Stewart (1998); Sharpe (1997)) for its linear stability properties. However, one should be aware that despite the simplicity of the model, its analysis via normal modes is complicated: postulation of the normal-mode solution of the governing equations linearized about a base steady-state solution leads to the two-point boundary-value problem for a system of ordinary differential equations (ODEs) with one boundary at the shock state and the other boundary in the equilibrium state at the end of the reaction zone, in which the system of ODEs becomes singular. Furthermore, due to the first-order of the reaction, the reaction zone length is infinite, such that one must truncate the computational domain at some point and formulate an approximate boundary condition at the “practical” end of the reaction zone, which demands

35 that perturbations are bounded (boundedness condition). After that one must integrate the system numerically, as it has variable coefficients due to the nonconstant base steady state, and obtain the eigenvalues of the problem (stability spectrum) with positive real parts indicating the growth of the amplitude of the perturbations with time and hence, instability. The eigenvalues are determined using the shooting method (Press et al., 2007) and there are two options here: start at the shock and integrate towards the chemical-equilibrium boundary ensuring that the boundedness condition is satisfied (an approach taken by Lee & Stewart (1990) and Short & Stewart (1998)), or start at the “practical” equilibrium boundary using asymptotic expansions to obtain the approximate values of the state variables as an initial condition for the shooting method and integrate towards the shock (approach used by Sharpe (1997)) ensuring that the Rankine–Hugoniot conditions are satisfied at the shock. Both approaches are complicated as they require a special treatment of the singularity. Our aim is to develop a new method of linear stability analysis and apply it to this model. One key advantage of the method over the normal-mode approach is that it avoids dealing with singular ODEs—instead, the stability problem is based on the linear nonsingular partial differential equations for the perturbations of the steady-state solutions. We fully explain in this chapter the method including the formulation of the governing equations in the shock-attached frame, formulation of the linear unsteady problem, numerical methods used to solve it, and the method of postprocessing the results of computations to extract, growth rates and frequencies (stability spectra) of the linear perturbations. Particularly, we show that the time series of the perturbation of detonation velocity are sufficient to extract stability spectra without the need to postprocess the full solution of the linear unsteady problem. Numerous stability results obtained by other researchers using the normal-mode approach allow us to verify the method.

36

2.2 Governing equations The one-dimensional reactive Euler equations for this system with one-step chemistry are the usual laws of conservation of mass, momentum, and energy, augmented with the balance law for the reaction progress variable:

ρ t ` + (ρu )x ` = 0, ³ ´ (ρu )t ` + ρu 2 + p ` = 0,

(2.2)

(ρe )t ` + (ρue + pu )x ` = 0,

(2.4)

x

(ρλ)t ` + (ρuλ)x ` = ρω,

(2.3)

(2.5)

where ρ , u , p , e , λ, and ω are density, flow velocity, pressure, total specific energy, reaction progress variable, and reaction rate, respectively. The reaction progress variable λ is defined such that the mass fraction of the reactant is 1 − λ and that of the product is λ (hence, λ varies from zero in the nonreacted state to unity in the fully reacted state). The specific total energy is e = e i + u 2 /2 with constitutive relations of a calorically perfect gas:

ei =

pv −Qλ, γ−1

pv = RT,

(2.6)

where e i is specific internal energy, γ the ratio of specific heats, v = 1/ρ specific volume, Q chemical heat release, R gas constant, T temperature. The reaction rate is taken here to follow simple-depletion first-order Arrhenius law: µ ¶ Eρ ω = k (1 − λ)exp − , p

(2.7)

where k is the rate constant and E the activation energy. To transform Eq. (2.2–2.5) to the shock-attached frame, we assume that in

37 the laboratory reference frame detonation wave moves from left to right with detonation speed D = D (t ` ) > 0 from the initial position x 0` = 0 at time t 0` = 0. Therefore, the transformation to the shock-attached frame is defined by `

`

`

x = x − xs = x −

t`

Z 0

D (τ)dτ,

t = t `,

such that the position of the shock in the shock-attached frame is x s = 0 for any time t ≥ 0. Consequently, by applying the chain rule to express partial derivatives with respect to t ` and x ` in terms of the derivatives with respect to t and x , the governing equations in the shock-attached frame are

ρ t + (ρ (u − D ))x = 0,

(2.8)

(ρu )t + (ρu (u − D ) + p )x = 0,

(2.9)

(ρe )t + (ρ (u − D ) e + pu )x = 0,

(2.10)

(ρλ)t + (ρ (u − D ) λ)x = ρω.

(2.11)

The numerical solution of the governing equations (2.8–2.11) requires an evolution equation for D as it appears explicitly in the equations. For this purpose, we use the equation derived in Taylor et al. (2009) by specializing it to one dimension: dM = s, dt

(2.12)

where M = −ρ a D is the normal mass flux into the shock, ρ a is the upstream density, and

s=

1 A0

(R s − A s ) ,

R s = Q (γ − 1)ρ s ωs ,

A0 =

2 γ+1

µ M va 3 +

γp a ρ a

A s = ρ s (c s2 −Us2 )u x |s ,

M2

¶ ,

(2.13) (2.14)

38 where c =

p

γp /ρ is the sound speed, U = u − D is the flow velocity in the shock-

attached frame, subscript “s” denotes the shock state and subscript “a” the ambient (upstream) state. Note that in Eq. (2.13-2.14) everything depends on the shock state through the Rankine-Hugoniot conditions except for the velocity gradient u x |s at the shock, which must be determined numerically by solving for the flow

downstream. The Rankine–Hugoniot conditions at x = 0 give all the state variables in terms of the mass flux (or, equivalently, the shock speed):

ps = −

γ−1 2v a M 2 pa + , γ+1 γ+1

vs =

γ−1 2 γp a va + , γ+1 γ+1 M2

Us =

M , ρs

λs = 0,

ρs =

1 vs

.

(2.15) Ambient state ahead of the shock is given by density ρ a , pressure p a , and velocity u a = 0, hence, the flow ahead of the shock is motionless. Additionally, we assume that the reaction rate ω = 0 in the ambient state such that the reaction is initiated by the leading shock wave. All equations are nondimensionalized using upstream scales (Erpenbeck scales), the density scale is ρ a , the pressure scale is p a , the velocity scale is

p

p a /ρ a , the

energy scale is RTa , the length scale is the length of the steady-state half-reaction zone, and the characteristic time scale is the length scale divided by the velocity scale. In these scales equations retain their form.

2.3 ZND solution The ZND solution is the steady-state traveling-wave solution of the equations (2.8–2.11) with the detonation wave travelling to the right with constant velocity ¯ . In the following, the steady-state quantities are denoted with the “bar” symbol. D Under the assumption of steadiness, the system (2.8–2.11) reduces to the algebraic-

39 differential problem:

¯ ρ¯U¯ = −ρ a D,

(2.16)

¯ 2, p¯ + ρ¯U¯ 2 = p a + ρ a D e¯ +

¯2 pa D p¯ U¯ 2 + = ea + + , ρ¯ 2 ρa 2 ¯ dλ¯ ω = . dx U¯

(2.17)

The algebraic expressions (2.16–2.17) can be manipulated in different ways such that ρ¯, U¯ , and p¯ depend on λ¯ , and it can be done in many different ways. Here we use the formulation from Kasimov & Stewart (2005). Introducing steady mass flux ¯ = −ρ a D¯ , momentum flux P¯ = p a + ρ a D¯ 2 , and energy flux H¯ = e a + p a /ρ a + D¯ 2 /2, M which are constant in the steady-state wave through the reaction zone, the ZND solution is given by:

v¯ =

γ P¯ (1 − δ), ¯2 γ+1 M

where

s δ=

1−

¯2¡ hM P¯ 2

¯ 2 v, p¯ = P¯ − M ¯

¢ ¯ , H¯ +Q λ

h=

¯ v, U¯ = M ¯

2(γ2 − 1) γ2

(2.18)

(2.19)

,

To obtain the solution in terms of x , the steady-state rate equation 2.17 is integrated starting at the shock with the initial condition λ¯ s = 0 and auxiliary relations (2.18) and (2.19). For self-sustained detonations with the Arrhenius kinetics (2.7), the flow becomes sonic relative to the shock at the end of the reaction zone where λ = 1. It can be shown (Kasimov 2004; Kasimov & Stewart 2005) that then δ = 0 at λ = 1, and it follows from Eq. (2.19) that D¯ = D CJ where

D CJ =

q

p c a2 + q + q,

q=

(γ2 − 1)Q , 2

c a2 =

γp a , ρa

(2.20)

40 and velocity D CJ is called the Chapman–Jouguet velocity, or equilibrium detonation velocity, as the sonic condition is satisfied at the point of chemical equilibrium λ = 1.

The pre-exponential factor k is defined such that the length scale is the steadystate half-reaction zone, that is, the width of the reaction zone in which λ¯ ∈ [0, 1/2]: Z 1/2 k=

0

¡

U¯ ¯ dλ. ¢ 1 − λ¯ exp(−E ρ¯/p¯)

(2.21)

2.4 Linear unsteady problem Next, we introduce the vector of state variables z = (ρ, u, p, λ)T and linearize it in addition to the detonation velocity as

z (x, t ) = z¯(x ) + ²z 0 (x, t ),

¯ + ²ψ0 (t ), D (t ) = D with 0 < ² ¿ 1, the bar denoting the steady-state, and the prime denoting the perturbation amplitudes. The reader will find all details of the linearization in the Appendix A. Upon linearization of the governing equations (2.8–2.11), we obtain

z 0t + A (z¯)z 0x + B (z¯)z 0 −

dz¯ 0 ψ = 0, dx

(2.22)

with  ¯ U   0 A=  0  

0

ρ¯

0



0

  U¯ 1/ρ¯ 0  ,  γp¯ U¯ 0   0 0 U¯



du¯ dx

    − 1 dp¯  ρ¯2 dx B =  C (ρ¯ω ¯ρ +ω ¯)    −ω ¯ρ



dρ¯ dx

0

du¯ dx dp¯ dx dλ¯ dx

0 γ ddux¯

+C ρ¯ω ¯p −ω ¯p

0 

  0   ,  C ρ¯ω ¯ λ    −ω ¯λ

(2.23)

41 where C = −(γ − 1)Q and subscripts ρ , p, and λ denote partial derivatives with respect to the corresponding steady quantities, with these derivatives being E ω ¯ ρ = k 1 − λ¯ exp − , p¯ ¡

¢

µ



ω ¯p =

kE ρ¯ ¡ p¯2

E ρ¯ 1 − λ¯ exp − , p¯ µ

¢



µ ¶ E ρ¯ ω ¯ λ = −k exp − . p¯

(2.24) The linearization of the shock-evolution equation (2.12) yields: dM 0 = s 0, dt

(2.25)

A 00 1 0 0 (R s − A s ) − 2 (R¯s − A¯ s ), s = A¯ 0 A¯

(2.26)

where 0

0

¯ 2 − γp a 0 2 3v a M M, ¯2 γ+1 M

¤¯ £¡ ¢ ρ¯ω ¯ρ +ω ¯ ρ 0 + ρ¯ω ¯ p p 0 ¯s ,

(2.27)

· ´ ¯ ³ ´ ¸¯¯ ¢ ¢ ³ 2 du¯ ¡ ¡ 0 0 0 2 du 0 2 2 ¯ ¯ ¯ A s = ρ¯ γ pv ¯ + v¯ p − 2UU + c¯ − U ρ + ρ¯ c¯ − U u x0 ¯¯ dx dx s

(2.28)

A 00 =

R s0 = Q (γ − 1)

0

with all the spatially-dependent quantities on the right-hand sides of R s0 and A 0s evaluated at the shock. The linearization of the Rankine–Hugoniot condi-

tions (2.15) results in following relations:

p s0 =

¯ 0 4v a M M, γ+1

v s0 = −

4γp a M 0, 3 ¯ (γ + 1) M

¯ v s0 + v¯s M 0 , Us0 = M

v0 ρ 0s = − 2s , v¯s

λ0s = 0. (2.29)

2.5 Numerical algorithms 2.5.1 Computational domain and integration of the ZND solution To determine stability properties of self-sustained detonations for any particular set of physical parameters Q , E , and γ, we solve an initial boundary value problem for (2.22–2.25) numerically to see if its solution grows or decays in time. In each simulation, after an initialization stage consisting of grid partitioning and the

42 Computational domain

Shock Ghost points

Biased stencils

Figure 2.1: Schematics of the numerical grid. evaluation of ZND quantities on the grid, the linear system is integrated in time with simultaneous recording of the time series of the perturbation of detonation velocity, ψ0 , and then the time series are postprocessed to extract a stability spectrum. Each of these stages is described in detail below. With the rate function (2.7), the reaction zone is formally infinite. For numerical simulations, this length must be truncated to some finite value L . We determine L by the closeness of the reaction progress variable of the ZND solution λ¯ to unity at the end of the reaction zone as measured by some tolerance tolλ . Once tolλ is specified, we compute L using L = dI e with I =

Z 1−tolλ 0

U¯ ¯ dλ, ¢ ¯ exp(−ρE k 1−λ ¯ /p¯) ¡

(2.30)

where k is computed from Eq. (2.21). The computational domain is then discretized using a uniform grid {x i }, i = 0, . . . , N with N = N1/2 L ; x 0 = −L and x N = 0 being left and right boundary points, respectively. Here N1/2 denotes the number of grid points per unit length (to remind, the unit length is the distance from the shock to the point where λ = 1/2, hence the index 1/2). The spatial step of the grid is ∆x = 1/N1/2 . Schematically, the grid is shown in Figure 2.1. The steady-state ZND profiles are computed on the grid by integrating ordinary differential equation (2.17) with auxiliary relations (2.18–2.19) via the numerical

43 integrator VODE (Brown et al., 1989) with the BDF method of the fifth order with absolute and relative tolerances set to 10−15 . Also, all other ZND quantities that enter matrices A (z¯) and B (z¯) in (2.23) are computed.

2.5.2 Integration of the linearized system After the initialization stage, an unsteady solution of the linearized system is computed by the method of lines. At each time step, the spatial derivatives are first approximated on the numerical grid thereby converting the system of PDEs (2.22–2.25) to a system of ODEs: dz i0 b(z, ¯ z 0 ), = −L dt ¡ ¢ dM 0 ¯ z0 , = sb z, dt

i = 1, . . . , N − 1,

(2.31) (2.32)

¯ z 0 ) and sb(z, ¯ z 0 ) are approximations of L (z, ¯ z 0 , z 0x ) = A (¯z )z 0x + B (z¯)z 0 − ddxz¯ ψ0 where Lb(z, ¯ z 0 , z 0x ) from (2.25), respectively. In these expressions, the from (2.22) and s (z,

spatial derivatives are evaluated using finite-difference formulas as given below, and subsequently the system is evolved in time over one time step ∆t . The spatial derivatives are calculated in the following four steps. 1. Approximation at points i = 0, . . . , N − 3. We compute left- and right-biased approximations of z x0 based on the upwind method of fifth order: z x0− ¯i ¯

=

−2z i0 −3 + 15z i0 −2 − 60z i0 −1 + 20z i0 + 30z i0 +1 − 3z i0 +2

60∆x 0 0 0 ¯ 3 z − 30 z − 20 z + 60z i0 +1 − 15z i0 +2 + 2z i0 +3 i i −2 i −1 0+ ¯ zx i = , 60∆x

,

(2.33) (2.34)

where subscripts “-” and “+” denote left- and right-biased approximations, respectively. 2. Approximation at i = N − 2. At this grid point, the Taylor expansion gives the

44 fifth-order formula 0 0 0 0 0 0 ¯ − 3z N + 30z N + 20z N − 60z N + 15z N −2z N −1 −2 −3 −4 −5 z x0 ¯N −2 = . 60∆x

(2.35)

3. Approximation at i = N −1. At this point, spatial derivatives are approximated using the fourth-order formula

z x0 ¯N −1 ¯

=

0 0 0 0 0 −z N + 6z N − 18z N + 10z N + 3z N −4 −3 −2 −1

12∆x

(2.36)

.

4. Approximation at i = N . At this point, only u x0 must be approximated because it appears in the shock-evolution equation (2.32). The approximation formula is

u x0 ¯N ¯

=

0 0 0 0 0 0 −12u N + 75u N − 200u N + 300u N − 300u N + 137u N −5 −4 −3 −2 −1

60∆x

.

(2.37)

Expressions (2.35)–(2.37) were first used for detonation simulations in Henrick et al. (2006) in their shock-fitting algorithm. Once the spatial derivatives are approximated over the grid, Lb is found using the global Lax–Friedrichs flux:

0

µ

0

b(z, ¯ z ) = L z, ¯ z, L

0+ z 0− x + zx

2



−α

0− z 0+ x − zx

2

,

(2.38)

where z x0− = z x0+ = z x0 for grid points i = {N − 2, N − 1, N } and α is the largest eigenvalue of A (z¯) in Eq. (2.23) over the numerical grid: ¯ u¯ − D, ¯ u¯ + c¯ − D}| ¯ x=xi . α = max {u¯ − c¯ − D, i =0,...,N

(2.39)

We found that a simple evaluation of L was not sufficient to guarantee numerical stability of simulations, hence the use of the Lax–Friedrichs flux in (2.38) that adds ¯ z 0) numerical diffusion, which stabilizes computations. After evaluation of Lb, sb(z,

45 is found using the approximation of u x |s given by (2.37). ¯ z 0 ) and sb(z, ¯ z 0 ) are evaluated on the whole grid, the system (2.31–2.32) Once Lb(z,

is evolved in time using the adaptive-step time integrator DOPRI5 (Hairer et al., 1993), which is an embedded pair of explicit Runge–Kutta methods of orders 5 and 4 used to estimate a local error. We set both absolute and relative tolerances of the time integrator to 10−14 . Such strict tolerances were chosen to minimize any influence of the errors due to integration in time on the solution. Simulation proceeds to the final time Tf , where typically Tf = 10 is used. While computing the solution, we record the time series of the perturbation of detonation velocity ψ0 (t ) at uniform time intervals ∆t = 0.005. When simulation reaches Tf , the solver computes the ratio of L 2 -norms of ψ0 (t ) for Tf /2 ≤ t ≤ Tf and 0 ≤ t ≤ Tf /2, and if this ratio is smaller than three, then Tf is increased to 100, as we found it necessary to record longer time series when the flow is not sufficiently unstable. The analysis of the resulting time series of the perturbation of detonation velocity allows us to extract growth rates and frequencies of unstable modes by a postprocessing method described below in Subsection 2.5.3. We also need to provide initial conditions for z 0 and M 0 . First, the perturbation of the detonation velocity by a small value A 0 =10−10 is specified. Subsequently, the initial value of M 0 is computed using M 0 = −ρ a A 0 . The Rankine–Hugoniot conditions are then used to find ρ 0s , u s0 , and p s0 from M 0 . Subsequently, initial conditions for perturbations over the whole grid are computed using

0

ρ 0 (x ) =

ρ 0s ρ¯s

ρ¯(x ),

0

u 0 (x ) =

u s0 u¯s

u¯ (x ),

0

p 0 (x ) =

p s0 p¯s

p¯(x ),

¯ ( x ). λ00 (x ) = A 0 λ

(2.40)

As can be seen from these formulas, the initial conditions for perturbations are specified as a multiple of the ZND solution. We found that such an initial condition minimizes the nonnormal transient behavior in the time series of ψ0 at early times.

46 The beginning of the time series is cut off before the post-processing in order to eliminate the initial numerical transients from the signal. Whenever the time integrator evaluates the right-hand side of system (2.31–2.32), boundary conditions are imposed as follows: at the shock (grid point i = N ), boundary conditions are specified using Eq. (2.29); at the downstream (left) end, no boundary condition should be specified for the system, however, boundary treatment is necessary for the finite-difference computations (Laney, 1998). This is computed by the zeroth-order extrapolation:

z i0 = z 00 for i ∈ {−3, −2, −1},

where negative indices are used for the ghost grid points to the left of the computational domain. For the fifth-order upwind method (2.33–2.34) three ghost points are needed. This boundary treatment formally is only first-order accurate and introduces spurious waves reflected from the downstream boundary. However, for sufficiently large computational domains, C + −characteristics (moving with speed ¯ + c¯) are almost vertical near this boundary, and, therefore, these spurious u¯ − D waves do not reach the shock and, hence, do not corrupt the time series of ψ0 , for the values of Tf used in the present work.

2.5.3 Postprocessing via dynamic mode decomposition In the postprocessing step, our goal is to extract growth rates and frequencies of unstable modes (discrete stability spectrum) from the time series of detonation velocity £ ¤T ψ0 (t ) = ψ00 , ψ01 , . . . , ψ0n ,

(2.41)

where n is the index of the final time step and values ψ0k are defined at times t k , k = 0, . . . , n such that t k − t k−1 = ∆t = const for all k = 1, . . . , n .

47 To process such a time series, an efficient and robust algorithm is required so that the stability spectrum is computed as accurately as possible. At the heart of our postprocessing algorithm is the dynamic mode decomposition (DMD), which was developed in recent years in the fluid-dynamics community for the analysis of the flows in numerical simulations and physical experiments. The DMD is a powerful method that is able to extract spectral information from complex data coming from different fluid-dynamic systems. It was introduced in Schmid & Sesterhenn (2008), Rowley et al. (2009), and Schmid (2010). The basic introduction to DMD in addition to several examples of its use in different fields can be found in Kutz et al. (2016a). In contrast with other methods for extraction of lowdimensional information from high-dimensional data, DMD modes, although not orthogonal like, for example, modes in proper orthogonal decomposition (Rowley & Dawson, 2017), carry information about a single frequency, similar to normal modes. One of the important features of the method is that, although DMD was developed and used originally by fluid dynamicists, the method is data-driven and is not coupled to the equations governing fluid flows. As such, it finds applications in many other fields besides fluid mechanics and can be used in an equation-free manner. Modern applications of this method include numerous fields that rely heavily on data analysis such as atmospheric research (Kutz et al., 2016b), neuroscience (Brunton et al., 2016), motion detection (Erichson & Donovan, 2016), biology (Proctor et al., 2016; Kunert et al., 2017), computational and experimental fluid mechanics (Schmid, 2011; Schmid et al., 2011; Tissot et al., 2014; Noack et al., 2016; Stankiewicz et al., 2016), etc. One can find also an application of DMD to detonation analysis (Massa et al., 2012), although in contrast with the present results, they did the analysis of nonlinear time series. There were also many attempts to analyze and extend DMD. Chen et al. (2012) provided a comparison between DMD and Fourier analysis and showed that

48 DMD of any data is unique. In Wynn et al. (2013), the authors generalized DMD to optimal mode decomposition and showed that for a noise-free linear system DMD correctly identifies system eigenvalues and, hence, is analogous to normal-mode analysis. Tu et al. (2014) gave a relatively simple description of the algorithm and showed that input data for DMD do not have to be sequential time series with a uniform time step between data snapshots. They also pointed out that snapshots must be sufficiently “rich” for the DMD to extract dynamic modes correctly. In Jovanovi´c et al. (2014), a sparsity-promoting version of the DMD algorithm was proposed to achieve a balance between the quality of the extracted dynamics and the number of modes required to represent the dynamics. Hemati et al. (2014) proposed a method to apply DMD “on-the-fly” for large sets of streaming data. Sayadi et al. (2015) developed a modification of DMD for the parameter-driven bifurcation analysis of thermoacoustically unstable systems, while Brunton et al. (2017) developed a variation of DMD algorithm to characterize chaotic dynamical systems. Several researchers considered sensitivity of DMD to the noise in input data. Duke et al. (2011) considered synthetic wave shapes and carried out a detailed study of errors in growth rates and frequencies under varying parameters. Bagheri (2014) demonstrated that output DMD eigenvalues are damped in the presence of noise, while Dawson et al. (2016) proposed four modifications for the DMD algorithm that alleviate noise effect. In Kou & Zhang (2017) and Zhang et al. (2017), criteria for choosing the rank of a low-order approximation are proposed, with Kou & Zhang (2017) suggesting to determine relevant modes using time-dependent amplitudes, while Zhang et al. (2017) were proposing to choose the mode eigenfunctions on the basis of how good they approximate the eigenfunctions of the corresponding Koopman operator (Koopman, 1931; Budiši´c et al., 2012; Mezi´c, 2013; Williams et al., 2015), which is an infinite-dimensional linear operator that

49 can be posed for a nonlinear dynamical system and whose eigenspectrum serves as an analog of normal modes of linear systems and, hence, is similar to the DMD spectrum. (Mezi´c, 2005; Rowley et al., 2009). We note however, that the DMD algorithm does not necessarily answer the question of what the rank of low-dimensional approximation should be, which to date is an unsolved problem (Zhang et al., 2017), especially when input data are contaminated with noise. Hence, we need a way to determine the best approximation with the ultimate goal being to find all modes that can also be found by normal-mode analysis. The remainder of this section proceeds as follows: first, for completeness, we give a brief generic description of DMD; then, we describe how we form input data for the DMD algorithm from the time series (2.41), and explain how we determine the best low-rank approximation; finally, we apply the algorithm to synthetic data to assess its performance.

2.5.3.1 Description of the DMD algorithm Suppose that we have a time series of the snapshots of the state of some dynamical system: [ x 0 , . . . , x n ]T ,

(2.42)

where x i ∈ Rm for all i = 0, . . . , n , and x i is a snapshot of the system state. Assume also that snapshots are taken at uniform time intervals ∆t . Our purpose is to find an eigendecomposition of a linear operator A such that

x k = Ax k−1

for all k = 1, . . . , n at least in the least-squares sense. The eigenvalues of A contain the information on the growth/decay rates and frequencies of oscillations that,

50 when combined, represent the time evolution of the dynamical system. To proceed, we define matrices X = [x 0 , . . . , x n−1 ]T ,

Y = [ x 1 , . . . , x n ]T

(2.43)

such that Y = AX and thus A = Y X † , where X † is a Moore–Penrose pseudoinverse of X . Although the problem is formally solved now, we cannot simply perform the eigendecomposition of A because A ∈ Rm×m with m generally being large. Instead, ˜ and its eigendecomposition are found. With a low-rank approximation matrix A given inputs (2.43), the DMD algorithm proceeds as follows: 1. Compute reduced singular value decomposition (SVD) of X : X = U ΣV T ,

(2.44)

with matrices U and V having orthonormal columns and Σ being diagonal with singular values σi , i = 1, . . . , min(m, n ), on the diagonal. 2. Truncate U , Σ, and V to rank r , which yields matrices U r , Σr , and V r . Truncation is necessary, for example, if all singular values σi in Σ with i > r are due to the noise in matrix X . Later on we explain how r is chosen in this work. ˜ ∈ Rr ×r : 3. Compute auxiliary matrix A ˜ = U Tr AU r = U Tr Y X †U r = U Tr Y V r Σ−r 1U Tr U r = U Tr Y V r Σ−r 1 . A

˜: 4. Find eigendecomposition of A

˜ W = W Λ, A

51 where Λ = diag(λ1 , . . . , λr ) is the matrix of eigenvalues and W the matrix of eigenvectors. ˜ are the eigen5. Find truncated eigendecomposition of A . All eigenvalues of A ˜ ∈ Rr ×r while values of A as well (the opposite is not necessarily true as A A ∈ Rm×m with r ≤ m ). The matrix of eigenvectors Φ of A can be found using

1 −1 Φ = Y V r Σ− r WΛ ,

with Φ being a rectangular matrix containing exact DMD modes as introduced in Tu et al. (2014). Knowing the eigenvalues Λ and eigenvectors Φ of A , we can reconstruct the time series (2.42): b k = Ax k−1 = A k x 0 = ΦΛk b, x

(2.45)

where b = Φ† x 0 is the vector of the initial amplitudes of DMD modes. Hence, DMD provides modes which are separated in terms of frequencies (that is, imaginary parts of λi ). To measure the quality of low-rank approximation, the relative error e DMD and residual error e resid are computed:

e DMD =

kb x − xk2 , kxk2

e resid = kY − ΦΛΦ† X k2 .

(2.46)

We also convert discrete-time eigenvalues λi that imply stability when they are inside the unit disc in C to the continuous-time eigenvalues αi that imply stability when they are in the left half-plane of C as they are commonly used in normal-mode analysis. The conversion formula is

αi =

log(λi ) ∆t

,

i = 1, . . . , r.

52

2.5.3.2 Application of DMD to the time series of detonation velocity The DMD, as introduced in Schmid (2010), assumes that each snapshot in a time series is a high-dimensional vector coming, for example, from a solution variable on a numerical grid at a given time step in simulations or from a large number of flow sensors in experiments. The time series of perturbation of detonation velocity ψ0 (2.41) that we record during our simulations are one-dimensional, that is, ψ0k ∈ R, k = 0, . . . , n . As pointed out in Tu et al. (2014, pp. 401–402), a one˜ which is 1 × 1, therefore it has only dimensional time series yields the matrix A one real eigenvalue, which can capture exponential growth or decay, but not oscillations. To overcome this rank-deficiency problem, we build an input Hankel matrix Z :



 0 0 ψ · · · ψ n−L   0    0  0 ψ · · · ψ  1 n−L+1   , Z =  . . . . . .  . . .      ψ0L−1 · · · ψ0n

(2.47)

which stacks time-shifted values of ψ0 . This approach is commonly used in other methods of time-series analysis such as Singular Spectral Analysis (Golyandina & Zhigljavsky 2013). Number L in (2.47) defines the upper bound on the number of modes that we can extract from the time series. In this work, L = 1000 is used. The input matrices X and Y (2.43) are built from Z by excluding the last and first columns of Z , respectively. As the time series of ψ0 inevitably contains noise induced by truncation and rounding errors of the numerical computations, truncation of matrices U , Σ, and ˜ ∈ Rr ×r whose eigenvalues V from Eq. (2.44) with a target rank r will yield a matrix A are all physical only if r is known a priori. However, if the rank choice is performed by a machine, there is always a possibility that the value of r will be larger than the rank of the ideal low-rank approximation, hence, spurious eigenvalues will

53 be present in the output of DMD. The other possibility is that r is smaller than the ideal rank, then, instead of a complex-conjugate pair of eigenvalues only one real eigenvalue could be found, which means, for a real-valued input, that instead of having oscillatory exponential mode, DMD finds a purely exponentially growing or decaying mode. As we aim to apply DMD in Section 2.6 for parametric studies of detonation stability, an algorithm is needed for choosing the target rank r automatically and as accurately as possible. The details of the algorithm that we

use are provided below. First, the beginning of the time series (2.41) must be cut off as it is necessary to exclude determining DMD modes from nonnormal behavior occurring due to the relaxation of the perturbation solution from volatile initial conditions (2.40) to a superposition of modes. This cutoff time depends on the length of the time series. Our simulations run by default to a final time Tf = 10, however, the solver can automatically increase Tf to 100 as described in Section 2.5.2. If the final time of the time series is larger than 10, which occurs when the solution is not strongly unstable, then the cutoff time is 10. This time was chosen empirically during experiments on obtaining the neutral stability curve for γ = 1.1 (see Subsection 2.6.3). If the final time is 10 or less, which happens for strongly unstable flows, then the cutoff time is chosen to be unity, also empirically, to exclude initial nonnormal behavior of the time series but to preserve less energetic modes. Once the cutoff time is chosen, the input matrix (2.47) is built and separated into X and Y . After that, candidates for a possible rank of X are determined by the

following procedure: we consider all singular values σk , k = 1, . . . , K of X normalized by σ1 with K being the smallest index corresponding to σK +1 < 10−10 ; then if 0.95 ≤ σk+1 /σk ≤ 1, then it is a strong indicator that σk and σk+1 correspond to the complex-conjugate pair of eigenvalues. Otherwise, if σk+1 /σk < 0.95 for some k , then k is added to the list of possible ranks as there is a sufficient gap between

54 σk and σk+1 , which indicates that σk+1 might be a spurious singular value due to

noise in the input data; for each rank from the list of all possible ranks we determine DMD decomposition and compute the fit and residual errors (2.46); then we consider two smallest fit errors e DMD,1 and e DMD,2 , and if 0.5 ≤ e DMD,1 /e DMD,2 ≤ 1, then we also take into account corresponding residual errors, so that the DMD decomposition corresponding to the smallest residual error is considered the best; otherwise, the best one corresponds to the decomposition with the fit error e DMD,1 . Once the best-rank DMD is determined, DMD eigenvalues that have the

real part larger than −1 and positive imaginary part are saved as an output of the postprocessing algorithm sorted by the increasing frequency and real part, so that the enumeration of modes is the same as in the normal-mode analysis. Saving eigenvalues with real parts between −1 and 0 is necessary to find neutral stability curves, particularly, in the cases when the fundamental mode is less unstable than other modes. For an example of a case like this, see Subsection 2.6.3.

2.5.3.3 Testing the method with synthetic data We apply the postprocessing algorithm to two synthetic time series chosen to represent typical unstable behavior of the solutions of the system (2.31–2.32). Example 1. The time series has a single exponentially growing and oscillating

mode, y 1 (t ) = A exp(αre t )sin(αim t ),

(2.48)

where α = αre + iαim = 3 + 2i, i being imaginary unity, t ∈ [0;51]. Example 2. The time series is a superposition of five exponentially growing and

oscillating modes of different growth rates and frequencies,

y 2 (t ) =

4 X i =0

A exp(αre,i t )sin(αim,i t ),

(2.49)

55 where α ∈ {0.7 + 0.1i, 0.8 + 1.57i, 0.6 + 2.76i, 0.5 + 3.88i, 0.01 + 15.62i}, t ∈ [0;21]. For both time series the amplitude of each mode is taken to be 10−10 and the time step is ∆t = 0.01. These time series are intentionally contaminated with a Gaussian noise of amplitude 10−13 , mean zero, and variance unity such that the noise is proportional to the signal y = y clean (1 + noise),

(2.50)

where y clean is a noiseless time series (either (2.48) or (2.49)). The added noise is to mimic the noise due to numerical simulation errors, where we assume that a simulation error is dominated by truncation error while round-off error is small. To assess the performance of the postprocessing algorithm, we measure relative errors of extracted growth rates and frequencies. For example, if one particular mode in a time series has an exact eigenvalue α while the postprocessing algorithm finds a corresponding approximate eigenvalue αˆ , then the errors are defined as e re =

|αre − α ˆ re | , αre

e im =

|αim − α ˆ im | . αim

(2.51)

In Table 2.1 we show the obtained modes along with the relative errors. Note that the algorithm was able to recover all the modes present in the time series. We find that in general e re is at least one order of magnitude larger than e im . Nevertheless, in the Example 2, the algorithm was able to recover the slowest growing mode (with growth rate 0.01), albeit with a larger relative error than for the other modes, even though the growth rate of this mode is 80 times smaller that the growth rate of the dominant mode α = 0.8 + 1.57i. These examples demonstrate that the postprocessing algorithm we use is able to extract stability spectra from time series consisting of several modes of very different growth rates and frequencies. That said, it must be clear that there are also limitations to the procedure: the quality of the extracted spectrum in general

56 Table 2.1: Results of applying the postprocessing algorithm to synthetic data (2.48) (Example 1) and (2.49) (Example 2). Columns 2 and 3 show found growth rates and frequencies, respectively; columns 4 and 5 show the relative errors of growth rates and frequencies, respectively, computed using (2.51). Example

α ˆ re

α ˆ im

e re

e im

1

0.30

0.20

5.37e-14

1.39e-12

2

0.70 0.80 0.60 0.50 0.01

0.10 1.57 2.76 3.88 15.62

1.16e-05 1.51e-08 3.93e-07 2.24e-07 4.30e-06

1.87e-05 1.94e-08 8.08e-08 5.12e-08 1.90e-09

depends on the quality of the time series at hand (its sampling rate, length, noise level, etc.). As the Fourier transform and power spectra are frequently used in analysis of oscillatory dynamics, here we contrast this traditional method with the present approach. Applying the fast Fourier transform (FFT) to the time series, one can determine approximate values of frequencies of the modes via the power spectrum and then use the nonlinear least-squares solver to find growth rates of the modes. However, note that for exponentially growing or decaying oscillations, FFT analysis may be too inaccurate. To illustrate, Figure 2.2 shows the time series (2.49) and its power spectra generated via FFT and DMD. While the DMD spectrum is sparse and is in excellent agreement with the time series, the FFT spectrum is essentially flat and dense in the sense that it is hard to find the particular modes present in the time series. The algorithm for extracting growth rates and frequencies of unstable modes presented in this section requires multiple application of DMD to find the best low-rank approximation. However, the most expensive operation in the algorithm (SVD of X ) is performed only once. Our postprocessing code is based on the Intel Math Kernel Library, which contains optimized routines for numerical linear

57 10−1

0.0025 0.0020

Amplitude

y2

0.0015 0.0010 0.0005 0.0000 −0.0005

10−3

(a)

(b)

10−5 10−7 10−9

10−11 0

5

10 t

15

20

0

2

4

6 8 10 12 14 16 Frequency

Figure 2.2: Comparison of FFT and DMD algorithms for the analysis of time series (2.49): (a) the time series; (b) the power spectra by FFT (solid line; here the frequencies are multiplied by 2π) and DMD (markers). algebra. The time that it takes to postprocess our data varies as it is a function of the length of the time series and the number of possible ranks. The typical time on a workstation with a Xeon processor is in the order of a minute. A reader familiar with the DMD may question why we apply the algorithm to the long time series (2.41) while we could save the full solution on a grid at each time step and use much shorter time series with tall and skinny input matrices X and Y . However, doing so does not alleviate the problems of choosing the best low-

rank approximation and deciding on the appropriate time length of the input data to properly capture low-frequency modes. Besides, as we process only onedimensional time series (2.41), without having to save the full solution on the numerical grid for each time step, there is a substantial saving in terms of disk space and corresponding input-output operations.

2.6 Results The stability of the steady-state self-sustained solutions for the model (2.2–2.5) depends on the values of E , Q , and γ (Erpenbeck 1964; Lee & Stewart 1990). Generally, increasing E has the effect of making the solution more unstable. The role of Q and γ is more complicated. For example, both very large and very small

58 Q correspond to stable detonations while there is an instability at intermediate

values of Q if E is large enough.

2.6.1 ZND and perturbation solutions The steady-state traveling-wave solutions of (2.8-2.11) are found by setting the time derivatives in the system to zero. In figure 2.3, we show such ZND solution profiles computed with fixed parameters γ = 1.2 and E = 30, and for Q = {0.5, 1, 10, 50, 100}. For this particular value of E , the case Q = 10 corresponds to

a strongly unstable solution, while all other values of Q correspond to stable or weakly unstable solutions. Note from the steady-state profiles that the strongly unstable case occurs for the ZND solution approaching a square-wave shape with a spiky heat-release zone, while all the other cases correspond to relatively smooth profiles. This effect of the steady-state solutions on their stability is well known. The limiting case of a true square-wave detonation with instantaneous heat release is known to possess a pathological spectrum with an infinite number of unstable eigenvalues with growth rates increasing unboundedly with the mode number (or its frequency), i. e., the square-wave stability problem is illposed. However, with a finite-rate chemistry, the problem remains well-posed even though the sharper the heat-release rate the larger the number of unstable modes and the larger their growth rates and frequencies. Computing the latter situations becomes a challenging problem. Now we analyze the time series of the perturbation of detonation velocity, ψ0 , for two cases that differ by the complexity of their unstable spectra. The parameters for these cases are γ = 1.2, Q = 50, and E = 35 or E = 40. Linear stability analysis (Sharpe 1997) predicts that both of these cases have several unstable modes. In addition, the fundamental mode for the second case is purely exponential and has two branches (namely, there are two eigenvalues with zero imaginary parts).

59 1.0

0.8

0.8

0.6

0.6 λ̄

ρ/̄ ρ̄s

1.0

0.4

0.4

0.2

0.2

0.0 −8

−7

−6

−5

−4 x

−3

−2

−1

0

0.0 −8 1.0

0.9

0.9

0.8

0.8

0.7

−5

−4 x

−3

−2

−1

0

−7

−6

−5

−4 x

−3

−2

−1

0

−7

−6

−5

−4 x

−3

−2

−1

0

0.7

0.6 0.5 −8

−6

u/̄ ūs

p/̄ p̄s

1.0

−7

0.6 −7

−6

−5

−4 x

−3

−2

−1

0

0.5 −8 1.0

2.5

0.8 0.6

T/̄ Ts̄

ω/̄ ω̄ max

2.0

0.4

1.5

1.0 −8

0.2 −7

−6

−5

−4 x

−3

−2

−1

0

0.0 −8

Figure 2.3: The steady solution with fixed γ = 1.2, E = 30, and various values of heat release: Q = 0.5 - solid, Q = 1.0 - dashed, Q = 10 - dotted, Q = 50 - solid with circles, Q = 100 - solid with stars. Density, pressure, velocity in the laboratory frame, and temperature are normalized by their respective values at the von Neumann (shock) state; the reaction rate is normalized by its maximum value.

60 Table 2.2: Unstable spectra for γ = 1.2, Q = 50, and E = 35 or E = 40 obtained from simulations with N1/2 = 1280. Here, i is the mode number; αre and αim are the growth rates and frequencies, respectively; e re and e im are the corresponding relative errors computed from the results with N1/2 = 640 and N1/2 = 1280. αre

i E = 35

0 1 E = 40 0 .0 0 .1

1 2 3

e re

αim

e im

0.42323 0.50528

3 × 10−11 1 × 10−11

0.17957 4.36602

2 × 10−10 1 × 10−13

0.18851 1.00847 1.01845 0.67850 0.18504

3 × 10−3 5 × 10−12 2 × 10−12 8 × 10−9 4 × 10−3

0.00000 0.00000 4.30457 7.90388 11.43653

2 × 10−16 2 × 10−16 3 × 10−13 3 × 10−10 6 × 10−5

Simulations are conducted with Tf = 50 for E = 35 and Tf = 40 for E = 40. Figure 2.4 displays the computed time series, with insets showing the solution at early times. The plots suggest that for E = 35 all modes are oscillatory, while for E = 40 an exponential nonoscillatory growth is present. Table 2.2 shows the unstable spectra for both of these cases obtained with N1/2 = 1280 as well as relative errors calculated by comparing the modes found with N1/2 = 640 and N1/2 = 1280. For the case E = 40, the number of modes found is the same as shown in figure 6 of Sharpe (1997). The fundamental mode is purely exponential and is split into two branches in agreement with Sharpe (1997), where the author notes that the split occurs for E ≥ 35.86. The growth rate of the upper branch of the fundamental mode is 1.008 in our calculations while it is 1.021 in Sharpe (1997) (extracted from figure 6 of Sharpe (1997) and rescaled to the Erpenbeck scales). During the integration of the linearized system (2.22–2.25) we obtain the spatial profiles of the perturbations that are superpositions of the eigenfunctions of the problem. Such profiles are shown in Figure 2.5 sampled at time t = 10 at E = 30 and at three values of Q , two of which are near the neutral curve and one somewhere

61 0.8 ψ′

0.4

6e-04 ψ′

0.6

(a)

2e-04 -2e-04

0.2

5

0.0 0

10

2000

ψ′

25

t

20

35

30

t

40

4e-06 ψ′

1500 1000

0

50 (b)

2e-06 0e+00

500 0

15

0

5

5 t 10

10 15 t

20

25

30

Figure 2.4: Computed time series of the perturbation of detonation velocity ψ0 for γ = 1.2, Q = 50, and (a) E = 35 or (b) E = 40. Insets show the evolution of ψ0 at early time of simulation. in the middle of the unstable domain. These are Q = 1 (stable), Q = 10 (strongly unstable), and Q = 50 (weakly unstable). Formulas to compute the perturbations of temperature and reaction rate are:

T 0 = T¯ρ ρ 0 + T¯p p 0 , ω0 = ω ¯ ρ ρ0 + ω ¯ p p0 + ω ¯ λ λ0 ,

where partial derivatives with respect to ρ , p , and λ are p¯ T¯ρ = − 2 , ρ¯

T¯p =

1 , ρ¯

and ω ¯ ρ, ω ¯ p , and ω ¯ λ are given by formulas (2.24).

62 0.8

1e−9

0.6

1.5

(a)

1.0

Q=1

0.4

0.0

0.0

−0.5

−0.2

−1.0

−0.4 −0.6 −5

0.1

−4

−3

−2

−1

0

(c)

−1.5 −5 2.0 1.5

−4

−3

−2

−1

0

(d)

1.0

0.0

0.5

Q = 10

0.0

−0.1

−0.5

−0.2

−1.0 −1.5

−0.3

−2.0

−0.4 −2.0 −1.5 1e−8 0.5 (e)

−1.0

−0.5

0.0

0.0

−2.5 −2.0 −1.5 1e−8 0.2 (f) 0.0

−1.0

−0.5

0.0

−0.2

Q = 50

−0.5

−0.4

−1.0

−0.6

−1.5 −2.0 −5

(b)

0.5

0.2

0.2

1e−9

−0.8 −4

−3

x

−2

−1

0

−1.0 −5

−4

−3

x

−2

−1

0

Figure 2.5: Perturbation profiles versus x for E = 30 and (a-b) Q = 1, (c-d) Q = 10, (e-f) Q = 50. Left column: ρ 0 (solid), u 0 (dashed), p 0 (dash-dotted); right column: λ0 (solid), T 0 (dashed), ω0 (dash-dotted). These profiles are plotted at t = 10.

63

2.6.2 Comparison with normal-mode results Here, we show a quantitative comparison with the normal-mode results found in the literature. As other researchers provide explicit values for unstable spectra for parameters γ = 1.2, Q = 50, we also take these parameters and vary the activation energy E . Table 2.3 shows the results of the comparison. Our results are in agreement with the normal-mode results to at least two significant digits for all considered cases. The normal-mode results are borrowed from the work of different researchers: Sharpe (1997) for E = 25.26 and E = 31.05, from Henrick et al. (2006) for E = 26, and from Lee & Stewart (1990) for E = 50. These researchers all used different scales, which are converted here to the same scales for comparison. We choose to convert to the Erpenbeck scales used in the present work. The conversion factors are as follows. Conversion factor from the scales of Sharpe (1997) to the Erpenbeck scales is nondimensional D¯ in the Erpenbeck scales, which is 6.809475 for γ = 1.2 and Q = 50. In Lee & Stewart (1990), the stability spectrum for E = 50 is given in scales of Abouseif & Toong (1982), and the following conversion

factor should be used: k , ¢−1 1 − λ¯ exp(ρE ¯ /p¯)dλ¯

R 1/2 ¡ 0

where k is the nondimensional Arrhenius rate constant computed using (2.21), and the integral in the denominator is the nondimensional Arrhenius rate constant in the scales of Abouseif & Toong (1982). For γ = 1.2, Q = 50, and E = 50, the numerical value of this factor is 0.938749. Conversion of the results of Henrick et al. (2006) is not needed because they provide the stability spectrum in the same scales as ours. The details of the conversion between different scales are provided in Appendix B. To obtain the results shown in Table 2.3, we run the simulations to final time Tf = 30 for E = 25.26 and E = 26 and to Tf = 10 for E = 31.05 and E = 50. The simu-

64

Heat release, Q

102 101 100

10−1

0

10

20 30 40 Activation energy, E

50

Figure 2.6: The neutral stability boundary computed with the present method (solid line) and its comparison with the normal-mode result of Lee and Stewart Lee & Stewart (1990, figure 7) (square markers). lations are conducted with resolutions N1/2 = 20, . . . , 1280 with refinement factor 2 to verify the convergence of the obtained stability spectra, which is observed as expected. Table 2.3 shows the spectra obtained with N1/2 = 1280. The errors on growth rates and frequencies are also shown in Table 2.3. They are computed using: ¯ ¯ ¯ αre,1 − αre,2 ¯ ¯, ¯ e re = ¯ ¯ αre,2

¯ ¯ ¯ αim,1 − αim,2 ¯ ¯, ¯ e im = ¯ ¯ αim,2

where subscripts 1 and 2 denote quantities obtained with N1/2 = 640 and N1/2 = 1280, respectively.

2.6.3 Neutral stability curves Next, we use the present method to generate neutral stability curves in (E ,Q ) plane for various γ. The corresponding algorithm is as follows: assuming that instability monotonically increases as E increases, the idea of bisection is used to find a critical value of E for each given Q such that the growth rate of a given mode is close to zero with tolerance 10−4 . The initial interval of search for E is taken [10;140] for all results shown in this section. Grid resolution used here is N1/2 = 40.

65

Table 2.3: Comparison of the unstable spectra between the present work and the normal-mode results for γ = 1.2 and Q = 50. Sources for the normal-mode results: ∗ Sharpe (1997, p. 2617); ∗∗ Henrick et al. (2006, p. 321); ∗∗∗ Lee & Stewart (1990, p. 127) except mode 0.0 which is extracted from (Short & Stewart, 1998, figure 6). Here, i is the mode number. For E = 50 the fundamental mode has two branches, which we denote with 0.0 and 0.1. Eigenvalues α are in the scales used in the present work. Notice that the errors are computed using grid refinement. Present work i

αre

e re

αim

Normal-mode analysis e im

αre

αim

E = 25.26∗ 0 −0.00017

1 × 10−4

0.53048

6 × 10−8

0.000

0.530

E = 26.0∗∗ 0 0.03709

5 × 10−9

0.52215

5 × 10−11

0.03710

0.52215

E = 31.05∗ 0 0.26756 1 −0.00060

1 × 10−10 2 × 10−8

0.40280 4.37774

3 × 10−10 1 × 10−12

N/A 0.00

N/A 4.38

E = 50.0∗∗∗ 0.0 0.09365 0.1 1.74458 1 1.76536 2 1.77399 3 1.67745 4 1.53541 5 1.34788 6 1.14096 7 0.91468 8 0.67788 9 0.43092 10 0.17764

3 × 10−7 2 × 10−11 6 × 10−11 6 × 10−12 1 × 10−9 1 × 10−8 4 × 10−8 1 × 10−7 5 × 10−7 1 × 10−6 4 × 10−6 2 × 10−5

0.00000 0.00000 4.10817 7.83005 11.42106 14.99313 18.55077 22.10524 25.65854 29.21139 32.76403 36.31692

2 × 10−16 2 × 10−16 1 × 10−12 2 × 10−11 1 × 10−10 4 × 10−10 6 × 10−10 7 × 10−10 2 × 10−10 2 × 10−9 5 × 10−9 1 × 10−8

N/A 1.743 1.764 1.772 1.676 1.534 1.346 1.140 0.913 0.677 0.431 0.177

N/A 0.000 4.104 7.823 11.42 14.98 18.53 22.09 25.64 29.19 32.73 36.28

66

Heat release, Q

(a)

101

102

(b)

Heat release, Q

102

γ = 1. 1 γ = 1. 2 γ = 1. 3 γ = 1. 4

101

100

100

0

20

40 60 80 100 120 140 Activation energy, E

0.0

0.2

0.4 0.6 0.8 1.0 Frequency, αim

1.2

1.4

Figure 2.7: Neutral stability for the fundamental mode at various γ: (a) neutral stability curves in (E ,Q ) plane, (b) frequency of oscillation along the neutral stability curves. Figure 2.6 displays the comparison between our results and those of Lee & Stewart (1990) for the fundamental mode and γ = 1.2. It can be seen that there is good agreement between the results. In Figure 2.7(a), the neutral stability curves of mode 0 in the (E ,Q ) plane are shown for various γ, while in Figure 2.7(b), the frequency of neutral oscillation is shown. The smaller values of γ are seen to extend the range of unstable E to smaller E when Q is larger than about 4, and to reduce the unstable range towards larger values of E at Q / 4. At the same time, the frequency of neutral oscillation is seen to decrease substantially with decreasing γ at large Q , while it is essentially independent of γ at small values of Q . However, for small Q , frequency exhibits nonmonotonic behavior for all γ except γ = 1.4.

For γ = 1.2, 1.3, and 1.4 neutral stability curves shown in Figure 2.7(a) represent stability boundaries for corresponding γ, which separate stable ZND solutions from unstable ones. However, the case γ = 1.1 is more complicated as, for 3.7 / Q / 36, mode 1 is more unstable than the fundamental mode. A typical time

series of the perturbation of detonation velocity ψ0 for such a situation is shown in Figure 2.8 for Q = 30.473 and E = 12.503, which corresponds to the mode 0 being neutrally stable (its growth rate is about −10−6 ) while mode 1 grows with a rate of

67 3

1e−10 Heat release, Q

2 ψ′

1

101

0 −1 −2 −3

0

10

20

t

30

40

10

50

15

20 25 30 Activation energy, E

35

Figure 2.8: Example time series of ψ0 for Figure 2.9: Neutral stability curves for γ = 1.1, in which mode 0 is neutrally sta- γ = 1.1 for modes 0 (solid line) and 1 ble while mode 1 is unstable. (dashed line). approximately 0.03. Figure 2.9 shows neutral stability curves for both mode 0 and 1 for γ = 1.1 on the same plot, where it can be clearly seen that the neutral stability boundary is determined by a composition of the neutral stability curves of these two modes. This intersection of neutral curves at γ = 1.1 has not been mentioned before in the literature, to the best of our knowledge.

2.6.4 Code verification by convergence study We verify the correctness of the unsteady linear solver described in Section 2.5 by studying self-convergence of the time series of ψ0 under grid refinement, as the agreement of the observed order of accuracy with the theoretical order of accuracy for solvers for partial differential equations is considered a demonstration of the correctness of the solver (Salari & Knupp, 2000; Roy et al., 2004; Roy, 2005). For this purpose we run simulations with parameters γ = 1.2, Q = 50, E = 50, and Tf = 20, collecting the sequence of time series {ψ0 }i , i = 1, . . . , 8 corresponding to resolutions N1/2 from 20 to 2560 with the refinement factor 2. Subsequently we measure the

relative errors between the time series

Ei =

kψ0i − ψ0i −1 k∞ kψ0i k∞

,

for i = 2, . . . 8,

68 and the observed order of accuracy r acc

r acc,i =

log(E i −1 /E i ) log2

for i = 3, . . . , 8.

Table 2.4 displays the obtained results; note that the observed order of accuracy r acc is five as expected. The final line of the table (for N1/2 = 2560) shows that r acc

decreases due to the accumulation of round-off errors. Table 2.4: Verification of the unsteady linear solver by studying self-convergence of the time series of detonation velocity ψ0 computed for γ = 1.2, Q = 50, E = 50, Tf = 20. N1/2 E r acc

20 40 N/A 3 × 10−2 N/A N/A

80 9 × 10−4 5.07

160 3 × 10−5 4.96

320 9 × 10−7 5.00

640 3 × 10−8 5.00

1280 1 × 10−9 4.81

2560 4 × 10−11 4.57

2.7 Conclusions In this chapter, we have introduced a novel method for linear stability analysis of detonation waves and applied it to the detonation model based on the reactive Euler equations with one-step chemistry. The method consists of two parts. First is a numerical solution of the time-dependent system linearized about the base ZND solution by a shock-fitting numerical algorithm of high order of accuracy. The result of such numerical calculation is a time series of the perturbation of detonation velocity for the linear problem. This solution is obtained only over relatively short time intervals which is computationally inexpensive. The second part of the method is an analysis of the obtained time series by using dynamic mode decomposition. As a result, we are able to calculate detailed stability properties of the model and compare our results to the normal-mode results obtaining excellent agreement, thus, validating the correctness of the method.

69 Being confident in the validity of the method, we extended the knowledge of the linear stability properties of the model by computing neutral stability boundaries for several values of specific-heats ratio γ and discovered that the neutral stability boundary for γ = 1.1 is determined by the composition of the first two modes, while for γ ∈ {1.2, 1.3, 1.4} it is determined only by the fundamental mode.

70

Chapter 3

Reactive Euler equations with two-step chemistry

3.1 Introduction In this chapter, we extend our work in Chapter 2 by applying the developed method for linear stability analysis to the detonation model based on the reactive Euler equations with a two-step reaction mechanism. The two-step model that we consider here consists of two consecutive reactions:

A → B → C,

(3.1)

where the first reaction is exothermic and the second endothermic. Again we assume that specific heats are equal for all species and there is no mole change in the mixture, such that all species have the same molar mass, and hence, gas constants. The reader may assume that this model is a simple extension of the one-step model of Chapter 2 by adding an extra reaction, and would be correct in thinking so if both reactions were exothermic. However, endothermicity in this model leads to qualitative changes of travelling-wave solutions compared to the solutions of the one-step model, as their existence is determined by the presence of the embedded sonic point inside the reaction zone (in contrast with the one-step model, where the sonic point is at the end of the reaction zone). The zone between the shock and the sonic point is determined uniquely, however, the flow between the sonic point and the equilibrium point is not, and can be

71 subsonic or supersonic, that is, the steady reaction zone is not actually steady. The detonation discussed in Chapter 2, is called equilibrium detonation (or CJ detonation) because the steady detonation velocity is determined through algebraic relations at the point of chemical equilibrium by the condition of the sonicity at that point and thus, it does not depend on the reaction rate but only on the overall heat release, which is maximum heat release due to the irreversible exothermic reaction. In contrast, the model considered in this chapter is often referred to as a “pathological” detonation or eigenvalue detonation (Fickett & Davis, 2000). It is quite different from the equilibrium detonation in the sense that the steady detonation velocity is determined by the reaction rates of both reactions and the maximum heat release obtained at the sonic point, that is, detonation velocity is an eigenvalue of a boundary-value problem for a system of ordinary differential equations with boundary conditions determined by the Rankine–Hugoniot conditions and the sonic-point conditions. The term “pathological detonation” (for a different model but with the same features) was coined by von Neumann (1942) and is based on the expectation that the detonation velocity is determined by the equilibrium heat release Q e , while it is actually determined by the maximum heat release Q max at the sonic point, hence, observed detonation velocity is larger than expected as Q e < Q max due to endothermicity. Due to the embedded sonic point, normal-mode analysis of this model is intricate (Sharpe, 1999) as it relies on the asymptotic expansions in the neighborhood of the sonic point to obtain the initial values for the integration with the shooting method towards the shock and towards the end of the reaction zone from the sonic point. Presence of the embedded sonic point thus makes normal-mode analysis difficult for this model. Therefore, our aim in this chapter is to apply the method to this model and investigate how the method performs when the base steady state contains an

72 embedded sonic point. We investigate both types of ZND solutions (subsonicsubsonic and subsonic-supersonic) possible for this model and compare our results qualitatively to the results of Sharpe (1999) for the subsonic-supersonic case. We also check the hypothesis that the stability properties are determined only by the region between the sonic point and the shock, albeit this analysis is challenging numerically as the subsonic-subsonic steady-state solutions are continuous but not differentiable at the sonic point, which leads to the linear hyperbolic system with discontinuous coefficients, hence, in addition, the perturbations develop discontinuities and fine grid resolution is required for the simulations. In spite of these difficulties, our method is capable of handling linear stability analysis of this model without any special treatment of the sonic point. Application of this method to this model is an important step in the assessment of the performance of the method for the detonation models with large kinetic mechanisms with possible nonmonotonic distribution of heat release in the reaction zone.

3.2 Governing equations Let the reactions (3.1) have rate constants k1 and k2 , heat releases Q 1 > 0 and Q 2 < 0, and activation energies, E 1 and E 2 , for the first (A → B) and second (B → C),

reactions, respectively. Then, if X A , X B , and X C denote the mass fractions of the species, the kinetic equations for the reactions follow the Arrhenius law: dX A E1 = −k 1 X A exp − , dt RT ¶ µ ¶ µ dX B E1 E2 = k 1 X A exp − − k 2 X B exp − , dt RT RT µ ¶ dX C E2 = k 2 X B exp − , dt RT µ



(3.2) (3.3) (3.4)

73 describing evolution of the mass fractions of the species in the absence of flow. Obviously, the sum of the mass fractions is unity: X A + X B + X C = 1, with initial mass fractions being X A = 1, X B = X C = 0. Instead of considering three equations for each species, we introduce the reaction progress variables for each reaction as λ1 = 1 − X A and λ2 = X C with λ1 and λ2 varying from zero at the shock to unity at

the end of the reaction zone with the reaction rates ¶ µ E1 ω1 = k 1 (1 − λ1 )exp − , RT µ ¶ E2 ω2 = k 2 (λ1 − λ2 )exp − . RT

(3.5) (3.6)

Then, transforming to the reference frame attached to the shock, the full system of the governing equations is

ρ t + (ρ (u − D ))x = 0,

(3.7)

(ρu )t + (ρu (u − D ) + p )x = 0, (ρe )t + (ρ (u − D ) e + pu )x = 0, (ρλ1 )t + (ρ (u − D )λ1 )x = ρω1 , (ρλ2 )t + (ρ (u − D )λ2 )x = ρω2 ,

(3.8)

with e = e i + u 2 /2 being the specific total energy, e i = pv /(γ − 1) − (Q 1 λ1 +Q 2 λ2 ) the specific internal energy, and v = 1/ρ the specific volume.

3.3 ZND solutions Again, as in Chapter 2, the ZND solution is defined by assuming steady-state solution of the system (3.7-3.8), that is, the wave moves with constant detonation velocity D¯ , which is a parameter to be determined. This leads to the algebraic-

74 differential problem :

¯ ρ¯U¯ = −ρ a D,

(3.9)

¯ 2, p¯ + ρ¯U¯ 2 = p a + ρ a D

(3.10)

e¯ + p¯ v¯ +

U¯ 2

2 ¯ dλ1

= ea + pava + =

¯2 D

,

2 ¯ 1 )exp(−ρE k 1 (1 − λ ¯ 1 /p¯) = , U¯ ¯ 1 − λ¯ 2 )exp(−ρE k 2 (λ ¯ 2 /p¯) = . U¯

ω ¯1

dx U¯ ¯ dλ2 ω ¯2 = dx U¯

(3.11) (3.12) (3.13)

¯ = −ρ a D¯ , P¯ = p a + ρ a D¯ 2 , H¯ = e a + p a v a + D¯ 2 /2 We introduce “big” quantities M and express U¯ and p¯ in terms of v¯: ¯ 2 v, p¯ = P¯ − M ¯

¯ v, U¯ = M ¯

and arrive at the quadratic equation for v¯ with the root

v¯ =

where

s δ = ± 1−

¯2¡ hM P¯ 2

γ P¯ (1 − δ), ¯2 γ+1 M

¢ H¯ + Q¯ ,

h=

2(γ2 − 1) γ2

,

¯ 1 +Q 2 λ¯ 2 . Q¯ = Q 1 λ

(3.14)

Also, it can be shown (Kasimov & Stewart, 2005) that ¶2 1 − M¯2 δ = , 1 + γM¯2 2

µ

where M¯ = U¯ /c¯ is the Mach number of the steady-state flow. Hence the sign of the square root in the expression for δ is determined by whether the flow is subsonic or supersonic. Also, δ = 0 when the flow is sonic. The regularity condition for the ZND solution to be continuous is determined

75 by the boundedness of the derivative of δ, hence ¯ 2 dQ dδ hM =∓ is finite, dx P¯ 2 δ dx when at the sonic point the maximum heat release is obtained:

δ∗ = 0 and

dQ ∗ =0 dx

(3.15)

for some x ∗ being the location of the sonic point. From the second of conditions (3.15) it follows using Eq. (3.12–3.13) that

δ∗ = 0 and Q 1 ω ¯ ∗1 +Q 2 ω ¯ ∗2 = 0,

(3.16)

which is known as the generalized Chapman–Jouguet condition (Kirkwood & Wood, 1954), that is, the detonation velocity is determined by the conditions (3.16) and generally calculated as follows (Lee, 2008; Higgins, 2012): for a guess D¯ , one solves steady reaction equations from the shock to the sonic point, x ∗ ; for arbitrary ¯ , conditions (3.16) will not be satisfied, and instead the steady-state solution D will blow up; however, for special D¯ , the solution will pass continuously (although not necessarily smoothly) through x ∗ to the equilibrium state at x → −∞, that is, ¯ is computed by shooting and determines the solution of the boundary-value D problem for system (3.9–3.13) with shock conditions and conditions (3.16) at the boundary points. This is the origin of the name “eigenvalue detonation” for this model. The general numerical procedure outlined above can be followed to determine ¯ . We, however, restrict our attention in this work to the case E 1 = E 2 and k1 = k2 , D for which equations simplify and detonation velocity can be found analytically. Considering only the reaction equations (3.12–3.13) and dividing the second-

76 reaction equation by the first-reaction equation, we arrive at dλ¯ 2 λ¯ 1 − λ¯ 2 = , dλ¯ 1 1 − λ¯ 1 which is solved as follows: denoting ψ = λ¯ 2 − λ¯ 1 , ξ = λ¯ 1 we obtain dψ −ψ dλ¯ 2 d(ξ + ψ) = 1+ = , = dξ dξ 1 − ξ dλ¯ 1 which gives a linear ODE of the first order for ψ: dψ ψ + +1 = 0 dξ 1 − ξ with the solution ψ = (1 − ξ) ln(1 − ξ). Transforming variables back, we obtain ¢ ¡ ¯ 2 = λ¯ 1 + (1 − λ¯ 1 ) ln 1 − λ¯ 1 , λ

(3.17)

and therefore, λ¯ 2 ≤ λ¯ 1 (the natural logarithm yields negative values when 0 ≤ ¯ 1 ≤ 1) with equality holding only for λ¯ 1 = 0 or λ¯ 1 = 1. This observation has λ the following consequence: both λ¯ 1 and λ¯ 2 monotonically increase towards the end of the reaction zone, as the reaction-rate expression of the first reaction is nonnegative, and the rate of the second reaction is also nonnengative as can be seen from the reaction-rate expressions (3.5–3.6) evaluated for the steady state. Hence, equilibrium heat release Q e = Q 1 +Q 2 is not a maximum value of Q in the reaction zone, which is illustrated on Figure 3.1. Now having the relation between λ¯ 2 and λ¯ 1 , we substitute them into (3.16) to obtain ¡ ¢ ¯ ∗ ) −Q 2 (1 − λ¯ ∗ ) ln 1 − λ¯ ∗ = 0, Q 1 (1 − λ 1 1 1

so either λ¯ ∗1 = 1 (which leads to an equilibrium solution as for the one-step model

77 45 40 35 30

Q

25 20 15 10 5 0 0.0

0.2

0.4

λ1

0.6

0.8

1.0

Figure 3.1: Total heat release Q¯ versus steady reaction progress variable λ¯ 1 . Notice that the maximum heat release Q¯ max is attained inside of the reaction zone. of Chapter 2) or λ¯ ∗1 = 1 − exp(Q 1 /Q 2 ). Both of these cases are considered below. ¯ ∗ = 1 must be discarded as it does not correspond Equilibrium solution. Case λ 1 to the condition δ = 0, moreover, the detonation velocity determined from δ = 0 for this case will be imaginary due to the negativity of the square-root expression for δ. Therefore, equilibrium solutions are impossible for this model. Solution with embedded sonic point. The sonic point is defined by condition

¯ ∗ = 1 − exp(Q 1 /Q 2 ). The flow is then divided into two parts by the sonic point. λ 1 Between the shock and the sonic point, the flow is always subsonic relative to the shock. Between the sonic point and the equilibrium point, the flow can be either subsonic or supersonic: flow is supersonic if expression (3.14) for δ is used with a “plus” sign, while subsonic if the sign is “minus”. Thus, two steady-state solutions with the embedded sonic point are possible for this model. Having determined the state λ¯ ∗1 , λ¯ ∗2 at the sonic point, we solve δ = 0, which ¡

¢

yields a similar result for steady detonation velocity D¯ as in the case of the one-step model:

s

¯= D

¢ γ2 − 1 ¡ ¯ ∗ +Q 2 λ¯ ∗ + c a2 + Q1λ 1 2

2

s

¢ γ2 − 1 ¡ ¯ ∗ +Q 2 λ¯ ∗ . Q1λ 1 2

2

(3.18)

78

3.4 Linear unsteady problem The linearized perturbation system is the same as it was in the previous chapter (see Eq. (2.22)), but with the vector of unknowns z = (ρ, u, p, λ1 , λ2 )T . The coefficients matrices are correspondingly redefined:  U¯   0   A= 0   0  

ρ¯

du¯ dx

   dp¯  − ρ¯12 dx    B = P C i (ρ¯ω ¯ iρ + ω ¯i)  i    −ω ¯ 1ρ   −ω ¯ 2ρ

0

1/ρ¯ 0

0 

0

0

0

0

0

du¯ dx

0 γ ddux¯ + ρ¯

dλ¯ 1 dx dλ¯ 2 dx

0



γp¯

dρ¯ dx

dp¯ dx

0



0





P

  

(3.19)

0 ,

0

  U¯ 0    0 U¯

0 0 ¯ip i Ci ω

ρ¯

P

¯ i λ1 i Ci ω

−ω ¯ 1p

−ω ¯ 1λ1

−ω ¯ 2p

−ω ¯ 2λ1

0



    0    P ρ¯ i C i ω ¯ i λ2     −ω ¯ 1λ2    −ω ¯ 2λ2

(3.20)

with C i = −(γ − 1)Q i for i = 1, 2. For the shock evolution equation we obtain

R¯s = (γ − 1)ρ¯s

X

Qi ω ¯ i s,

(3.21)

i

R s0 = (γ − 1)

X i

Qi

£¡ ¢ ¤¯ ρ¯ω ¯ iρ + ω ¯ i ρ 0 + ρ¯ω ¯ i p p 0 ¯s .

(3.22)

All the other expressions remain as before. The equations are rescaled using the Erpenbeck scales such that ρ a = p a = 1, velocity is scaled by

p

p a /ρ a , spatial scale is that of the half-reaction length. With

these scales, equations maintain their form. Activation energy and heat release are rescaled by the ambient temperature multiplied by R .

79 Numerical solution of this problem and postprocessing are performed as described in Subsection 2.5.

3.5 Results The model with two reactions represents an example of complex chemistry that brings in additional complications, which are absent in the one-step model. First, we now have two traveling-wave solutions for the same set of parameters of the problem, and therefore we must investigate the stability of both. Second, the solution has an embedded sonic point, that is the sonic condition is reached at a finite distance from the lead shock. This divides the post-shock flow into two regions that play distinct roles in the detonation dynamics. Following prior research on the role of sonic points in unsteady detonation dynamics (Kasimov & Stewart, 2004; Stewart & Kasimov, 2005), we expect that the dynamics is affected only by the solution between the shock and the sonic point. How exactly this idea plays out in the linear stability calculations is not obvious, and this will be explored in this section.

3.5.1 ZND and perturbation solutions Now the model has several additional parameters coming from the more complex chemistry. Comprehensive investigation of the entire parameter space would be a difficult task. Instead, we focus here on particular questions that have certain physical interest. Namely, we investigate how the level of endothermicity affects stability, and in addition, we compare the stability behavior of the two steady-state solutions that co-exist at the same parameter set. First, we study the behavior of ZND solutions when Q 2 is varied while all other parameters are kept fixed. Figure 3.2 shows profiles of ZND variables ρ¯, u¯ , p¯, T¯ , λ¯ 1 , and λ¯ 2 for the subsonic-supersonic case with γ = 1.2, Q 1 = 50 and E = 25 for several

80 values of Q 2 , obtained with N1/2 = 40. As can be seen from the figure, the lower the values of Q 2 , the larger the gradients of the profiles become. The maximum temperature decreases with increasing magnitudes of Q 2 , as expected. Figure 3.3 displays the same results as in Figure 3.2, but for the subsonicsubsonic solution, which is obtained now with N1/2 = 640 to accurately resolve the weak discontinuity in the solution. It can be seen that as the magnitude of Q 2 increases, and the sonic point moves toward the shock, the kink in the profiles of ρ¯, u¯ and p¯ becomes more pronounced.

Now we consider how the perturbation profiles z 0 (x, t ) = (ρ 0 , u 0 , p 0 , λ01 , λ02 )T behave taken at a particular time with Tf = 4. We take γ = 1.2, E = 15 and Q 1 ∈ {10, 20, 50} with Q 2 = −0.75Q 1 . Figure 3.4 shows the perturbations for the subsonic-

supersonic case, while Figure 3.5 shows the subsonic-subsonic case. The most prominent distinction between the two figures is that, while profiles for the subsonic-supersonic case are smooth, profiles for the subsonic-subsonic case develop discontinuities at the sonic point caused by discontinuous coefficient matrices A and B in (3.19–3.20): ZND profiles themselves are continuous, while their derivatives with respect to x have jump discontinuities at the sonic point, hence, discontinuous perturbations. In addition, we noticed that the fifth-order upwind method (2.33–2.34) with Lax–Friedrichs flux led to severe Gibbs oscillations on both sides of the jump in the profiles of the subsonic-subsonic case. For this reason, we replaced the upwind method with the WENO method (Jiang & Peng, 2000) for all simulations based on the subsonic-subsonic case. Another important observation from comparing figures 3.4 and 3.5 is that perturbations are nearly the same in the region between the sonic point and the leading detonation shock: for example, the relative difference between density perturbations is only 1% in this region and is likely caused by the smearing of the discontinuities in the neighborhood of the sonic point.

81 1.0

0.8

0.8

0.6

0.6 u/̄ ūs

ρ/̄ ρ̄s

1.0

0.4

0.4

0.2 0.0 −8

0.2 −7

−6

−5

−4 x

−3

−2

−1

0

0.0 −8

−7

−6

−5

−4 x

−3

−2

−1

0

−7

−6

−5

−4 x

−3

−2

−1

0

−7

−6

−5

−4 x

−3

−2

−1

0

1.0 2.5

0.8 0.6 p/̄ p̄s

T/̄ Ts̄

2.0

0.4

1.5

0.2 0.0 −8

−7

−6

−5

−4 x

−3

−2

−1

0

1.0 −8

0.8

0.8

0.6

0.6 λ2̄

1.0

λ1̄

1.0

0.4

0.4

0.2

0.2

0.0 −8

−7

−6

−5

−4 x

−3

−2

−1

0

0.0 −8

Figure 3.2: ZND subsonic-supersonic profiles of density ρ¯, flow velocity u¯ , pressure ¯ 1 and λ¯ 2 for γ = 1.2, Q 1 = 50, E = 25, and p¯, temperature T¯ , and progress variables λ different values of Q 2 . Solid line corresponds to Q 2 = −1; dashed line to Q 2 = −10; dotted line to Q 2 = −20; solid line with circle markers to Q 2 = −30. Density, flow velocity, pressure, and temperature are normalized by their values on the shock.

82 1.0

0.8

0.8

0.6

0.6 u/̄ ūs

ρ/̄ ρ̄s

1.0

0.4

0.4

0.2 0.0 −8

0.2 −7

−6

−5

−4 x

−3

−2

−1

0

0.0 −8

−7

−6

−5

−4 x

−3

−2

−1

0

−7

−6

−5

−4 x

−3

−2

−1

0

−7

−6

−5

−4 x

−3

−2

−1

0

1.0 2.5

0.8 0.6 p/̄ p̄s

T/̄ Ts̄

2.0

0.4

1.5

0.2 0.0 −8

−7

−6

−5

−4 x

−3

−2

−1

0

1.0 −8

0.8

0.8

0.6

0.6 λ2̄

1.0

λ1̄

1.0

0.4

0.4

0.2

0.2

0.0 −8

−7

−6

−5

−4 x

−3

−2

−1

0

0.0 −8

Figure 3.3: ZND subsonic-subsonic profiles of density ρ¯, flow velocity u¯ , pressure ¯ 1 and λ¯ 2 for γ = 1.2, Q 1 = 50, E = 25, and p¯, temperature T¯ , and progress variables λ different values of Q 2 . Solid line corresponds to Q 2 = −1; dashed line to Q 2 = −10; dotted line to Q 2 = −20; solid line with circle markers to Q 2 = −30. Density, flow velocity, pressure, and temperature are normalized by their values on the shock.

83 1e−9

6 5

0.5

4

0.0

3

−0.5

2

−1.0

1

−1.5

0

Q1 = 10

1.0

−2.0 −5 −4 1e−9

−3

−2

−1

0

−1 −5 −4 1e−10 6 5

1

4

0

3

−1

2

−2

1

−3

0

Q1 = 20

2

−4 −5 −4 1e−9

−3

−2

−1

0

1e−10

−1 −5 −4 1e−11 3

−3

−2

−1

0

−3

−2

−1

0

−2

−1

0

2

0.0

1 −0.5 Q1 = 50

0 −1

−1.0

−2 −3

−1.5 −2.0 −5

−4 −4

−3

x

−2

−1

0

−5 −5

−4

−3

x

Figure 3.4: Perturbation profiles for the subsonic-supersonic case for γ = 1.2, E = 15, Q 1 ∈ {10, 20, 50} with Q 2 = −0.75Q 1 . On the left: solid line ρ 0 , dashed line u 0 , dash-dotted line p 0 . On the right: solid line λ01 , dashed line λ02 . Perturbations are not scaled and are shown as they are at Tf = 4. Notice, that the left limits of the x -axes are taken to be 5 for clarity and are much smaller than the reaction lengths.

84 1e−9

6 5

0.5

4

0.0

3

−0.5

2

−1.0

1

−1.5

0

Q1 = 10

1.0

−2.0 −5 −4 1e−9

−3

−2

−1

0

−1 −5 −4 1e−10 6 5

1

4

0

3

−1

2

−2

1

−3

0

Q1 = 20

2

−4 −5 −4 1e−9

−3

−2

−1

0

1e−10

−1 −5 −4 1e−11 3

−3

−2

−1

0

−3

−2

−1

0

−2

−1

0

2

0.0

1 −0.5 Q1 = 50

0 −1

−1.0

−2 −3

−1.5 −2.0 −5

−4 −4

−3

x

−2

−1

0

−5 −5

−4

−3

x

Figure 3.5: Perturbation profiles for the subsonic-subsonic case for γ = 1.2, E = 15, Q 1 ∈ {10, 20, 50} with Q 2 = −0.75Q 1 . On the left: solid line ρ 0 , dashed line u 0 , dashdotted line p 0 . On the right: solid line λ01 , dashed line λ02 . Perturbations are not scaled and are shown as they are at Tf = 4. Notice, that the left limits of the x -axes are taken to be 5 for clarity and are much smaller than the reaction lengths.

85

3.5.2 Migration of spectra with increasing endothermicity Here we investigate how the stability behavior changes as the endothermicity increases, that is, as Q 2 increases in magnitude while all the other parameters are kept fixed. First, consider the subsonic-supersonic case. We choose N1/2 = 40, γ = 1.2, E = 25, Q 1 = 50 and let Q 2 vary in the range −45 ≤ Q 2 ≤ 0 with step ∆Q 2 = 0.0025

for Q 2 ∈ [−0.7; −0.8] ∪ [−13.9; −14] ∪ [−19.9; −20] ∪ [−25.2; −25.3] and step ∆Q 2 = 0.1 otherwise. Figure 3.6 shows how growth rates change as Q 2 decreases. Initially, at Q 2 = 0 detonation corresponds to the equilibrium CJ case and is stable for the γ, E

and Q 1 that we use (critical value of activation energy is 25.26). As Q 2 decreases, the growth rate of the fundamental mode (mode 0) increases almost linearly such that at Q 2 ≈ −0.79 detonation becomes unstable, and at Q 2 ≈ −25.2225 mode 0 splits into two purely exponential branches, with growth rates increasing and decreasing for the top and bottom branches, respectively. At the same time at Q 2 ≈ −13.1 the stable mode 1 appears in the spectrum, which becomes unstable at Q 2 ≈ −13.9425. As the growth rate of this mode increases faster than that of mode

0, at Q 2 ≈ −19.92 it becomes the dominant mode and remains so all the way down to Q 2 = −45. In total, up to six modes appear in the spectrum as Q 2 → −45. Here we point out that taking |Q 2 | larger than 45 makes computations prohibitively expensive as the effective reaction lengths become hundreds of thousands of nondimensional units, because as |Q 2 | → |Q 1 | the overall heat release becomes effectively zero and detonation degenerates to an inert weak shock. Next, consider the subsonic-subsonic case. The physical parameters are the same as in the previous case. However, due to the discontinuous coefficients in the linearized system (2.31), now we must use drastically higher resolutions (here N1/2 = 1280) than for the subsonic-supersonic case, otherwise, we are unable

to extract low-energy modes. We also emphasize that it is significantly more

86 1.2 1.0 0.8 αre

αim

0.6

18 16 5 14 4 12 3 10 2 8 6 1 4 2 0 0 −45 −40 −35 −30 −25 −20 −15 −10 −5 Q2

0.4 0.2

5 4 3 2 1 0 0.0 −45 −40 −35 −30 −25 −20 −15 −10 −5 0 Q2

0

Figure 3.6: Growth rate versus Q 2 for −45 ≤ Q 2 ≤ 0 with other parameters γ = 1.2, E = 25, Q 1 = 50 for subsonic-supersonic case. 1.2 1.0 0.8 αre

0.6 0.4 0.2

5 4 3 2 1 0 0.0 −45 −40 −35 −30 −25 −20 −15 −10 −5 0 Q2

Figure 3.7: Growth rate versus Q 2 for −45 ≤ Q 2 ≤ 0 with γ = 1.2, E = 25, Q 1 = 50. Solid lines: subsonic-supersonic case, dashed lines: subsonic-subsonic case. Note that in most cases the solid and dashed curves overlap and are thus indistinguishable. difficult to postprocess the time series of perturbation of detonation velocity for the subsonic-subsonic case. The reason is that the time series contain substantial amount of numerical noise caused by the numerical smearing of the jump in the solution in the neighborhood of the sonic locus. Figure 3.7 shows the spectra for subsonic-supersonic and subsonic-subsonic cases together. It can be seen that the spectra are close to each other. In fact, they are essentially indistinguishable for the oscillatory part of mode 0, upper branch of exponential part of mode 0, and for modes 1, 2, and 3. However, for the subsonic-subsonic case, the lower branch of exponential part of mode 0 shows significant numerical noise as Q 2 decreases. The same behavior is seen with

87 mode 4: the reason being that the discontinuity in the solutions for the subsonicsubsonic case generates noise that contaminates the time series of ψ0 . As the lower branches of modes 0 and 4 contribute much less energy to the time series than other modes, they are more difficult to distinguish from the noise, which is an issue that can probably be fixed by using an improved algorithm for solving hyperbolic problems with discontinuous coefficients. Nevertheless, the computed values are near their counterparts for the subsonic-supersonic case. Despite some differences in the spectra, Figure 3.7 serves as a demonstration that the stability of ZND solutions of the two-step model with endothermicity does not depend on whether the underlying ZND solution is of the subsonic-supersonic or the subsonic-subsonic case. That is, the detonation stability is apparently determined solely by the dynamics of the reaction zone between the sonic point and the leading detonation shock. As a final note on this section, we give the reader an idea of the computational requirements for these calculations. For both cases, subsonic-supersonic and subsonic-subsonic, the number of values of Q 2 we took was approximately 500. Then, as the resolution necessary for the subsonic-supersonic case is only N1/2 = 40, we were able to compute Figure 3.6 overnight on a workstation with 16 cores. In contrast, the subsonic-subsonic case requires much higher resolution of N1/2 = 1280, so these computations required a day with 512 parallel workers on a parallel computer, where a worker simulates and postprocesses results for a given Q 2 independently of other workers.

3.5.3 Neutral stability curve In this section, we study how endothermicity affects the neutral stability boundary in comparison with the one-step model. The comparison is based on the idea that the detonation dynamics in both cases is determined by the energy released

88 between the shock and sonic point. Therefore, the one-step model with heat release Q is compared with the two-step model with maximum heat release Q ∗ equal to Q . Figure 3.8 shows the neutral stability curves in the (E ,Q ) parametric plane. We use γ = 1.2 and therefore the neutral stability curve for the one-step model is the same as in Figure 2.6. The neutral curve for the two-step model is obtained for the ZND solution corresponding to the subsonic-supersonic case. The range of Q 1 is 10−0.35 –102.35 with 256 logspaced points while Q 2 is varied in proportion to Q 1 : Q 2 = −0.75Q 1 . For each Q 1 we find the critical value of E using the algorithm

described in Subsection 2.6.3. For values of 10−0.35 ≤ Q 1 / 0.8 (approximately 20 points) the algorithm did not converge and we omit these values from the figure. The curves shown demonstrate that endothermicity consistently shifts the neutral stability boundary toward smaller values of E over the whole range of Q that we investigated. The difference, ∆E , between critical activation energies for the two models behaves with varying Q as follows: the minimum value ∆E ≈ 1.5 is attained for 5 / Q / 15; as Q increases, ∆E increases monotonically and reaches ∆E ≈ 4.5 as Q → 100; for low values of Q , ∆E increases quickly such that ∆E → ∞ as the detonation wave becomes an inert shock in the limit Q → 0. Additionally, while for the one-step model all ZND solutions are stable at E ≤ 14.2 for the full range of Q considered here, for the two-reaction model all ZND solutions are stable only

for E ≤ 12.6. Thus, the two-step exothermic-endothermic heat release leads to an increased instability compared to the equivalent one-step heat release. In order to provide some qualitative explanation for this observation, we look at how the steady-state thermicity behaves for the same values of E and Q ∗ for the two models. The thermicity for the two-step model is given by σ¯ 1 ω ¯ 1 + σ¯ 2 ω ¯ 2 , where σ¯ i = (γ − 1)Q i /c¯2 , i = 1, 2. Figure 3.9 shows the thermicity profiles for γ = 1.2, Q ∗ = 10, and E = 15

89 1.0 0.8 101

Thermicity

Max heat release, Q

102

100

10

−1

0.6 0.4 0.2 0.0

0

10

20 30 40 Activation energy, E

50

Figure 3.8: Neutral stability curves for two-step model with subsonicsupersonic steady-state solution (solid line) and one-step model (dashed line).

−0.2 −10

−8

−6

x

−4

−2

0

Figure 3.9: Thermicity profiles for Q ∗ = 10, E = 15 for two-step model with subsonic-supersonic steady-state solution (solid line) and one-step model (dashed line).

which correspond to the unstable steady-state solutions from Figure 3.8. It is clear that for the two-step model the thermicity has significantly a larger peak value which in addition is located closer to the lead shock than the value for the one-step model. Thus, a possible reason for the increased instability is this sharper profile of the thermicity function that makes for a more effective resonant cavity for the perturbations between the shock and the sonic point. Recall that the thermicity is a source term in the equation for the right-going characteristics p˙ + ρc u˙ = ρc 2 (σ1 ω1 + σ2 ω2 )

with dots denoting the derivatives along the trajectories. As such the thermicity could be looked at as a quantity showing how much the chemical reactions amplify the compression waves that influence the lead shock. The larger the magnitude of thermicity, the stronger the effect of amplification. The sharpness of the thermicity function reflects in the sensitivity of this influence on the changes in the local state in the reaction zone, the stronger sensitivity implying more instability. This is a qualitative argument explaining the shift of the stability boundary.

90

3.6 Conclusions In this chapter we applied the method developed in this work for linear stability analysis to the detonation model based on the reactive Euler equations with the two-step kinetic mechanism with one reaction being endothermic which gives rise to the equilibrium detonations with embedded sonic point. We demonstrated that the method is capable of analysis of this model without any special treatment for the sonic point. However, numerical solution of the linear unsteady problem for this model is challenging when the underlying base steady-state solution is of the subsonic-subsonic type as the linearized system has discontinuous solutions in this case, hence, fine grid resolution is required for the extraction of stability information from the perturbations of detonation velocity. We also compared migration of the linear spectrum when endothermicity is varied and found that the both steady-state solution possess the same spectrum with minor differences in less energetic modes due to the numerical noise generated around the discontinuity in the solution. Nevertheless, we demonstrated that the linear spectrum (and hence detonation dynamics) is determined by the region of the reaction zone between the sonic point and the shock, at least in the linear case.

91

Chapter 4

Fickett’s detonation analog

A toy problem can be useful in learning about a complicated property of fluid motion if the simple toy problem describes that property, even if it does not describe most other properties of the motion. P. G. Drazin (2002, p. 12)

4.1 Introduction In this chapter we consider a model that is much simpler than conventional detonation models are, yet is able to reproduce the unsteady dynamics of these models in terms of the possibility of stable and unstable ZND solutions. There are several such models, for example, Fickett’s model (Fickett, 1979), Majda’s model (Majda, 1981), or the model by Kasimov et al. (2013). However, Majda’s model contains a diffusive term as it is based on viscous Burgers’ equation (Burgers, 1948), and thus the model is not hyperbolic, and the model by Kasimov et al. (2013) contains nonlocal source term, which is a very unusual feature for the conventional detonation models. Therefore, we choose Fickett’s model (also called Fickett’s detonation analog in the following text) as it most closely resembles mathematically hyperbolic systems considered in the previous chapters. Fickett’s detonation analog was proposed by Fickett (1979). This model is not rigorously derived from equations of reactive fluid mechanics using asymptotic

92 analysis. Quoting the author of the model himself (Fickett, 1985b, section 2.1), “It is perhaps best regarded as a mathematical construct, invented to act as a prototype for the Euler equations usually used to describe the physical system.” Fickett created his model by coupling inviscid Burgers’ equation with an equation of one-step reaction (here the flux in the first equation is taken to be in the form used by Fickett (1984)):

ut +

´ 1³ 2 u + qλ = 0, x 2

(4.1)

λt = ω,

(4.2)

where u is the quantity analogous to the velocity or density in physical systems (in (Fickett, 1979) the symbol ρ is used, whereas we use u to emphasize that it is a one-dimensional vector), λ is the reaction progress variable, q is the total heat release, and ω is the reaction rate. Later, Fickett studied his model with different expressions for reaction rate (Fickett, 1985a), however, for all cases he considered the ZND solution was stable. In Clarke et al. (1993a) and Clarke et al. (1993b) Fickett’s model with a very simple reaction rate expression and analytic stable solution was used for assessment of numerical methods applicable for reactive flows. It was understood that the instability in this model can be obtained by changing chemistry. Radulescu & Tang (2011) considered Fickett’s detonation analog in which the flow following the shock is spatially divided into two distinct zones: an induction zone with no reaction and whose length depends on the strength of the shock, and a reaction zone itself with exothermic reaction. By employing this model they obtained nonlinear dynamics with limit cycles and period-doubling bifurcations with eventual onset of chaos similar to the dynamics observed in the reactive Euler equations. We chose to use the one-step chemistry with an expression for reaction rate

93 appearing in the derivation of the asymptotic model of the reactive Euler equations by Faria (2014); Faria et al. (2015). It was noted that their expression for the reaction rate being supplied to the Fickett’s detonation analog gives rise to the spectrum of solutions including stable and unstable ZND solutions, similar to observed in the detonation models with one-step chemistry based on the Euler equations. So, in this chapter we aim to study the linear stability properties of this model using the method developed in this work, and also nonlinear stability properties of this model.

4.2 Governing equations In this section we define the governing equations, ambient and Rankine–Hugoniot conditions, and the shock-evolution equation.

4.2.1 Equations in laboratory and shock-attached frames Governing equations for the Fickett’s detonation analog consist of equations of motion and reaction:

ut +

´ 1³ 2 u + qλ = 0, x 2 λt = ω (u, λ) ,

(4.3a) (4.3b)

where t is the time, x the coordinate, u the velocity, λ is the progress variable showing the degree of completeness of the chemical reaction (λ = 0 points to no reaction state, λ = 1 points to complete reaction state), q the heat release, and ω the reaction rate which is a function of state variables u and λ. Reaction rate is defined by the expression (Faria, 2014; Faria et al., 2015) ¡ ¡p ¢¢ ω(u, λ) = k (1 − λ)exp θ qu + qλ ,

(4.4)

94 where k is the preexponent (constant defined in Section 4.3), and θ the activation energy, which is also constant. We also define the thermicity coefficient

σ=

q

2

,

which is used in the following text. All the quantities in this model are assumed to be nondimensional. As in previous chapters, we transform the governing equations to the shockattached frame by assuming that the detonation wave moves to the right with velocity D (t ) > 0. Then the equations in a conservative form in the shock-attached frame become

ut +

1 ((u − 2D ) u + qλ)x = 0, 2

(4.5a)

λt − Dλx = ω.

(4.5b)

and in nonconservative form, which will be more convenient for linearization in Section 4.4, they are

u t + (u − D ) u x + σλx = 0,

(4.6a)

λt − Dλx = ω.

(4.6b)

4.2.2 Ambient and Rankine–Hugoniot conditions As ambient conditions (flow state ahead of the shock) we take a trivial constant state: u a = 0,

λa = 0.

(4.7)

The Rankine–Hugoniot condition for λ follows from the physical consideration.

95 We know that ahead of the shock λ = 0. The shock per se is very thin, such that no chemical reaction occurs inside it. Therefore

λs = 0.

(4.8)

To derive the Rankine–Hugoniot condition for u we rewrite Eq. (4.6a) in a conservative form: ut +

³³ u

´ ´ − D u + σλ = 0,

2

x

and integrate over domain [x 1 , x 2 ] containing shock: x 1 < x s < x 2 : ∂ ∂t

xs

µZ

x1

Z u dx +

x2 xs

¶ Z u dx = −

x 2 ³³ u

2

x1

´ ´ − D u + σλ d x. x

(4.9)

Note that x 1 , x 2 , and x s do not depend on t . Therefore, Eq. (4.9) can be simplified to: xs

Z

x1

Z ut d x +

x2 xs

ut d x = −

³³ u

2

´ ´¯x2 ¯ − D u + σλ ¯ . x1

Now we let x 1 → x s− and x 2 → x s+ , such that integrals on the left-hand side vanish, which yields ³u

2

¯x 2 ´ ¯ − D u + σλ¯ = 0. x1

Using ambient conditions (4.7) and condition (4.8) we arrive at ³u

s

2

´

− D u s = 0,

and the Rankine–Hugoniot condition for u is

u s = 2D.

(4.10)

96

4.2.3 Shock-evolution equation In equations (4.5) detonation speed D appears explicitly, therefore, we need a way to compute D . In this subsection we derive the shock-evolution equation for this purpose. We want to find the rate of change of u on the shock path:

u| ˙ s = [u t + Du x ] |s .

Using the equation (4.3a):

u t = −uu x − σλx ,

we arrive at u| ˙ s = [(D − u )u x − σλx ] |s .

By noting that no chemical reaction occurs right behind the shock and using (4.3b) we replace λx : ˙ s = λt + Dλx = 0, λ| therefore, λ x |s = −

ω ¯¯ ¯ . D s

Using the relation u| ˙ s=

du s dD =2 , dt dt

which connects the evolution of u on the shock path with function u s defined by the Rankine–Hugoniot condition (4.10), and find that the shock-evolution equation is dD 1 h σ i¯¯ = (D − u )u x + ω ¯ . dt 2 D s

(4.11)

97 The shock-evolution equation was first derived by Fickett (1985a) (with the name “shock-change relation”). The equation allows us to know detonation speed D at any time if the initial condition for D is provided.

4.3 ZND solution System (4.5) admits the ZND solution—one-dimensional travelling-wave solution with constant detonation velocity such that an observer moving at the same velocity sees wave structure that does not change in time. To denote steady quantities we use the “bar” symbol here and further. So, for constant speed D¯ and vanishing time derivatives, equations (4.5) reduce to q d u¯ 2 ¯ u¯ + λ¯ = 0, −D dx 2 2 dλ¯ ω ¯ =− . ¯ dx D ¶

µ

(4.12) (4.13)

Using ambient conditions 4.7 we obtain from Eq. (4.12) an algebraic relation u¯ 2

2

¯ u¯ + −D

q

2

¯ = 0, λ

which is a quadratic equation for u (λ) with roots: ¯± u¯ = D

q

¯ ¯ 2 − q λ, D

and we take the root that satisfies Rankine–Hugoniot conditions (4.8, 4.10):

¯+ u¯ = D

q

¯ ¯ 2 − q λ. D

(4.14)

Therefore, steady structure is defined by substituting (4.14) into (4.13) and

98 integrating with Rankine–Hugoniot conditions as initial conditions starting at the shock: ¯ u¯s = 2D,

¯s = 0 λ

(4.15)

towards the chemical equilibrium. From Eq. (4.14) it follows that the Chapman–Jouguet detonation velocity is defined by D CJ =

p

q.

(4.16)

We define preexponent k such that characteristic length scale is half-reaction length l 1/2 of the steady reaction zone: Z 1/2 ¯ Z l 1/2 D dλ¯ = d x = 1, ω ¯ 0 0

therefore, Z 1/2 k=

0

¯ D ¯ ¢ ¡ ¡p ¢¢ dλ, ¯ ¯ 1 − λ exp θ q u¯ + q λ

¡

(4.17)

where u¯ is defined by (4.14). For further work we will need expressions for du¯ /dx and dλ¯ /dx : ω ¯ dλ¯ =− , ¯ dx D

(4.18) dλ¯

du¯ 1 −q dx 1 qω ¯ = q = . q dx 2 ¯ 2 ¯ 2 D¯ D¯ 2 − q λ¯ D − qλ

(4.19)

4.4 Linear unsteady problem In this section we linearize equations (4.6) about the steady state found in Section 4.3. Let z = (u, λ)T denote the vector of state variables. Then, expanding z and ¯ + ψ0 (t ), we arrive at the linearized D about the steady state: z = z¯ (x ) + z 0 (x, t ), D = D

99 governing equations z 0t = −A (z¯)z 0x − B (z¯)z 0 +

where



¯

u¯ − D A=

0





σ  , ¯ u¯ − D

dz¯ 0 ψ, dx

du¯ dx

(4.20)



0   B = , −ωu −ωλ

in which the steady quantities are known, with partial derivatives of the reaction p rate being ω ¯ u = θ qω ¯ , ωλ = k exp θ u¯ + q λ¯

¢¢ ¡ ¢ ¯ ) − 1 , and the perturbations θq (1 − λ

¡ ¡

u 0 , λ0 and ψ0 are to be found. The linearized Rankine–Hugoniot conditions are

λ0s = 0,

u s0 = 2ψ0 .

(4.21)

Lastly, the linearized shock-evolution equation is · ¸ dψ0 1 qk exp(2θ D¯ )(θ D¯ − 1) 0 0 ¯ = ψ − Du x |s . ¯2 dt 2 D

(4.22)

Linear unsteady computations start with the evaluation of the following parameters: value of self-sustained detonation velocity D¯ , which is found using Eq. (4.16), value of the pre-exponential factor k (by (4.17)), and the reaction zone length L which is found by integrating Eq. (4.13) up to λ¯ = λ∗ using separation of variables: ∗

L (λ ) =

Z

λ∗

0

¯ D

¯ ¢ dλ, ω ¯ u, ¯ λ¯ ¡

(4.23)

where u¯ = u¯ λ¯ is given by relation (4.14). Notice that in Eq. (4.23) λ∗ < 1 be¡ ¢

cause the integral on the right-hand side is divergent due to the factor 1 − λ¯ ¡

¢

in the denominator of the integrand. Therefore, for practical computations we should not integrate the ZND solution to λ¯ = 1 but to λ∗ = 1 − τλ , where τλ is a tolerance describing the deviation of the practical final value of λ¯ from the chemical-equilibrium value λ¯ = 1. Therefore, the length of the computational

100 domain L is defined by numerical evaluation of the integral (4.23) with λ∗ = 1 −τλ . Once these parameters are computed, all steady-state quantities are computed on the grid. As initial conditions we use the conditions

u 0 (x, 0) = 2 A 0

u¯ (x ) , u¯s

¯ ( x ), λ0 (x, 0) = A 0 λ

ψ0 (0) = A 0 ,

(4.24)

where A 0 is the amplitude of the initial perturbation. Notice that the initial condition for ψ0 is chosen such that initial and the Rankine–Hugoniot conditions match. Then unsteady linear simulations are conducted using the algorithm described in Subsection 2.5.

4.5 Nonlinear simulations by Godunov method 4.5.1 Description of the method To solve the nonlinear system (4.5), we use a second-order MUSCL scheme with the minmod flux limiter (LeVeque, 2002), which we describe below along with the algorithm for the solution of the Riemann problem required for this method. MUSCL scheme is a conservative Godunov-type method for a hyperbolic system: z t + f ( z )x = 0

with z ∈ Rn being the vector of unknowns and f : Rn → Rn being the flux function. Partitioning computational domain into cells [x i −1/2 , x i +1/2 ], we obtain a semidiscretized scheme: d¯z i ∆t =− (F i +1/2 − F i −1/2 ) , dt ∆x

101 where z¯ i is the average value of z for the computational cell centered at x i , and the fluxes at the cell boundaries F i +1/2 are approximated by reconstructing z i +1/2 . The reconstruction is done through the solution of the Riemann problem with the initial conditions z i +1/2 =

   z

L

  z , R

x < x i +1/2 , x > x i +1/2 ,

where z L and z R are given immediately on the left and the right of x i +1/2 , respectively. For a second-order scheme, they are found by assuming linear distribution of z i (x ) in the cell [x i −1/2 ; x i +1/2 ]: z i (x ) = z¯ i + σi (x − x i ),

where σi is the slope of the linear reconstruction, therefore, z L = z¯ i −1 + σi −1 z R = z¯ i − σi

∆x

2

∆x

2

,

.

For the second-order MUSCL schemes, σi = σi (¯z i − z¯ i −1 , z¯ i +1 − z¯ i ). We compute σi using the minmod flux limiter LeVeque (2002). The formula for its computation is given componentwise below:

σi (a, b ) =

     min(a, b )/∆x,    

max(a, b )/∆x,        0 ,

if a > 0, b > 0, if a < 0, b < 0, if ab < 0.

For the model under consideration, solution z i +1/2 of the Riemann problem corresponds to state (u i +1/2 , λi +1/2 )T , which is found by analysis of the waves that

102 can occur in the Riemann problem and is given in the next subsection. This solution propagates along the line (x i +1/2 , t )T in the phase plane (x, t ). As each Riemann problem assumes x i +1/2 = 0 in a local reference frame, then this solution is time-independent. Once all the local Riemann problems are solved over the whole grid, fluxes F i +1/2 at the edges x i +1/2 can be computed. On the left part of the domain we

introduce a ghost point with extrapolation u −1 = u 0 such that the left-most flux F −1/2 can be computed. On the right boundary the flux F N +1/2 is computed by

means of the Rankine–Hugoniot conditions u N +1/2 = 2D , λN +1/2 = 0, which leads to the fluxes F N +1/2 = 0 both for u and λ for the chosen ambient conditions.

4.5.2 Riemann problem for the Fickett’s model Below we provide the solution of the Riemann problem for our system, which in the characteristic form is dx dp = 0 on = u − D, dt dt dλ dx = 0 on = −D, dt dt

(4.25) (4.26)

where p=

1 ((u − 2D ) u + qλ) 2

is the flux and, from the characteristic equations themselves, we can see that the first wave is nonlinear as it depends on u (and correspondingly on λ due to coupling), while the second wave is linear and, therefore, the second wave is always the left-going wave with speed −D . Case 1. u L < u R , λL < λR . In this case from the second characteristic equation

it follows that λ experiences a jump along the line x = −D t , hence this line is the trajectory of the contact discontinuity moving to the left. Then, for M-region,

103 which is to the right of the contact, we have λM = λR . The invariant p of the first characteristic equation does not change, hence u must jump also:

pM = pL



´ 1³ ´ 1³ 2 u M + qλM = u L2 + qλL , 2 2

(4.27)

from which it follows

uM =

q q u L2 + q (λL − λM ) = u L2 + q (λL − λR ) < u L < u R .

(4.28)

As u M < u R , it follows that the first characteristic equation leads to the centered rarefaction wave, with constraining characteristics

x HEAD = (u R − D ) t ,

x TAIL = (u M − D ) t ,

with the solution inside the rarefaction wave being u = x /t + D , which is found from the integration of the trajectory of the first family of characteristics with initial conditions x = 0, t = 0. Case 2. u L < u R , λL ≥ λR . In this case λM = λR and u M is found exactly as before

but now u L ≤ u M < u R . Hence, again one wave is a contact discontinuity and another wave is a rarefaction wave, that is, this case is the same as case 1. Case 3. u L ≥ u R , λL ≥ λR . In this case λM = λR again and

uM =

q u L2 + q (λL − λR ) ≥ λL ≥ λR ,

that is, C + -characteristics from M- and R-regions collide and form the shock wave that moves with the speed s=

uR + uM

2

− D.

Case 4. u L ≥ u R , λL < λR . In this case λM = λR and u M is found according to the

104 above formula. Following this, the first wave is either a centered rarefaction or shock wave depending on whether u M < u R or u M ≥ u R and the properties of these waves a found according to the above formulas. For the Godunov method, we are interested in the solution of the Riemann problem along the trajectory x = 0 only. Summarizing the cases considered above, we can conclude that if u R > D (subsonic flow), then u i +1/2 = u M , λi +1/2 = λM . However, while developing the nonlinear unsteady solver it was noticed that occasionally during a simulation the flow becomes supersonic with u R < D , particularly near the left boundary, and then not only the contact moves to the left but the rarefaction or shock wave as well. In these cases the solution of the Riemann problem is taken to be z i +1/2 = u R , λi +1/2 = λR . One possible explanation is that, although theoretical chemical equilibrium λ = 1 is attained at x → −∞, in numerical simulations it occurs at a finite distance from the leading shock due to numerical errors. To illustrate the complete solutions to the Riemann problem and characterstic trajectories, we solve two such problems. For both problems we take q = 9, hence, D=

p

q = 3, λL = 1.0, λR = 0.0, and final time t = 1. Figure 4.1 shows the solution

of the Riemann problem for u L = 4, u R = 3, along with the characteristic plane. As can be seen, the solution consists of three distinct regions: in the region to the left of the contact discontinuity is (u L , λL ), between the contact and the shock (u M = 5, λR ), and to the right of the shock (u R , λR ). As u R = D , the characteristics in the rightmost region are vertical. The shock speed in this case is s = 0.5. Figure 4.2 shows the solution of the Riemann problem for u L = 4 and u R = 6, along with the characteristic plane. In this case, the solution consists of four distinct regions: to the left of the contact discontinuity x = −D t solution is (u L , λL ), between the contact and the tail of the rarefaction wave solution is (u M = 5, λR ), inside the rarefaction wave it is (u = x /t + D, λR ), and to the right of the head of the

105 rarefaction solution solution is (u R , λR ). The head and tail of the rarefaction are determined by the lines x = (u R − D ) t and x = (u M − D ) t , respectively.

4.5.3 Computation of detonation velocity Note that we need to update D in the beginning of the time step to carry out the solution of the Riemann problems and the boundary conditions. To update detonation velocity D , we use a procedure from Kasimov (2004) and Kasimov & Stewart (2004), slightly modifying it as for the Godunov method the grid is staggered, that is, the von Neumann state of the lead shock is located at x N +1/2 = 0. We will denote quantities with index N + 1/2 using subscript “s” in this subsection for clarity of notation. The equation of the C + -characteristic is

u

du dλ dx +σ − σω(u, λ) = 0 on = u − D. dt dt dt j

j

At the time step j + 1 we are given u i , λi and D j for i = 0, . . . , N , and aim to find D j +1 . We assume that during the time step ∆t from t j to t j +1 characteristics behave as straight lines and then approximate the trajectory of the characteristic that arrives at the shock at time t j +1 using the forward Euler method (LeVeque, 2007): x s − x∗ = j

³

j u∗ − D j

´

∆t ,

j

where u ∗ = u (x ∗ , t j ). To find u ∗ , which is the flow velocity at time step j at location x ∗ , we interpolate linearly using the value at the shock and the average of u in the

rightmost computational cell: j

j

us − u N xs − x N

j

=

j

u∗ − u N x∗ − x N

.

106

5.0

u

4.5 4.0 3.5 3.0 −4

−3

−2

−4

−3

−2

−4

−3

−2

1.0

x

−1

0

1

−1

0

1

−1

0

1

0.8 λ

0.6 0.4 0.2 0.0

1.0

x

0.8

t

0.6 0.4 0.2 0.0

x

Figure 4.1: Solution of the Riemann problem and characteristics trajectories for q = 9, D = 3, (u L , λL ) = (4, 1), (u R , λR ) = (3, 0) up to the final time t = 1. Bottom plot shows characteristics of the first family. Characteristics of the second family are not shown except for the trajectory of the contact discontinuity x = −D t (all other characteristics of this family are parallel to it).

107

6.0

u

5.5 5.0 4.5 4.0 −4

−3

−2

−1

0 x

1

2

3

4

−4

−3

−2

−1

0 x

1

2

3

4

−4

−3

−2

−1

0 x

1

2

3

4

1.0 0.8 λ

0.6 0.4 0.2 0.0

1.0 0.8

t

0.6 0.4 0.2 0.0

Figure 4.2: Solution of the Riemann problem and characteristics trajectories for q = 9, D = 3, (u L , λL ) = (4, 1), (u R , λR ) = (6, 0) up to the final time t = 1. Bottom plot shows characteristics of the first family. Characteristics of the second family are not shown except for the trajectory of the contact discontinuity x = −D t (all other characteristics of this family are parallel to it).

108 Using relations x s − x N = ∆x /2, x N = −∆x /2, we arrive at the following formula for j

u∗ :

³

j

j

u∗ = us + j

j

´

j

2 us − u N

x∗ ,

∆x

j

where u s = 2D j and u N are known quantities. Similarly we obtain the formula for j

λ∗ :

³

j

j

λ∗ = λs +

j

j

2 λs − λN

´ x∗

∆x

j

with λ∗ = 0. Approximating the trajectory equation with the forward Euler method j +1

xs

j − x∗

=

³

j u∗ − D j

´

∆t

j

and by substituting the expression for u ∗ we find that x ∗ can be computed explicitly: j

x∗ =

D j ∆t ³

−1 − j

j

j

2 u s −u N

´.

∆x

j

After computing x ∗ , u ∗ , and λ∗ , we can compute D j +1 by solving a nonlinear equation arising from the discretization of the C + -equation: j +1

u∗

us

j

j +1

j

³ ´ − u∗ λ − λ∗ j +1 j +1 j +1 +σ s − σω u s , λs = 0, where u s = 2D j +1 , ∆t ∆t

j +1

λs

= 0.

The time integrator is an adaptive-step Runge–Kutta integrator DOPRI5 (Hairer et al., 1993) as used for linear unsteady simulations. The reader may object that strong-stability-preserving Runge–Kutta integrators (Gottlieb et al., 2009) are more appropriate for the problem at hand as we expect internal shock being generated inside the reaction zone. However, for the reactive flow, error estimation and consequently adaptive choice of a time step are crucial components of a time integrator, as one cannot rely on the Courant–Friedrichs–Lewy condition (Courant

109 et al., 1928) in the computations of the time step due to chemical reactions occurring at much smaller time scales than wave propagation. The presence of the adaptive-step DOPRI5 integrator in the numerical environment that we use (Scipy library of the scientific Python stack), hence, determined our choice of the time integrator. Relative tolerance and absolute tolerance of the DOPRI5 integrator are set to 10−8 and 10−6 , respectively, for the nonlinear simulations.

4.6 Results If not stated explicitly, all the results described here are obtained using the following parameters: q = 4, initial perturbation is A 0 = 10−10 , and the length of the computational domain is taken to be such that λ¯ ∈ [0, 1 − tolλ ] with tolλ = 10−6 .

4.6.1 ZND solutions Figure 4.3 shows the profiles of velocity u¯ and reaction progress variable λ¯ versus spatial coordinate x in the ZND solution for fixed value of heat release q = 4 and varying activation energy θ = {0.5, 1, 2, 5}. The figure demonstrates that, as θ increases, the profiles develop a steep slope, and hence, the ZND solution becomes analogous to “square-wave” detonations, which are known to be unstable (Fickett (1985b); Short (1997); He (1999)). Therefore, we expect that, for any fixed q , with the increase of θ , the ZND solution loses stability for some critical value of activation energy.

4.6.2 Demonstration of the stable and unstable solutions Figure 4.4 shows the time series of the perturbation of detonation velocity ψ0 for stable and unstable solutions for the same value of q with θ = 0.92 and θ = 0.95. In Figure 4.4the stable case is shown, in which initial perturbation of the detonation

110

Figure 4.3: ZND profiles for q = 4: a) velocity u¯ , b) reaction progress variable λ¯ as θ is varied: θ = 0.5 (solid line), θ = 1 (dashed line), θ = 2 (dashed-dotted line), θ = 5 (dotted line). As θ increases, ZND profiles become steeper and approach “square-wave” structure, in which reaction effectively occurs in a narrow region near the lead shock at x = 0.

1e−10

1e−8

ψ′

ψ′

2 0 −2

0.0 −0.5

0

50 t

100

0

50 t

100

Figure 4.4: Time series of the perturbation of detonation velocity ψ0 for q = 4: a) stable solution with θ = 0.92; b) unstable solution with θ = 0.95. velocity decays such that for large times D = D CJ . In contrast, Figure 4.4b shows the unstable case, in which detonation velocity grows away from the CJ value. In summary, Figure 4.4 demonstrates that Hopf bifurcation occurs in the system (4.3) with production rate (4.4). Figure 4.5 shows corresponding perturbation profiles that are sampled at time t = 50.

111

1e−10

1e−9

u′

1

u′

1

0

0 −1

−20

0

−20 0.0

−2

λ′

λ′

0

1e−11

−10 x

1e−9

−10 x

0

−10 x

0

−0.5 −1.0

−20

−10 x

0

−20

Figure 4.5: Perturbation profiles for q = 4 recorded at time t = 50: a-b) perturbations of velocity u 0 , c-d) perturbations of reaction progress variable λ0 . Subfigures a) and c) are for θ = 0.92 and subfigures b) and d) are for θ = 0.95.

112

4.6.3 Comparison of linear and nonlinear solutions In this subsection, we compare the solutions of the linear and nonlinear problems for verification purposes. If parameters of the problem are such that the ZND solution is unstable and the initial amplitude of the perturbation of the ZND solution is small relative to the ZND solution itself, then for early times both linear and nonlinear solutions should agree with each other. Figure 4.6 shows detonation velocity with respect to time obtained both from linear and nonlinear solvers for heat release q = 4 and activation energy θ = 0.95, with resolutions N1/2 = 40 and N1/2 = 320 for linear and nonlinear solvers, correspondingly, and initial amplitude A 0 = 10−4 . This choice of resolutions and initial amplitude is dictated by the difference in the order of accuracy of the solvers (fifth for linear and second for nonlinear). In Figure 4.6a the detonation velocity is shown for early times, when the amplitude of perturbation is small. As we can see, the solutions of the linear and nonlinear problems are indistinguishable from each other: they both exhibit exponential growth with the same frequency of oscillations. In Figure 4.6b the detonation velocity is shown at later times when the amplitude of perturbation becomes comparable with D CJ . It can be seen that the linear solution continues to grow exponentially while the nonlinear solution “saturates” and sets on a limit cycle. An agreement of the two solutions at early times additionally verifies the correctness of both linear and nonlinear solvers.

4.6.4 Migration of linear spectrum as activation energy is varied Figure 4.7 shows how the linear spectrum changes for q = 4 as activation energy θ increases in a range [0.90;1.15] with step ∆θ = 0.001. Only one mode in the

spectrum was found, and its behavior is as follows. Growth rate increases from −0.081 to 0.493 effectively linearly with a bifurcation to instability at θcrit = 0.937 ±

0.001. Frequency exhibits a more complicated behavior: it initially increases from

2.02

2.2

2.00

2.0

1.98

D

D

113

50

75

100 t

125

150

1.8

150

175

200 t

225

250

Figure 4.6: Comparison of solutions of linearized (solid line) and nonlinear (dashed line) problems: a) at early times when the amplitude of perturbation is small comparing to D CJ = 2 and two solutions are indistinguishable; b) at later times when the amplitude of perturbation is large and two solutions diverge: linear solution continues to grow exponentially while nonlinear solution exhibits a limit cycle. 0.864 at θ = 0.90 and reaches maximum value 0.870 at θ = 0.95, then decreases to 0.738 at θ = 1.15.

4.6.5 Neutral stability Now we turn our attention to the determination of the neutral stability boundary for a wide range of parameters θ and q , that is, the curve in a (θ, q ) plane that separates stable ZND solutions from unstable ones. We generate 256 linearly spaced values of q in a range [0.81, 16] (corresponding to D¯ in a range [0.9, 4]) and for each q we find the critical value of θ such that the growth rate of instability is close to zero with a tolerance 10−3 . Finding θcrit is based on the idea of bisection where a range of θ is recursively divided into subranges unless the above condition on the growth rate is satisfied. Initial interval of search for θ is taken [0.2;5]. Grid resolution used here is N1/2 = 40. Figure 4.8 shows obtained neutral stability curve in the (θ, q ) plane and the frequency of oscillation αim along the curve. The algorithm for finding the critical θ did not converge for three values of q and the outliers were removed from the

Frequency αim

Growth rate αre

114

0.50 0.25 0.00 0.90

0.95 1.00 1.05 1.10 Activation energy θ

1.15

0.90

0.95 1.00 1.05 1.10 Activation energy θ

1.15

0.85 0.80 0.75

Figure 4.7: Migration of the linear spectrum for q = 4 as activation energy θ is varied. Each curve consists of 251 points. Only one mode in the spectrum was found.

Heat release q

115

15 10 5 0

1 2 3 4 Activation energy θ

0.5

1.0 1.5 Frequency αim

Figure 4.8: Neutral stability for the Fickett’s model: a) neutral stability curve in (θ, q ) plane; b) frequency of oscillation αim along the neutral stability curve. Each curve consists of 256 linearly spaced points for q ∈ [0.81;16]. Table 4.1: Critical values of activation energy θcrit and the frequency of oscillation αim on the neutral stability boundary for several values of heat release q . i

q

θcrit

αim

1 2 3 4 5 6

0.81 1.00 2.00 4.00 9.00 16.00

4.625000 3.746094 1.873438 0.937109 0.416504 0.234204

0.391473 0.434969 0.615148 0.869979 1.304970 1.739916

figure. Data demonstrate that when q increases, θcrit decreases while αim increases. For the lowest q considered here (q = 0.81), θcrit = 4.625 and αim = 0.391, and for the largest q (q = 16), θcrit ≈ 0.234 and αim = 1.74. We also provide information about the critical values of activation energy θcrit and the frequency of oscillation αim on the neutral stability boundary for several values of q in Table 4.1.

4.6.6 Comparison with normal-mode analysis In this subsection we compare the results obtained by unsteady linear analysis with the results obtained with normal-mode analysis. We conduct normal-mode

116 analysis using the approach by Lee & Stewart (1990). In this approach the governing equations are linearized about the ZND solution and the boundary-value problem is posed for a complex eigenvalue α with Rankine–Hugoniot conditions on the right boundary (on the shock) and complex-valued boundedness condition H on the left boundary (at the chemical equilibrium). The boundedness condition

expresses the fact that eigenfunctions of the problem should not “blow up” at the equilibrium point of the reaction zone. Full details of the approach applied to the Fickett’s detonation analog are described in Appendix C. We apply the normal-mode analysis for the problem with parameters: heat release q = 4 and activation energy θ = 0.95, which is known from unsteady linear computations described above to be an unstable case. In the first part of the normal-mode analysis we do a “carpet search” to find a minimum of the boundedness function H . First, we plot the carpet by sampling the boundedness function in a range αre ∈ [0, 0.5], αim ∈ [0, 1] with resolution 101 points along each axis. The carpet is shown as a contour plot on Figure 4.9. As can be seen from the contour plot, there is a minimum in the neighborhood of α = 0.03 + 0.85i. Using the obtained values of α as an initial guess, we now solve the boundary value problem using the Newton method for the system of equations

Hre (α) = 0,

Him (α) = 0

in the domain [0;1 − 10−3 ]. We also attempted to work in the domain [0;1 − 10−4 ] with no success due to integration from the shock towards the equilibrium being ill-conditioned, which is a known issue in application of normal-mode analysis to detonation stability (Sharpe, 1997). Figure 4.10 shows a ZND solution versus the spatial coordinate along with the perturbation eigenfunctions u 0 and λ0 . Figure shows that the eigenfunctions

117 log(1 + abs(H))

1 0.9 0.8 0.7

ℑ α

0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

ℜ α

Figure 4.9: Carpet plot showing an approximate location of the eigenvalue α. are very smooth, which is expected as activation energy θ = 0.95 is only slightly larger than the critical value θcrit = 0.937, and vanish as the equilibrium part of the reaction zone is approached, similar to what is found in the analysis of the reactive Euler equations (Lee & Stewart, 1990; Sharpe, 1997; Short & Stewart, 1998). Table 4.2 shows the comparison of growth rates and frequencies of perturbations obtained via unsteady linear computations and normal-mode computations. The results of the two different approaches match to two significant digits for growth rates and to three digits for frequencies, which verifies the correctness of linear unsteady computations.

4.6.7 Nonlinear dynamics of the model Now we turn our attention to the nonlinear dynamics of the model and solve system (4.5) directly using the algorithm described in Subsection 4.5. We aim to

118

ZND solution

4 u λ ω

3 2 1 0 -10

-9

-8

-7

-6

-5

-4

-3

-2

-1

0

-4

-3

-2

-1

0

-4

-3

-2

-1

0

Real part

0.5

0

-0.5

-1 -10

Re u Re λ -9

-8

-7

-6

-5

Imaginary part

2

1

0

-1 -10

Im u Im λ -9

-8

-7

-6

-5 x

Figure 4.10: ZND solution and perturbations obtained via normal-mode analysis with parameters q = 4, θ = 0.95.

Approach

αre

αim

Unsteady linear solver Normal-modes

0.02909 0.02926

0.87041 0.87073

Table 4.2: Comparison of the growth rates and frequencies of perturbations obtained via unsteady linear computations and normal-mode computations for q = 4, θ = 0.95.

119 study the different nonlinear regimes that this model is capable of and plot the bifurcation diagram for this model. We run simulations with the following parameters: N1/2 = 640, final Tf = 1000, q = 4, and θ in a range [0.95, 1.15] with step 0.001 (251 simulations in total). The

outcome of the simulations is the time series of detonation velocity D that are postprocessed by extracting local minima in the time window t ∈ [800;1000]. To overcome the noise in the time series due to numerical artifacts, point i of the time series is considered as a local minimum if it is not larger than the 100 neighboring points on both sides of it. Plotting of the minima against θ demonstrate the occurring bifurcations as θ is increased. Figure 4.11 (top) shows bifurcation diagrams plotted over each other for resolutions N1/2 = 320 and N1/2 = 640. Figure suggests that the bifurcation data are effectively the same in the region θ / 1.07 for which the solution is either stable or on a stable limit cycle. For larger θ the bifurcation diagrams significantly differ from each other: for N1/2 = 320 there is a wide window of period-6 behavior in the neighborhood of θ = 1.10, while it occurs only for θ = 1.088 for bifurcation data obtained for N1/2 = 640. Furthermore, additional minima in the neighborhood of D = 1.90 appear for N1/2 = 640 as θ → 1.15.

Figure 4.11 (bottom) shows the obtained bifurcation diagram for simulations with N1/2 = 640 and N1/2 = 1280. It can be seen that the bifurcation data are effectively the same in the region θ / 1.077 for which the solution is either stable or on a stable limit cycle. For larger θ , where chaotic regimes are realized, the bifurcation diagrams differ but only in minor details, which is expected as chaotic dynamics is very sensitive to every parameter of a simulation, especially a significant parameter such as grid resolution. For both N1/2 = 640 and N1/2 = 1280, period-6 dynamics is observed for θ = 1.089. In conclusion, it can be noted that the nonlinear dynamics of the model is

120 Table 4.3: Bifurcation points and corresponding approximations of the first Feigenbaum constant. n

θn

δn

0 1 2 3 4

0.936 1.001 1.055 1.069 1.072

N/A N/A 1.204 3.857 4.667

quite rich but only after proper study of the sensitivity of the simulations to the grid resolution can one be sure about the obtained results. Feigenbaum (1978; 1979) discovered that many dynamical systems of the form x i +1 = f (x i ; a ), where a is a parameter of the system, exhibit a universality in the

structure of their bifurcation diagrams, namely, that the rate of bifurcations is independent of the form of f . The first Feigenbaum constant δ is defined as a limit of the ratios of the differences of values of a bifurcation parameter in successive period doublings in the logistic map (Strogatz, 2015):

δ = lim δn = 4.669201609103, n→∞

δn =

a n−1 − a n−2 . a n − a n−1

We check whether a sequences of δn determined for the present model (with parameter a being activation energy θ ) from the bifurcation diagram 4.11 converges to δ. Table 4.3 shows the sequences of θn and corresponding δn . Values of θn are determined visually by zooming into the bifurcation diagram and system-

atically choosing the left values of θ in the neighborhoods of bifurcation points. The uncertainty of θn is 0.001. It is difficult to determine θn for n > 4, as density of local minima becomes very dense for θ ≥ 1.072. Nevertheless, the closeness of δ4 = 4.667 to the true value of δ is remarkable.

121

Figure 4.11: Bifurcation diagrams for the Fickett’s model for q = 4 and θ in a range [0.95;1.15]. Two bifurcation diagrams are plotted over each other for two grid resolutions on each subfigure. Top: N1/2 = 320 (black circles), N1/2 = 640 (blue squares). Bottom: N1/2 = 640 (black circles) and N1/2 = 1280 (blue squares).

122

4.6.8 Analysis of stable limit cycles We study time series for θ = {0.95, 1, 1.004, 1.055, 1.065, 1.089} that correspond to stable limit cycles with different numbers of periods as evident from the bifurcation diagram Figure 4.11 (bottom). We also plot phase portraits for each case. Phase portraits require knowledge of the acceleration of the front d D /d t which is approximated with central finite differences: d D i D i +1 − D i −1 = for i = 1, 2, . . . , dt t i +1 − t i −1

(4.29)

and initial acceleration is taken to be zero: d D 0 /d t = 0. To regularize numerical differentiation, time series of D were smoothed using a simple moving average algorithm with window size 11, and then d D /d t computed from (4.29) was smoothed using the same algorithm with window size 5. Figures 4.12–4.13 show nonlinear dynamics for θ = 0.95 and θ = 1, which is a stable limit cycle in these two cases. Figures 4.14–4.15 demonstrate period-2 limit cycles for θ = 1.004 and θ = 1.055, respectively. From these figures we can see that as θ increases, distance between relative maxima and minima increases as well: for θ = 1.004 the ratio of relative maxima is ≈ 1.11, while for θ = 1.055 this ratio is ≈ 1.84. Figure 4.16 shows period-4 nonlinear oscillations for θ = 1.065 and Figure 4.17 shows period-6 nonlinear oscillations for θ = 1.089. Note that as θ increases, the range of detonation acceleration d D /d t increases as well: for weakly unstable cases d D /d t is in the order of unity, while for strongly unstable cases it is in the order of 1000, as can be seen from the Figures 4.15–4.17.

4.6.9 Code verification To confirm the correctness of the obtained results, we verify the solvers used to conduct the numerical simulations in this chapter. It is widely accepted (Salari &

123

Figure 4.12: Analysis of nonlinear dynamics for q = 4, θ = 0.95. Top: time series of detonation velocity D for t ∈ [800;1000]. Bottom: numerically generated phase portrait d D /d t versus D for t ∈ [0;1000].

124

Figure 4.13: Analysis of nonlinear dynamics for q = 4, θ = 1. Top: time series of detonation velocity D for t ∈ [800;1000]. Bottom: numerically generated phase portrait d D /d t versus D for t ∈ [0;1000].

125

Figure 4.14: Analysis of nonlinear dynamics for q = 4, θ = 1.004. Top: time series of detonation velocity D for t ∈ [800;1000]. Bottom: numerically generated phase portrait d D /d t versus D for t ∈ [0;1000].

126

Figure 4.15: Analysis of nonlinear dynamics for q = 4, θ = 1.055. Top: time series of detonation velocity D for t ∈ [800;1000]. Bottom: numerically generated phase portrait d D /d t versus D for t ∈ [0;1000].

127

Figure 4.16: Analysis of nonlinear dynamics for q = 4, θ = 1.065. Top: time series of detonation velocity D for t ∈ [800;1000]. Bottom: numerically generated phase portrait d D /d t versus D for t ∈ [0;1000].

128

Figure 4.17: Analysis of nonlinear dynamics for q = 4, θ = 1.089. Top: time series of detonation velocity D for t ∈ [800;1000]. Bottom: numerically generated phase portrait d D /d t versus D for t ∈ [0;1000].

129 Knupp, 2000; Roy, 2005) that the most stringent test case is a comparison of the observable order of accuracy with the theoretical order of accuracy. We run simulations for q = 4 and θ = 1 with several grid resolutions and compute relative errors: Ei =

kψ0i − ψ0i −1 k kψ0i k

,

i = 2, . . . N ,

(4.30)

where ψ0i is the time series of the perturbation of detonation velocity obtained with i th grid resolution, N is the total number of resolutions, and norms are L 1 , L 2 , and L ∞ . Then we compute the observed order of accuracy:

ri =

log(E i /E i −1 ) , log(∆x i /∆x i −1 )

i = 3, . . . , N ,

where ∆x i is a spatial step for the i th grid resolution. If the observed order of accuracy matches the theoretical fifth order, then one can conclude that the implementation of the solver is correct. Table 4.4 shows the obtained errors and orders of accuracy for the linear solvers where errors were computed with the different norms listed above. It can be seen that the observed order of accuracy is five, as expected, hence, the implementation of the linear solver is correct. The sudden drop of the order of accuracy for the resolution N1/2 = 1280 is due to the domination of the rounding error of the floating-point arithmetic over the truncation errors of the numerical schemes.

4.7 Conclusions In this chapter we applied our approach to linear stability analysis of detonations to the Fickett’s detonation analog with one-step kinetics. The analog exhibits instability to longitudinal perturbations similar to the one-dimensional reactive Euler equations. The results of our computations demonstrate that for sufficiently

130 Table 4.4: Convergence study for the linear solver with the Fickett’s model and parameters q = 4, θ = 1. N1/2 is the resolution per half-reaction zone, E 1 , E 2 , E ∞ are relative errors computed by (4.30) in L 1 , L 2 , and L ∞ norms, respectively; r 1 , r 2 , r ∞ are corresponding observed order of accuracy. Negative values for orders of accuracy are due to the domination of the round-off error of floating-point arithmetic over truncation error of numerical methods used. N1/2

E1

r1

E2

r2

E∞

r∞

20 40 80 160 320 640 1280

N/A 9 × 10−5 3 × 10−6 8 × 10−8 2 × 10−9 6 × 10−11 6 × 10−11

N/A N/A 5.20 5.10 5.06 5.42 -0.10

N/A 9 × 10−5 3 × 10−6 8 × 10−8 2 × 10−9 7 × 10−11 7 × 10−11

N/A N/A 5.20 5.10 5.06 5.25 -0.05

N/A 1 × 10−4 3 × 10−6 8 × 10−8 3 × 10−9 3 × 10−10 4 × 10−10

N/A N/A 5.19 5.10 5.06 3.13 -0.49

high values of activation energy θ the system switches from a stable ZND solution to an unstable one through what appears to be a Hopf bifurcation. We also study the migration of the linear spectrum as θ increases and determine that even far from the critical value of θ , the linear spectrum consists of only one mode with a growth rate increasing linearly as θ increases. We determine the neutral stability boundary that separates stable ZND solutions from unstable ones and found that as heat release increases, critical activation energy decreases while frequency of oscillations increases. We verify our approach by comparing solutions of the linearized system with the solutions obtained from numerical simulations of the nonlinear system and found that both linear and nonlinear solutions agree with each other for early times when the amplitude of perturbation is small relative to the ZND solution. Also, we compare the growth rates and frequencies of perturbations obtained from the results of unsteady linear computations with the results of more conventional normal-mode analysis and find an agreement between the two results. We also investigated the nonlinear dynamics of the Fickett’s model and found

131 that it exhibits rich nonlinear dynamics with stable limit cycles with period-1, period-2, period-4, and period-6 oscillations, and chaotic behavior for regimes far from the stability. Also, we verified the computations by convergence study and found an agreement between observed and theoretical orders of accuracy. In summary, the application of our approach to the Fickett’s model shows that it is a viable alternative to analysis by normal modes that allows avoidance of imposing the radiation condition on the downstream boundary of the reaction zone, while giving an advantage of using highly accurate numerical method to the solution of the linearized system. Moreover, this model possesses a rich set of qualitatively different nonlinear solutions, and can be used as an example of the simple 2 × 2 hyperbolic system of balance laws with unstable steady-state solutions and, also for verification of numerical solvers for such systems.

132

Chapter 5

Conclusions

5.1 Summary In this thesis, we developed and implemented a new method to study the stability of steady solutions of detonation models. Our approach is based on the linearization of the governing equations written in a shock-attached frame and numerical solution of the linearized problem. Resultant time series of detonation velocity are postprocessed using a novel technique of dynamic mode decomposition, which is aimed at identification of modes of the system that grow or decay exponentially in time with fixed frequencies, and these modes are identical to the normal modes in linearized systems. We implement and test our algorithm using a detonation model based on the reactive Euler equations with one irreversible reaction. Comparing our stability results with the results of linear stability analysis via normal-mode approach, which are widely present in the literature for this model, we have found excellent agreement both for growth rates and frequencies of the perturbations of the base flow for all cases considered. Computation of the neutral stability boundary for a wide range of activation energy E and heat release Q for γ = 1.2 also have an agreement with the results of the normal-mode approach. Additionally, we have studied neutral stability boundaries for other values of γ (1.1, 1.3, and 1.4) and have found that for low values of Q the critical value of E decreases as γ increases, while for, large values of Q , the critical values of E increase as γ increases. Also,

133 we have studied the dependence of the frequencies along the neutral stability curves on γ and have found that for low values of Q frequency is insensitive to γ, while, for large values of Q , frequency increases as γ increases. One particularly interesting result was that although for γ ≥ 1.2 neutral stability curves of the fundamental mode are neutral stability boundaries that separate stable base flows from unstable ones, for γ = 1.1 that was not the case: for mid-range of Q the first overtone is more unstable than the fundamental mode, hence, the neutral stability boundary is determined by a union of neutral stability curves for these two modes. The advantage of the present approach over the traditional normal-mode approach is that analytic complications of formulating boundary-value problem, namely, formulation of the boundary condition at the end of the reaction zone, are replaced by straightforward linearization with simple outflow boundary condition. However, difficulty arises in the postprocessing step as time series of detonation velocity inevitably contain noise due to numerical approximation, hence, it may become challenging to identify all normal modes correctly. Having developed and tested this approach to linear stability analysis, we have considered an extension of the model based on the reactive Euler equations with one-step kinetics, to the model that includes two reactions: one is exothermic and one is endothermic. Analysis of the ZND solutions for this model showed that there are two types of admissible ZND solutions, which are the same for the region between the embedded sonic locus and the leading shock, and are different between the end of the reaction zone and the sonic locus. We have studied the effects of endothermicity on stability and have shown that, as endothermicity increases, more and more unstable modes are excited in the stability spectrum, rendering the base ZND flow unstable. We have also shown that stability spectrum is independent of the postsonic part of the ZND solution, which serves as a

134 demonstration of the postulate that intrinsic stability of detonations is determined solely by flow between the sonic locus and the detonation front. Using the particular reaction-rate expression obtained in the asymptotic theory of detonations, we analyze the Fickett’s detonation analogue both for linear and nonlinear stability. We have shown that this model demonstrates rich nonlinear dynamics with multiple bifurcations and chaotic behavior.

5.2 Future Research Work The work presented in this thesis can be extended in the following directions. 1. Application of the present method of linear stability analysis to more complicated models based on reactive Euler equations coupled with thermodynamic relations of calorically perfect gases with multistep reaction mechanisms (Short & Dold, 1996; Short & Sharpe, 2003; Gorchkov et al., 2007). 2. Application of the method to more complicated models based on reactive Euler equations coupled with thermodynamic relations of ideal, but not calorically perfect, gases with multistep reaction mechanisms (Aslam & Powers, 2009; Romick et al., 2012). 3. Application of the method to detonation models for explosives with nonideal equation of state (Short et al., 2006). 4. Application of the method to analysis of detonation waves in different geometries, such as cylindrical detonations (Wescott et al., 2004; Kasimov & Korneev, 2014). 5. Application of the method to detonation models with losses (Semenko et al., 2016; Sow et al., 2017).

135 6. Extension of the method to include perturbations in the direction transverse to the lead shock (Short & Stewart, 1998; Stewart & Kasimov, 2006). 7. Thorough investigation of nonlinear dynamics of the Fickett’s model with the reaction-rate expression (4.4) by using different numerical schemes for approximation of spatial derivatives (including characteristic decomposition) and different time integrators (strong-stability-preserving Runge–Kutta schemes (Ketcheson & Robinson, 2005; Gottlieb et al., 2009; Hadjimichael et al., 2013; Bresten et al., 2016) and implicit-explicit Runge–Kutta schemes (Pareschi & Russo, 2005; Zharovsky et al., 2015; Boscarino & Pareschi, 2017)) to analyze sensitivity of the nonlinear dynamics to numerical methods and to determine which schemes are the most robust in terms of preservation of the solution invariants (u > 0 and 0 ≤ λ ≤ 1) and handling of the secondary shock waves.

136

APPENDICES

137

Appendix A

Details of the linearization of the reactive Euler equations

In this appendix, we give the details of the linearization of the governing equations for the model considered in Chapter 2. Nonconservative form. For the linearization purposes it is more convenient to

use the nonconservative form of the equations (2.8–2.11):

ρ t + (u − D ) u x + ρu x = 0,

(A.1)

1

u t + (u − D ) u x + p x = 0 , ρ µ ¶ ¾ ½µ ¶ 1 1 + (u − D ) = 0, e t + (u − D ) e x + p ρ t ρ x λt + (u − D ) λx = ω.

(A.2) (A.3)

By using Eq. (2.6) for e i with e = e i + u 2 /2 and using Eq. (A.1) and (A.3), we simplify system (A.1–A.3) to the following form:

ρ t + (u − D ) u x + ρu x = 0, u t + (u − D ) u x +

1 ρ

p x = 0,

p t + (u − D ) p x + γpu x = Q (γ − 1)ρω, λt + (u − D ) λx = ω.

(A.4) (A.5) (A.6) (A.7)

138

Linearization We expand all dependent variables about the ZND solution:

¯ + ²ψ0 (t ), D (t ) = D ρ (x, t ) = ρ¯(x ) + ²ρ 0 (x, t ), u (x, t ) = u¯ (x ) + ²u 0 (x, t ), p (x, t ) = p¯(x ) + ²p 0 (x, t ),

¯ (x ) + ²λ0 (x, t ), λ(x, t ) = λ and the dependent function ω = ω(ρ, p, λ) is expanded using Taylor series as ω ¯+ ¡ ¢ ²ω0 = ω ¯ +² ω ¯ ρ ρ0 + ω ¯ p p0 + ω ¯ λ λ0 .

To linearize, we substitute the expansions into the equations (A.4–A.7) and leave only O (²) terms. Continuity equation. Linearized equation (A.4) is ¡ ¢ du¯ 0 dρ¯ 0 ¯ ρ 0x + ρu ρ 0t + u¯ − D ¯ x0 + ρ + (u − ψ 0 ) = 0 . dx dx Momentum equation. Linearized equation (A.5) is ¡ ¢ 1 1 dp¯ 0 du¯ 0 ¯ u x0 + p x0 − u t0 + u¯ − D ρ + (u − ψ 0 ) = 0 . 2 ρ¯ dx ρ¯ dx Energy equation. Linearized equation (A.6) is

p t0

+ γpu ¯ x0

µ

¯ p x0 −C ρ¯ω + u¯ − D ¯ρ +ω ¯ ρ0 + γ ¡

¢

¡

¢

where we introduced C = Q (γ − 1).

du¯ dp¯ 0 −C ρ¯ω ¯ p p 0 −C ρ¯ω ¯ λ λ0 + (u − ψ 0 ) = 0 , dx dx ¶

139 Reaction equation. Linearized equation (A.7) is ¡ ¢ dλ¯ 0 ¯ λ0x − r¯ρ ρ 0 − r¯p p 0 − r¯λ λ0 + λ0t + u¯ − D (u − ψ 0 ) = 0 . dx

Subsequently the vector form (2.22) of these linearized equation is easily obtained. ¯ + ²M 0 (t ), we obtain the following Adding expansion for the mass flux: M (t ) = M linearized Rankine–Hugoniot conditions (2.15) through the Taylor-series expansion:

p s0 =

¯ 0 4v a M M, γ+1

Us0 = −

4γp a M 0, ¯3 (γ + 1) M

¯ v s0 + v¯s M 0 . v s0 = M

Then, expanding ρ = 1/v in Taylor series, we obtain: 0

v s0

ρs = − 2 . v¯s

Perturbation of the detonation velocity ψ0 is related to the perturbation of mass flux M 0 by ψ0 = −v a M 0 . For λ we have: λ¯ + ²λ0 = 0, thus, λ0 = 0.

140

Appendix B

Conversion between different scales used in the detonation stability research

Different researchers working in detonation stability research rescale the governing equations in different ways with the most widely used being the upstream Erpenbeck scales (Erpenbeck, 1964), postshock scales of Lee & Stewart (1990), and scales of Sharpe (1997). Hence we need to convert between these scales to compare our results (obtained with the Erpenbeck scales) to theirs for the comparison in Chapter 2. Table B.1 gives the definitions of nondimensional quantities in scales of Erpenbeck, Sharpe, and Lee and Stewart. In this table we use the following notation: the “tilde” symbol over a quantity denotes that the quantity is dimensional; subscript “a” is a short form of “ambient” and denotes the quantity ahead of the detonation wave; subscript “c” is for “characteristic”; D˜¯ is dimensional ZND (constant) detonation velocity; c˜a is dimensional speed of sound in the medium ahead of the detonation wave; R˜ is a gas constant which is given by a ratio of the universal gas constant to the molar mass of the species A or B (recall that we consider a reaction in which the molar mass of both species is the same, has gas constant is the same); dimensional quantities E˜act and Q˜ are intensive thermodynamic quantities per unit of mass. Notice that in all scales characteristic length is defined as the length at which the mass fraction of product B increases from 0 to 1/2 (half-reaction length). Hence, it follows from this definition that the nondimensional preexponential

141 Table B.1: Scales used in the detonation stability research. Quantity

Erpenbeck

Density, ρ

ρ = ρ˜a

Pressure, p

p = p˜a

Velocity, u

u = c˜ /u˜pγ

ρ˜

ρ˜

p=

˜

Q=

Length, x

c¯s

˜ R˜ T¯s

˜

a

E = ˜E2 ¯ D

Q˜ R˜ T˜a

Q = ˜2 ¯ D

E act = E act ˜



Q=

Q˜ ˜ R T˜¯s

x = ˜x˜

`c

x = ˜x˜

t˜ t˜c

t˜ t˜c

t = t˜t˜

t=

k = t˜c k˜ R 1/2 U˜ 0

`c c

k = t˜c k˜ R 1/2 U˜

¯

˜¯ dλ ω

k = t˜c k˜

¯

˜¯ dλ ω

0

¡ p ¢ `˜c / c˜a / γ

Characteristic time, t˜c

u = u˜˜

`c

t=

Characteristic length, l˜c



p= ˜ p¯s

x = ˜x˜

Time, t Preexponential factor, k



˜¯ 2 ρ˜a D ¯ D

E = R˜ET˜

Heat release, Q

ρ= ˜ ρ¯s

u = u˜˜

a

Lee & Stewart ρ˜

ρ = ρ˜a



Activation energy, E

Sharpe

˜¯ `˜c /D

R 1/2 U˜ 0

˜¯ ω

dλ¯

`˜c /c˜¯s

factor k is computed using formula Z 1/2 k=

0



¯ /p¯) (1 − λ)exp(−ρE

dλ¯

in all three scales, although, obviously, the numerical value of k is different in different scales due to the different velocity and time scales. Below we use subscripts “E”, “S”, and “LS” to denote nondimensional quantities in Erpenbeck, Sharpe, and Lee&Stewart scales, respectively.

Note on the Erpenbeck scales Computing quantities in Erpenbeck scales is performed as follows. Ambient quantities ρ a and p a are both unity in this scale. Then the nondimensional sound p speed is c a = γ. If Q is given in the Erpenbeck scales, then it is easy to compute D¯

and all other ZND quantities.

142 Table B.2: Conversion factors used to convert from the Erpenbeck scales to the scales of Sharpe (1997) and Lee & Stewart (1990). Conversion factor

Quantity

to Sharpe scales

to Lee & Stewart scales

Density

1

ρ a /ρ¯s

Pressure

¡ 2¢ ¯ c a2 / γD ¡p ¢ ¯ c a / γD ¢ ¡ ¯2 c a2 / γD ¡ 2¢ ¯ c a2 / γD

p a /p¯s ¡p ¢ c a / γc¯s

Velocity Activation energy Heat release Length Time Preexponential factor

c a2 /c¯s2 c a2 /c¯s2

1

1

p ¯ γD /c a ¡p ¢ ¯ c a / γD

¡p ¢ γc¯s /c a ¡p ¢ c a / γc¯s

Conversion from the Erpenbeck scales to other scales Using Table B.1 it is easy to find conversion factors from the Erpenbeck scales to the scales of Sharpe and Lee&Stewart. Conversion factors are ratios of the characteristic scales of the corresponding quantity and are shown in Table B.2. The conversion factors are computed by using quantities in the Erpenbeck scales p

(for example, c a in the table is γ in the Erpenbeck scales). We also need a conversion factor from the scales of Abouseif & Toong (1982) to the Erpenbeck scales, because Lee & Stewart (1990) provide the eigenvalues in table 1 of their paper in these scales. As eigenvalues have dimension of 1/time, then to convert from Erpenbeck scales to Abouseif & Toong, the following conversion factor should be used:

R 1/2 1 ¯ 0 kT ¯ /p¯) dλ (1−λ)exp(−ρE =R , 1/2 U¯ kE dλ¯ 0

¯ /p¯) (1−λ)exp(−ρE

where subscript “T” denotes a quantity in the scales of Abouseif and Toong.

143

Appendix C

Details of the normal-mode analysis for the Fickett’s detonation analog

C.1 Linearization with normal-modes We linearize the Fickett’s detonation analog using normal-modes expansions ¯ t + ψ ( t ), x s` = D ¯ +ψ ˙ ( t ), D (t ) = D u (x, t ) = u¯ (x ) + u 0 (x )exp(αt ) ,

¯ (x ) + λ0 (x )exp(αt ), λ(x, t ) = λ ψ(t ) = exp(αt ),

˙ (t ) = α exp(αt ), ψ

which leads to the system of equations

¯ αz 0 + A

dz 0 ¯ 0 + C z − αb¯ = 0, dx

(C.1)

where   0 u  0 z =  , λ0



¯ u¯ − D

¯ = A

0



σ  , ¯ −D



du¯ dx



0   C¯ =  , −ω ¯ u −ω ¯λ







du¯  dx  =  dλ¯ dx

and it is understood that matrices A and C , as well as vector b , implicitly depend on spatial coordinate x .

144 Then the system can be written in the conventional form for the system of first-order ODEs: £ ¡ ¢ ¤ dz 0 ¯ −1 − αI + C¯ z 0 + αb¯ , =A dx

where ¯ A

−1

=

1 ¯ 2 − u¯ D¯ D

 ¯ −D 

0

(C.2)

 −σ  . ¯ u¯ − D

Linearized Rankine–Hugoniot conditions are:

u 0 = 2α,

λ0 = 0.

(C.3)

The system (C.2) must be integrated using the initial value (C.3) towards the equilibrium of the reaction zone where the condition derived in Section C.3 is imposed. This condition is called boundedness condition and it can be satisfied only for specific values of α. As parameter α in the system is not known in advance, we must take a guess on it and then iterate until the boundedness condition is satisfied. That is, the problem is a boundary-value problem for a complex-valued eigenvalue α, and the method of the solution described above, is a shooting method.

C.2 Formulation of stability problem in reaction variable The formulation described in previous section is not very useful if one wants to apply adaptive ODE solver to the problem, otherwise the ZND solution must be integrated from the shock for each internal step of the ODE solver. Here, we transform the spatial coordinate to the coordinate of the reaction progress variable ¯ (x ) of the ZND solution. Then the ZND solution in this case is reduced to s=λ the evaluation of several algebraic expressions entering the matrices A and C in

145 addition to vector b . Precisely, by using transformation ds ω ¯ =− ¯ dx D



d ω ¯ d =− , ¯ ds dx D

we have ¯ d dz 0 D =− , ds ω ¯ dx and hence ¯ A

¯ dz 0 ¯£ ¡ ¢ ¤ dz 0 D D ¯ =− A =− − αI + C¯ z 0 + αb¯ ds ω ¯ dx ω ¯

and then system C.2 becomes ¯ −1 £ ¡ ¢ ¤ dz 0 D ¯ =− A − αI + C¯ z 0 + αb¯ . ds ω ¯

(C.4)

C.3 Derivation of the boundedness condition In Stewart & Kasimov (2005) the method of derivation of the boundedness condition which is simpler than the method of Lee & Stewart (1990) is provided. We apply this method of derivation here to derive the boundedness condition for the normal-mode analysis of the Fickett’s detonation analog. Precisely, we take the governing equations:

u t + uu x + σλx = 0, λt = ω,

and, multiplying the second equation by σ/u , we add these equations to obtain

146 the characteristic equations:

(u t + uu x ) +

σ σ (λt + uλx ) = ω, u u λt = ω,

that is dx du σ dλ σ + = ω on = u, dt u dt u dt dλ dx =ω on = 0. dt dt These equation can be simplified by multiplying the first equation by u :

u

dλ dx du +σ = σω on = u, dt dt dt dλ dx =ω on = 0. dt dt

(C.5a)

Boundedness condition can be derived from the first characteristic equation by linearizing it. We consider that the flow at the point of chemical equilibrium can be represented as a small perturbation of ZND flow with parameter ² defining the magnitude of the perturbation:

x ∗ = x¯∗ + ²x 0 (t ), z ∗ = z¯ (x ∗ ) + ²z 0 (x ∗ , t ),

where “overbar” denotes ZND quantities and “prime” denotes perturbations, subscript “*” denotes evaluation at the point of chemical equilibrium. Due to the perturbation of the location of the equilibrium point, z¯ (x ∗ ) can be expanded into

147 the O (1) and O (²) parts: d¯z ¯¯ z¯ (x ∗ ) = z¯ x¯∗ + ²x (t ) = z¯ (x¯∗ ) + ² x 0 ( t ). dx ¯x=x¯∗ 0

¡

¯

¢

Now we linearize the first characteristic equation keeping only O (²) terms: µ ¶ ¡ ¢ dλ0 du¯ dλ¯ dx 0 du 0 +σ + u¯ +σ =σ ω ¯ u u0 + ω ¯ λ λ0 , u¯ dt dt dx dx dt

where all quantities are evaluated at the point x¯∗ . Now recall that due to the expression for ω, it follows that ω ¯ u,

du¯ dλ¯ , go to zero as x → x¯∗ . dx dx

Therefore, the linearized equation simplifies to du 0 dλ0 u¯ +σ = σω ¯ λ λ0 , dt dt and, finally, we assume the normal-modes form of the perturbations

z 0 = z˜ (x )exp(αt ) ,

to obtain the boundedness condition ¡ ¢ ˜ − σω α u¯ u˜ + σλ ¯ λ λ˜ = 0.

C.4 Formulation of the stability problem Combining material described in the previous sections of this Appendix together, we formulate the full stability problem for the Fickett’s detonation analog. The full problem is a boundary-value problem for complex-valued unknowns z 0 and

148 eigenvalue α for the system of ODEs: ¯ −1 £ ¡ ¢ ¤ dz 0 D ¯ =− A − αI + C¯ z 0 + αb¯ for 0 ≤ s < 1, ds ω ¯

(C.6)

subject to the shock boundary conditions  2α z 0 =   at s = 0 

(C.7)

0

and the boundedness boundary conditions

¡ ¢ ¯ λ z 20 → 0 as s → 1. H (α) = α u¯ z 10 + σ z 20 − σω

(C.8)

To solve the problem numerically we use the shooting method. We define the computational domain as [0;1 − tols ], where tols is a tolerance determining how close to the chemical equilibrium the reaction is. As α is not known in advance, for any given set of problem parameters q and θ we specify a grid of α in some range, solve the initial-value problem (C.6)–(C.7) and compute H (α) (C.8). By plotting resulting values of log(1 + |H (α)|) on this grid, we can find the approximate location α∗ of the minimum of H (α). Then using α∗ as an initial guess, we solve the optimization problem aiming to zero |H (α)|. We implement the computations using MATLAB which allows us to separate real and imaginary parts of complex-valued expressions automatically. We use solver ode45 with relative and absolute tolerances set to 10−10 to solve the initial-value problem for system of ODEs and fsolve with tolerance 10−15 on the optimization step to find roots of |H (α)|.

149

Appendix D

Software and computational environment

Most of the computational work in this thesis was conducted using the foundational scientific libraries for the Python1 programming language, which are NumPy2 , SciPy3 , and Matplotlib4 . For both linear and nonlinear solvers subroutines approximating spatial derivatives and subroutines solving the Riemann problem for the Fickett’s model for nonlinear solver, were written in C and wrapped by Cython5 using GNU C compiler6 for compilation. It was crucial to use a distribution for scientific Python which was compiled against high-performance linear-algebra libraries, such as Intel Math Kernel Library7 , because our postprocessing algorithm relies on the linear-algebra subroutines. Our distribution was Anaconda Python Distribution8 of versions 4.2 and 5.0, while GCC compiler version was 6.2. The normal-mode analysis for the Fickett’s model was conducted using MATLAB9 . As a computational environment, we used two workstations working under Ubuntu10 operating system. Multiparametric studies, such as computation of 1 https://www.python.org/ 2 http://www.numpy.org/ 3 http://www.scipy.org/ 4 https://matplotlib.org 5 http://cython.org/ 6 https://gcc.gnu.org/ 7 https://software.intel.com/en-us/mkl 8 https://www.anaconda.com/ 9 https://www.mathworks.com/products/matlab.html 10 https://www.ubuntu.com/

150 neutral stability boundaries, were conducted using Shaheen 2 supercomputer of the KAUST Supercomputing Lab. The reproducibility package for the results from Chapters 2 and 3 is available at https://github.com/dmitry-kabanov/euler1d-reproducibility. The reproducibility package for the results from Chapter 4 is available at https://github.com/dmitry-kabanov/fickettmodel-reproducibility.

151

REFERENCES

Abouseif, G. E. & Toong, T. Y., 1982. Theory of unstable one-dimensional detonations. Combustion and Flame 45, 67–94. doi: 10.1016/0010-2180(82)90034-7. Aslam, T. D. & Powers, J. M., 2009. The dynamics of unsteady detonation in ozone. 47th AIAA Aerospace Sciences Meeting pages 5–8. doi: 10.2514/6.2009-632. Bagheri, S., 2014. Effects of weak noise on oscillating flows: Linking quality factor, Floquet modes, and Koopman spectrum. Physics of Fluids 26, 094104. doi: 10.1063/1.4895898. Berthelot, M. & Vieille, P., 1881. Sur la vitesse des propagation des phenomenes explosifs dans le gaz [on the velocity of propagation of explosive processes in gases]. Comptes rendus de l’Académie des Sciences 93, 18–21. Boscarino, S. & Pareschi, L., 2017. On the asymptotic properties of IMEX Runge–Kutta schemes for hyperbolic balance laws. Journal of Computational and Applied Mathematics 316, 60–73. doi: 10.1016/j.cam.2016.08.027. Bourlioux, A., Majda, A. J. & Roytburd, V., 1991. Theoretical and numerical structure for unstable one-dimensional detonations. SIAM Journal on Applied Mathematics 51, 303–343. doi: 10.1137/0151016. Bresten, C., Gottlieb, S., Grant, Z., Higgs, D., Ketcheson, D. I. & Németh, A., 2016. Explicit strong stability preserving multistep Runge–Kutta methods. Mathematics of Computation 86, 747–769. doi: 10.1090/mcom/3115. Brown, P. N., Byrne, G. D. & Hindmarsh, A. C., 1989. VODE: A variable-coefficient ODE solver. SIAM Journal on Scientific and Statistical Computing 10, 1038–1051. doi: 10.1137/0910062. Brunton, B. W., Johnson, L. A., Ojemann, J. G. & Kutz, J. N., 2016. Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition. Journal of Neuroscience Methods 258, 1–15. doi: 10.1016/j.jneumeth.2015.10.010.

152 Brunton, S. L., Brunton, B. W., Proctor, J. L., Kaiser, E. & Kutz, J. N., 2017. Chaos as an intermittently forced linear system. Nature Communications 8. doi: 10.1038/s41467-017-00030-8. Buckmaster, J. & Neves, J., 1988. One-dimensional detonation stability: The spectrum for infinite activation energy. Physics of Fluids 31, 3571–3576. doi: 10.1063/1.866874. Budiši´c, M., Mohr, R. & Mezi´c, I., 2012. Applied koopmanism. Chaos: An Interdisciplinary Journal of Nonlinear Science 22, 047510. doi: 10.1063/1.4772195. Burgers, J. M., 1948. A mathematical model illustrating the theory of turbulence. Advances in Applied Mechanics 1, 171–199. doi: 10.1016/S0065-2156(08)701005. Cai, W., Oh, W. & Zhu, Y., 1996. Direct numerical calculations of a neutral stability curve for one-dimensional detonations. SIAM Journal on Scientific Computing 17, 814–829. doi: 10.1137/0917053. Chapman, D. L., 1899. On the rate of explosion in gases. Philosophical Magazine Series 5 47, 90–104. doi: 10.1080/14786449908621243. Chen, K. K., Tu, J. H. & Rowley, C. W., 2012. Variants of dynamic mode decomposition: Boundary condition, Koopman, and Fourier analyses. Journal of Nonlinear Science 22, 887–915. doi: 10.1007/s00332-012-9130-9. Clarke, J. F., Karni, S., Quirk, J. J., Roe, P. L., Simmonds, L. G. & Toro, E. F., 1993a. Numerical computation of two-dimensional unsteady detonation waves in high energy solids. Journal of Computational Physics 106, 215–233. doi: 10.1016/S00219991(83)71104-6. Clarke, J. F., Roe, P. L., Simmonds, L. G. & Toro, E. F., 1993b. Numerical studies of a detonation analogue. Journal of Energetic Materials 7, 265–296. doi: 10.1080/07370658908014900. Colella, P., Majda, A. & Roytburd, V., 1986. Theoretical and numerical structure for reacting shock waves. SIAM Journal on Scientific and Statistical Computing 7, 1059–1080. doi: 10.1137/0907073. Colella, P. & Woodward, P. R., 1984. The piecewise parabolic method (PPM) for gas-dynamical simulations. Journal of Computational Physics 54, 174–201. doi: 10.1016/0021-9991(84)90143-8.

153 Courant, R., Friedrichs, K. & Lewy, H., 1928. Über die partiellen Differenzengleichungen der mathematischen Physik. Mathematische Annalen 100, 32–74. doi: 10.1007/BF01448839. Courant, R. & Friedrichs, K. O., 1948. Supersonic Flow and Shock Waves. Interscience. Dawson, S. T. M., Hemati, M. S., Williams, M. O. & Rowley, C. W., 2016. Characterizing and correcting for the effect of sensor noise in the dynamic mode decomposition. Experiments in Fluids 57, 42. doi: 10.1007/s00348-016-2127-7. Döring, W., 1943. Über Detonationsvorgang in Gasen [On detonation processes in gases]. Annalen der Physik 43, 421–436. doi: 10.1002/andp.19434350605. Drazin, P. G., 2002. Introduction to Hydrodynamic Stability. Cambridge University Press. doi: 10.1017/CBO9780511809064. Duke, D., Soria, J. & Honnery, D., 2011. An error analysis of the dynamic mode decomposition. Experiments in Fluids 52, 529–542. doi: 10.1007/s00348-0111235-7. Erichson, N. B. & Donovan, C., 2016. Randomized low-rank dynamic mode decomposition for motion detection. Computer Vision and Image Understanding 146, 40–50. doi: 10.1016/j.cviu.2016.02.005. Erpenbeck, J. J., 1962. Stability of steady-state equilibrium detonations. Physics of Fluids (1958–1988) 5, 604–614. doi: 10.1063/1.1706664. Erpenbeck, J. J., 1964. Stability of idealized one-reaction detonations. Physics of Fluids (1958–1988) 7, 684–696. doi: 10.1063/1.1711269. Faria, L. M. Qualitative and asymptotic theory of detonations. Ph.D. thesis, King Abdullah University of Science and Technology, 2014. http://hdl.handle.net/ 10754/335798. Faria, L. M., Kasimov, A. R. & Rosales, R. R., 2015. Theory of weakly nonlinear self-sustained detonations. Journal of Fluid Mechanics 784, 163–198. doi: 10.1017/jfm.2015.577. Fedkiw, R. P., Merriman, B. & Osher, S., 1997. High accuracy numerical methods for thermally perfect gas flows with chemistry. Journal of Computational Physics 132, 175–190. doi: 10.1006/jcph.1996.5622.

154 Feigenbaum, M. J., 1978. Quantitative universality for a class of nonlinear transformations. Journal of Statistical Physics 19, 25–52. doi: 10.1007/bf01020332. Feigenbaum, M. J., 1979. The universal metric properties of nonlinear transformations. Journal of Statistical Physics 21, 669–706. doi: 10.1007/bf01107909. Fickett, W., 1979. Detonation in miniature. American Journal of Physics 47, 1050– 1059. doi: 10.1119/1.11973. Fickett, W., 1984. Shock initiation of detonation in a dilute explosive. Physics of Fluids 27, 94–105. doi: 10.1063/1.864493. Fickett, W., 1985a. Introduction to Detonation Theory. University of California Press. Fickett, W., 1985b. Stability of the square-wave detonation in a model system. Physica D: Nonlinear Phenomena 16, 358–370. doi: 10.1016/01672789(85)90014-4. Fickett, W. & Davis, W. C., 2000. Detonation: Theory and Experiment. Dover Publications. [Reprint of the University of California Press, 1979]. Fickett, W., Jacobson, J. D. & Schott, G. L., 1972. Calculated pulsating onedimensional detonations with induction-zone kinetics. AIAA Journal 10, 514– 516. doi: 10.2514/3.50130. Fickett, W. & Wood, W. W., 1966. Flow calculations for pulsating one-dimensional detonations. Physics of Fluids (1958-1988) 9, 903–916. doi: 10.1063/1.1761791. Godlewski, E. & Raviart, P.-A., 1999. The linearized stability of solutions of nonlinear hyperbolic systems of conservation laws: A general numerical approach. Mathematics and Computers in Simulation 50, 77–95. doi: 10.1016/S03784754(99)00062-2. Golyandina, N. & Zhigljavsky, A., 2013. Singular spectrum analysis for time series. Springer. doi: 10.1007/978-3-642-34913-3. Gorchkov, V., Kiyanda, C. B., Short, M. & Quirk, J. J., 2007. A detonation stability formulation for arbitrary equations of state and multi-step reaction mechanisms. Proceedings of the Combustion Institute 31, 2397–2405. doi: 10.1016/j.proci.2006.07.219.

155 Gottlieb, S., Ketcheson, D. I. & Shu, C.-W., 2009. High order strong stability preserving time discretizations. Journal of Scientific Computing 38, 251–289. doi: 10.1007/s10915-008-9239-z. Hadjimichael, Y., Macdonald, C. B., Ketcheson, D. I. & Verner, J. H., 2013. Strong stability preserving explicit Runge–Kutta methods of maximal effective order. SIAM Journal on Numerical Analysis 51, 2149–2165. doi: 10.1137/120884201. Hairer, E., Wanner, G. & Nørsett, S. P., 1993. Solving Ordinary Differential Equations I: Nonstiff Problems. Springer. doi: 10.1007/978-3-540-78862-1. Harten, A., Engquist, B., Osher, S. & Chakravarthy, S. R., 1987. Uniformly high order accurate essentially non-oscillatory schemes, III. Journal of Computational Physics 71, 231–303. doi: 10.1016/0021-9991(87)90031-3. He, L., 1999. On the two-dimensional instability of square-wave detonations. Combustion Theory and Modelling pages 37–41. doi: 10.1088/1364-7830/3/2/306. Hemati, M. S., Williams, M. O. & Rowley, C. W., 2014. Dynamic mode decomposition for large and streaming datasets. Physics of Fluids 26, 111701. doi: 10.1063/1.4901016. Henrick, A. K. Shock-fitted Numerical Solutions of One- and Two-dimensional Detonation. Phd thesis, University of Notre Dame, 2008. http://etd.nd.edu/ ETD-db/theses/available/etd-04172008-123816/. Henrick, A. K., Aslam, T. D. & Powers, J. M., 2006. Simulations of pulsating onedimensional detonations with true fifth order accuracy. Journal of Computational Physics 213, 311–329. doi: 10.1016/j.jcp.2005.08.013. Higgins, A., 2012. Steady one-dimensional detonations. In Zhang, F., ed., Shock Waves Science and Technology Reference Library, Volume 6: Detonation Dynamics, pages 33–105. Springer. doi: 10.1007/978-3-642-22967-1_2. Hildebrand, F. B., 1987. Introduction to numerical analysis. Dover Publications. Hwang, P., Fedkiw, R. P., Merriman, B., Aslam, T. D., Karagozian, A. R. & Osher, S. J., 2000. Numerical resolution of pulsating detonation waves. Combustion Theory and Modelling 4, 217–240. doi: 10.1088/1364-7830/4/3/301.

156 Jiang, G.-S. & Peng, D., 2000. Weighted ENO schemes for Hamilton–Jacobi equations. SIAM Journal on Scientific Computing 21, 2126–2143. doi: 10.1137/S106482759732455X. Jiang, G.-S. & Shu, C.-W., 1996. Efficient implementation of weighted ENO schemes. Journal of Computational Physics 126, 202–228. doi: 10.1006/jcph.1996.0130. Jouguet, E., 1905. Sur la propagation des réactions chimiques dans les gaz [On the propagation of chemical reactions in gases]. Journal des Mathématiques Pures et Appliquées 6, 347–425. Jovanovi´c, M. R., Schmid, P. J. & Nichols, J. W., 2014. Sparsity-promoting dynamic mode decomposition. Physics of Fluids 26, no. 2, 024103. doi: 10.1063/1.4863670. Kasimov, A. R. Theory of instability and nonlinear evolution of self-sustained detonation waves. Ph.D. thesis, University of Illinois at Urbana-Champaign, 2004. http://hdl.handle.net/2142/87720. Kasimov, A. R., Faria, L. M. & Rosales, R. R., 2013. Model for shock wave chaos. Physical Review Letters 110, 104104. doi: 10.1103/PhysRevLett.110.104104. Kasimov, A. R. & Korneev, S. V., 2014. Detonation in supersonic radial outflow. Journal of Fluid Mechanics 760, 313–341. doi: 10.1017/jfm.2014.598. Kasimov, A. R. & Stewart, D. S., 2002. Spinning instability of gaseous detonations. Journal of Fluid Mechanics 466, 179–203. doi: 10.1017/S0022112002001192. Kasimov, A. R. & Stewart, D. S., 2004. On the dynamics of self-sustained onedimensional detonations: A numerical study in the shock-attached frame. Physics of Fluids 16, 3566–3578. doi: 10.1063/1.1776531. Kasimov, A. R. & Stewart, D. S., 2005. Asymptotic theory of evolution and failure of self-sustained detonations. Journal of Fluid Mechanics 525, 161–192. doi: 10.1017/S0022112004002599. Ketcheson, D. I. & Robinson, A. C., 2005. On the practical importance of the SSP property for Runge-Kutta time integrators for some common Godunov-type schemes. International Journal for Numerical Methods in Fluids 48, 271–303. doi: 10.1002/fld.837.

157 Kirkwood, J. G. & Wood, W. W., 1954. Structure of a steady-state plane detonation wave with finite reaction rate. The Journal of Chemical Physics 22, 1915–1919. doi: 10.1063/1.1739939. Koopman, B. O., 1931. Hamiltonian systems and transformation in Hilbert space. Proceedings of the National Academy of Sciences 17, 315–318. http://www.pnas. org/content/17/5/315. Kou, J. & Zhang, W., 2017. An improved criterion to select dominant modes from dynamic mode decomposition. European Journal of Mechanics - B/Fluids 62, 109–129. doi: 10.1016/j.euromechflu.2016.11.015. Kunert, J. M., Proctor, J. L., Brunton, S. L. & Kutz, J. N., 2017. Spatiotemporal feedback and network structure drive and encode Caenorhabditis elegans locomotion. PLOS Computational Biology 13, e1005303. doi: 10.1371/journal.pcbi.1005303. Kutz, J. N., Brunton, S. L., Brunton, B. W. & Proctor, J. L., 2016a. Dynamic Mode Decomposition. Society for Industrial and Applied Mathematics. doi: 10.1137/1.9781611974508. Kutz, J. N., Fu, X. & Brunton, S. L., 2016b. Multiresolution dynamic mode decomposition. SIAM Journal on Applied Dynamical Systems 15, 713–735. doi: 10.1137/15M1023543. Laney, C. B., 1998. Computational Gasdynamics. Cambridge University Press. Lee, H. I. & Stewart, D. S., 1990. Calculation of linear detonation instability: Onedimensional instability of plane detonation. Journal of Fluid Mechanics 216, 103–132. doi: 10.1017/S0022112090000362. Lee, J. H. S., 2008. The Detonation Phenomenon. Cambridge University Press. doi: 10.1017/CBO9780511754708. LeVeque, R. J., 2002. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press. doi: 10.1017/CBO9780511791253. LeVeque, R. J., 2007. Finite Difference Methods for Ordinary and Partial Differential Equations. Society for Industrial and Applied Mathematics. doi: 10.1137/1.9780898717839.

158 Majda, A., 1981. A qualitative model for dynamic combustion. SIAM Journal on Applied Mathematics 41, 70–93. doi: 10.1137/0141006. Mallard, E. & Chatelier, H. L., 1881. Sur la vitesse de propagation de l’inflammation dans les melanges explosifs. Comptes rendus de l’Académie des Sciences 93, 145– 148. Massa, L., Kumar, R. & Ravindran, P., 2012. Dynamic mode decomposition analysis of detonation waves. Physics of Fluids 24, 066101. doi: 10.1063/1.4727715. Mezi´c, I., 2005. Spectral properties of dynamical systems, model reduction and decompositions. Nonlinear Dynamics 41, 309–325. doi: 10.1007/s11071-0052824-x. Mezi´c, I., 2013. Analysis of fluid flows via spectral properties of the Koopman operator. Annual Review of Fluid Mechanics 45, 357–378. doi: 10.1146/annurevfluid-011212-140652. Ng, H. D., Higgins, A. J., Kiyanda, C. B., Radulescu, M. I., Lee, J. H. S., Bates, K. R. & Nikiforakis, N., 2005a. Nonlinear dynamics and chaos analysis of onedimensional pulsating detonations. Combustion Theory and Modelling 9, 159– 170. doi: 10.1080/13647830500098357. Ng, H. D., Radulescu, M. I., Higgins, A. J., Nikiforakis, N. & Lee, J. H. S., 2005b. Numerical investigation of the instability for one-dimensional Chapman–Jouguet detonations with chain-branching kinetics. Combustion Theory and Modelling 9, 385–401. doi: 10.1080/13647830500307758. Ng, H. D. & Zhang, F., 2012. Detonation instability. In Zhang, F., ed., Shock Waves Science and Technology Reference Library, Volume 6: Detonation Dynamics, pages 107–212. Springer. doi: 10.1007/978-3-642-22967-1_3. Noack, B. R., Stankiewicz, W., Morzynski, ´ M. & Schmid, P. J., 2016. Recursive dynamic mode decomposition of transient and post-transient wake flows. Journal of Fluid Mechanics 809, 843–872. doi: 10.1017/jfm.2016.678. Pareschi, L. & Russo, G., 2005. Implicit–explicit Runge–Kutta schemes and applications to hyperbolic systems with relaxation. Journal of Scientific Computing 25, 129–155. doi: 10.1007/BF02728986.

159 Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P., 2007. Numerical Recipes: The Art of Scientific Computing. Cambridge University Press, 3 edition. Proctor, J. L., Brunton, S. L. & Kutz, J. N., 2016. Dynamic mode decomposition with control. SIAM Journal on Applied Dynamical Systems 15, 142–161. doi: 10.1137/15M1013857. Radulescu, M. I. & Tang, J., 2011. Nonlinear dynamics of self-sustained supersonic reaction waves: Fickett’s detonation analogue. Physical Review Letters 107, 164503. doi: 10.1103/PhysRevLett.107.164503. Romick, C. M., Dame, N. & Powers, J. M., 2012. Dynamics of unsteady inviscid and viscous detonations in hydrogen-air. 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition pages 1–11. doi: 10.2514/6.2012-981. Rowley, C. W. & Dawson, S. T., 2017. Model reduction for flow analysis and control. Annual Review of Fluid Mechanics 49, 387–417. doi: 10.1146/annurev-fluid010816-060042. Rowley, C. W., Mezi´c, I., Bagheri, S., Schlatter, P. & Henningson, D. S., 2009. Spectral analysis of nonlinear flows. Journal of Fluid Mechanics 641, 115. doi: 10.1017/s0022112009992059. Roy, C. J., 2005. Review of code and solution verification procedures for computational simulation. Journal of Computational Physics 205, 131–156. doi: 10.1016/j.jcp.2004.10.036. Roy, C. J., Nelson, C. C., Smith, T. M. & Ober, C. C., 2004. Verification of Euler/Navier–Stokes codes using the method of manufactured solutions. International Journal for Numerical Methods in Fluids 44, 599–620. doi: 10.1002/fld.660. Salari, K. & Knupp, P. Code verification by the method of manufactured solutions. Technical report, Sandia National Labs, United States, 2000. doi: 10.2172/759450. Samtaney, R., 2009. A method to simulate linear stability of impulsively accelerated density interfaces in ideal-MHD and gas dynamics. Journal of Computational Physics 228, 6773–6783. doi: 10.1016/j.jcp.2009.05.042.

160 Sayadi, T., Schmid, P. J., Richecoeur, F. & Durox, D., 2015. Parametrized data-driven decomposition for bifurcation analysis, with application to thermo-acoustically unstable systems. Physics of Fluids 27, 037102. doi: 10.1063/1.4913868. Schmid, P. & Sesterhenn, J. Dynamic Mode Decomposition of numerical and experimental data. In APS Division of Fluid Dynamics Meeting Abstracts, page MR.007. 2008. Schmid, P. J., 2010. Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics 656, 5–28. doi: 10.1017/S0022112010001217. Schmid, P. J., 2011. Application of the dynamic mode decomposition to experimental data. Experiments in Fluids 50, 1123–1130. doi: 10.1007/s00348-010-0911-3. Schmid, P. J., Li, L., Juniper, M. P. & Pust, O., 2011. Applications of the dynamic mode decomposition. Theoretical and Computational Fluid Dynamics 25, 249– 259. doi: 10.1007/s00162-010-0203-9. Semenko, R., Faria, L. M., Kasimov, A. R. & Ermolaev, B. S., 2016. Set-valued solutions for non-ideal detonation. Shock Waves 26, no. 2, 141–160. doi: 10.1007/s00193-015-0610-3. Sharpe, G. J., 1997. Linear stability of idealized detonations. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 453, 2603–2625. doi: 10.1098/rspa.1997.0139. Sharpe, G. J., 1999. Linear stability of pathological detonations. Journal of Fluid Mechanics 401, 311–338. doi: 10.1017/S0022112099006655. Sharpe, G. J., 2001. Numerical simulations of pulsating detonations: II. Piston initiated detonations. Combustion Theory and Modelling 5, 623–638. doi: 10.1088/1364-7830/5/4/307. Sharpe, G. J. & Falle, S. A. E. G., 1999. One-dimensional numerical simulations of idealized detonations. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 455, 1203–1214. doi: 10.1098/rspa.1999.0355. Sharpe, G. J. & Falle, S. A. E. G., 2000. Numerical simulations of pulsating detonations: I. Nonlinear stability of steady detonations. Combustion Theory and Modelling 4, 557–574. doi: 10.1088/1364-7830/4/4/310.

161 Shchelkin, K. I., 1959. Two cases of unstable combustion. Zhurnal Experimentalnoi i Teoreticheskoi Fiziki 36, 416–420. http://www.jetp.ac.ru/cgi-bin/r/ index/e/9/2/p416?a=list. Shepherd, J. E., 2009. Detonation in gases. Proceedings of the Combustion Institute 32, 83–98. doi: 10.1016/j.proci.2008.08.006. Short, M., 1997. Multidimensional linear stability of a detonation wave at high activation energy. SIAM Journal on Applied Mathematics 57, 307–326. doi: 10.1137/S0036139995288101. Short, M., Bdzil, J. B. & Anguelova, I. I., 2006. Stability of chapman–jouguet detonations for a stiffened-gas model of condensed-phase explosives. Journal of Fluid Mechanics 552, 299–309. doi: 10.1017/s0022112005008347. Short, M. & Dold, J. W., 1996. Linear stability of a detonation wave with a model three-step chain-branching reaction. Mathematical and Computer Modelling 24, 115–123. doi: 10.1016/0895-7177(96)00144-6. Short, M., Kapila, A. K. & Quirk, J. J., 1999. The chemical-gas dynamic mechanisms of pulsating detonation wave instability. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 357, 3621–3637. doi: 10.1098/rsta.1999.0513. Short, M. & Sharpe, G. J., 2003. Pulsating instability of detonations with a two-step chain-branching reaction model: Theory and numerics. Combustion Theory and Modelling 7, 401–416. doi: 10.1088/1364-7830/7/2/311. Short, M. & Stewart, D. S., 1998. Cellular detonation stability. Part 1. A normal-mode linear analysis. Journal of Fluid Mechanics 368, 229–262. doi: 10.1017/S0022112098001682. Sow, A., Semenko, R. E. & Kasimov, A. R., 2017. On a stabilization mechanism for low-velocity detonations. Journal of Fluid Mechanics 816, 539–553. doi: 10.1017/jfm.2017.70. Stankiewicz, W., Morzynski, ´ M., Kotecki, K., Roszak, R. & Nowak, M., 2016. Modal decomposition-based global stability analysis for reduced order modeling of 2D and 3D wake flows. International Journal for Numerical Methods in Fluids 81, 178–191. doi: 10.1002/fld.4181.

162 Stewart, D. S. & Kasimov, A. R., 2005. Theory of detonation with an embedded sonic locus. SIAM Journal on Applied Mathematics 66, 384–407. doi: 10.1137/040616930. Stewart, D. S. & Kasimov, A. R., 2006. State of detonation stability theory and its application to propulsion. Journal of Propulsion and Power 22, 1230–1244. doi: 10.2514/1.21586. Strogatz, S. H., 2015. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. Westview Press, Boulder, CO, 2 edition. Taylor, B. D. Instability of steady and quasi-steady detonations. Dissertation, University of Illinois at Urbana-Champaign, 2010. http://hdl.handle.net/ 2142/18646. Taylor, B. D., Kasimov, A. R. & Stewart, D. S., 2009. Mode selection in weakly unstable two-dimensional detonations. Combustion Theory and Modelling 13, 973–992. doi: 10.1080/13647830903324186. Tissot, G., Cordier, L., Benard, N. & Noack, B. R., 2014. Model reduction using dynamic mode decomposition. Comptes Rendus Mécanique 342, 410–416. doi: 10.1016/j.crme.2013.12.011. Tu, J. H., Rowley, C. W., Luchtenburg, D. M., Brunton, S. L. & Kutz, J. N., 2014. On dynamic mode decomposition: Theory and applications. Journal of Computational Dynamics 1, 391–421. doi: 10.3934/jcd.2014.1.391. Tumin, A., 2007a. Initial-value problem for small disturbances in an idealized onedimensional detonation. Physics of Fluids 19, 106105. doi: 10.1063/1.2793156. Tumin, A., 2007b. Multidomain spectral collocation method for stability analysis of detonations. AIAA Journal 45, 2356–2359. doi: 10.2514/1.29722. Tumin, A. Stability of idealized one-reaction detonations revisited. In 45th AIAA Aerospace Sciences Meeting and Exhibit. 2007c. doi: doi:10.2514/6.2007-987. von Neumann, J. Theory of detonation waves. Technical report 549, National Defense Research Committee of the Office of Scientific Research and Development, 1942.

163 Voytsekhovskiy, B. V., Mitrofanov, V. V. & Topchiyan, M. E., 1963. Struktura Fronta Detonatsii v Gazakh. Sibirskoye Otdeleniye Akademii Nauk SSSR. [English translation: The Structure of a Detonation Front in Gasses. Wright-Patterson Air Force Base, Report FDT-MT-64-527 (AD 633-812)]. Wescott, B. L., Stewart, D. S. & Bdzil, J. B., 2004. On self-similarity of detonation diffraction. Physics of Fluids 16, 373–384. doi: 10.1063/1.1633552. White, D. R., 1961. Turbulent structure of gaseous detonation. Physics of Fluids 4, 465–480. doi: 10.1063/1.1706350. Williams, M. O., Kevrekidis, I. G. & Rowley, C. W., 2015. A data-driven approximation of the Koopman operator: Extending dynamic mode decomposition. Journal of Nonlinear Science 25, 1307–1346. doi: 10.1007/s00332-015-9258-5. Wynn, A., Pearson, D. S., Ganapathisubramani, B. & Goulart, P. J., 2013. Optimal mode decomposition for unsteady flows. Journal of Fluid Mechanics 733, 473– 503. doi: 10.1017/jfm.2013.426. Xu, S., Aslam, T. & Stewart, D. S., 1997. High resolution numerical simulation of ideal and non-ideal compressible reacting flows with embedded internal boundaries. Combustion Theory and Modelling 1, 113–142. doi: 10.1080/713665233. Zel’dovich, Y. B., 1940. K teorii rasprostranenia detonatsii v gazoobraznykh sistemakh. Zhurnal Experimentalnoi i Teoreticheskoi Fiziki 10, 542–568. [English translation: 1950. “The Theory of the Propagation of Detonation in Gaseous Systems”. National Advisory Committee for Aeronautics Technical Memorandum 1261]. Zel’dovich, Y. B. & Kompaneets, A. S., 1955. Theory of detonation. Gostekhizdat, Moscow. [English translation: 1960. “Theory of Detonations”, Academic Press, New York]. Zhang, H., Dawson, S. T. M., Rowley, C. W., Deem, E. A. & Cattafesta, L. N., 2017. Evaluating the accuracy of the dynamic mode decomposition. arXiv preprint . Zharovsky, E., Sandu, A. & Zhang, H., 2015. A class of implicit-explicit two-step Runge–Kutta methods. SIAM Journal on Numerical Analysis 53, 321–341. doi: 10.1137/130937883.