VALIDATION OF THE
KURGANOV-TADMOR/KURGANOV-NOELLE-PETROVASCHEME FOR ROTATING DETONATION ENGINESIMULATIONS USING OPENFOAM

Chibuikem Uchenna Ajaero

A thesisin
The DepartmentofMechanical, Industrial and Aeropace Engineering

Presented in Partial Fulfillment of the Requirements
For the Degree of Master of Applied Science
Concordia University
Montréal, Québec, Canada

January 2019
c Chibuikem Uchenna Ajaero, 2019Concordia University
School of Graduate Studies

This is to certify that the thesis prepared
By:
Entitled:

Chibuikem Uchenna Ajaero
Validation of the Kurganov-Tadmor/Kurganov-NoellePetrova Scheme for Rotating Detonation Engine Simulations using OpenFOAM

and submitted in partial fulﬁllment of the requirements for the degree ofMaster of Applied Science
complies with the regulations of this University and meets the accepted standards with respect to originality and quality.
Signed by the ﬁnal examining commitee:
Dr. Mehdi Hojjati
Dr. Hoi Dick Ng
Dr. Xiaolei Wang
Dr. Charles Basenga Kiyanda

Chair
Examiner
Examiner
Supervisor

Approved
Chair of Department or Graduate Program Director

20
Dr. Amir Asif, Ph.D, PEng Dean
Gina Cody School of Engineering and Computer ScienceAbstract
Validation of the Kurganov-Tadmor/Kurganov-Noelle-Petrova Scheme for
Rotating Detonation Engine Simulations using OpenFOAM
Chibuikem Uchenna Ajaero

Detonation waves are a challenging ﬁeld of study given the short time and length scales
involved in the phenomenon. Such waves exhibit a complex structure consisting of a leadshock and shock waves travelling transversely to a detonation’s normal propagation direction. The interaction between the shocks and the rapid chemical reactions they trigger
results in the emergence of a natural length scale, the detonation cell size. There are still
no complete theory or model that can accurately predict the cell size or a detonation waves
initiation, propagation and failure dynamics. The numerical simulation of detonation waves
is also challenging, due to the rapid reaction rates encountered. The Open source Field
Operation And Manipulation (OpenFOAM) framework, commonly referred to as OpenFOAM, Computational Fluid Mechanics (CFD) package is increasingly used and referenced.
One drawback of the stock OpenFOAM package is that the only ﬁnite volume numerical
scheme available for the solution of the Euler equations in conservative form is the KuganovTadmor (KT)/Kuganov-Noelle-Petrova (KNP) numerical scheme. Moreover, combustion is
not implemented, hence which needs to be modiﬁed to simulate detonation waves; a coupling of compressible ﬂows and reaction. This particular scheme is 2nd order accurate insmooth region based on the idea behind the Lax-Friedrichs scheme and which does not
involve the solution to a Riemann problem in order to evaluate the intercell ﬂuxes. Thisis unlike the methods currently used in detonation research, which nearly always consistof Godunov-type schemes with an approximate Riemann Solver such as Harten-Lax-van
Leer-Contact (HLLC). OpenFOAM, with the KNP scheme, was recently used to simulate
the two-dimensional structure of detonation waves despite having not been fully validated
for the detonation simulation. Eﬀorts to get access to the codes used proved abortive. In
this work, we create a custom solver named rhoCentralFoamreac which we used to evaluate(validate) the KNP scheme for detonation cases, by simulating a standard 1D detonationcase that usually results in pulsating propagation with a single mode. Metrics for detailed
comparison and convergence studies are the oscillation peak pressure and period. Using the
KNP scheme, we then examine OpenFoam as a CFD tool for the simulation of a detonationiii

based engine, where an initiated wave propagates circumferentially in a combustion chamber, commonly referred to as a rotating detonation engine. We study the eﬀect of diﬀerent
ignition methods, and initiation ﬂow ﬁelds (subsonic and supersonic) on the formation ofthese rotating detonation waves.

iv

Acknowledgments
I cannot thank my parents enough for the role you have constantly played in my life, foreven making this adventure possible in the ﬁrst place. I owe you a world of gratitude, butthat in itself is not enough. Thank you for support and for always being there. I wish to
thank my supervisor, Professor Charles Kiyanda, for his immense assistance and guidance
during the course of this work. He kept the door open every time I needed to walk throughit for clariﬁcations. Considering that I did use the door a lot, I just want to say thank you.
Finally, thank you Lord for being my strength.

v

To a most dear friend, Emmanuel Ugbaja. You will forever be in my heart.

viContents
List of Figuresx

List of Tables

xiv

1 INTRODUCTION
1.1 DETONATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 PULSE DETONATION ENGINE . . . . . . . . . . . . . . . . . . . .
1.1.2 STANDING DETONATION WAVE ENGINE . . . . . . . . . . . . .
1.1.3 ROTATING DETONATION ENGINE . . . . . . . . . . . . . . . . .
1.2 DETONATION AND DEFLAGRATION . . . . . . . . . . . . . . . . . . . .
1.2.1 WHAT IS A DEFLAGRATION WAVE . . . . . . . . . . . . . . . . .
1.2.2 WHAT IS A DETONATION WAVE . . . . . . . . . . . . . . . . . .
1.2.3 STRUCTURE OF A DETONATION WAVE . . . . . . . . . . . . .
1.2.4 ESSENTIAL GEOMETRIC RELATIONS . . . . . . . . . . . . . . .
1.2.5 THERMAL EFFICIENCY OF A DETONATION WAVE . . . . . . .
1.2.6 ONE-DIMENSIONAL DETONATION DYNAMICS . . . . . . . . .
1.3 THESIS OVERVIEW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
4
5
5
6
6
7
7
8
9
10
12

2
2.1
2.2
2.3

2.4

THE KURGANOV-TADMOR/KURGANOV-NOELLE-PETROVASCHEME AND OpenFOAM
14
REVIEW OF NUMERICAL DISCRETIZATION . . . . . . . . . . . . . . . 14
GOVERNING EQUATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.1 CONSERVATION LAWS IN INTEGRAL FORM . . . . . . . . . . . 16
TYPES OF NUMERICAL SCHEMES . . . . . . . . . . . . . . . . . . . . . 17
2.3.1 GODUNOV’S METHOD (FIRST-ORDER UPWIND METHOD) . . 17
2.3.2 MUSCL-HANCOCK APPROACH FOR HIGHER RESOLUTION . . 17
CENTERED SCHEMES (LAX-FRIEDRICHS FAMILY OF SCHEMES . . . 18

vii

2.5

2.6

2.4.1 LAX-FRIEDRICHS SCHEME . . . . . . . . . . . . . . . . . . . . . .
2.4.2 NESSYAHU-TADMOR (NT) SCHEME . . . . . . . . . . . . . . . .
2.4.3 DESCRIPTION OF THE KNP NUMERICAL SCHEME . . . . . . .
OVERVIEW OF OpenFOAM . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.1 WHAT IS OpenFOAM . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5.2 FILE STRUCTURE OF OpenFOAM . . . . . . . . . . . . . . . . . .
2.5.3 POST PROCESSING . . . . . . . . . . . . . . . . . . . . . . . . . .
CREATION OF CUSTOM SOLVER FOR DETONATION SIMULATIONS(rhoCentralFoamreac) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18
19
20
23
23
24
27
27

3 VALIDATION OF THE KT/KNP SCHEME FOR DETONATIONSIMULATION
30
3.1 NON-REACTIVE VALIDATION . . . . . . . . . . . . . . . . . . . . . . . . 31
3.1.1 CASE DESCRIPTION AND SET-UP . . . . . . . . . . . . . . . . . 31
3.1.2 BOUNDARY AND INITIAL CONDITIONS . . . . . . . . . . . . . . 31
3.1.3 RESULTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2 REACTIVE VALIDATION . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2.1 CASE SETUP DESCRIPTION . . . . . . . . . . . . . . . . . . . . . 37
3.2.2 INITIAL CONDITIONS . . . . . . . . . . . . . . . . . . . . . . . . 38
3.2.3 SOLUTION METHOD . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.3 RESULTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.4 ANALYSIS AND CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . 40
4 SIMULATION OF COLD ENGINE START
4.1 PROBLEM DESCRIPTION/CASE SET UP . . . . . . . . . . . . . . . . . .
4.1.1 MESH GENERATION . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 GEOMETRY SETUP . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.3 ACTUAL ENGINE GEOMETRY . . . . . . . . . . . . . . . . . . .
4.1.4 AREA RATIO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 BOUNDARY CONDITIONS . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 FIXED VALUE BOUNDARY CONDITION . . . . . . . . . . . . .
4.2.2 INLET-OUTLET BOUNDARY CONDITION . . . . . . . . . . . .
4.2.3 TOTAL PRESSURE/TEMPERATURE BOUNDARY CONDITION
4.2.4 PLENUM PRESSURE BOUNDARY CONDITION . . . . . . . . . .
4.3 CHOICE OF FLOW CONDITIONS . . . . . . . . . . . . . . . . . . . . . .

viii

46
46
47
51
52
52
53
53
54
55
56
56

4.3.1

SUPERSONIC FLOW FIELDS: BACK PRESSURE OF 124kPa and
100kPa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 SUBSONIC FLOW FIELDS: CASES WITH BACK PRESSURE OF
619 kPa and 650 kPa . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.3 MACH NUMBER . . . . . . . . . . . . . . . . . . . . . . . . . . . .
INVESTIGATION OF MASS FLOW RATE . . . . . . . . . . . . . . . . . .
4.4.1 RESULTS AND ANALYSIS . . . . . . . . . . . . . . . . . . . . . . .
SELECTED FLOW FIELDS FOR INITIATION . . . . . . . . . . . . . . . .
4.5.1 FLOW FIELD FOR SUBSONIC INITIALIZATION . . . . . . . . .
4.5.2 FLOW FIELD FOR SUPERSONIC INITIALIZATION . . . . . . . .

60
60
60
62
66
66
66

5 SIMULATION OF ROTATING DETONATION ENGINE INITIATION
5.1 PROBLEM DESCRIPTION/CASE SET UP . . . . . . . . . . . . . . . . . .
5.1.1 MESH DESCRIPTION . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.2 BOUNDARY CONDITION . . . . . . . . . . . . . . . . . . . . . . .
5.1.3 INITIAL CONDITIONS . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.4 DIFFERENT TYPES OF INITIATION . . . . . . . . . . . . . . . .
5.2 RESULTS WITH LOW RESOLUTION MESH . . . . . . . . . . . . . . . .
5.2.1 SUBSONIC FLOW FIELD INITIATION . . . . . . . . . . . . . . . .
5.2.2 INITIATION WITH SUPERSONIC FLOW FIELD . . . . . . . . .
5.3 EFFECT OF HIGHER RESOLUTION ON DETONATION INITIATION .
5.3.1 SUBSONIC CASE . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.2 SUPERSONIC CASE . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.3 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67
67
67
68
68
68
68
70
72
77
77
78
79

6 CONCLUSION
6.1 RECOMMENDATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

85
86

4.4
4.5

Appendicesix

60

List of Figures
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

Images are of detonation waves in a mixture of CH4 + 2O2 from [KH13]
showing coupling eﬀect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Rotating Detonation Engine from [LB14] . . . . . . . . . . . . . . . . . . . .
Detonation Cell Size with complex ﬁsh scale structure [VJ14] . . . . . . . . .
Control volume for 1-D dynamics [Tat15] . . . . . . . . . . . . . . . . . . . .
Rayleigh-Hugoniot curve showing diﬀerent regions . . . . . . . . . . . . . . .
Representation of the scheme in ﬁnite volume formulation . . . . . . . . . .
Overview of the structure of OpenFOAM [Ope16] . . . . . . . . . . . . . . .
OpenFOAM case formulation tree [Ope16] . . . . . . . . . . . . . . . . . . .
Pressure ﬁle at time 0 (P) . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Inclusion of Transport Equation to rhoCentralFoam. Fig10h shows the advection of ”Y” after some time . . . . . . . . . . . . . . . . . . . . . . . . . .
Images are from simulation of detonation waves in [MTE17b] . . . . . . . . .
Case set-up for non-reactive validation (shock tube problem) . . . . . . . . .
Initial State (t = 0.000s) . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Simulated State (t = 0.007s) . . . . . . . . . . . . . . . . . . . . . . . . . . .
Initial condition for P , U and T (t = 0.000 s). . . . . . . . . . . . . . . . . .
End state showing P,U,T and ρ (t = 0.007 s) . . . . . . . . . . . . . . . . . .
Pressure vs position for the shock tube problem using the minmod limiter
CFL = 0.09 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Flow simulations at diﬀerent Courant-Friedrichs-Levy number (CFL) numbers
Flow simulations using superbee limiter, indicating intense amplitude of oscillation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Disturbance amplitude DP (Pa) vs CFL number for the minmod, van Alabada,van Leer and superBee limiters for a pressure ratio of 10. . . . . . . . . . . .
Disturbance amplitude DP (Pa) vs CFL number for the minmod, van Alabada,van Leer and superBee limiters for a pressure ratio of 1.5. . . . . . . . . . . .

x

4
6
8
10
12
22
23
24
25
28
31
32
32
32
33
33
34
35
35
36
36

22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

Disturbance amplitude DP vs CFL number for the minmod, van Alabada,van Leer and superBee limiters for a pressure ratio of 2. . . . . . . . . . . . . 37
Disturbance amplitude vs CFL number for the minmod, van Alabada, vanLeer and superBee limiters for a pressure ratio of 5. . . . . . . . . . . . . . . 38
Case set-up for reactive validation (piston problem) . . . . . . . . . . . . . . 39
Result from octave codes for Zeldovich-Von Neumann-Döring (ZND) model
41
Shock peak pressure vs x axis for Δx = half reaction zone length (L1/2 ) /10
and Δx = L1/2 /20 for a domain 200 L1/2 long . . . . . . . . . . . . . . . . . 42
Pshock /P0 vs position for Δx = L1/2 /10 for ZND initiated case (600 L1/2 long). 42
Pshock /P0 vs time for Δx = L1/2 /10 for ZND initiated case (600 L1/2 long). . 42
Pshock /P0 vs position for Δx = L1/2 /20. (ZND)). . . . . . . . . . . . . . . . . 43
Pshock /P0 vs time for Δx = L1/2 /20. (ZND). . . . . . . . . . . . . . . . . . . 43
Pshock /P0 vs position for Δx = L1/2 /40. (ZND) . . . . . . . . . . . . . . . . . 43
Pshock /P0 vs time for Δx = L1/2 /40. (ZND). . . . . . . . . . . . . . . . . . . 43
Pshock /P0 vs position for Δx = L1/2 /60. (ZND) . . . . . . . . . . . . . . . . . 43
Pshock /P0 vs time for Δx = L1/2 /60. (ZND). . . . . . . . . . . . . . . . . . . 43
Pshock /P0 vs position for Δx = L1/2 /80. (ZND) . . . . . . . . . . . . . . . . . 44
Pshock /P0 vs time for Δx = L1/2 /80. (ZND). . . . . . . . . . . . . . . . . . . 44
Pshock /P0 vs position for Δx = L1/2 /100. (ZND) . . . . . . . . . . . . . . . . 44
Pshock /P0 vs time for Δx = L1/2 /100. (ZND). . . . . . . . . . . . . . . . . . 44
Pshock /P0 vs position for Δx = L1/2 /80 for an extended domain. (ZND) . . . 44
Pshock /P0 vs time for Δx = L1/2 /80 for an extended domain. (ZND). . . . . 44
Image shows single mode pulsating detonation [NNL00]. . . . . . . . . . . . 45
Pshock /P0 vs Number of points per L 1 [NNL00]. . . . . . . . . . . . . . . . . 45
2
Rotating Detonation Engine [LB14] . . . . . . . . . . . . . . . . . . . . . . . 46
Representaion of an RDE [Sha15] . . . . . . . . . . . . . . . . . . . . . . . . 46
Image shows unwrapping of an RDE (3D to 2D). . . . . . . . . . . . . . . . 47
2D RDE representaion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2D RDE representaion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2D Mesh with one obstacle . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
Sketch of domain with 2 obstacles . . . . . . . . . . . . . . . . . . . . . . . . 48
2D mesh generation with two obstacles . . . . . . . . . . . . . . . . . . . . . 49
Initial pressure for the case . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2D mesh with three obstacles . . . . . . . . . . . . . . . . . . . . . . . . . . 49
Image showing mesh blocks(11 in this case) . . . . . . . . . . . . . . . . . . 49
Vertices for 2 obstacles using octave code. . . . . . . . . . . . . . . . . . . . 50
xi

55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

Image of mesh with 10 obstacles (code generated) . . . . . . . . . . . . . . .
Image of mesh with 100 Obstacles (code generated) . . . . . . . . . . . . . .
A resulting ﬂow ﬁeld using engine A. . . . . . . . . . . . . . . . . . . . . . .
Pressure ﬂow ﬁelds in the ”slip” BC test with 1 obstacle . . . . . . . . . . .
Pressure ﬂow ﬁelds in the ”slip” BC test with 2 obstacle . . . . . . . . . . .
InletOutlet boundary condition [Cfd] . . . . . . . . . . . . . . . . . . . . . .
Plots showing ﬁxed boundary condition with ﬁeld variables indicating unphysical behaviour in the boundary . . . . . . . . . . . . . . . . . . . . . . . . . .
Plots obtained using total pressure and temperature boundary condition . .
Plots showing plenum boundary condition with ﬁeld variables . . . . . . . .
Steady state subsonic ﬂow ﬁeld using a back pressure of 650 kPa . . . . . . .
Mach number for subsonic and supersonic initiation ﬂow ﬁelds . . . . . . . .
Mass ﬂow rate across boundaries for Pb = 100 kPa. . . . . . . . . . . . . . .

50
50
51
54
54
55
57
58
59
61
62
63
63
64
65
69
70
71

Mass ﬂow rates across boundaries for Pb = 250 kPa. . . . . . . . . . . . . . .
Mass ﬂow rate across at boundaries for Pb = 500 kPa . . . . . . . . . . . . .
Mass ﬂow rate against time across boundaries for Pb = 650 kPa . . . . . . .
Subsonic ﬂow ﬁeld using a back pressure of 619 kPa . . . . . . . . . . . . . .
Various initiation techniques considered . . . . . . . . . . . . . . . . . . . .
Detonation initiated with a high temperature using a subsonic ﬂow ﬁeld . . .
Detonation initiated with a high pressure-temperature zone using a subsonic
ﬂow ﬁeld . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
Detonation initiated further downstream with a high temperature zone using
a subsonic ﬂow ﬁeld . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
Detonation initiated with a high pressure and temperature zone using a subsonic Flow ﬁeld . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
Detonation initiated with a high temperature zone using a supersonic ﬂow ﬁeld 76
Detonation initiated with a high temperature and pressure zone using a supersonic ﬂow ﬁeld . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
Temperature based initiation for supersonic ﬂow with a pressure ratio of 1.2857 81
Pressure-Temperature based initiation for subsonic ﬂow with a pressure ratioof 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
Pressure-Temperature based initiation for supersonic ﬂow with a pressure ratio of 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
Initial unstable detonation wave (Subsonic case) . . . . . . . . . . . . . . . 82
Image indicates decoupling of shock pressure-reaction(Unsustained Detonation for subsonic case) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
xii

83
84

Pressure-reaction coupling indicating a detonation wave (supersonic case) . .
Coupling remains as wave travels further (supersonic case) . . . . . . . . . .

xiii

83
84

List of Tables
1
2
3
4

Table of CFL and amplitude of pressure oscillations for the 4 chosen limiters
for a pressure ratio of 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Table of CFL and amplitude of pressure oscillations for the 4 chosen limiters
for pressure ratio of 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Table of CFL and amplitude of pressure oscillations for the 4 chosen limiters
using pressure ratio of 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Table of CFL and amplitude of pressure oscillations for the 4 chosen limiters
using pressure ratio of 1.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xiv

Abbreviations
CFD Computational Fluid Mechanics
CFL Courant-Friedrichs-Levy number
CJ Chapman-Jouguet
DCJ detonation velocity at Chapman-Jouguet (CJ) point
DDT Deﬂagration to Detonation Transition
FVM Finite Volume Methods
HLL Harten-Lax-van Leer
HLLC Harten-Lax-van Leer-Contact
KNP Kuganov-Noelle-Petrova
KT Kuganov-Tadmor
ODW Oblique Detonation Wave
OpenFOAM Open source Field Operation And Manipulation
OWDE Oblique Wave Detonation Engine
PDE Pulse Detonation Engine
RDE Rotating Detonation Engine
SDE Standing Detonation Engine
WENO Weighted Essentially Non-oscillatory Schemes
ZND Zeldovich-Von Neumann-Döring

Nomenclaturec unwrapped circumference. m
ΔR total length of the reaction zone. m
L length of engine. m
1

Ea activation energy. kJ/kg
l injection width. mṁ mass ﬂow rate. kg/s
n number of obstacles
γ adiabatic index
Hf heat of formation. kJ/kg
h height of ﬂuid. me internal energy per unit mass. kJ/kg
k pre-exponential constant 1/s
L1/2 half reaction zone length. m
M moleculer mass of mixture. kg/kmol
Q heat release. kJ/kg
u velocity. m/s
ts obstacle thickness. m

2

Chapter 1
INTRODUCTION
1.1

SHOCK INTERACTIONS AND COMBUSTIONAS A PREREQUISITE FOR DETONATION

It is common practice to generate power from low-speed (subsonic), constant pressure combustion processes. Of growing interest however has been the idea of realizing this highly
craved demand for power generation through a rather more spontaneous heat release process, as is the case with detonation waves [EGL91]. The interaction of shock waves in
a compressible medium, and exothermic chemical reactions, both occurring on short timescales, gives rise to detonation waves, wherein we have a coupling of hydrodynamic andchemical processes. The self-propagating waves similar to that shown in ﬁgure 1 occur on
short time and length scales, with propagation speeds in the order of 2–3 km/s in reactive
gases [FD79]. These waves generate pressures 6–8 times the initial reactant pressure, andhigh temperatures (above 2000 K). High and rapid energy conversion is a signiﬁcant component of a detonation process which is quite promising particularly as it relates to the ﬁeldof propulsion. The period between 1940s/50s saw an increase in activities in the detonationﬁeld, aimed at exploiting its promise of low weight to high power ratio. These eﬀorts mainly
focused on understanding the phenomena, with huge interest shown to the development ofwaves from slow moving combustion waves to fast travelling (detonation), as well as assessing the actual practicability of the process. Three interesting areas of application have been
examined in the course of utilizing detonation waves for propulsion. These include the Pulse
Detonation Engine (PDE), the Oblique Wave Detonation Engine (OWDE) and the Rotating
Detonation Engine (RDE). We examine these 3 brieﬂy.

3

Figure 1: Images are of detonation waves in a mixture of CH4 + 2O2 from [KH13] showing couplingeﬀect.

1.1.1

PULSE DETONATION ENGINEPDE utilize a detonation wave for propulsion through the use of an intermittent detonation
cycle process, where a detonation wave is initiated at one end of a mixture ﬁlled chamber[EGL91]. They are based on axially propagating detonation waves mainly in the downstream
direction. A lot of studies have been done on this engine. Thomas [Bus96] received a patentfor his work on the PDE coupled to a fuel and air source in 1996. Nichols et al. looked at
the practicability of a reaction device functioning on intermittent gaseous detonation waves[Nic57] by examining the thrust, fuel-air ﬂow among others. PDE research continues tothis day. We note here that a challenge facing the PDE is the fact that its operation isintermittent, needing re-initiations of high frequency, which could be energy demanding.
Add that to the low thrust it generates, which in itself results from its mode (pulsed) ofoperation, necessitated by the need to ”reﬁll” the engine, during which time there is zerothrust generation [ZWW]. The intermittent nature of the device also makes it impossible todesign a practical nozzle. Hence the focus on PDEs these days gears towards deﬂagaration todetonation transition studies; using obstacles to achieve ﬂame acceleration thereby reducing
energy demand [Tat15]. There has been a test ﬂight on an engine modiﬁed to run on pulsed
detonation [MMM98].

4

1.1.2

STANDING DETONATION WAVE ENGINEIn this detonation-based type of engine, a detonation wave is stabilized within a chamber,
oblique or normal to a wedge [ZWW], with fuel injected and mixed with supersonic air ﬂow,
ahead of the wave. The mixture is then compressed by the resulting shock wave after which
detonation occurs, and the combustion products used to generate thrust. The Standing
Detonation Engine (SDE) requires high inﬂow velocities and can only work with the inﬂow
of high mach numbers [Kai09] hence even though the idea seems quite promising from a
conceptual stand point, the SDE is faced with certain diﬃculties that challenge its practicalimplementation. Studies show that it has been quite challenging to sustain detonation
waves in the chamber for long [ZWW]. Normal waves are more diﬃcult to stabilize andthey usually result in higher stagnation losses in pressure than with oblique waves, hence
the desire to rather utilize oblique shock waves for detonation. Braun [BLWC13] and others
describe oblique detonation as a combustion process initiated by an oblique shock wave.
Powers et al. [PJBJ94] also explains an oblique detonation wave as one wherein the pressureﬁeld induced by combustion behind the shock wave inﬂuences the wave itself, resulting ina rotation of the oblique wave to a normal orientation, the limiting case being an oblique
Chapman-Jouguet detonation, in which there is a mach 1 normal velocity component behind
the wave. An increase in heat addition ultimately means a decoupling from the anchoring
point, thereby making it rotate to a more normal orientation [PJBJ94]. The idea that byutilizing an oblique shock to initiate a detonation, we can achieve the stabilization of adetonation wave in a supersonic ﬂow ﬁeld has escalated interest in its use for hypersonic
ﬂight [PJBJ94]. Carrier et al. [CFM+ 92] looked at utilizing lasers to initiate such waves byanalyzing the practicability of a supersonic combustor which utilizes a stabilized,conicallyconﬁgured Oblique Detonation Wave (ODW). We note that in spite of eﬀorts like this, andmany years of research, there remains the challenge of stabilizing a steady oblique detonationin high speed, combustible mixtures.

1.1.3

ROTATING DETONATION ENGINEThis represents an approach to ensuring pressure gain combustion and has been receiving
increased attention in the last 10 years, mainly because of its ability to provide signiﬁcant
improvements in thermal eﬃciency, and the promise of engine design simpliﬁcation. Far backas in the early 1960s, Voitsekhovskii et al. carried out experiments on continuously rotating
detonations [VM63]. Tobita et al. [TFW10] made an application for both Japanese andUS based patents for their rotating detonation based device. The RDE mode of operation
presents a certain advantage over the PDE concept. The RDE requires a single initial
5

Figure 2: Rotating Detonation Engine from [LB14]

activation/ignition energy to have its four main processes occurring at very high frequencies(in order of 103 Hz) [PFHS15] during operation (introduction of reactants, mixture blowdown in the chamber, purging of combustion products) since the detonation wave propagatescircumferentially in an annular chamber while consuming combustible products. Hence it
is always generating thrust and it doesn’t require Deﬂagration to Detonation Transition(DDT) devices. Figure 2 shows an overview of the RDE design for rocket mode operation.
Here we basically see an annular chamber of a given diameter, through which we have ﬂuid
injection into the combustion chamber. Depending on the conﬁguration, it could be premixed
or separately injected, in which case the oxidizer and the fuel enter through diﬀerent ﬂow
channels. Engine dynamics and geometry give rise to a large number of governing parameters[BZV06]. The goal is that a RDE functions as one, by having at least one stable detonation
wave travelling circumferentially in the annular chamber.

1.2

DETONATION AND DEFLAGRATION

Generally speaking, there are two types of self-propagating combustion waves: deﬂagarationand detonations [Lee08]. We provide some insight on these modes of combustion.

1.2.1

WHAT IS A DEFLAGRATION WAVE

Deﬂagration is a low speed (subsonic) mode of combustion involving a mixture of a fuel andan oxygen source, that results in heat release and rise in temperature. It is common to thinkof it as a constant pressure process since it usually occurs with little changes in pressure.
Common applications include the internal combustion engine and in gas turbines. They areused to create useful work and generate power when initiated, using some spark mechanism,
and allowing the ﬂame front to propagate through the fuel-Air mixture usually at speeds

6close to 1 m/s [Tat15].

1.2.2

WHAT IS A DETONATION WAVE

Detonation on the other hand can be described as a supersonic burning process involving a
combustible mixture, wherein a shock wave triggers a chemical reaction, and is sustained bythis exothermic reaction. This enables it to propagate at speeds in the order of 2000 m/s[Tat15]. The wave thickness and whether it can be considered as a shock wave accompanied
by a chemical reaction depends on its structure [FD79]. We note therefore that the supersonicwave with a shock structure raises the temperature and pressure of the mixture, [Tat15]
hence the reason devices utilizing a detonation wave are often referred to as self-pressurizing.
Generally speaking, detonation velocities obtained in gaseous mixtures are usually about 4–
7 times higher than the corresponding sound speed in their reactant media [GSMK12] andthe velocity depends on the mixture composition, its temperature and pressure. The use of
detonation as a combustion process is motivated by the expectation of improved performance,with great expectations for jet propulsion, considering that chemical energy release occurs
in a very short time over a small region in space, in comparison to the deﬂagration mode.

1.2.3

STRUCTURE OF A DETONATION WAVE

Generally speaking, detonations are unstable, possessing a cellular frontal structure [Lee08]
with cell boundaries deﬁned by the connecting points of transversely propagating shock waves
fed by energy release in the front that is non uniform. This implies, in contradiction to the CJ
theory, that there exists an interactive relationship between the wave front and the ongoingenergy-release happening behind it [Lee08]. Figure 3 shows a typical representation of thisdetonation wave structure with the intersection between the normal, lead shocks and thetransverse shock waves denoted, as they are commonly called, as triple points. The averagewidth of successive triple points is the cell width, λ. This length scale λ, in geometries
where the overall scale is large compared to the detonation cell width, is a characteristic ofthe detonable mixture (its constituents and stoichiometry) and of its initial thermodynamic
state (mainly its initial pressure, but also its initial temperature). In the one-dimensional
analysis of detonation waves (the ZND model), the only emerging length scale is the steadydetonation thickness, i.e. the distance between the lead front and the end of reaction. Since,for some reaction models, the detonation thickness can be inﬁnite, it is customary to use thehalf-reaction zone length, L1/2 as a reference detonation length scale. L1/2 corresponds to thephysical extent between the lead shock wave and the point, within the detonation structure,
where half of the chemical potential energy has been released as heat. The detonation cell7

Figure 3: Detonation Cell Size with complex ﬁsh scale structure [VJ14]

width can be broadly correlated to the chemical length scale, either the overall detonation
thickness or L1/2 [ST65].

1.2.4

ESSENTIAL GEOMETRIC RELATIONSAs stated earlier, there needs to be a steady renewal of the chemical reactants , and hence, thelayer of combustible mixture, to sustain a rotating/continuous spin detonation wave[BZV06].
Bykovskii et al. in [BZV06] highlighted that this layer, h, should not be any smaller thancertain critical value h∗ . They showed that a change in the length of engine (L) of an annularcylindrical chamber does not strongly inﬂuence the detonation process [BZV06]; making L
large enough results in development of boundary layer by lowering eﬀective cross sectional
area which in turn aids momentum heat losses [BZV06]. It can be seen from experiments that
L should be between 1.5 to 2 times higher than h for obtaining stable transverse detonation
waves (TDWs). The critical value the of height of ﬂuid (h), is related to the total lengthof the reaction zone (ΔR ), in the detonation wave as shown in [BZV06]. Therefore theexpression below is presented;h ≈ h∗ ≈ (17 ± 7)ΔR .
(1)
In the absence of physical processes (premixed gaseous mixture) [BZV06] or in cases where
their duration can be neglected,
ΔR ≈ 0.7λ(2)8where λ the width of the self-excited cell of the detonation front, the detonation cell size,
readily measured and known for many mixtures. Combining 1 and 2 as in [BZV06] givesh ≈ (12 ± 5)λ.

(3)which agrees well with experimental values.
Another useful geometric parameter [BZV06] is given as
K = l/h = πdc /nh.

(4)which has been shown to be roughly the same for all annular cylindrical chambers using agaseous oxidizer, given as K = 7±2. The parameter n is the number of transverse detonation
waves rotating in the chamber. Thus with K known, for n = 1, from 4, we obtain that for aprescribed mean pressure in a chamber, the minimum chamber diameter [BZV06] must be
such that(dc )min = hK/π.
(5)For chambers operating on a gaseous oxidizer, we have (dc )min = 40ΔR with the value ofthe minimum diameter known, the minimum and optimal length of the chambers can be
determined approximately as:
Lmin ≈ 2h ≈ dc /n.

(6)Lopt ≥ 4h ≈ 0.71 ≈ 2dc /n.

(7)These equations are important since they serve as benchmarks for engine design. The also
highlight the usefulness of knowing the cell size of a mixture.
Thus given a cell size which is a function of a given mixture composition. We can determine
the geometry needed such that an engine designed to run on such mixture can support a
rotating detonation wave. We can also simulate such an engine numerically since we know
geometric parameters and can specify an equivalence ratio. The minimum diameter andchannel width have been shown to be approximately d ≤ 30ΔR and δ ≤ 2.5ΔR , respectively[LB14]. These implies that the RDE can be scaled up or down depending on how reactive
the mixture components are.

1.2.5

THERMAL EFFICIENCY OF A DETONATION WAVE

Cycles based on detonations are more thermodynamically eﬃcient. Detonations, which aremodelled as constant volume combustion (Humphrey and Fickett-Jacobs thermodynamic
9

Figure 4: Control volume for 1-D dynamics [Tat15]

cycles) processes produce a higher thermal eﬃciency (1.3 -1.5 times more) [Tat15] thanthat of a constant pressure combustion cycle at an equivalent pressure ratio, meaning that
this can result in a similar increase in fuel eﬃciency provided other mechanical and related
eﬃciencies can be maintained [BP94]. The formation and propagation of a detonation wave
compresses the gas ahead of it causing a dramatic increase in pressure and temperature after
the combustion process. This can be described by the one dimensional Chapman-Jouguet
theory and the Zeldovich-Von Neumann-Döring (ZND) model.

1.2.6

ONE-DIMENSIONAL DETONATION DYNAMICS

The conservation equations are presented for a combustion wave with a constant area ﬂow.
Figure 4 describes a control volume for analysis. The following four assumptions are pertinent:
1. the detonation is steady,2. the ﬂow is one-dimensional and planar,
3. the products approach a state of equilibrium at a distance behind the detonation front,
4. the detonation velocity obtained is the minimum permitted by the conversion conditions.
The mass, momentum and energy balances across the cross sectional area suggest that a
relationship exists between the states on two sides of a shock [Dah11]. This results in thefollowing expressions based on the assumptions made earlier:10ρa V a = ρb V b(8)

Pa + ρa Va2 = Pb + ρb Vb2
V3
V3
ρ a V a h a + ρa a = ρ b V b h b + ρb b
2
2

(9)
(10)

These are known as the ”energy and mass balances”. They express a relationship between
thermodynamic quantities before and after a shock, with the left and right values representing
low and high pressure regions in this case. A combination of these equations result inγ − 1 Pb
P a Pb − P a 1
1( )−( + ) = −(hb − ha )
−
γ
ρb
ρa
2
ρb ρa

(11)

It is pertinent to state that these relations have been derived without considering the speciﬁc
processes within the control volume. Thus the distinguishing factor between a detonation
and a shock is that whereas for a shock wave the chemical properties remain the same, inthe case of a detonation, there is an energy release due to changes in chemical composition.
Thus, the main diﬀerence in analysis would be the enthalpy behind the detonation, i.e. at
the exit of the control volume, hb [Dah11]. For a given value of Pa and change in enthalpy(−(hb − ha )), a representative plot of Pb against 1/ρb is known as the Hugoniot curve, andpoint (Pa , 1/ρa ) refers to the unreacted state. It represents the connection of all points
representing a solution of the conservation equations as restated in equation 11. Two points
known as the Chapman-Jouguet points exist by drawing tangents through the Hugoniot
origin. The mass and momentum equations can be combined to obtain the expression:
ρ2a va2 =Pb − P a
1
− ρ1b
ρa

(12)

This yields the Rayleigh relation. A combination of equations 11 and 12 results in theRayleigh-Hugoniot relation which shows the regions within which we can have a detonation
or deﬂagaration. Figure 5 shows the various regions of the Raleigh-Hugoniot curve. We notethat regions 1 and 2 are possible areas for strong and weak detonation, while regions 3 and
4 are those for strong and weak deﬂagaration. However, looking at region 5, we can see that1− ρ1b results in a negative number, and Pa − Pb gives a positive result. These two answers
ρaare thermodynamically impossible; therefore a solution to the Rayleigh-Hugoniot line fordetonation and deﬂagaration at these points does not exist [Dah11]. Chapman [Cha99] andJouguet [Jou05] independently stated that detonations propagate at a given velocity, which is
the minimum for the branch that can possibly result in a detonation. This velocity is known

11

Figure 5: Rayleigh-Hugoniot curve showing diﬀerent regions[hys16]as the Chapman-Jouguet (CJ) velocity. Further manipulations of equation 11 show that themach number behind the shock wave is unity (Mb = 1) at the Chapman-Jouguet points.
For emphasis, the CJ theory is arrived at with complete disregard for the detonation structure. It is essentially a consideration of the possible solutions of the steady one-dimensional
conservation equations that link the upstream and downstream equilibrium states of thereactants and products, respectively. In the absence of some description for the structure
however, the propagation mechanism of the detonation wave cannot be known. Althoughshock-induced ignition was known to the early pioneers who discovered the phenomenon, it
was Zeldovich (1940), Von-Neumann (1942), and Döring (1943) [Zel40, VN42, Dör43] who
independently described the model of the detonation structure as comprising a lead shock
closely followed by a ﬁnite chemical reaction zone. Their model is known as the ZND modeland is presented in chapter 3.

1.3

THESIS OVERVIEW

Detonation as a concept is introduced, discussing various eﬀorts at utilizing it for power
branch generation in chapter 1. Chapter 2 discusses some select numerical schemes andtheir formulations and presents results obtained from in-house implementation of the schemes
using Octave. It also presents the implementation of the Kuganov-Tadmor (KT)/KuganovNoelle-Petrova (KNP) scheme in the Open source Field Operation And Manipulation (OpenFOAM). Chapter 3 evaluates the KNP scheme with a view to determine the conditions underwhich the scheme is suitable for the simulation of detonation waves. This was done through12

a reactive and non reactive validation process using a one-dimensional shock tube case.
Chapter 4 is concerned with the simulation of the cold engine starting process for an RDE.
Chapter 5 shows diﬀerent methods of initiation of an RDE. Chapter 6 presents conclusion
and recommendations based on ﬁndings made, and also provides some of the challenges during the course of this work. Certain tables/codes generated during the course of this work
are included in the appendix section.

13

Chapter 2
THE KURGANOVTADMOR/KURGANOV-NOELLEPETROVA SCHEME AND
OpenFOAM
2.1

REVIEW OF NUMERICAL DISCRETIZATION

For partial diﬀerential equations, it is common to employ numerical solution schemes thatconsist of discretizing the spatial extent, discretizing time, and ﬁnally expressing derivative
quantities approximately on those discretized grids. The approximate integrated functionsof the diﬀerential system are then obtained at discrete time values and either (a) at discrete
points in space, in which case they are called ﬁnite diﬀerence or ﬁnite element methods,or (b) over a certain region of space (i.e. as an average value over a given spatial cell),in which case the are referred to as ﬁnite volume methods. For problems involving shock
wave and other discontinuities, ﬁnite volume methods are the almost exclusively appropriate
methods. Finite diﬀerence schemes are, except to some degree for 1st order accurate schemes,
incapable of capturing discontinuities. Finite volume methods discretize the integrated formof a system of PDE, making them able to capture discontinuities in function values, also
called weak solutions.
In this chapter, an introduction to the Euler equations that describe the compressible ﬂow isdone. Next is a review of the Godunov-type schemes, most commonly employed for compressible ﬂow problems, as well as the Lax-Friedrichs family of numerical schemes. A summaryof the KT and KNP numerical schemes for compressible ﬂow, which are an extension of14

the Lax-Friedrichs scheme is presented. Finally, the OpenFOAM numerical framework isreviewed.

2.2

GOVERNING EQUATIONSThe Euler equations are the balance of mass, momentum and energy, and govern the dynamics of inviscid, compressible ﬂow.
ρt + (ρu)x + (ρv)y + (ρw)z = 0,(13)(ρu)t + (ρu2 + p)x + (ρuv)y + (ρuw)z = 0,(14)(ρv)t + (ρuv)x + (ρv + p)y + (ρvw)z = 0,

(15)(ρw)t + (ρuw)x + (ρvw)y + (ρw + p)z = 0,

(16)

Et + [u(E + p)]x + [v(E + p)]y + [w(E + p)]z = 0.

(17)

In equations 13–17, the Euler equations are expressed in terms of the primitive or physicalvariables density, ρ(x, y, z, t), pressure P (x, y, z, t), velocity V (x, y, z, t) = [u, v, w] and totalenergy E(x, y, z, t), whereV2
E =ρ e+.
(18)
2
The primitive variables are more convenient to describe the problem’s initial and boundary
conditions, but, as shown by Toro [Tor] and Hou et al. [HL94], numerical methods that
discretize the primitive-based, diﬀerential form of the equations are unable to accurately
capture discontinuities. It is necessary to integrate the PDE system in its conservative form
where the integrated functions are the conserved variables⎤ρ⎢ ⎥⎢ ρu ⎥⎢ ⎥⎥U =⎢⎢ ρv ⎥
⎢ ⎥
⎣ρw⎦
E

⎡

⎡

F

⎤
ρu
⎢ 2
⎥
⎢ ρu + p ⎥⎢⎥⎥=⎢⎢ ρuv ⎥⎢⎥
⎣ ρuw ⎦u(E + p)

⎡

⎤
ρv
⎢ 2
⎥
⎢ ρu + p ⎥
⎢
⎥
⎥
G=⎢
⎢ ρvu ⎥⎢⎥
⎣ ρvw ⎦
v(E + p)

⎡

H

⎤
ρw
⎢
⎥
⎢ ρw2 + p ⎥⎢⎥⎥=⎢⎢ ρwu ⎥⎢⎥
⎣ ρwv ⎦
w(E + p)(19)

and the system can be expressed as
Ut + F (U )x + G(U )y + H(U )z = 0,(20)

Only one equation is still required to close the set, the equation of state that relates e =e(P, ρ). The ideal gas equation of state is utilized for analysis (i.e ρ = P/RT ).
15

2.2.1CONSERVATION LAWS IN INTEGRAL FORM

A useful form of expressing conservation is its integral form. In one-dimension, the Euler
equation is:
Ut + F (U )x = 0(21)where subscript t indicates time dependence. This continuous mathematical expression canbe conﬁned within a ﬁnite extent by discretization in time and space. In terms of volume
considerations, a domain can split this into m computing cells (volumes) Ii = [xi− 1 , xi+ 1 ]22
with each of width Δx = xi+ 1 − xi− 1 , where subscript ”i” ranges from 1 to m. A control
2
2
volume Vc = Ii × [tn , tn+1 ] is also deﬁned [ZM10]. By evaluating the spatial integral ofequation 21 above given by the expression,
d
dt

xi+ 1

2

U (x, t)dx = F (U (xi− 1 ), t) − F (U (xi+ 1 , t))
2

xi− 1

(22)

2

2

as well as the temporal integration between tn and tn < tn + 1 represented by

xi+ 1
2

xi− 1

2

U (x, t

n+1

xi+ 1

2

)dx =

tn+1

n

U (x, t )dx+
tn

xi− 1

tn+1

F (U (xi− 1 , t)dt−
2tn

F (U (xi+ 1 , t))dt (23)2

2

It is possible therefore to represent equation 21 in integral conservation form [ZM10]. By
xi+ 1
1
n
2 u(x, t ) dt
deﬁning the variable in a cell Uin as equal to some integral average within the cell Δx
x 1
i− 2

and the ﬂux through an interface, Fxi+ 12equation can rewritten 21 as:1
as given by the expression Δt

Uin+1 = Uin +Δt(F 1 − Fi+ 1 )2
Δx i− 2

tn+1
F [U(xi +1/2,t) ]dt,
tn

(24)

The above mathematical expression has been derived purely from calculus, without any
applied constraints imposed so far. It has been assumed that knowledge of the analytical
functions U (x, t) and F (x, t) exists. [ZM10] . However what is of interest to us is a numericalformulation, which is bounded by some ﬁnite domain. Hence there is need to determine theconstituents of the expression in equation 24 above by some form of estimation. Numerical methods diﬀer in the way these constituents are evaluated. Bearing this in mind, apresentation of some select numerical schemes studied in the course of this work is done.

16

2.3

TYPES OF NUMERICAL SCHEMESA description of the two types of numerical methods mentioned above is presented: Godunov
based methods and central schemes, with some examples.

2.3.1

GODUNOV’S METHOD (FIRST-ORDER UPWIND METHOD)

We note the expression in ﬁgure 24 evolves from a mere mathematical expression to a numerical scheme [ZM10] when we make approximations for the evaluation of the numerical ﬂuxes
and cell averages, Ui . Fundamental to Godunov’s (First order) approach is the idea that at
a certain time step, n, the control volume has a piecewise constant grid data [God59] of the
xi+ 1
1
n
2 u(x, t ) dx. Doing this results in a loss of some level of cell data accuracy,
form uni = Δx
x 1
i− 2

which gets improved upon by utilizing higher order reconstructions, of which a number ofmethods exist. Sergei Godunov [God59] is acknowledged to be the ﬁrst to successfully carryout a conservative extension of the ﬁrst–order upwind scheme of Courant, Isaacson and Rees
(CIR) [CIR52] to nonlinear systems of hyperbolic conservation laws. We note that the main
highlight of upwind methods is that the discretization depends on the sign, and hence direction of wave propagation [ZM10] . Godunov’s ﬁrst–order upwind method is a conservative
method of the form of equation 21, and numerical ﬂuxes (Fi− 1 , Fi+ 1 ) are evaluated by
2
2
solving a local Riemann problem.

2.3.2

MUSCL-HANCOCK APPROACH FOR HIGHER RESOLUTION

Van Leer [VL76] , is credited with the idea of modifying the piecewise constant proposed by
Godunov as a necessary step in achieving high resolution schemes (second order and above).
Thus one can for instance perform a linear interpolation of cell data to the interface. His
idea is known today as the MUSCL (Monotone Upstream-Centered Schemes for conservationlaws) or variable extrapolation approach [VL76] , for conservation Laws. Godunov’s orderbarrier theorem implies that one cannot have a scheme higher than ﬁrst order that is free fromspurious oscillations [God57]. Therefore, to enjoy the beneﬁt of higher resolution methods,
they must be used together with a limiter. A piecewise linear, local reconstruction of cellΔii)data, uni is achieved using the linear expression ui (x) = uni + (x−xΔi , and the ratio ΔxrefersΔx1
to a slope of ui (x) in cell Ii . At the extreme ends , uLi = ui (0) = uni − 2 Δi , uR
i = ui (Δx) =1nui + 2 Δi are usually called boundary extrapolated values. Toro [Tor] states that mostapproaches neglect the generalized Riemann problem, utilizing the boundary extrapolated

17

values uLi , uR
i for each function of ui (x). Therefore, the result is a constant data Riemannproblem;⎧⎨ uR ,
if x < 0.
iu(x, 0) =(25)⎩uL , if x > 0.
i+1
for the cells. The Muscl-Hancock approach requires that these extrapolated values are
evolved by half a time step according to the relations below
1 Δt[F (uLi ) − F (uR
i )]
2 Δx
1 Δt
R
[F (uLi ) − F (uRūRi = ui +i )]2 ΔxūLi = uLi +

(26)
(27)

after which the solution is evolved to the next time step by evaluating the piece-wise constantdata Riemann problem. Simulation of test cases for the MUSCL scheme for the Burgers
equation using the same discretized domain as with the Godunov ﬁrst order test cases toevaluate the eﬀectiveness of diﬀerent numerical schemes were also conducted and have notbeen included in this work. Octave codes have been uploaded to GitHub and are available
on request.

2.4

CENTERED SCHEMES (LAX-FRIEDRICHS FAMILY OF SCHEMES

A characteristic trait of Godunov-type schemes is that the evaluation of ﬂux quantities, anecessary step in solving ﬁnite volume based problems, involves the solution to a Riemannproblem. This demands the characteristic decomposition to enable the evolution of the
problem in space and time. However, the centered schemes achieve the solution of a quantity
at the next time step without solving a Riemann problem, rather they do this by integratingover the Riemann fan. This integration implies that both the left and right travelling waves
are catered for in this method. Hence there is no need to solve a Riemann problem at every
cell interface, thus these methods could be computationally inexpensive. A presentation ofsome of the schemes belonging to this family and their numerical formulations.

2.4.1

LAX-FRIEDRICHS SCHEME

The Lax-Friedrichs scheme applied to the linear advection equation [Tor] is given as
=
un+1
i

1 + CF L n
CF L − 1 n(ui−1 ) +(ui+1 )
2CF L
2CF L
18

(28)

where CFL is the Courant-Friedrichs-Levy number (CFL). Given the 1D conservative equation for linear advection,
ut + a(u)x = 0(29)
the ﬂux term is evaluated as
n+ 12

Fi

=

1 + CF L
CF L − 1
F (uni ) +F (uni+1 )
2CF L
2CF L

(30)

and its solution at a cell interface, i can be obtain. Looking at equation 28, it is seen that
it is a weighted integral average of a corresponding Riemann problem with diﬀerent left andright states at time t = 12 Δt. For a nonlinear system of conservation laws
Ut + F (U )x = 0(31)

considering the scheme to be an integral average within cells, a conservative ﬂux generalization for the LaxFriedrichs scheme yields :1
1 Δx nn
Fi+ 1 = (Fin + Fi+1 ) +(U − Ui+1)
2
2
2 Δt i(32)

and hence obtain the general solution of the Lax-Friedrichs scheme evolved in time as
1 n
1 Δt nnn
Uin+1 = (Ui−1(Fi−1 − Fi+1
+ Ui+1
)+
)
2
2 Δx

(33)

Equation 32 is the numerical ﬂux for the Lax-Friedrichs scheme when applied to systems ofconservation laws. It can be easily noted that the ﬂux evaluation and the evolution of thesolution in time is done without solving a Riemann problem.

2.4.2

NESSYAHU-TADMOR (NT) SCHEME

Consider again the conservation equation in 1D given below:dU∂+(F (U )) = 0
dt
∂x

(34)

As stated earlier, Godunov type schemes are upwind in nature. This can be demonstrated
as follows: if the characteristic speeds within all the neighbouring ﬁnite control volumes are

19

positive (or negative), an integration of equation 34 reduces touj (t + Δt) = uj (t) − λ[F (uj (t)) − F (uj−1 (t))](35)

uj (t + Δt) = uj (t) − λ[F (uj+1 (t)) − F (uj (t))](36)

where the second equation corresponds to the case with negative speed [NT90]. In thisdt
convention of NT scheme, λ = dx
. For cases where there are existing waves with bothsigns for the characteristic speeds (right and left going waves) in the control volume, thisbecomes complicated and the ﬂux evaluation demands the ”direction of the wind” is taken
into account [NT90] so as to diﬀerentiate both waves within the Riemann fan, hence theexact solution to the Riemann problem in this case becomes challenging. Nessyahu et al.
[NT90] proposed the NT scheme as a way of by-passing the need to solve multiple Riemannproblems for such cases. Starting from the staggered version of the Lax-Friedrichs scheme
given as:1
uj+ 1 (t + Δt) = [vj + uj+1 ]λ[F (uj+1 (t)) − F (uj (t))](37)
22they developed a central diﬀerence scheme that integrates over the Riemann fan. We notethat such integration will cater to waves in both directions. Hence, it eliminates the needto solve a Riemann problem. However since the Lax-Friedrichs scheme is quite diﬀusive andﬁrst order, they introduced MUSCL type interpolants thereby making the scheme of higher
order(second order). The reader is referred to [NT90], however a summary of the process
results in the NT scheme described by a ﬁrst order predictor step:uj (t +Δt1
) = uj (t) − λFj
2
2

(38)

and the second order corrector step
1
1
Δt
Δt
uj+ 1 (t + Δt) = [uj (t) + uj+1 (t)] + [uj − uj+1 ] − λ[F (uj+1 (t +) − F (uj (t +)))] (39)22822where the prime, ” ” indicates a derivative. The non-oscillatory properties of the scheme
have been demonstrated in [NT90] . We note that there are two versions of the scheme: thestaggered and its non-staggered forms. Codes were implemented for both versions in Octave.

2.4.3

DESCRIPTION OF THE KNP NUMERICAL SCHEME

The reader is referred to the works of Kurganov et al. in [KBS15, KNP01, KT00], where adetailed description of the mathematical formulation of the second order numerical scheme
and its implementation in OpenFOAM has been described. While the KT is a central
20scheme based on the Lax-Friedrichs scheme, the KNP, its generalization, is central- upwind;central since it involves integration over the Riemann fan with no decomposition required,and upwind because it takes into account the wave direction of travel by measuring the”one-sided local” speeds of propagation which for a one dimensional case results in the
eigenvalue of the ﬂux Jacobian, as opposed to the spectral radius for KT, which is lessprecise [KNP01]. Our simulations used the KNP scheme. As mentioned earlier, conservative
formulation of the Euler equation is of particular interest,in which case it is mandatory thatsome ﬂux evaluation technique is used to evolve the equation in time. The ﬁnite volume
formulation for ﬂux evaluation through a cell interface for KNP as presented in [KBS15] istherefore highlighted in this section. The system solved is usually the reactive inviscid Euler
equations in conservation form (here in its 1D, single reaction form):∂U∂F+= S(G)
∂t
∂x
where



U = ρ, ρu, ρ(e + u2 /2 + yQ), ρy

(40)

(41)

is the vector of conserved variables, with the primitive variables being density, (ρ), velocity(u), internal energy per unit mass (e), and extent of reaction y. The vectors

and



F = ρu, ρu2 , ρ(e + P/ρ + u2 /2 + yQ), ρuy

(42)



s(U) = 0, 0, 0, r

(43)

are the ﬂuxes, and the source terms associated with the one-step Arrhenius chemical reaction
rate, r. The system is most often simulated employing ﬁnite volume methods, with numerical
schemes capable of handling large, rapid gradient changes, i.e., shock waves. The implementation of the algorithm has been discussed above. The KT/KNP scheme is a Riemann-free
ﬂux evaluation scheme where the ﬂux, F (U ) of the quantity U , is computed as
P P
P min
N N
P min
F (U ) = ψR φR = ψRP ((αRφR + α R
αR ) + ψRN (αR
φR − α R
αR )).

(44)

where the above is in the original KNP notation from [KBS15] as implemented in OpenFOAM. The subscript R implies that a quantity is computed at the cell interface, and thesuperscripts, ”P ” and ”N ”, indicate two diﬀerent cells, from which the cell values are interpolated to the cell face, i.e. the two sides of a cell face as represented in ﬁgure 6.ψ is thevector of conserved quantities, U, S is the cell face normal area, and phi is ṁ/ρ = V · S, the21

volumetric ﬂow rate. In a more familiar notation, 44 can be rewritten as:
P
P min
N
P min
F (U ) = URP ((αR(VRP · S)PR + αRαR ) + URN (αR(VRP · S)N
R − αR αR )).

(45)

P min
The quantity (αRαR ) can be deﬁned as ωf , a product of a weighting and the volume ﬂux,
P
max
max
min
N
min
max
min
where αR = αR /(αR+ αR), and αR
= αR
/(αR
+ αR
). We note that for the KT
Pvariant, the value of αRand ωf is 0.5, making the scheme central.

Figure 6: Representation of the scheme in ﬁnite volume formulation

The wave speeds are evaluated as follows:max
αR= max(φPR + cPR |SR |, φPR + cN
R |SR |)

(46)

min
N
= −min(φPR − cPR |SR |, φNαR
R − cR |SR |)(47)

or in more familiar notation in the detonation community,max
αR= max(VRP · |SR | + cPR |SR |, VRN · |SR | + cN
R |SR |)

(48)min
= −min(VRP · |SR | − cPR |SR |, VRN · |SR | − cN
αR
R |SR |)

(49)

where V represents the macroscopic or particle velocity within a given cell, c is the thermodynamic sound speed, (cPR = γR TRP , cN
γR TRN ) and SR is the surface area normal ofR =the interface. The details of the scheme are available in [KBS15]. The conditions at the cell
faces are evaluated using a linear, slope-limited interpolation. For non-reactive cases, theinviscid Euler equations are solved, while for the reactive cases, a single step, irreversible
reaction with heat release (Q) and a rate r = ke[−Ea/RT ] is coupled to the system and solved.

22

2.5

OVERVIEW OF OpenFOAM

2.5.1

WHAT IS OpenFOAM

Simulations for this work have been based on the OpenFOAM framework, commonly referred
to as OpenFOAM. It is a library of over 100 C++ codes [Ope16] that can be grouped together
as needed, to create executables, referred to in OpenFOAM as applications. OpenFOAM
has over 250 of them pre-built [Ope16]. These can be categorized into two: solvers heavily
employed for problem formulation and analysis in continuum mechanics, and utilities that
execute some form of data manipulation. OpenFOAM is such that users are allowed thefreedom to design their own custom solvers and utilities that relate better to their desiredcase or speciﬁc problem. This is one of the major reasons it is gaining traction in the CFD
community [WTJF98]. A representation of the basic structure of OpenFOAM is shownin ﬁgure 7, where it can be seen that there are three phases to an OpenFOAM problem:
pre-processing, solving and post-processing. OpenFOAM, when installed, comes with a

Figure 7: Overview of the structure of OpenFOAM [Ope16]

case set-up and result analysis section, referred to as the pre-processing and post-processing
environments, which in actual terms employ the use of these utilities in their operation. For
instance, to deﬁne a hexagonal mesh, one will have to specify the co-ordinates in a text ﬁlenamed blockMeshDict, which can be called by the utility ”blockMesh” to generate the meshdata. One can then use, for instance, the setFields utility to specify initial conditions within
two regions of the domain as with a basic multiphase ﬂow problem for instance. The wholeessence of this is to guarantee some form of consistency particularly with regards to datautilization across all environments [Ope16].

23

Figure 8: OpenFOAM case formulation tree [Ope16]

2.5.2

FILE STRUCTURE OF OpenFOAM

To run an OpenFOAM case, a basic ﬁle description is presented in ﬁgure 8. Included here
as well are some utilities to help facilitate understanding, particularly as they relate to thisthesis. We see that the we can deﬁne an OpenFOAM case or problem by populating the 3
major ﬁle folders or directories included in ﬁgure 8. These will be discussed, making reference
to a shock tube case that has been built in OpenFOAM. We go over these 3 folders described
below.
0 FOLDER(TIME DIRECTORY)
The time directory contains the initial and boundary conditions of ﬁeld variables we are
seeking solutions to, or those they depend on. e.g. we may not need to initialize density since
we can obtain it by specifying temperature and gas constant. Initial ﬁelds in OpenFOAM
simulations must be speciﬁed as a requirement even for cases that are within the domainof steady ﬂow. We note that the directory is called a 0 folder as we ordinarily begin our
simulations from a ”0” time. However as we advance in time, OpenFOAM writes out othertime folders based on a user speciﬁed output time and interval. Referencing our shock tube
case, we need to evolve the velocity, pressure and temperature as we advance in time. Hencein our zero (0) folder, we have speciﬁed the U, P and T ﬁles which contain the initial andboundary conditions for the case. Figure 9 shows a sample pressure ﬁle, known as a P ﬁle,24

Figure 9: Pressure ﬁle at time 0 (P)for our case. A number of things are striking. We see an array speciﬁed as dimensions. In
OpenFOAM, dimensions are speciﬁed by considering the S.I unit. The ﬁrst element of thearray is always the unit of mass, in kg, and the next being length in meters. The rest of theelements are time (s), temperature (K), number of moles (mol), current (A) and luminous
intensity (cd). Thus, derived quantities are a combination of two or more of these as seenin the unit of pressure kg/ms2 . OpenFOAM data is speciﬁed using dictionaries: an entity
containing data entry(ies) that can be called upon by utilizing keywords [Ope16] . Most
OpenFOAM data ﬁles are themselves dictionaries, and can contain one or more data entries.
Thus the key word ”InternalField” is a dictionary containing one keyword, and is assigned
the value 0 to represent the pressure values within the domain initially. The boundary ﬁeld
is a dictionary containing the boundary condition for the ﬁeld variable (P in this case).
Here we see that a ”zeroGradient” boundary condition has been assigned to the ﬁeld at aregion deﬁned as ”sides” and a boundary named empty is assigned an ”empty” boundary
condition. There are a host of boundary conditions available in OpenFOAM. In the casediscussed above, the internal ﬁeld is constant or uniform. OpenFOAM also allows for nonuniform internal ﬁelds. Thus we can use a non-uniform ﬁeld obtained from another software(e.g. octave) , or from an existing simulation, as initial data for another case.
CONSTANT FOLDER
A basic constant folder will contain 3 elements: A thermophysical properties (thermophysicalProperties in OpenFOAM) ﬁle, a turbulence properties (turbulenceProperties in OpenFOAM) ﬁle , and a folder known as the PolyMesh. A thermophysical properties ﬁle describes

25

the thermophysical models utilized in the OpenFOAM case in question, where the thermalenergy, compressibility and/or mass transfer is of particular interest [Ope16] . These properties can be constant, or can be a function of thermodynamic properties such as T,P. The
thermophysicalProperties dictionary is called upon by any solver that utilizes the thermophysical model library [Ope16] . OpenFOAM creates these models as pressure -Temperaturesystem which in turn can be used to create other models. A lot of these can be found in
OpenFOAM and are readily described in the user guide. The thermophysical properties(keyword thermophysicalProperties) begins with a ”thermotype”, which deﬁnes the basic
equation of state and models for ﬂuid transport and mixture modelling among others. The
turbulenceProperties ﬁle contains the turbulenceProperties dictionary wherein we can specify the turbulence model we wish to use. For our non-turbulent cases, the ”simulationType”
dictionary is always set to laminar, and inviscid ﬂow is enforced by specifying a ”0” value ofviscosity in the chemical properties ﬁle. The polyMesh folder is home to the mesh. Simple
OpenFOAM meshes are generated by deﬁning a blockMeshDict ﬁle (usually found in thesystem folder) which, when called by the blockMesh utility, generates the mesh as speciﬁedin the ﬁle. Here, we can specify the co-ordinates of the domain, and indicate the meshblocks and desired number of cells within each block. Once we run the ”blockMesh” utility,
OpenFOAM populates the PolyMesh folder with ”points”, ”Faces”, ”owner” ”neighbour”
and ”boundary” ﬁles, where Points refers to a vector of cell vertices, faces being a list of
all the existing faces, with each face being made up a list of vertices sourced from ”points”.
Owner and Neighbour represent owner and neighbour cells respectively. The boundary folder
describes a list of patches, that is being a list of faces grouped together to act as a boundary.
E.g, an inlet boundary is a group of faces that have been formed into a patch.
SYSTEM FOLDER
The system folder contains utilities that are used to carry out mesh manipulation and control simulation cases. Some dictionary entries found here include the blockMeshDict, whichhouses the geometrical data, the controlDict where we can specify the time step, CFL criteria among others, and the fvSchemes where we can describe the choice of discretization,interpolation and reconstruction methods for ﬁeld variables. A useful dictionary to highlight
as well is the decomposeParDict which enables one to decompose a domain into a numberof units which can then be run in parallel on diﬀerent processors. This is a useful tool forrunning OpenFOAM cases in computer clusters.

26

2.5.3

POST PROCESSING

OpenFOAM comes with an open source visualization tool known as paraFoam, (a fork ofparaview) which enables visualization of simulation results. Usually we tell OpenFOAM thetime it should output our results; for instance we can say every 1 × 106 s of run time. We canthen visualize these results and produce animations based on them. We can also perform
numerous operations on the results with the help of various tools it comes with. We can also
compute derived ﬁelds which is a function of two or more output ﬁelds.

2.6

CREATION OF CUSTOM SOLVER FOR DETONATION SIMULATIONS (rhoCentralFoamreac)

OpenFOAM consists of a signiﬁcant amount of ”solvers” which refers to a collection of C++
codes written and compiled to solve a given type of problem. These include among others,
combustion solver, compressible and imcompressible ﬂow solvers, etc. Detonation waves
involve a complex coupling of shocks and reactions, hence any detonation solver must be able
to capture shocks eﬀectively and also have the ability to solve for chemical reactions. Thereexist compressible ﬂow solvers e.g., rhoCentralFoam, sonicFoam, etc. The details of thesesolvers are described in [Ope16]. In OpenFOAM5, however, there is not a ﬂux evaluating
solver for compressible ﬂows with chemical reactions that can capture this phenomenon
eﬀectively. Of all the compressible ﬂow schemes available, only the rhoCentralFoam solver
evaluates the ﬂuxes in a true ﬁnite volume manner. The rhoCentralFoam solver is writtensuch that it can be employed for cases whose ﬂow behaviour can be described by the Euler
equations. However a solver was needed for the reactive Euler equations which describe thebehaviour of detonation waves. Hence a new solver was created for that purpose which wasnamed rhoCentralFoamreac. This needed to be compiled in all servers and or clusters where
it was utilized. Starting oﬀ with the rhoCentralFoam solver already in existence, a copy of
the C++ based codes that included the various header ﬁles necessary to successfully executethe application was created. Considering it was not intended to solve for full chemistry
as the focus was in one-step governed reaction kinetics, there was need to modify the code
executable to account for this. First the transport equation was included into our code which
is given by the equation below;
∂ρuY
∂ρY
+
=0
(50)
∂t
∂x
where Y accounts for ﬂuid transport. To evaluate this, a shock tube case was run with the
newly modiﬁed solver to see how it behaves with the inclusion of the transport term Y. This

27

Initial States(a) P at t = 0(b) U at t = 0(c) T at t = 0(d) Y at t = 0

State at t = 0.007(e) P at t = 0.007(f) U at t = 0.007(g) T at t = 0.007(h) Y at t = 0.007

Figure 10: Inclusion of Transport Equation to rhoCentralFoam. Fig10h shows the advection of”Y” after some time

was done by modifying the source code to enable us solve for an imaginary property being
convected through our domain, deﬁned by the transport equation, and then compiling a newsolver know as rhoCentralFoamnor using the openFoam ”wmake” command. It can be seenthat for a given transport quantity when solved using this solver, based on certain initialcondition, the property gets advected through the domain, and the simulated result gives
the expected behaviour as shown in ﬁgure 10. This inclusion was done such that the user
has the option of employing the KNP or the KT scheme (its slightly more diﬀusive version)as implemented in OpenFOAM. The ﬂux terms of the transport equation were evaluatedusing the scheme as well after interpolating cell values to faces following the implementation
described in [KBS15]. With the solver modiﬁed to handle ﬂuid transport, the reactive
source terms were then included. Simulations were based on the one step Arrhenius kinetics
equation, hence the source code was modiﬁed to account for that. The procedure as well asa brief overview of how the solver works has been enumerated in 7:
1. Create copy of rhoCentralFoam solver.
2. Create Y in the zero folder.
3. Rename source code and modify to include transport terms.
4. Modify createFields.H to include volScalerField Y to be obtained from 0 folder.
5. Modify ”createFields.H” to include conserved volScalerField rhoY, (ρ × Y ).
6. Include concentration change and integrate rate law equations rhoY ∗(e(−(dt/rkA)∗e1) and correct boundary condition
28(−Ea∗psi) )−7. Compile solver by running ”wmake” within the directory using a terminal(Ubuntu inthis case) .
The new solver, rhoCentralFoamreac was then used to validate the KNP/KT scheme fordetonation simulation, and to perform RDE simulations, as discussed in the following chapters.

29

Chapter 3
VALIDATION OF THE KT/KNP
SCHEME FOR DETONATIONSIMULATIONThe open source Computational Fluid Mechanics (CFD) package, OpenFOAM, as it is commonly known, is increasingly used and referenced. One drawback of the stock OpenFoam
package is that the only ﬁnite volume numerical scheme available for the solution of the Euler equations in conservative form is the KT, and its variant, KNP, numerical scheme. This
family of scheme is 2nd order accurate, and is based on the idea behind the Lax-Friedrichs
scheme which does not involve the solution to a Riemann problem in order to evaluate theintercell ﬂuxes. This is unlike the methods currently used in detonation research, which aremainly Godunov-type schemes with an approximate Riemann solver [Tor] at each cell interface. Common modern ﬂux evaluation methods encountered include the class of Weighted
Essentially Non-oscillatory Schemes (WENO), and schemes derived from the Harten-Lax-van
Leer (HLL) method, such as the Harten-Lax-van Leer-Contact (HLLC). OpenFoam, with
the KNP scheme, was recently used to simulate the two-dimensional structure of detonation
waves in [MTE17b] despite having not been fully validated for the detonation simulation.
Figure 11a and 11b show the 1D simulated detonation structure and the 2D simulated detonation cells obtained in that work. They created a solver known as rhoCentralRfFoam using
the OpenFOAM framework which they used to simulate planar detonations in H2-O2-Ar,
and also [MTE17a] simulated a two-dimensional detonation cellular structure in H2-O2Ar mixtures with OpenFOAM employing full chemistry. While qualitatively correct, thoseresults are based on a scheme that was not validated for detonation simulations. The validation of the KT/KNP scheme for detonation waves is done by simulating a 1D problem
whose expected result is a pulsating detonation with a single oscillation mode. This is done30(a) Concentration vs space for a H2 -02 -Ar Mixture(b) Simulated Detonation cell sizes

Figure 11: Images are from simulation of detonation waves in [MTE17b]

by carrying out reactive and non-reactive validation studies. Useful metrics for detailed
comparison and convergence studies are the oscillation peak pressure and period. Firstly,
an examination of how it handles non-reactive ﬂows with shock waves is performed, andthereafter, the 1-D dynamics of a propagating detonation wave is investigated.

3.1

NON-REACTIVE VALIDATION

3.1.1

CASE DESCRIPTION AND SET-UP

This approach was ﬁrst implemented using a one-dimensional shock tube problem implemented within the rhoCentralFoam solver in OpenFoam. This is a compressible non-Riemannbased solver. The case was set up building from the already existing shock tube case in the
OpenFoam tutorial. The domain was set to be 0.1m long. This was meshed with 10000 cellswith uniform cell distribution across the domain. A sketch of the domain with boundary
conditions is shown in ﬁgure 12

3.1.2

BOUNDARY AND INITIAL CONDITIONS

A transmissive (zero gradient) boundary condition was applied on the sides for pressure,temperature and velocity, while the other sections of the domain, e.g., front face (relevant for2D and 3D) were speciﬁed as empty to enforce a one-dimensional simulation in x. The density
is evaluated using the perfect gas equation of state. Two regions of distinct, thermodynamic
31

Figure 12: Case set-up for non-reactive validation (shock tube problem)

Figure 14: Simulated State (t = 0.007s)

Figure 13: Initial State (t = 0.000s)

properties were speciﬁed as initial condition: a driver and a driven section within the tube,
separated by a contact surface. The initial pressure was 100 kPa and 10 kPa, initial velocity
is 0 everywhere, initial temperature was set to 348.432 K and 278.746 K for the driver anddriven ends respectively. The mixture inside the tube is treated as an ideal gas with zeroviscosity, a moleculer mass of mixture (M ) of 28.98 g/mol and a heat of formation (Hf ) of
2.544e+06 kJ/kmol. The performance of four diﬀerent limiters: Minmod, Van Leer, Van
Albada and the superBee limiters in OpenFOAM, are evaluated. The simulation was run
for 0.007s and the initial and end states presented in ﬁg 13 and ﬁg 14).

3.1.3

RESULTS

Figure 16 presents images obtained from paraview, prior to and at the end of simulation
cases for non-reactive validation. This indicates the four diﬀerent states obtained in thesolution to the shock tube problem. where ﬁgures 16a, 16b, 16c and 16d are the resulting
solutions obtained by specifying P , U and T initially as shown in ﬁgures 15a, 15b and 15c.
ρ is obtained from the ideal gas equation of state as stated earlier.
For each limiter, the maximum CFL number allowed in a computation was gradually increased, starting with an initial value of 0.05. For low enough CFL numbers, the solution
32

(a) Pressure(b) Temperature(c) Velocity

Figure 15: Initial condition for P , U and T (t = 0.000 s).

(a) Pressure(b) Temperature(c) Velocity(d) Density with patch names

Figure 16: End state showing P,U,T and ρ (t = 0.007 s)

33

Figure 17: Pressure vs position for the shock tube problem using the minmod limiter CFL = 0.09

behaves as expected, with four apparently constant states: the low pressure state on the
far right, the high pressure state on the far left and the two intermediate states around thecontact surface. In ﬁgure 17, the distribution of pressure against space is shown at the endof the simulation, for the minmod limiter with a CFL number of 0.09. As the CFL number
is increased further, non-physical oscillations are readily visible and degrade the solution,
as shown in ﬁgures 18a and 18b, for the minmod limiter with CFL of 0.23 and 0.33. It is
important to determine under which conditions those unphysical oscillations appear, since it
is important to understand how the KT/KNP scheme performs when simulating unreactive
ﬂows, before it can be utilized for simulations involving detonation wave . Careful examination of the solutions obtained show these unphysical oscillations originate from the contact
surface and propagate both towards the shock front and through the expansion fan. It canbe noted that there are measurable oscillations even for cases that appear oscillation free,such as that shown in ﬁgure 17. The maximum amplitude of the oscillations, in the pressureﬁeld, after a simulation has been run, are taken as a metric for the level of degradation of
the solution. As all the computed cases are at the same initial conditions, and for the same
duration, Δt, the absolute disturbance maximum amplitude is used in dimensional scale. An
example of a solution dominated by unphysical oscillations is shown for the case of superbee
limiter in ﬁgure 19a and ﬁgure 19b, with CFL = 0.05 and CFL of 0.16. It is evident that34

(a) Pressure vs cell position, (CFL = 0.23)(minmod)(b) Pressure vs cell position, (CFL = 0.33)(minmod)

Figure 18: Flow simulations at diﬀerent CFL numbers(a) Pressure (Pa) vs cell position, (CFL = 0.05)(superbeee)(b) Pressure (Pa) vs cell position, (CFL = 0.16)(superbeee)

Figure 19: Flow simulations using superbee limiter, indicating intense amplitude of oscillation

the choice of limiter can be signiﬁcant as even for very low CFL number as shown in ﬁg 19a
there exists large amplitude oscillations. In the case shown in ﬁg 19b, the unphysical oscillations have even reached the right propagating shock and have started to inﬂuence the highpressure state on the left side of the domain. Thus, an examination of these oscillations andthe corresponding eﬀect changes in pressure ratios between the driver and driven pressures
could have on them was performed. In ﬁgure 20, a presentation of the pressure disturbance
amplitude as a function of CFL number for a pressure ratio of 10 for the 4 limiters considered is shown. For the test conditions evaluated, the minmod limiter is the best performer
while the superbee limiter fails to accurately compute the solution at almost all the CFL
numbers considered. There are easily recognizable outlier points, for example at CFL = 0.11
and 0.17 for the minmod limiter. Those correspond to cases where a numerical instability
in the platform resulted in a failure of the code to complete the simulation. The sourceof this error has not been successfully identiﬁed and may not originate with the numericalsolver itself. The van Leer and van Albada limiters perform only slightly worse than theminmod limiter, except above limiting CFL numbers of 0.12 and 0.22 respectively. Abovethese critical values, the algorithm fails to compute a solution to the problem that is not35

dominated by those spurious oscillations. It can be seen that, even for the minmod limiter,
to maintain an accuracy of 1% in the states on either side of the contact surface, a CFL
number below a value of about 0.15 should be used. This critical CFL is very low. With
more commonly used solvers in the detonation ﬁeld, the CFL number is most often in therange of 0.5–0.8. The process was repeated for various pressure ratios, each time evaluating

Figure 20: Disturbance amplitude DP (Pa) vs CFL number for the minmod, van Alabada, vanLeer and superBee limiters for a pressure ratio of 10.

the amplitude of oscillation resulting from the numerical method. The results for a pressure
ratio of 1.5 and 2 are shown in ﬁgure 21 and ﬁgure 22, and that for a pressure ratio of 5presented in ﬁgure 23.

Figure 21: Disturbance amplitude DP (Pa) vs CFL number for the minmod, van Alabada, vanLeer and superBee limiters for a pressure ratio of 1.5.

36Figure 22: Disturbance amplitude DP vs CFL number for the minmod, van Alabada, van Leerand superBee limiters for a pressure ratio of 2.

3.2

REACTIVE VALIDATION

Having determined the conditions under which the KT/KNP scheme can reproduce the shocktube validation case, the simulation of reactive compressible ﬂow can be evaluated. This isdone by simulating a standard detonation case often used in the literature as a validationcase, corresponding to a one-dimensionally unstable, over-driven detonation. The chemicalrate parameters are Q = Ea = 50RTo , adiabatic index (γ) = 1.2, detonation overdrive2
factor f = DD2 s = 1.6, where DCJ is the detonation velocity at CJ point (DCJ ), and k isCJ
the pre-exponential constant 1/s (k), which is arbitrary. It was chosen here to correspondto the combustion of methane and oxygen, such that the half reaction zone length (L1/2 ) =1.9 × 10− 4 m. As described in chapter 2, the 1-step Arrhenius kinetics is used to model thereaction.

3.2.1

CASE SETUP DESCRIPTION

An over-driven, unstable detonation is simulated by considering the case of a piston pushingon the gas at a velocity higher than the CJ particle speed. The ﬂow is simulated in theframe of reference of the piston. Two initial conditions were considered:
1. A small, high pressure and temperature initiation region is initialized near the pistonand the dynamics of the initiation process are simulated;
2. The ZND structure is initiated near the piston and the wave instability is allowed todevelop.
37

Figure 23: Disturbance amplitude vs CFL number for the minmod, van Alabada, van Leer andsuperBee limiters for a pressure ratio of 5.

The domain size was set to be a multiple of the half-reaction zone length, most often 200L1/2(0.038 m in physical terms) to 600L1/2 (0.101410 m).
Boundary Conditions
The piston face is treated as a solid wall. The opposite end of the domain is an inlet providing
fresh reactants ﬂowing in at the imposed particle velocity, U = -1811.6 m/s

3.2.2

INITIAL CONDITIONS

Fig 24 shows the overall initial conditions. In a small region near the piston, a driver section
was initiated.
INITIATING PISTON
The driver section was initially static (U = 0), with a high pressure of 100 kPa, and atemperature of 298 K. The heat release term, Y , was set to 4.2 × 106 KJ/Kg everywherein the domain. The velocity elsewhere was set at the same imposed particle velocity U =-1811.6 m/s.

38

Figure 24: Case set-up for reactive validation (piston problem)

1-D ZND STRUCTURE
An alternative initiation method that bypasses the initiation transient associated with thelocalized high pressure and temperature driver section consists in initiating, near the piston,
the 1D ZND proﬁle of the steady wave and letting instabilities grow. The details of the ZND
model can be found in [FD79]. In summary, the gas is compressed by a strong shock to a highpressure, the von Neumann spike, followed by chemical reaction. The reaction rate occurs
within a ﬁnite time frame and over a certain physical extent, instead of in an instantaneous
process. For a steady wave, a reaction progress variable can be used to describe the state
points completely along the Rayleigh-Hugoniot curve. Thus, at each point in the reaction,
the value of pressure, P , speciﬁc volume, v, and velocity, U can be obtained. An Octave code
was written to integrate the 1D Euler equations for a steady detonation wave and calculate
the physical extent between the front of the shock and the point where half of the reaction iscomplete, L1/2 . This integration is performed assuming a single Arrhenius chemical reaction
rate with speciﬁed values of heat release Q and activation energy (Ea). This wave proﬁle isthen used as the initial condition within the driver section near the piston. The ZND wave
proﬁle was calculated with the same spatial grid extent as the computation considered, i.e.
L1/2 /Δx = 10, 20, 40, etc. Sample results of the ZND structure are shown in ﬁgure 25 andthe Octave code used to generate the ZND structure is reproduced in the appendices section.

3.2.3

SOLUTION METHOD

The 1D CFD Euler calculations were performed with varying grid spacing corresponding to
fractions of the half-reaction zone length, i.e. Δx = L1/2 /C, where C = 10, 20, 40, 80, etc up

39

to 100. The OpenFOAM postprocessing utility sampleDict was used to produce 1D results
data sets (remember that OpenFOAM always, by default, produces 3D data sets) that werethemselves post-processed using a python script to extract the shock pressure as a functionof time and space.

3.3

RESULTSThe results obtained from the piston set up and the ZND initiated cases are presented. The
domain was also extended to 600L1/2 (0.101418 m) from the initial 200L1/2 used for thepiston problem. In each simulated case, the mesh size was adjusted to reﬂect the desired
number of cells per half reaction zone length, e.g 12000 cells in x gives 20points per L1/2 forthe mesh 600L1/2 long. 6000 cells corresponds to 10 points, etc. The peak shock pressure,(Pshock ) is plotted against time and space. It can be seen that the scheme reproduces theoscillatory structure of a detonation wave as shown in Figures 27, 29, 31, 33, 35, and 37
which shows results obtained for 10, 20, 40, 60, 80 and 100 points per half reaction zone
length initiated using results from the ZND code with a CFL number of 0.1. After a certain
time, t, the detonation wave hits the domain boundary, at which point it no longer is adetonation wave, hence the reason for the artifact seen towards the right end of the domain
as in 34 and 38 for plots against time.
The domain length was doubled thereby allowing the detonation wave travel longer. This wasdone to see if the propagating detonation wave would result to a single mode of oscillation
at a distance further away. This however was still not the case as can be seen in ﬁgures 39
and 40 which show dimensionless pressure proﬁles plotted against space and time.

3.4

ANALYSIS AND CONCLUSION

For the simple, overdriven detonation problem, the KNP scheme recovers the oscillating
solution. Findings from a number of linear stability analysis [BM95, NNL00] show that theresulting ZND solution gives a single mode instability for the chosen values of heat release(Q) and Ea as shown in ﬁgure 41. However, a single mode of oscillation is not obtained, andeven after increasing the resolution up to 100 cells per half reaction zone length it still fails
to converge to a single mode. Instead, a proﬁle that is rather irregular with multiple modesis the resulting solution. Figure 42 shows the comparison of the oscillation peak pressure(Pshock ) for various other schemes as seen in [NNL00]. A comparison of the KT/KNP scheme,with other schemes is made. It can be noted that the peak pressure values obtained for eachcase varied with space, with the maximum values peaking above 100 kPa (≈ 108 kPa), quite
40

(a) Pressure, P against cell position in space(b) Temperature,T against cell position in space(c) Particle velocity, Up against cell position in space(d) Reaction rate against cell position in space

Figure 25: Result from octave codes for ZND model41

Detonation Pressure proﬁle for Δx = L1/2 /10 and Δx = L1/2 /20 for a domain 200 L1/2 long

Figure 26: Shock peak pressure vs x axis for Δx = L1/2 /10 and Δx = L1/2 /20 for a domain 200
L1/2 long

Figure 27: Pshock /P0 vs position for Δx = Figure 28: Pshock /P0 vs time for Δx = L1/2 /10
L1/2 /10
for
ZND initiated case (600 L1/2 long).
for ZND initiated case (600 L1/2 long).

42Figure 29: Pshock /P0 vs position for Δx = Figure 30: Pshock /P0 vs time for Δx = L1/2 /20.
L1/2 /20. (ZND)).
(ZND).

Figure 31: Pshock /P0 vs position for Δx = Figure 32: Pshock /P0 vs time for Δx = L1/2 /40.
L1/2 /40. (ZND) .
(ZND).

Figure 33: Pshock /P0 vs position for Δx = Figure 34: Pshock /P0 vs time for Δx = L1/2 /60.
L1/2 /60. (ZND) .
(ZND).

43

Figure 35: Pshock /P0 vs position for Δx = Figure 36: Pshock /P0 vs time for Δx = L1/2 /80.
L1/2 /80. (ZND) .
(ZND).

Figure 37: Pshock /P0 vs position for Δx = Figure 38: Pshock /P0 vs time for Δx =L1/2 /100. (ZND) .
L1/2 /100. (ZND).

Figure 39: Pshock /P0 vs position for Δx = Figure 40: Pshock /P0 vs time for Δx = L1/2 /80
L1/2 /80 for an extended domain. (ZND) .
for an extended domain. (ZND).

44

Figure 41: Image shows single mode pulsating detonation [NNL00].

Figure 42: Pshock /P0 vs Number of points per L 1 [NNL00].
2

a high value considering that most of the schemes discussed in [NNL00] show a convergence
to a peak pressure value of ≈ 98.6 kPa. Out of curiosity, a reduction in the CFL number forsome of the simulation cases to 0.05 was done, and and test cases carried out, to see if there
would be a diﬀerent result; that however was not the case. In conclusion, the ZND solutionof a detonation wave solved employing the KNP scheme cannot recover the regular, single
mode, periodic oscillation for an overdriven pulsating detonation with an overdrive factor of
1.6, which is usually the case with numerical schemes employed in problems that deal withhydrodynamic-chemical reaction problems, used for the detonations studies, e.g., the SLIC
scheme [NNL00].

45

Chapter 4
SIMULATION OF COLD ENGINE
START
4.1

PROBLEM DESCRIPTION/CASE SET UP

The evaluation process is continued by studying the RDE in this section. We see that it ispossible to obtain a circumferentially propagating detonation wave within an annular combustor by ﬂowing combustive mixtures through the chamber as demonstrated in [BLWC13]and [Sha15]. Figures 43 presents a visual representation of the concept as applicable torocket propulsion, and ﬁgure 44 shows a rotating detonation wave propagating circumferentially along the chamber. Though the KT/KNP scheme has not performed as we would
expect it to in the 1-D validation, it is possible that the detonation structure details areirrelevant to 2D RDE simulations. The suitability of OpenFOAM is also evaluated for thesesimulations.

Figure 43: Rotating Detonation Engine [LB14]
Figure 44: Representaion of an RDE [Sha15]

46

Figure 45: Image shows unwrapping of an RDE (3D to 2D).

Figure 46: 2D RDE representaion

Figure 47: 2D RDE representaion[TSS+ 14]

[KS13]For ease of analysis, it is common to see simpliﬁed, 2-D simulations of the RDE phenomenon
in literature. Towery et al. [TSS+ 14] for instance, show an unwrapped 2D model of the detonation wave. This approach is elaborated upon by presenting ﬁg 45, where the unwrapping
process has been demonstrated. This approach minimizes the resource requirements needed
for such simulations, and simpliﬁes the analysis involved. Hence a 3D combustion channel
can be reduced to an ”unwrapped” 2D version for analysis. Figures 46 and 47 show similar
approaches as employed in RDE simulations in [TSS+ 14] and [KS13].

4.1.1

MESH GENERATIONThe RDE mesh was developed from that used for the validation processes. This adaptation
was undertaken because it is desired to simulate the injection of an RDE, which in ourmodel involves a premixed ﬂuid being injected into a combustion chamber from a plenum.
To achieve this, the domain was taken as it were and modiﬁed the blockMeshDict found in47

Figure 48: 2D Mesh with one obstacle

Figure 49: Sketch of domain with 2 obstacles

the system folder of our built case, by adjusting the vertex points to include a square-shaped
void 0.05 m wide, (referred to as obstacles going forward) placed at a distance a little further
away from the upstream section of our computational mesh. The resulting domain is shownin ﬁg 48. The obstacle is to serve as the solid boundary within the injection cross sectionalarea, and the spaces below and above the void serving as the injection ﬂow part (oriﬁce
in this case). Using this mesh, and by employing similar approach as described above, thevertex points were further adjusted to include two and three obstacles as shown in ﬁgures 49,
50 and 52, such that they are placed above each other, with the ﬂuid ﬂow part in-between
the obstacles as desired. It can be seen that this has the eﬀect of increasing the number ofblocks and faces in the resulting mesh data written into the constant folder in OpenFoam,
when the ”blockMesh” command is run through a terminal. It can be seen that the numberof blocks generated when the ”blockMesh” command is run on OpenFoam increases as thenumber obstacles are increased. This results because the mesh is such that there are strips ofblocks extending from the inlet region to an obstacle, and then from behind it to the end of
the domain (the outlet) for each case. The OpenFoam command ”paraFoam-block” was run
to generate images of the number of mesh points/vertices as well as the faces to help visualize48

Figure 50: 2D mesh generation with two obstacles

Figure 51: Initial pressure for the casethese points. The obstacles were extended from the boundaries for the case with 3 walls,
hence the reason it looks a little diﬀerent. Details of these are found in ﬁgures 52 and 53.
The mesh vertices were generated, for the 3 obstacle cases above, by hand, so that a program

Figure 52: 2D mesh with three obstacles

Figure 53: Image showing mesh blocks(11 in thiscase)

can be written that determines the obstacles width, for a choice of domain conﬁguration,
by taking in as inputs the number of desired obstacles, the void dimensions, and the oriﬁce
diameter(or width). Certain trends were established through this. For instance, we see thatthe total number of points needed to generate a mesh expressed as a function of numberof obstacles (n) is 4 ∗ (2 ∗ n + 2), where n is the desired number of obstacles (solid walls
between oriﬁces). It can also be seen that the total number of blocks needed in the mesh is
given as 3 ∗ (2 ∗ n + 1) − n. A block in OpenFoam is described by eight(8) points. Startingfrom the origin, a numbering system was chosen that respected OpenFoam’s format forface description (e.g outward pointing normal for faces) and generated an 8 column matrix,
which was then populated using a code written in Octave, depending on the desired numberof obstacles. Each face therefore is a resulting combination of four numbers, from the ﬁrst
four elements of the matrix (representing the point where z = 0) or from the last four column
elements in each row representing the z direction (z = 0). Using this, the mesh input ﬁle,blockMeshDict, can be generated to be used with the blockMesh command. Figure 54 shows
the generated meshes for the case with 3 obstacles in octave. The code was improved upon49

Figure 54: Vertices for 2 obstacles using octave code.

Figure 55: Image of mesh with 10 obstacles (codeFigure 56: Image of mesh with 100 Obstacles
generated)(code generated)by making it capable of grouping the faces into Hexagonal blocks, and also generate a text
ﬁle with the same format as the custom blockMeshDict, OpenFoam’s text ﬁle for mesh data,
which we can use to generate the PolyMesh folder, for a simulation case. The text ﬁle should
be copied into the system folder in the case directory. The updated code behaved as expected
and in ﬁgures 55 and 56, images showing cases with 10 and 100 obstacles generated using
it are presented. The obstacle width is given as obstacle thickness (ts ) = (c/n) - injection
width (l), where n is the chosen number of obstacles , c is the unwrapped circumference (c)and l is the injection width. To achieve periodicity of the unwrapped 2D domain geometry,it was speciﬁcally ensured that for each mesh gotten using this code, the injection width at
the top-most and bottom-most parts of the domain sum up to l (one injection width) (i.el = 12 l + 12 l). The sum of all the injection widths (oriﬁces), l, give the total oriﬁce area of
the RDE, in 2D.

50

(a) Steady state pressure proﬁle(b) Image showing shock stuck at inlet(c) Pressure ﬂow ﬁeld(d) Velocity

Figure 57: A resulting ﬂow ﬁeld using engine A.

4.1.2

GEOMETRY SETUP

The simulations in this section are based studies on two engine conﬁgurations, an imaginary
engine with equal width and span, which has been named ”engine A”, and then an actual
conﬁguration from the literature. Using the assumed engine conﬁguration, preliminary simulations were conducted to determine the eﬀect ramping up the injection velocity (and hence
mass ﬂow rate) will have on the ﬂow ﬁeld of an operating RDE, beginning with 100 m/s insteps of 100, up to 700 m/s. These were run on the Compute Canada Cedar and Grahamcluster, using a 10 points per L1/2 reﬁned mesh. The initiation of a detonation wave wasattempted in some cases. Two observations are worthy of note using the assumed engine:(1)It can be seen that there is a pressure rise in the inlet section of the domain resulting froma reﬂected shock travelling back to the inlet, and ”sticking” there, which is completely unphysical. This was the case even after ramping up the velocity, and initializing the domain
with lower pressures set at the discharge end with the hope that this shock gets ”sucked” inat some point.(2) It is important to stabilize the ﬂow ﬁeld before initiating a detonation. A
ﬁxed value inlet with a zero gradient outlet condition for the ﬁeld variables was imposed asboundary condition. Considering that the exact reason for the unphysical behaviour was notknown, it was also thought necessary to examine the boundary condition. Some obtained
results for a case with injection velocity of 700 m/s are presented in Figure 57, where it isseen that the pressure and velocity proﬁles indicate no detonation wave(s) present.

51

4.1.3ACTUAL ENGINE GEOMETRYThis case is based on the works of J. Kindracki and P. Wolanski in [KWG11]. They carriedout an experimental study of the rotating detonation engine with most necessary details
reported. It has an internal diameter of 38 mm. The length of the channel was varied from
23-50 mm; the stagnation pressure and temperature are 700 kPa and 293 K respectively.
Fuel was injected through an oriﬁce of diameter between 0.7-1 mm, while oxygen came inthrough a 0.5-1 mm diameter oriﬁce. For the purpose of our simulations which are based uponpremixed reactants within the combustion chamber, using the same stagnation conditions,
it was assumed that the total injection diameter was 8mm (l), and the span was 50 mm.
Using this l value, and making use of our mesh-generator code, it is possible to generate amesh for n number of obstacles of our choice. The simulations performed were based on 10
injectors given a total injection area of 0.08 × h m2 , and neglecting thickness (for 2-D), h
as in our case, 0.08 m. Each injector area has a throat of 0.008 m. We will refer to thisengine as the Wolanski engine going forward. By assuming a 1-D converging diverging nozzle
conﬁguration problem, we can make useful assumptions that aide engine modeling.

4.1.4

AREA RATIOHaving successfully devised a means to generate n number of obstacles, one may want toknow the ratio between the injection area and the total area of the inlet, including the
obstacles. This is a very useful metric in compressible ﬂow problems. It can be seen thatfor 2D cases, this becomes a length ratio. Of interest usually is the mass ﬂow rate, ṁ forany condition of ﬂow, including situations where the ﬂow is chocked at the oriﬁce. For suchcases, we can deﬁne a throat ,t, and its area At , known as throat area. We can determine
the throat area using 1D isentropic ﬂow. In actual terms, since we are considering 2D cases,
we are most interested in the throat width. Given the width of an oriﬁce or nozzle, andthe total number available, we can determine the total injection width. The area of theunwrappd RDE with witdh (circumference) c, for a given annular thickness δ, is given as
ARDE × π × c × δ. Hence the area ratio becomes,φ=AtlAreaof T hroat==
Areaof RDEAnulus
ARDE
c

(51)

This is important since by considering that mass ﬂow rate, ṁ = ρ × U × A, where A refers
to the cross sectional area of ﬂow, we can compute the injection velocity, U for chocked ﬂow,
which we can specify as input in the initial condition folder, (named ”0” ) for OpenFoam.
We can also determine the area ratio needed to support a given mass ﬂow rate.

52

4.2

BOUNDARY CONDITIONS

Our geometry consists of 5 main regions where we can expect ﬂuid-boundary interactions;
Top, bottom, obstacle, inlet and outlet. These interactions are usually modelled by imposing
boundary conditions which reﬂect the physical behaviour of the region, and hence must makerealistic choices. The mode of operation of our engine is thus; Fluid ﬂows into the chamber
through the inlet boundary at certain conditions, undergoes an expansion process driven bya detonation wave within the chamber, and then exits the chamber via its outlet boundary
end. It is required that these two regions are the only pathways for entrance or exit of ﬂuid
particles through the domain. Of mention again at this point is the fact that OpenFoam
meshes by default are in 3D. Hence one must deﬁne two other boundaries, front and backwhich represent the front and rear faces of our domain. However a boundary condition known
as ”empty”, in OpenFoam lingu, is imposed which enforces a wall boundary condition on thefront and back faces, and also guaranties a 2d simulation. A periodic boundary condition,(known in OpenFoam terms as cyclic boundary condition) is imposed at the top and bottom
boundary, such that whatever ﬂow in the Y direction ”reappears” in the domain from theother end. We want the obstacles to act as solid walls with no ﬂuid ﬂow through, as wellas have the ability to reﬂect oﬀ shocks that bounce at them, hence the slip conditions
are applied on those. In using OpenFoam and the custom solver, rhoCentralFoamreac, forthese simulations, it is required to specify boundary conditions for the following four (4)thermodynamic properties (pressure, P , temperature, T , velocity, U and heat release, Y ).
It is also assumed that the engine is constantly fed with combustible reactants, hence thereason for imposing the ﬁxed value BC . The outﬂow BC is a linear extrapolation with zerodζ
gradient dx= 0, where ζ = {P, T, U, Y } These boundary conditions were tested with asimple case involving a shock travelling into the domain which encounters an obstacle (orobstacles) within it. Figures 58 and 59 show the pressure proﬁle of a reﬂected shock oﬀ
for test cases having one (1) and two (2) obstacles, with slip boundary condition imposed
which is expected. However this shock travels back to the inlet and ”sticks” to it, leading
to a constant pressure rise in the region. This result, together with that from the assumed
engine conﬁguration above, implied that it became necessary to explore other conditions,particularly for the inlet boundary. Some of the cases considered are described below.

4.2.1

FIXED VALUE BOUNDARY CONDITION

This implements the Dirichlet condition where face values are deﬁned according to a speciﬁedreference value: φf = φref . The syntax is deﬁned in [Ope16] and was widely used for the53

(a) t = 0.00017 s(b) t = 0.00029 s(c) t = 0.00043 s(d) t = 0.00046 s

Figure 58: Pressure ﬂow ﬁelds in the ”slip” BC test with 1 obstacle(a) t = 0.00006 s(b) t = 0.00018 s(c) t = 0.00033 s(d) t = 0.00038 s

Figure 59: Pressure ﬂow ﬁelds in the ”slip” BC test with 2 obstacle

validation process. It is assumed that the mixture is injected into the engine at a certain(ﬁxed) thermodynamic state, hence the reason the ﬁxed value BC was imposed.
RESULT
Using this BC results in a unphysical behaviour in which a reﬂected shock from the obstacles
travels to the inlet patch and remains there as shown in ﬁgures 57c and 61b. This raises thepressure, and aﬀects the temperature and velocity at this boundary patch during simulation.

4.2.2

INLET-OUTLET BOUNDARY CONDITION

This is usually applied in situations where one wishes to account for back ﬂows. It applies azero gradient condition to the forward moving ﬂow and then applies a speciﬁed (ﬁxed value)
constant to the backward ﬂow. Figure 60 gives a representation of this. This was applied at
the inlet for the pressure ﬁeld, and for other ﬁelds as well.
RESULT
The inlet-outlet performed more like the ﬁxed value boundary condition, failing to provide
any outﬂow in the inlet boundary. Thus the reﬂected shock oﬀ the obstacle travels backtowards the inlet and remains there.

54

Figure 60: InletOutlet boundary condition [Cfd]4.2.3

TOTAL PRESSURE/TEMPERATURE BOUNDARY CONDITION

The total pressure (total temperature) boundary condition is an extension of the ﬁxed value
boundary condition in OpenFoam, used to specify the static pressure at the boundary basedon the user deﬁned total or stagnation pressure and total temperature as well. It is usefulfor modelling plenum conditions. Variants for compressible and incompressible ﬂows exist
in OpenFoam. For the cases involving supersonic ﬂow, for instance, it evaluates the staticTo
conditions by using the isentropic relations. Thus, Tp = 1+ γ−1, Where Tp is deﬁned as the|u|2γ

temperature of the patch (static), and γ is the adiabatic index. Figures 62 show the results
obtained for this case.
RESULTS
The total pressure condition performed well for our tests. A back pressure value was speciﬁed
for the domain based on isentropic ﬂow calculations. Also considered was a case with zero
gradient boundary condition at the outlet. In all cases, we see that we do not obtain theunphysical solution as before, in the inlet boundary, and depending on the back pressure
55

value, we are able to obtain a steady solution. Figure 62 presents plots for the pressure,temperature and velocity.

4.2.4

PLENUM PRESSURE BOUNDARY CONDITION

OpenFOAM uses this boundary conditions to set a plenum inlet condition. It models a zero
dimensional volume of gas enclosed ahead of the inlet patch and the static pressure exerted
on the inlet depends on the thermodynamic state within the enclosed volume. This model
neglects momentum and advances the density and temperature of the plenum with run-time.
The smaller the volume of the plenum, the more responsive it is to changes in mass ﬂow
rate and the further the model approximates a ﬁxed mass ﬂow system. If we signiﬁcantly
increase the volume, we obtain a ﬁxed mass ﬂow rate situation.
RESULTS
The plenum boundary condition also eliminates the jump in pressure observed at the inlet.
As implied from our deﬁnition earlier, when it is assumed that a large enough volume isinvolved, this condition approximates to a ﬁxed total pressure condition. Since a plenum
that signiﬁcantly maintains the stagnation conditions is desired, and in order not to introduce
extra dynamics into our system, the total pressure and temperature boundary conditions was
utilized for our simulations. Figure 63 represents ﬁeld parameters obtained using the plenum
condition.

4.3

CHOICE OF FLOW CONDITIONSUsing the total pressure and temperature boundary conditions, four diﬀerent cases were set
up, reﬂecting four possible solutions for an isentropic converging-diverging nozzle problem,by varying the back pressure in each case to determine the ﬂow ﬁeld. Using the stagnation
conditions, the limiting back pressure values were determined for the subsonic (619 kPa) andsupersonic (124 kPa) ﬂow solutions with shock at the exit. Two other cases (100 kPa and650 kPa) were chosen to enable us to see the results for other scenarios. The resulting ﬂow
ﬁelds are discussed in more detail.

56

(a) Initial pressure state(b) pressure plotted against distance, x.

(c) Velocity using ﬁxed boundary condition(d) Temperature plotted against distance, x

Figure 61: Plots showing ﬁxed boundary condition with ﬁeld variables indicating unphysical behaviour in the boundary57

(a) Pressure plot showing shock and expansion
wave.

(b) Pressure after 25 time steps(c) Velocity after 25 times steps(d) Temperature

Figure 62: Plots obtained using total pressure and temperature boundary condition

58

(a) Initial pressure state(b) Pressure proﬁle using plenum boundary condition(c) Velocity after 25 time outputs(d) Temperature after 25 time outputs

Figure 63: Plots showing plenum boundary condition with ﬁeld variables

594.3.1

SUPERSONIC FLOW FIELDS: BACK PRESSURE OF 124kPa
and 100kPa

These back pressures result in a supersonic ﬂow at the nozzle exit, as shown in ﬁgure 62.
However we do not get to a stable solution partly because the solver cannot handle thevortices that develop as the simulation is advanced in time. The case with back pressure of
124 kPa is the boundary case for a shock at the exit of the diverging end, while the casewith a back pressure of 100 kPa should in reality result in oblique shocks at the exit. We
encounter unphysical inﬂows into the domain from the outlet after the simulation runs fora long time, which eventually makes it fail. Hence we were not able to achieve a steadysolution of the supersonic cases. We highlight that in view of these, we set up another casewith an expected supersonic solution having a back pressure of 250 kPa. This was able torun for an extended period of time, without failing, until we ended the simulation. Despite
reaching a metastable condition, disturbances eventually grew from the back boundary andpolluted the solution.

4.3.2

SUBSONIC FLOW FIELDS: CASES WITH BACK PRESSURE OF 619 kPa and 650 kPa

These result in subsonic ﬂow ﬁelds, with 619 kPa being the limit based on 1-D isentropic
ﬂow calculations. Steady solutions of these cases were obtained by allowing the simulations
to run for a long time . The corresponding pressure and velocity proﬁles for the case withback pressure of 650 kPa is presented in ﬁgure 64.

4.3.3

MACH NUMBER

Here the Mach number for the selected ﬂow ﬁelds for both the subsonic and supersonic cases
are presented as shown in ﬁgures 65a and 65b.

4.4

INVESTIGATION OF MASS FLOW RATE

From continuity, the mass ﬂow rate at the outlet should equal that at the inlet, for a steady
ﬂow. Some anomalies were noticed when using the post processing utility of determining
mass ﬂow rate inbuilt in OpenFOAM, at the preliminary stage of simulation, where diﬀerent
values for the inlet and outlet patches were obtained. We want to be able to examine theimpact variations in supply mass ﬂow rate could have on the operations of an RDE, hence
the reason this is of importance. This lead us to investigate the accuracy of the method of60

(a) Pressure proﬁle vs x axis(b) Temperature proﬁle vs x axis(c) Velocity ﬁeld(d) Velocity plot vs x

Figure 64: Steady state subsonic ﬂow ﬁeld using a back pressure of 650 kPa

61

(a) Mach number for subsonic initiation ﬂow ﬁeld(b) Mach number for supersonic initiation ﬂow ﬁeld

Figure 65: Mach number for subsonic and supersonic initiation ﬂow ﬁelds

calculation implemented in OpenFoam. To do this, the value at the inlet and outlet of ourdomain were evaluated for 6 back pressure test cases: 100 kPa, 124 kPa, 250 kPa, 650 kPa,
500 kPa and 619 kPa. Plots of plots mass ﬂow rate (ṁ) against time for the 3 cases: 250 kPa,
650 kPa and 500 kPa are presented. The ”integrate patch” utility was utilized to determine
the surface integral over the patch faces in question (inlet and outlet) as a function of time.
A comparison of these values with that obtained from 1-D isentropic ﬂow analysis. The
mass ﬂow rate into a control volume is given by the product of the density, velocity and thesurface normal to it. In OpenFoam, a negative surface normal indicates an inlet area.

4.4.1

RESULTS AND ANALYSIS

The mass ﬂow rate results obtained for the 4 cases in ﬁgures 66 to 69 are shown. The plotsin blue indicate the inlet mass ﬂow rate, while that for the outlet boundary are in red.
The supersonic cases do not achieve a steady solution and does not result in a convergenceof mass ﬂow rate values at both inlet and outlet patches even after a long time. They startoﬀ with a large amplitude of oscillation which decreases with time, until the simulation
eventually fails, due to the development of intense vortices which cannot be handled by theoutlet boundary condition, as shown in ﬁgures 66a. However it can be seen that at sometime steps before the solution goes unstable, the same values are obtained approximately,
implying a reasonable degree of convergence at that point. From the outlet ﬂow proﬁle, it62

(a) Mass ﬂow rate against time at inlet boundary for Pb =100 kPa.

(b) Mass ﬂow rate against time at outlet boundary for Pb =100 kPa

Figure 66: Mass ﬂow rate across boundaries for Pb = 100 kPa.

(a) Mass ﬂow rate against time at inlet boundary for Pb =250 kPa.

(b) Mass ﬂow rate against time at outlet boundary for Pb =250 kPa(c) Inlet and outlet ﬂow ﬁelds plotted together
showing some oscillations

Figure 67: Mass ﬂow rates across boundaries for Pb = 250 kPa.

63

(a) Mass ﬂow rate against time at inlet boundary for Pb =500 kPa.

(b) Mass ﬂow rate against time at outlet boundary for Pb =500 kPa(c) Inlet and outlet ﬂow ﬁelds plotted together
showing ﬂow is steady

Figure 68: Mass ﬂow rate across at boundaries for Pb = 500 kPa

can be seen that the mass ﬂow rate dips to a negative value which suggests an inﬂow into
the boundary at some point, after which the simulation fails. The case with a 250 kPa backpressure as shown in ﬁgure 67 gives a more stable result. The velocity goes supersonic at
the nozzle exit, but even then, the ﬂow swirls after some time, and a ﬁnal steady solution
is not achieved. The calculated mass ﬂow rates for the supersonic ﬂow with a back pressureof 100 and 250 kPa is 1.225800 × 10−1 kg/s and 1.124639 × 10−1 kg/s respectively, while forthe subsonic cases with a back pressure of 619 and 650 kPa, the resulting mass ﬂow rate is3.397538 × 10−2 kg/s and 2.128058 × 10−2 kg/s respectively. The cases with back pressuresof 250 kPa, 500 kPa and 650 kPa are shown in ﬁgures 67, 68 and 69. The subsonic cases
resulted in a convergence of inlet and outlet values to a constant value, after running for along time (over 3700 times steps).

64(a) Mass ﬂow rate against time at inlet boundary for Pb =650 kPa.

(b) Mass ﬂow rate against time at outlet boundary for Pb =650 kPa(c) Inlet and outlet ﬂow ﬁelds plotted together
showing ﬂow is steady

Figure 69: Mass ﬂow rate against time across boundaries for Pb = 650 kPa

65

4.5

SELECTED FLOW FIELDS FOR INITIATION

4.5.1

FLOW FIELD FOR SUBSONIC INITIALIZATIONThe case with the back pressure set at 619 kPa was selected, and its steady state solution
using a coarse mesh obtained. The CFL of 0.1 was maintained, and the calculated mass ﬂow
rate outputted at every time step. This served as our stability and convergence criteria. The
corresponding pressure proﬁle is shown in ﬁgure 70a.

4.5.2

FLOW FIELD FOR SUPERSONIC INITIALIZATIONThe selected ﬂow ﬁeld to be tested for supersonic initiation is that obtained using the casewith a back pressure of 124 kPa. The same process as in the subsonic case above wasrepeated. The ﬂow simulation fails eventually due to vortices and oscillations in the exit,
and also because our boundary condition does not account for the solution of a ConvergingDiverging nozzle with an expansion or oblique shock at the exit. The proﬁle obtained issimilar to the one shown in ﬁgure 62. To conclude this section, at this point two ﬂow ﬁelds,
one subsonic and one supersonic have been chosen. These become the initial condition for
the initiation of a rotating detonation wave, discussed in the proceeding chapter, where adescription of some of the various initiation techniques attempted are presented.

66

Chapter 5
SIMULATION OF ROTATINGDETONATION ENGINEINITIATION
5.1

PROBLEM DESCRIPTION/CASE SET UP

It has been established that the ﬂow ﬁeld for the Wolanski engine, with a certain numberof injectors, can be modelled as a converging-diverging nozzle problem from gas dynamics,
using the equations of 1-D isentropic nozzle. Solutions for subsonic and supersonic cold start
conditions for this engine depending on the choice of back pressure have been obtained. In
this chapter, it is desired to initiate the RDE, and evaluate diﬀerent initiation patterns,
considering a case with high temperature, another with high pressure and temperature, andone with burnt mixture within the domain.

5.1.1

MESH DESCRIPTION

The same domain as with the cold start has been retained; i.e., the unwrapped 2-D domain.
However the mesh was reﬁned to have 0.7 points per L1/2 , resulting in 86850 cells, and thecase with 10 points per half reaction zone length, resulted in 17 million cells (17042520). Thiswas run on Compute Canada high performance clusters Graham and Cedar. OpenFOAM
allows for parallel decomposition of mesh for computing, hence the domain was decomposed
into 32 units for the same number of processors, each having over 500,000 computational
cells. It was split along the x-axis.

675.1.2BOUNDARY CONDITION

The same boundary conditions as with the cold start case was used to generate the mesh,i.e., a total pressure and temperature condition with a back pressure value depending on if
the initiating ﬂow is to be subsonic or supersonic. First the mesh was generated with theseand for the actual initiation ﬂow simulation, after which the outlet boundary condition waschanged from the back pressure condition (ﬁxed value pressure) to a boundary extrapolated
ﬂow condition (zero gradient). This was done to prevent expansion waves at the divergent
section, which degrade the strength of the initiation zone, thereby preventing a transition toa detonation wave, particularly for the supersonic case.

5.1.3

INITIAL CONDITIONS

For the purpose of simulating the RDE, two initiation ﬂow ﬁelds were considered; subsonicand supersonic cases. The selected ﬂow ﬁeld results in chapter 4 were utilized as our initial
ﬂow ﬁelds. For both cases, the results were mapped from a steady time step obtained fromthe coarse mesh cold ﬂow simulations, using the ”mapFields” utility, which had 16200 cells,corresponding to 0.2 points per half reaction zone length, to the ﬁner meshes considered. It
should be noted that both meshes are ”consistent”, i.e., of similar geometry and boundaries,
hence there was no need to create a ”mapFieldsDict” in the system folder for these cases.

5.1.4

DIFFERENT TYPES OF INITIATION

Four diﬀerent initiation techniques were considered as shown in ﬁg 71; a purely temperature
based initiation in which case the temperature within a small box or zone is set to an elevatedvalue as in ﬁg 71a, and an initiation pattern with both elevated temperature and pressure
values as in ﬁg 71b, another with burnt gases in the combustion chamber as in ﬁq 71d, andﬁnally a case where the initiation was done further downstream of the domain as in ﬁg 71c.
These were tried for both the supersonic and subsonic ﬂow conditions, for both the high andlow resolution cases that were simulated.

5.2

RESULTS WITH LOW RESOLUTION MESH

Results for the low resolution mesh with 0.7 points per L1/2 simulated on a local device arepresented and discussed below.

68(a) Pressure proﬁle(b) Temperature proﬁle(c) Velocity ﬂow ﬁeld(d) velocity proﬁleFigure 70: Subsonic ﬂow ﬁeld using a back pressure of 619 kPa

69

(a) Temperature based initiation(b) Pressure-Temperature based initiation(c) Pressure-Temperature based initiation far away from thedomain(d) Pressure-Temperature initiation with burnt gases

Figure 71: Various initiation techniques considered

5.2.1

SUBSONIC FLOW FIELD INITIATION

The ”setFields” utility was used to set a high temperature region of 2000 K within a few
cells contained in a small box, by using the ”boxTocell” utility, after which the simulation
was advanced in time. The hot spot(small box) was set to be 0.005 m away from theinjectors and was of same magnitude in length. It makes sense to do this before decomposing
the domain using the ”decomposePar” command, otherwise one will have to map the cells
from the individual processors to those contained in the ﬁner mesh. Detonation waves that
propagated in both the axial and circumferential directions in the combustion chamber were
initiated. These waves are symmetric. The von Neumann spike pressure, as outputed fromthe simulation after 9 time outputs is 2.69168 × 107 kPa , while the pressure behind thedetonation wave is about 8 × 106 kPa, roughly 10 times the initial value, which correlateswith what is obtainable in literature. The shock pressure is 2.77533 × 107 kPa just before
the simulation ends. We note that the simulation fails as the exists a detonation wave
travelling back axially to the inlet boundary, through the injector holes, at a very highpressure. The pressure, temperature and velocity plots are shown in ﬁgure 72. Results forpressure-temperature based initiations are presented in ﬁgure 73 for a hot spot introduced

70

(a) Initiation spot (T)(b) Pressure(c) Temperature(d) Velocity(e) Pressure ﬂow ﬁeld in x direction(f) Temperature proﬁle in in x direction(g) Velocity temperature proﬁle in x direction(h) Heat release plot in x direction

Figure 72: Detonation initiated with a high temperature using a subsonic ﬂow ﬁeld

71

in a region close to the nozzle outlet. The same temperature was maintained as with thetemperature based initiation, of 2000 K. The pressure in the region was however set to
900 kPa, about 1.3 (1.2857) times higher than the stagnation pressure. It can be seen thatthe initiation of a rotating detonation wave was achieved, which had the same von Neumannpeak pressure (2.69168 × 107 kPa) value after nine (9) time outputs, as with the temperature
based initiation, with the pressure behind the detonation also remaining fairly the samevalue. This suggests that the temperature seems to be the main driver of the shock strength
for this choice of pressure ratio.
A temperature based initiation with the hot spot placed further away from the domain
was performed. Plots in the y direction showing two detonation waves travelling with thesame peak pressures are presented. Figure 74 presents results for the ﬁeld variables as thedetonations travels in space. Of signiﬁcance is the observed shock-reaction coupling eﬀect
which is the hallmark of a detonation. Initiation cases in which the reactive mixture fromcertain regions of the domain are cut oﬀ were considered to see if that results to obtain
one initial detonation wave travelling in the y-direction. In practice, this can be viewed as
igniting a detonation chamber ﬁlled with a mixture of burnt and un-burnt gases. Figure 75
presents results obtained, and we see from the image for the pressure shown in ﬁgure 75b,that this approach results in one detonation wave travelling faster than the other.

5.2.2

INITIATION WITH SUPERSONIC FLOW FIELD

Using the case with a back pressure of 124 kPa, we initialized the supersonic cases to run
with a CFL of 0.1, and with Δt set as 1×106 s. The supersonic cases never achieved a steady
ﬂow solution as stated earlier, however, by making a comparison between the inlet and outlet
mass ﬂow rates, and also visualizing the ﬂow ﬁeld using Paraview, it was possible to pick out
a time step that appeared stable to some extent. Hence the simulation results were mappedfrom that time step into our ﬁner mesh. The process of initiation using the temperature,and pressure-temperature based approaches were consided. The setFields utility, was used
to set up hot spots as described earlier. The same temperature of 2000 K was maintained
for this case. Figure 76 presents results obtained for a temperature based initiation for thesupersonic ﬂow ﬁeld. We see that a symmetric detonation wave propagating in the x and y
directions is achieved, and one also recovers the oscillating behaviour of a detonation wave.
It can be seen also that the shock pressure for this case was of the same order of magnitudewith the activation energy, and the propagating detonation is oscillatory in nature. The
peak shock pressure after 10 seconds is 5.7013 × 106 kPa, with the highest particle velocity
being 1305.01 m/s Results for pressure-temperature based initiations for this ﬂow ﬁeld are72

(a) Initiation spot (P)(b) Pressure proﬁle(c) Temperature proﬁle(d) Velocity(e) Pressure plotted over a line in y direc- (f) Temperature plotted over a line in y (g) Velocity plotted over a line in y direction
directiontion(h) Heat release plotted over a line in ydirection

Figure 73: Detonation initiated with a high pressure-temperature zone using a subsonic ﬂow ﬁeld

73

(a) Initiation spot (T)(b) Two travelling detonation waves(c) Temperature(d) Velocity(e) Heat release(f) Shock pressure- reaction couplingin y direction(g) Detonation Temperature in y direction(h) Velocity in y direction

Figure 74: Detonation initiated further downstream with a high temperature zone using a subsonic
ﬂow ﬁeld74

(a) Initiation spot (T)(b) Velocity(c) Temperature(d) Heat release(f) Velocity in x direction(e) Pressure in x direction(g) Shock-reaction coupling in x direction

Figure 75: Detonation initiated with a high pressure and temperature zone using a subsonic Flow
ﬁeld

75

(a) Initiation spot (T)(b) Temperature(c) Velocity(d) Reaction(e) Pressure proﬁle in y direction ﬂow ﬁeld(f) Temperature in y direction(g) Velocity proﬁle in Y direction(h) Shock pressure-reaction coupling in y direction

Figure 76: Detonation initiated with a high temperature zone using a supersonic ﬂow ﬁeld

76

shown in Figure 77. The same ratio as with the subsonic cases was utilized. The peakpressure value after 10 time outputs was 5.7013 × 106 kPa, slightly lower than that for thetemperature based initiation. We see that two symmetric detonation waves propagating in
the y direction is the resulting solution. Another mode of initiation for the supersonic ﬂow
mixture was considered, in this case, a certain section of the domain is ﬁlled with burnt
gases. This was done by setting a zero value to the activation energy at the exit end of thedomain. The pressure and temperature values within the hot spot was set to the same value
as with those for the pressure-temperature based initiation, by using the setFields utility.
Results for this initiation procedure are presented in ﬁgure 78

5.3

EFFECT OF HIGHER RESOLUTION ON DETONATION INITIATION

Using both ﬂow ﬁelds (subsonic and supersonic), simulations were performed on a highresolution mesh, with 10 points per L1/2 , in Compute Canada, resulting in over 17 million
cells. The domain was decomposed into 32 strips, using the ”decomposePar” utility in
OpenFoam, with the mesh split in x-direction, each having over 500,000 cells. Again attempt
was made at initiating a detonation wave using temperature and pressure-temperature based
initiation techniques as with the low resolution cases. The results obtained for each mode of
initiation for both the subsonic and supersonic ﬂow ﬁelds considered are discussed below.

5.3.1

SUBSONIC CASE

First the subsonic high resolution mesh was initiated with a hot spot of temperature of2000 K. However, a detonation wave was not achieved. It is important to restate here thatthe hallmark of a detonation wave is a coupling between the shock and reaction, thus one
should expect to see this coupling when a visualization of simulation results is performed in
the post processing software used for this work, Paraview. It is worthy of mention that eventhough one gets an initial jump to a high pressure value, and the reactants are burnt out after
ignition, the ﬂow does not transition into a detonation wave. Another initiation approachfor the RDE, done by setting a high pressure-temperature region within the combustion
chamber just after the nozzle as in ﬁgure 71b, using the setFields utility as above. Four(4) pressure ratios were considered: 1.3 (as above), 1.5, 2 and 5, each case maintaining thesame temperature. A sustained detonation wave for the subsonic ﬂow ﬁeld for all the ratios
considered was not achieved. There temporarily existed a detonation wave propagating
for a few time steps, however this quickly decouples. It is important to highlight that77

even though combustion continues, the engine doesn’t function as one that is detonation
based. Eventually, the simulation fails after some time as there is an axially back-travellingdetonation wave, with a shock pressure value higher than the stagnation pressure, and hence
the inlet boundary values, which reaches the inlet, thereby preventing further injection.
Results for the subsonic case with a pressure ratio of 5 are presented, with a look at ﬁgures 79a
and 79c showing the decoupling between the shock and chemical reaction behind it, indicatingthat a sustained propagating detonation wave is not achieved. A presentation of 1-D slice
images as shown in ﬁgures 81 and 82 highlight this further, where by a coupling of thepressure and reactive terms exist initially, decoupling after some time.

5.3.2

SUPERSONIC CASE

For this case, the high resolution mesh was initially initiated with a temperature of 2000 K.
A detonation wave for the supersonic case with this technique was also not initiated. It isa struggle to see even a signiﬁcant jump in pressure after ignition, as is usually the casewhen a detonation is formed, and there exists a coupling between the shock pressure andreaction. One would also expect to see this coupling when a visualization of the solution isdone in Paraview. It is seen that even though this results to a jump in pressure value, andthe reactants are burnt out after ignition, it does not transition into a detonation wave. The
RDE was then initiated by setting a high pressure-temperature region within the combustion
chamber using the setFields utility as above. Four (4) pressure ratios were also considered:
1.2857(as above), 1.5, 2 and 5, each time maintaining the same temperature of 2000 K. It
was observed that this fails to achieve a detonation (or sustained) wave for both the subsonicand supersonic ﬂow ﬁelds for all the ratios except for the supersonic case with a pressure
ratio of 5. This case results in a sustained detonation wave circumferentially propagatingin our domain. However the simulation fails as well after some time as we have an axial
back-travelling detonation wave with a pressure higher than the inlet boundary pressure,which hits the inlet thereby preventing further injection. Results for the supersonic ﬂow
initiations using a pressure ratio 5 are presented in ﬁgures 83 and 84, where it can be seenthat the propagating pressure wave remains coupled with the chemical reaction occurring
behind it. A rotating detonation wave is achieved and sustained for the supersonic caseusing this pressure ratio, with a Von-Neumann peak pressure of 1.1 × 107 kPa after 14 time
outputs.

78

5.3.3

CONCLUSION

In conclusion, circumferentially propagating detonation waves have been initiated by using
diﬀerent initiation techniques; a high temperature spot of 2000 K, and a high pressuretemperature initiation zone as well. For the low resolution mesh, the eﬀect of utilizing ahigh pressure-temperature initiation instead of only a high temperature zone didn’t seem tobe obvious, but became so when the high resolution (10 points per L1/2 ) mesh was utilized,
and for the same ignition pressure ratio (of 5) and temperature, a transition to detonation
waves was achieved using a supersonic ﬂow ﬁeld as the initial ﬂow ﬁeld.

79

(a) Pressure(b) Temperature(c) Velocity(d) Heat release(e) Pressure proﬁle in y direction(f) Temperature proﬁle in y direction(g) Velocity proﬁle in y direction(h) Shock-Reaction coupling in y direction

Figure 77: Detonation initiated with a high temperature and pressure zone using a supersonic ﬂow
ﬁeld

80

(a) Pressure(b) Velocity(c) Temperature(d) Heat release

Figure 78: Temperature based initiation for supersonic ﬂow with a pressure ratio of 1.2857(a) Pressure(b) Velocity(c) Temperature(d) Heat release

Figure 79: Pressure-Temperature based initiation for subsonic ﬂow with a pressure ratio of 581

(a) Pressure(b) Velocity(c) Temperature(d) Heat release

Figure 80: Pressure-Temperature based initiation for supersonic ﬂow with a pressure ratio of 5

Figure 81: Initial unstable detonation wave (Subsonic case)82

Figure 82: Image indicates decoupling of shock pressure-reaction(Unsustained Detonation for subsonic case)

Figure 83: Pressure-reaction coupling indicating a detonation wave (supersonic case)83

Figure 84: Coupling remains as wave travels further (supersonic case)84

Chapter 6
CONCLUSION
A validation study has been conducted for the KT/KNP scheme. We were able to showthat the minmod limiter produces the least disturbance amplitude for the various limitersconsidered, with superbee performing the worst. Critical CFL before a signiﬁcant jumpin disturbance amplitude is observed increases with increasing pressure ratio. For pressure
ratios up to 5, a maximum CFL number of 0.15 gives a 1% degree of accuracy using theminmod limiter, with lower ratios allowing for higher accuracy. CFL ≤ 0.1 agrees with rulesof thumb disseminated throughout the community. For the simple, over-driven detonation
problem, KNP scheme recovers the oscillating solution of a detonation wave. Findings from asigniﬁcant number of linear suitability analysis show however that the resulting ZND solution
gives a single mode instability [NNL00] for the chosen choice of parameters for a pulsating
detonation wave, with an overdrive factor of 1.6. However, this mode of oscillation is notseen in this case, even after increasing the resolution up to 100 cells per L1/2 . Instead, we
obtain a proﬁle that is rather irregular in pattern. We note that the peak pressure values
obtained for each case was above 100 kPa (about 108 kPa), quite a high value consideringthat most of the schemes discussed in [NNL00] show a peak pressure value below 100 kPa,
for the various cells considered.
Cold start conditions for initiation using the Wolanski engine have been determined, andthe initiation of a propagating detonation wave demonstrated. It has been established thatregardless of the initial ﬂow ﬁeld-be it supersonic or subsonic-,a rotating detonation wave canbe achieve. For the same stagnation pressures however, one would require a higher pressure
ratio to initiate a detonation using a subsonic ﬂow ﬁeld, as opposed to a corresponding
supersonic ﬂow ﬁeld. The signiﬁcance of a highly resolved mesh in detonation simulations
has been demonstrated, as this is crucial to determining physically acceptable results.
It was thus determined that under certain conditions, the KT/KNP can accurately solve thestandard shock tube problem. What has so far been an anecdotal rule-of-thumb has been85

formalized. It was also determined that a reactive solver based on the KT/KNP scheme
is unable to accurately capture the 1-D dynamics of detonation waves. Detonations were
initiated in a premixed RDE setup, but always resulted in a detonation travelling backinto the plenum section. A conditional reactive solver (for e.g one that suppresses reaction
for spatial positions upstream of the injection holes) could be implemented to mimic nonpremixed RDE operation.

6.1

RECOMMENDATION

Higher resolution cases (with more points per L1/2 ) need to be tested for the initiation
process, to see how the scheme performs under such cases, as well as to compare peakpressure and velocity values. Also, the inlet boundary condition needs to be modiﬁed toenable it adapt to on-coming ﬂows to it, since a detonation wave will always cause a rise inpressure, which in turn will seek to drive the ﬂow upstream of it backwards. This will entail
creating a custom boundary condition, within the OpenFoam framework, to that eﬀect.
Having been constrained by the validation results to a CFL number of 0.1, it becomes
necessary to implement a numerical scheme that allows one run at even higher numbers
than that, without spurious oscillations. Overall, we do not recommend using the KT/KNP
scheme for detonation simulation.

86

AppendicesTable 1: Table of CFL and amplitude of pressure oscillations for the 4 chosen limiters for a pressure
ratio of 10.

(*)represents CFL numbers that guarantee 1% error

CFL
0.05
0.06
0.07
0.08
0.09
0.1
0.11
0.12
0.13
0.14
0.15
0.16
0.17
0.18
0.19
0.2
0.21
0.22
0.23
0.24
0.25
0.26

SuperBee(dp)77237
267.2(*)
54742.06
36153.5
12538069167.2
84169.2
Ftr
69481.2
91258.35
72417.5
86432.71
1163.290264.49
94596.69
83276.64
26541.2
106738110998.5
93775.585220.3
62853.4

VanAlbada(dp)
115.9
180.1
230.8
Ftr
306.8(*)
407.8
469.4
512.3
561.5
703.8
612.3
716.1
826.7
842.4
19256.8
1000.21163.51555.115613
15831.4
17712.5
20305

MinMod (dp)
83
104.3
127.8
152.9
184.3
204
20716.5
297.3
317.6(*)
364.9
419.4
484.5
826.7
581.1
615.3
710.8
767.5
785.5
883.4
915.9
1047.21121.7

Vanleer(dP)
73
81.7
117.9
12219.8
229
228.6
297.7(*)
406
1666.112845.6
14749.8
17085.9
17223.3
18106.3
19600.72057.325534.1
22977.927202.3
25785.4
30827.430570.2

Table 2: Table of CFL and amplitude of pressure oscillations for the 4 chosen limiters for pressure
ratio of 2.

(*)represents CFL numbers that guarantee 0.1% error

CFL
0.05
0.06
0.07
0.08
0.09
0.10
0.11
0.12
0.13
0.14
0.15
0.16
0.17
0.18
0.19
0.20
0.21
0.22
0.23
0.24
0.25
0.26

VanAlbada (dp)
0.3
0.4
0.4
703
0.2
0.2
0.9
0.1
6.5(*)
34.4
106.1
253.8
1278.84971.931784.9
19509.97001.48341.17961.28774.58739.99045.1SuperBee (dp)34643.01
30547.937165.3535529.2143249.41
31232.47
10756.2
19732.37
14130.24
14826.14
16620.48
34643.01
17212.33
15320.68
14339.22
15363.5
14913.28
16516.96
17131.22
16803.63
16969.06
17732.7

Vanleer (dp)
4.4
4.3
5.7
1419.97.3
7.5(*)
25.3
133.1
5408.86794.16322.37206.58632.19014.19914.710941.7
10251.1
12105.3
12138.38
12444.513510.27
13384.94

MinMod (dP)
0.5
0.5
0.6
0.7
0.8
0.8
294.5
1.4
4.3
7
6.8
8.5
8.8(*)
30.8
18.4
36.7
54.2
192.5
3491.73666.65429.95741.3

Table 1 presents the various Finite Volume Methods (FVM) test cases evaluated, their CFL numbers, and pressure disturbance amplitude values(dp)for each, with a number of cases failing to run (Ftr). We also note that although all cases were set to end at the same time, (0.007),some didn’t,
particularly with the superbee limiter. We mark with ∗, the critical CFL number for the 4 limiters.

Table 3: Table of CFL and amplitude of pressure oscillations for the 4 chosen limiters using pressure
ratio of 5(*)represents CFL numbers that guarantee 1% error

CFL
0.05
0.06
0.07
0.08
0.09
0.1
0.11
0.12
0.13
0.14
0.15
0.16
0.17
0.18
0.19
0.2
0.21
0.23
0.24
0.25
0.26

Vanleer(dp)
12.8
14.7
17.8
1250
28.9
32.1
40.3
69.8
114.5
213.3(*)
384.6
4404.510573.611310.1
13127.914894.3
15727.9
15950.7
18912.922197.921926.6

Minmod(dp)
21.5
25.6
29.5
21.5
21.5
21.5
44.6
1609.8
64.3
69.3
77.9
86.2
96.8
107.3
124.2
134.8
141.5
225.6
194.6
222.7(*)
263.1

VanAlbada (dp)
12.6
20.7
30.5
2212.4
52.5
69.9
78.9
103
104.7
142.4
155.5
162.5
208.1
195.8(*)
2776.93472.4
325.2
392.4
459.8
732.7
732.7

SuperBee(dP)
59117
67428
52553
40368.5
57997
61737
64700
23176
63612.766364.8
62244
72266
66927
76708
45177.896736.7
51001
53346
62656
81719
73943

Table 4: Table of CFL and amplitude of pressure oscillations for the 4 chosen limiters using pressure
ratio of 1.5(*)represents CFL numbers that guarantee 0.1% error.

CFL
0.05
0.06
0.07
0.08
0.09
0.1
0.11
0.12
0.13
0.14
0.15
0.16
0.17
0.18
0.19
0.2
0.21
0.22
0.23
0.24
0.25
0.26

Vanleer(dp)
0.3
0.5
0.7
1423.5
1(*)
45.5
877.7
2545.54209.55618.816252.366537.87295.028195.428274.458972.979958.369825.24
11021.36
11436.22
12286.66
12167.49

VanAlbada (dp)
1.2
0.9
0.8
424.3
1.6
0.7
0.8
1.2
10.3(*)
38.8
294.7
1641.24411.25256.17
18755.5427324.18
6914.776892.927387.477988.18155.868291.08

Minmod (dp)
5.4
4.7
3.4
2.9
2.2
2.4
223
2.3
2.3
2.8
3.7
18
11.7(*)
41.2
108.1
117
3386.44181.45482.855624.996345.65
6166.28

SuperBee(dP)
29367.75
21181.44
24962.17
18244.87
26378.03
28265.94
29249.06
18437.71
10251.05
10174.27
10212.92
10523.67
10766.38
10482.85
10448.48
11591.58
11888.02
12837.94
12001.94
12581.73
12890.19
14205.06

Listing 1: Codes utilized for ZND solution
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

gamma = 1 . 2 ;P 0 = 1 0 0 ; #kPa
T 0 = 298;
#K
#D = 2 0 0 0
f = 1 . 6 ; #o v e r d r i v e = Dˆ 2 / ( D c j ) ˆ2
#k = 0 . 6 9 3 1 5
k = ( 1 / 2 . 5 e −9)
R = 8 . 3 1 4 ; #KJ/KMolK
M = 2 8 . 9 6 ; #Kg/KMol
Rsp = R / M; # k J /KgKr h o 0 = P 0 / ( Rsp ∗ T 0 ) #Kg/mˆ 3 ;c 0 = s q r t (gamma ∗ Rsp ∗1000 ∗ T 0 ) ;lambda = 1 ;Q = 50;Ea= 5 0 ;q = Q∗ Rsp ∗1000∗ T 0 ; # KJ/Kg
e a= Ea∗ Rsp ∗1000∗ T 0 ; # KJ/Kg[ D c j , P c j , r h o c j , V c j , C c j , M cj ,D = sqrt ( f ) ∗ D cj ;T cj ,

u cj ] = znd polyfn (P 0 ,r h o 0 , T 0 , Q, gamma, Rsp , c 0 ) ;tspan = [ 0 0.5 0 . 9 9 9 9 ]o p t = o d e s e t ( ” R e l T o l ” , 1 e −12 , ” I n i t i a l S t e p ” , 1 e −1 , ” MaxStep ” , 1 e −1) ;# , ” A b s T o l ” , 1 e −8) ;[ lambda , x t ] = ode45 (@( lambda , x t ) z n d i n t e g r a t e ( lambda , xt , P 0 , r h o 0 , T 0 , Q, gamma, Rsp , c 0 , D, Ea , k ) ,tspan , [ 0 , 0 ] ) ;t12 = xt ( 2 , 2 )x12 = x t ( 2 , 1 )xr = xt ( 3 , 1 )N12 = 4 0 ; #Number o f p o i n t s p e r h a l f −r e a c t i o n z o n e l e n g t h
dx = x12 /N12 ; function [ dL , d t ] = z n d i n t e g r a t e x ( x , y , P 0 , r h o 0 , T 0 , Q, gamma, Rsp , c 0 , D, Ea , k )[ P , U, c , rho , V, x ,M, T ] = o v e r d r i v e n D ( P 0 , r h o 0 , T 0 , Q, gamma, Rsp , c 0 , D, y ( 1 ) ) ;dL = ((1 − y ( 1 ) ) . ∗ k . ∗ e .ˆ( − Ea∗ T 0 . / ( T) ) ) / (U.−D) ;d t = 1 . 0 / (U.−D) ;endfunction
xspan = 0 : dx : x r ;#N= f l o o r ( x r / d x ) +1;#x s p a n = l i n s p a c e ( 0 , x r , N) ;o p t = o d e s e t ( ” R e l T o l ” , 1 e −12 , ” I n i t i a l S t e p ” , 1 e −1 , ” MaxStep ” , 1 e −1) ;# , ” A b s T o l ” , 1 e −8) ;[ x s o l , l t ] = ode45 (@( x , l t ) z n d i n t e g r a t e x ( x , l t , P 0 , r h o 0 , T 0 , Q, gamma, Rsp , c 0 , D, Ea , k ) , xspan ,[0 ,0]) ;x = xsol ;lambda = l t ( : , 1 ) ;t= lt (: ,2) ;P = 0∗ lambda ;U = 0∗ lambda ;c = 0∗ lambda ;r h o = 0∗ lambda ;V = 0∗ lambda ;xFD = 0∗ lambda ;M = 0∗ lambda ;T = 0∗ lambda ;j = length ( lambda ) ;for i = 1 : j ;[ P( i ) , U( i ) , c ( i ) , r h o ( i ) ,V( i ) ,xFD( i ) ,M( i ) ,T( i ) ] = o v e r d r i v e n D ( P 0 ,lambda ( i ) ) ;endr a t e = k∗(1− lambda ) . ∗ e .ˆ( − Ea∗ T 0 . / T) ;# f i g u r e (1)# p l o t ( x , P/ P 0 , ” o r ” )# grid# t i t l e ( ’ Pressure vs x ’ )# f i g u r e (2)# p l o t ( x , U/ c 0 )# grid# t i t l e ( ’ P a r t i c l e Velocity vs x ’ )# f i g u r e (3)# p l o t ( t / t 1 2 , T/ T 0 , ” r ” )# grid# t i t l e ( ’ Temperature vs t ’ )# f i g u r e (4)r h o 0 , T 0 , Q, gamma, Rsp , c 0 , D,

70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#

p l o t ( x , r h o / r h o 0 , ”+ b k ” )gridt i t l e ( ’ Density vs x ’ )f i g u r e (5)p l o t ( x , lambda , ”+ b k ” )gridt i t l e ( ’ lambda v s x ’ )f i g u r e (6)p l o t ( t / t 1 2 , lambda , ”+ b k ” )gridt i t l e ( ’ lambda v s t ’ )f i g u r e (7)p l o t ( x / x12 , T/ T 0 , ” r ” )gridt i t l e ( ’ Temperature vs x ’ )f i g u r e (8)p l o t ( t / t12 , r a t e )t i t l e ( ’ Rate v s time ’ )

#
#
#
#

fid p
fid U
fid Y
fid T

#
#
#
#

fskiplfskiplfskiplfskipl====fopen (” p ” ,f o p e n ( ”U” ,f o p e n ( ”Y” ,f o p e n ( ”T” ,” r+t ” )” r+t ” )” r+t ” )” r+t ” )( fid p ,22)
( fid T ,22)
( fid Y ,22)
( fid U ,22)

# m = l e n g t h (P) ;# f o r i = m: − 1 : 1%f p u t s ( f i d , ’% f \n ’ , P ( 1 : 4 ) )#f p r i n t f ( f i d p , ’ % 3 0 . 1 5 f ’ , P( i ) ∗ 1 0 0 0 )%f p u t s ( f i d , ’% f \n ’ , P ( 1 : 5 ) )#fskipl ( fid p )##f p r i n t f ( fid T ,f s k i p l ( fid T )’ % 3 0 . 1 5 f ’ , T( i ) )##f p r i n t f ( fid Y ,f s k i p l ( fid Y )’%30.15 f ’ ,# f p r i n t f ( fid U ,#f s k i p l ( fid U )q ∗(1− l a m b d a ( i ) ) )’ ( % 3 0 . 1 5 f 0 0 ) ’ , U( i )−U( end ) )# endfor#f c l o s e ( f i d p )# fclose ( fid U )# fclose ( fid Y )#fclose ( fid T )%o p e n
fid P
fid U
fid T
fid Y

the f i l e s for reading= fopen ( ’ p ’ , ’ r t ’ ) ;= fopen ( ’U ’ , ’ r t ’ ) ;= fopen ( ’T ’ , ’ r t ’ ) ;= fopen ( ’Y ’ , ’ r t ’ ) ;only%r e a d s t h e c o n t e n t s o f t h e f i l e s i n t o v a r i a b l e s . Those v a r i a b l e s%c e l l array v a r i a b l e s
P f i l e d a t a = t e x t s c a n ( f i d P , ’%s ’ , ’ D e l i m i t e r ’ , ’ \n ’ ) ;U f i l e d a t a = t e x t s c a n ( f i d U , ’%s ’ , ’ D e l i m i t e r ’ , ’ \n ’ ) ;T f i l e d a t a = t e x t s c a n ( f i d T , ’%s ’ , ’ D e l i m i t e r ’ , ’ \n ’ ) ;Y f i l e d a t a = t e x t s c a n ( f i d Y , ’%s ’ , ’ D e l i m i t e r ’ , ’ \n ’ ) ;%T h i s i s t o e l i m i n a t e one l a y e r%n e e d e d , b u t i t i s .
P f i l e d a t a = P f i l e d a t a {1};U f i l e d a t a = U f i l e d a t a {1};T f i l e d a t a = T f i l e d a t a {1};Y f i l e d a t a = Y f i l e d a t a {1};

of

the

cell

array .

areI ’m n o t 100% s u r e why i t ’ s

144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217

%c l o s e t h e o r i g i n a lfclose ( fid P ) ;fclose ( fid U ) ;fclose ( fid T ) ;fclose ( fid Y ) ;files%s e e how any p o i n t s
m = length (P) ;

are

in

the

znd

solution%r e p l a c e t h e v a l u e s w i t h t h e z n d d a t a i n t h e p r o p e r l i n e s
f o r i = m: − 1 : 1
P f i l e d a t a {21+m−i } = s p r i n t f ( ’%f ’ , P( i ) ∗ 1 0 0 0 ) ;U f i l e d a t a {21+m−i } = s p r i n t f ( ’ (% f 0 0 ) ’ ,U( i )−U( end ) ) ;T f i l e d a t a {21+m−i } = s p r i n t f ( ’%f ’ , T( i ) ) ;Y f i l e d a t a {21+m−i } = s p r i n t f ( ’%f ’ , ( ( ( 1 − lambda ( i ) ) ∗q ) ) ) ;end%o p e n f i l e s f o r w r i t i n g . I f you u s e t h e same f i l e n a m e s ,%t h e c o n t e n t s o f t h e f i l e w i l l b e o v e r w r i t t e n
f i d P = fopen ( ’ p ’ , ’ wt ’ ) ;f i d U = fopen ( ’U ’ , ’ wt ’ ) ;f i d T = fopen ( ’T ’ , ’ wt ’ ) ;f i d Y = fopen ( ’Y ’ , ’ wt ’ ) ;i . e . p , T , Y, U

%w r i t e t h e f i l e d a t a t o t h e f i l e
f p r i n t f ( f i d P , ’%s \n ’ , P f i l e d a t a { : } ) ;f p r i n t f ( f i d U , ’%s \n ’ , U f i l e d a t a { : } ) ;f p r i n t f ( f i d T , ’%s \n ’ , T f i l e d a t a { : } ) ;f p r i n t f ( f i d Y , ’%s \n ’ , Y f i l e d a t a { : } ) ;%c l o s e t h e f i l e s
fclose ( fid P ) ;fclose ( fid U ) ;fclose ( fid T ) ;fclose ( fid Y ) ;that#p r i n t −d p n g

p l o t 1 . png#p r i n t −d p n g

p l o t 1 . png

you a r e

outputtingto%f u n c t i o n e v a l u a t e s CJ p r o p e r t i e sfunction [ D c j , P c j , r h o c j , V c j , C c j , M cj , T c j , u c j ] = z n d p o l y f n ( P 0 ,#M c j = ( 2 ∗ q /R∗ T 0 )p h i = ( 2 ∗ Q∗ ( ( (gammaˆ 2 ) −1)/ gamma) ) ˆ 0 . 5 ;#p ( x ) = mˆ2 − 1 − p h i ∗ m.
p = [ 1 , −phi , − 1 ] ;a = roots ( p ) ;M cj = max( a ( 1 ) , a ( 2 ) ) ;#D c j = M c j ∗ C c j ;P c j = ( 1 + (gamma ∗ ( M cj ˆ 2 ) ) ) / ( 1 + gamma) ;P cj = P cj ∗ P 0 ;V c j = ( 1 + (gamma ∗ ( M cj ˆ 2 ) ) ) / ( M cj ˆ2 ∗ ( 1 + gamma) ) ;V cj = V cj / rho 0 ;T c j = ( P c j ∗ V c j ) / Rsp ;r h o c j = 1/ V c j ;C c j = s q r t (gamma ∗ Rsp ∗1000 ∗ T c j ) ;D c j = M cj ∗ c 0 ;u c j = D c j ∗(1− r h o 0 / r h o c j ) ;r h o 0 , T 0 , Q, gamma, Rsp , c 0 )%f u n c t i o n e v a l u a t e s r e a c t i o n p r o g r e s s a s a f u n c t i o n o f s p a c e .
function [ dL , d t ] = z n d i n t e g r a t e x ( x , y , P 0 , r h o 0 , T 0 , Q, gamma, Rsp , c 0 , D, Ea , k )[ P , U, c , rho , V, x ,M, T ] = o v e r d r i v e n D ( P 0 , r h o 0 , T 0 , Q, gamma, Rsp , c 0 , D, y ( 1 ) ) ;dL = ((1 − y ( 1 ) ) . ∗ k . ∗ e .ˆ( − Ea∗ T 0 . / ( T) ) ) / (U.−D) ;d t = 1 . 0 / (U.−D) ;endfunction
function [ dx , d t ] = z n d i n t e g r a t e ( lambda , x , P 0 , r h o 0 , T 0 , Q, gamma, Rsp , c 0 , D, Ea , k )[ P , U, c , rho , V, x ,M, T ] = o v e r d r i v e n D ( P 0 , r h o 0 , T 0 , Q, gamma, Rsp , c 0 , D, lambda ) ;

218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234

dx = (U.−D) /((1 − lambda ) . ∗ k . ∗ e .ˆ( − Ea∗ T 0 . / ( T) ) ) ;d t = 1 . 0 / ( ( 1 − lambda ) . ∗ k . ∗ e .ˆ( − Ea∗ T 0 . / ( T) ) ) ;endfunction%f u n c t i o n e v a l u t e s p a r a m e t e r s o f an o v e r d r i v e n f u n c t i o n .
function [ P , U, c , rho , V, x ,M, T ] = o v e r d r i v e n D ( P 0 , r h o 0 , T 0 , Q, gamma, Rsp , c0 , D, lambda )M = D/ c 0 ;x = ( ( (M − M.ˆ −1) ˆ2 − 2 . ∗ (gamma. ˆ 2 − 1 ) . ∗ ( lambda . ∗ Q / gamma) ) . / ( (gamma + 1 ) . ˆ 2 . ∗ M. ˆ 2 ) ) . ˆ 0 . 5 ;V = ( (gamma + (Mˆ−2) ) / (gamma + 1 ) ) − x ;V = V/ r h o 0 ;P = ( ( (gamma∗Mˆ 2 ) + 1 ) / (gamma + 1 ) ) + gamma ∗ Mˆ2 ∗ x ;P = P ∗ P 0;
U = ( (M − (Mˆ−1) ) / (gamma + 1 ) ) + M∗x ;U = U∗ c 0 ;r h o = 1 . /V ;T = P . ∗ V. / Rsp ;c = s q r t (gamma. ∗ Rsp . ∗ 1 0 9 0 . ∗T) ;Listing 2: Shortcuts for ComputeCanada and Linux
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

cp −r chebzy@graham . computecanada . c a : ˜ / s c r a t c h / shockTube o u t p u t l i m%c o p y f i l e s f r o m c o m p u t e c a n a d a t odirs c p −r chebzy@graham . computecanada . c a : ˜ / s c r a t c h / d e t o n a t i o n 1 0 p t s g r a h a m t e s t c o m p l e t e%c o p y f i l e s f r o mcomputecanada t o d i r
s c p −r . / shockTube chebzy@graham . computecanada . c a : ˜ / shockTubexs c p −r ˜/ g r a h a m t e s t c o m p l e t e / d e t o n a t i o n 1 0 p t s chebzy@graham . computecanada . c a : ˜ / // t o s e n d t o r e m o t e/home/ c h e b z y / s c r a t c h / d e t o n a t i o n 2 0 p t s /systemrm −r ∗−05 remove m u l t i p l e f i l e s e n d i n g w i t h −05# s w a p o f f −v / s w a p f i l e# rm / s w a p f i l e%d e l e t e s w a p f i l e f o r m e d when you
savingrm − r f /home/ c h e b z y / g r a h a m t e s t c o m p l e t e / d e t o n a t i o n 1 p t s 2 d / p r o c e s s o r ∗pwd p r i n t w o r k i n g d i r e c t o r y#! / b i n / b a s h#SBATCH −−t i m e = 0 0 : 0 1 : 0 0#SBATCHa c c o u n t =d e f −c k i y a n dscancel%Use command w i t h j o b ID t o c a n c e l a j o b#! / b i n / b a s h#SBATCH −−t i m e = 0 0 : 3 0 : 0 0#SBATCH −−a c c o u n t=d e f −c k i y a n dmodule p u r g e
module load openfoam / 5 . 0
cp −r ˜/ shockTube
scratch /cd s c r a t c h / shockTube
cp −r 0 . o r i g 0/home/ c h e b z y / s c r a t c h / d e t o n a t i o n 1 0 p t s o r r e 0 5 /system
blockMesh
setFields
rhoCentralFoam
p o s t P r o c e s s −f u n c s a m p l e D i c t%p o s t p r o c e s s s a m p l i n g p o i n t s f o r d a t a

edit

a

file

and

exit

without

Flow v i s u a l i z t i o n o v e r a s e r v e r u s i n g P a r a v i e w%V i s u a l i z e y o u r s i m u l a t i o n s
1 . F i r s t , i n s t a l l on y o u r d e v i c e t h e same ParaView version a s t h e one a v a i l a b l e on t h e c l u s t e r you w i l lu s i n g ; l o g i n t o Cedar o r Graham and s t a r t a s e r i a l GPU i n t e r a c t i v e j o b .
s a l l o c −−t i m e = 1 : 0 0 : 0 −−n t a s k s =1m e m −p e r−cpu =2024 −−g r e s=gpu : 1 −−a c c o u n t=
make s u r e t h e p a r a v i e w you a r e r u n n i g i s t h e same a s t h a t i n t h e c l u s t e r . Load t h e v a r i a b l e s g r a d u a l l y ,one a t a t i m e .
s a l l o c −−t i m e = 0 : 3 0 : 0 −−n t a s k s =16 −−cpus−per−t a s k =2 −−mem−p e r−cpu =4024
module load p a r a v i e w−o f f s c r e e n / 5 . 3 . 0
srun pvservermesa
s s h c h e b z y @ c e d a r . compuetecanada . c a −L 1 1 1 1 1 : c d r 7 6 8 : 1 1 1 1 1a c c o u n t =d e f −c k i y a n d

fore f f e c t i v e u s e o f t h e p r o c e s s o r , run y o u r case i n t h e decomposed foam . I n p a r a v i e w . T h i s g i v e s
p r o c e s s i n g speed as in v i s u a l i z a t i o n speeds a l l o c −−t i m e = 1 : 3 0 : 0 −−n t a s k s =32 −−cpus−per−t a s k =3 −−mem−per−cpu =3024 −−a c c o u n t=d e f −c k i y a n d .
module load p a r a v i e w−o f f s c r e e n / 5 . 3 . 0
s r u n p v s e r v e r −−mesa
cp −r / d i r / d i r%u s e d

this

to

copy

be

files

cp −r /home/ c h e b z y / s c r a t c h / d e t o n a t i o n 8 0 p p t s / p o s t P r o c e s s i n g / s a m p l e D i c t / p p s c r i p t p l o t . py s c r a t c h /case2wallsoolong16 / postProcessing / sampleDict

the

best

47
48
49
50
51
52
53
54
55
56

57
58
59
60
61
62
63
64
65
66
67
68

module load python
python p p s c r i p t p l o t . py

% run p y t h o n%l o a d p y t h o n c o m p u t e c a n a d as c r i p t f o r making p l o t s o f p s h o c kagainsts p a c e and t i m e

s c p t o o k l o n g enough . F i r s t case was a t 1 1 : 4 6 l a s t was 1 6 : 4 8
p o s t P r o c e s s −f u n c ’ p a t c h i n t e g r a t e ( name=o u t l e t , mdot ) ’ %E v a l u a t e mass f l o w r a t e a c r o s s a s u r f a c e
t h i s case )t o u c h ” $ {PWD##∗/}.foam ” % t o c r e a t e . foam f i l e n e e d e d f o r v i s u a l i z a t i o n u s i n g p a r a v i e w( outletin

When b u i l d i n g y o u r case t o t e s t v a r i o u s v e l o c i t y s a m p l e s , make s u r e t h a t t h e i n l e t and t h e u n i f o r mv e l o c i t y v a l u e s a r e t h e same i f n o t you
w i l end up h a v i n g a f l o a t i n g p o i n t e x c e p t i o n e r r o r which w i l lmake l e a d t o a f o r k e r r o r i n computecanada .
cd g l o b u s c o n n e c t p e r s o n a l − 2 . 3 . 5 /./ globusconnectpersonal

Setting f i e l d default values−−> FOAM Warning :From function b o o l s e t C e l l F i e l d T y p e ( c o n s t Foam : : word &, c o n s t Foam : : fvMesh &, c o n s t l a b e l L i s t &, Foam : :I s t r e a m &) [ w i t h Type = Foam : : V e c t o r<d o u b l e >; Foam : : l a b e l L i s t = Foam : : L i s t <i n t >]C o n t r o l D i c t problem

ensiureCodes a r e

a t Github

available

that:starttimecorrespondschebzy thesis

to 0 ordesiredv a l u e ( n mapped c a s e s )Bibliography[BLWC13] Eric M Braun, Frank K Lu, Donald R Wilson, and José Camberos. Airbreathingrotating detonation wave engine cycle analysis. Aerospace Science and Technology, 27(1):201–208, 2013.
[BM95]

Anne Bourlioux and Andrew J. Majda. Theoretical and numerical structure of
unstable detonations. Philosophical Transactions: Physical Sciences and Engineering, 350(1692):29–68, 1995.

[BP94]T Bussing and George Pappas. An introduction to pulse detonation engines. In
32nd Aerospace Sciences Meeting and Exhibit, page 263, 1994.

[Bus96]

Thomas RA Bussing. Rotary valve multiple combustor pulse detonation engine,
May 7 1996. US Patent 5,513,489.

[BZV06]

Fedor A Bykovskii, Sergey A Zhdan, and Evgenii F Vedernikov. Continuous spin
detonations. Journal of Propulsion and Power, 22(6):1204–1216, 2006.

[Cfd]

Cfdsupport. Inletoutlet boundary condition. https://www.cfdsupport.com/
OpenFOAM-Training-by-CFD-Support/node114.html. Accessed: 2018-11-25.

[CFM+ 92] G Carrier, F Fendell, R McGregor, S Cook, and M Vazirani. Laser-initiatedconical detonation wave for supersonic combustion. Journal of Propulsion andPower, 8(2):472–480, 1992.
[Cha99]David Leonard Chapman. Vi. on the rate of explosion in gases. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 47(284):90–
104, 1899.

[CIR52]

Richard Courant, Eugene Isaacson, and Mina Rees. On the solution of nonlinear
hyperbolic diﬀerential equations by ﬁnite diﬀerences. Communications on Pureand Applied Mathematics, 5(3):243–255, 1952.

[Dah11]

AE Dahoe. Tutorial on deﬂagagration and detonations. Joint European summer school on fuel cell and hydrogen technology, University of Ulster, United
Kingdom, 8 2011.

[Dör43]

Werner Döring. Über den detonationsvorgang in gasen. Annalen der Physik,435(6-7):421–436, 1943.

[EGL91]

Shmuel Eidelman, W Grossmann, and I Lottati. Review of propulsion applications and numerical simulations of the pulsed detonation engine concept. Journalof Propulsion and Power, 7(6):857–865, 1991.

[FD79]

Wildon Ficket and William C Davis. Detonation: theory and experiment. Dover
Publications, New York, 1979.

[God57]Sergei Konstantinovich Godunov. Diﬀerence method of computation of shockwaves. Uspekhi matematicheskikh nauk, 12(1):176–177, 1957.

[God59]Sergei Konstantinovich Godunov. A diﬀerence method for numerical calculationof discontinuous solutions of the equations of hydrodynamics. Matematicheskii
Sbornik, 89(3):271–306, 1959.

[GSMK12] Boris E Gelfand, Mikhail V Silnikov, Sergey P Medvedev, and Sergey V Khomik.
Thermo-gas dynamics of hydrogen combustion and explosion. Springer Science& Business Media, 2012.
[HL94]

Thomas Y Hou and Philippe G LeFloch. Why nonconservative schemes converge
to wrong solutions: error analysis. Mathematics of computation, 62(206):497–530,
1994.

[hys16]

hysafe.org2016. Detonations and shockwaves. detonations and shock-waves - fundamentals of hydrogen safety, course materials., 08 2016.

[Jou05]

E Jouguet. Sur la propagation des réactions chemiques dans les gaz. Journal deMathématiques Pures et Appliquées, 1:347–425, 1905.

[Kai09]K Kailasanath. Research on pulse detonation combustion systems: a status
report. In 47th AIAA aerospace sciences meeting including the new horizons
forum and aerospace exposition, page 631, 2009.

[KBS15]

Matvey Kraposhin, Arina Bovtrikova, and Sergei Strijhak. Adaptation of
kurganov-tadmor numerical scheme for applying in combination with the piso

method in numerical simulation of ﬂows in a wide range of mach numbers. Procedia Computer Science, 66:43–52, 2015.
[KH13]CB Kiyanda and AJ Higgins. Photographic investigation into the mechanism of
combustion in irregular detonation waves. Shock Waves, 23(2):115–130, 2013.

[KNP01]

Alexander Kurganov, Sebastian Noelle, and Guergana Petrova. Semidiscrete
central-upwind schemes for hyperbolic conservation laws and hamilton–jacobi
equations. SIAM Journal on Scientiﬁc Computing, 23(3):707–740, 2001.

[KS13]

Kailas Kailasanath and DA Schwer. Rotating detonation engine research at nrl.
Technical report, 2013.

[KT00]

Alexander Kurganov and Eitan Tadmor. New high-resolution central schemes
for nonlinear conservation laws and convection–diﬀusion equations. Journal ofComputational Physics, 160(1):241–282, 2000.

[KWG11] Jan Kindracki, Piotr Wolański, and Zbigniew Gut. Experimental research on therotating detonation in gaseous fuels–oxygen mixtures. Shock waves, 21(2):75–84,
2011.
[LB14]

Frank K Lu and Eric M Braun. Rotating detonation wave propulsion: experimental challenges, modeling, and engine concepts. Journal of Propulsion andPower, 30(5):1125–1142, 2014.

[Lee08]John H. S. Lee. The Detonation Phenomenon. Cambridge University Press, 2008.

[MMM98] Franklin Mead, Jr, Leik Myrabo, and Donald Messitt. Flight and ground testsof a laser-booted vehicle. In 34th AIAA/ASME/SAE/ASEE Joint Propulsion
Conference and Exhibit, page 3735, 1998.
[MTE17a] LF Gutiérrez Marcantoni, J Tamagno, and S Elaskar. Two-dimensional numericalR
simulations of detonation cellular structures in h2o2ar mixtures with openfoam.
International Journal of Hydrogen Energy, 42(41):26102–26113, 2017.
[MTE17b] LF Gutiérrez Marcantoni, J Tamagno, and Sergio Elaskar. rhocentralrﬀoam:
An openfoam solver for high speed chemically active ﬂows–simulation of planar
detonations–. Computer Physics Communications, 219:209–222, 2017.
[Nic57]

R. B Morrison JA Nicholls, HR Wilkinson. Intermittent detonation as a thrustproducing mechanism. Journal of jet propulsion, 27(5):534–541, 1957.

[NNL00]

H.D. Ng, N. Nikiforakis, and J.H.S Lee. Assessment of high resolution centeredscheme for detonation modelling. Technical report, McGill University, Montreal,
Canada, 2000.

[NT90]

Haim Nessyahu and Eitan Tadmor. Non-oscillatory central diﬀerencing for hyperbolic conservation laws. Journal of computational physics, 87(2):408–463, 1990.

[Ope16]OpenCFD. OpenFOAM - The Open Source CFD Toolbox - User’s Guide.
OpenCFD Ltd., United Kingdom, 5 edition, 7 April 2016.

[PFHS15] Daniel E Paxson, Matthew Fotia, John Hoke, and Frederick Schauer. Comparisonof numerically simulated and experimentally measured performance of a rotating
detonation engin. In 53rd AIAA Aerospace Sciences Meeting, page 1101, 2015.
[PJBJ94]

Kumar A. (eds) Powers J.M. Buckmaster J., Jackson T.L. Oblique detonations:
Theory and propulsion applications. In Combustion in High-Speed Flows. ICASE/LaRC Interdisciplinary Series in Science and Engineering, vol 1. Springer,
Dordrecht, pages 345–347, 1994.

[Sha15]

Nihar N Shah. Computer modeling of a rotating detonation engine in a rocket
conﬁguration. Technical report, Air force Institute of Technology WrightPatterson AFB OH Graduate School of Engineering and Management, 2015.

[ST65]KI Shchelkin and Ya K Troshin. Gas dynamics of detonations. Mono Book Corp.,
Baltimore, 10, 1965.

[Tat15]Christopher Tate. Investigation of pulse detonation engines: the eﬀect of variableblockage ratio on the deﬂagration to detonation transition. 2015.

[TFW10]Akihiro Tobita, Toshitaka Fujiwara, and Piotr Wolanski. Detonation engine andﬂying object provided therewith, August 31 2010. US Patent 7,784,267.

[Tor]

Eleuterio F Toro. Riemann solvers and numerical methods for ﬂuid dynamics: apractical introduction.

[TSS+ 14] Colin A Towery, Katherine M Smith, Prateek Shrestha, Peter E Hamlington,
and Marthinus Van Schoor. Examination of turbulent ﬂow eﬀects in rotating
detonation engines. In 44th AIAA Fluid Dynamics Conference, page 3031, 2014.
[VJ14]

D Mahaboob Valli and DRTK Jindal. Pulse detonation engine: Parametersaﬀecting performance. International Journal of Innovative Research in Science,
Engineering and Technology, 3(4), 2014.

[VL76]

B Van Leer. Muscl, a new approach to numerical gas dynamics. In Computing
in Plasma Physics and Astrophysics, 1976.

[VM63]BV Voitsekhovskii and VV Mitrofanov. M. e˜ topchiyaa, structure of the detonation front in gases [ia russian], izd. Sibirsk. Otd. Akad. Nauk SSSR, Novosibirsk,
1963.

[VN42]

John Von Neuman. Theory of detonation waves. Technical report, Institute forAdvanced Study Princeton NJ, 1 Einstein Drive Princeton, New Jersey 08540
USA, 1942.

[WTJF98] Henry G Weller, Gavin Tabor, Hrvoje Jasak, and Christer Fureby. A tensorial approach to computational continuum mechanics using object-oriented techniques.
Computers in physics, 12(6):620–631, 1998.
[Zel40]Ya B Zeldovich. On the theory of detonation propagation in gaseous systems.
Zh. Eksp. Teor. Fiz, 10(5):542–568, 1940.

[ZM10]

Olindo Zanotti and Gian Mario Manca. A very short introduction to godunov methods. Albert Einstein Institute, Max-Planck-Institute for Gravitational
Physics, Potsdam, Germany, 2010.

[ZWW]

Rui Zhou, Dan Wu, and Jianping Wang. Progress of continuously rotating detonation engines. Chinese Journal of Aeronautics, 29.