Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction
(USC Thesis Other)
Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Direct Numerical Simulation of Mixing and Combustion
under Canonical Shock Turbulence Interaction
by
Xiangyu Gao
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulllment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(Mechanical Engineering)
August 2020
Copyright 2020 Xiangyu Gao
Acknowledgements
First and foremost, I would like to thank my advisor, Prof. Ivan Bermejo-Moreno for all the
patient and helpful instructions. Working with Ivan as his rst student was quite rewarding.
I really appreciate that the door of his oce was always open to me, and I was able to
stop by any time and discuss academic and life problems with him, even at 1am or 2am. I
also want to thank him for taking me to the 2018 Center for Turbulence Research Summer
Program. I was able to live in Stanford campus and experience the life at Palo Alto (where
I call the center of the universe) for a month, which was quite a special experience in my
life. Finally, I want to thank him for allowing me to explore opportunities in industry and
other elds, including supporting my master's degree in computer science.
I would like to thank Prof. Julian Domaradzki, Fokion Egolfopoulos, Carlos Pantano-
Rubino, Aiichiro Nakano for accepting to be in my qualifying exam and defense committee,
and providing academic insights to polish my research. Special thanks goes to Prof. Nakano
and his wife. I want to thank them for their hospitality and inviting me to their house for
tasty food and Maotai. I also want to thank my master thesis advisor, Prof. Jun Chen and
my undergrad research advisor Prof. Xiaohui Su for getting me interested in the research of
uid mechanics and turbulence. I also want to thank Prof. Johan Larsson for sharing his
original version of Hybrid code with me and collaborate with me on my papers.
The greatest thanks belongs to my parents. I want to thank them for their endless
support while I pursue my new life and academic dream in the US, especially my father,
who taught me the importance for diligence.
I would like to thank my labmates Naili, Alex, Jon and Jonas for all the fun we had out
of the oce and research. I want to thank Dr. Subodh Tiwari for all the tasty drinks we
had together. I feel so lucky to meet him at Argonne National Lab and explore so many
nice restaurants in KTown with him. I also want to thank Dr. Zaki Hasnain for teaching
me so much English. Thanks also goes to Dr. Yongqian Ma. I cannot secure a job before
the outbreak of COVID-19 without his help in stochastic calculus. I also want to thank
ii
other of my friends, Dr. Qian Wan, Wenlong, Xiaobin, Qingquan, all my roommates, and
all the friends I met at DUT alumni association. You made my PhD journey more than an
academic challenge. I worked ve years in RRB 206. Even though it is about to be torn
down, my memory of this building and with all the sta in the building will live on.
This work is impossible without the computational resources from USC HPCC (High
Performance Computing Center), XSEDE (Extreme Science and Engineering Discovery En-
vironment), and Discretionary and INCITE program of ALCF (Argonne Leading Computa-
tional Facility).
iii
Table of Contents
Acknowledgements ii
List Of Tables vi
List Of Figures vii
Abstract xv
Chapter 1: Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.1 Shock Induced Mixing . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.1.1 Theoretical work . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.1.2 Experimental work . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.1.3 Numerical work . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.1.4 Eects of Compressibility on Turbulent Mixing . . . . . . . 10
1.2.2 Shock-Induced Ignition . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.2.2.1 Experimental Work . . . . . . . . . . . . . . . . . . . . . . . 11
1.2.2.2 Numerical Work . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3 Outline of present work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Chapter 2: Methodology 16
2.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.1.1 Constitutive Relations for Passive Scalar Mixing . . . . . . . . . . . . 18
2.1.2 Constitutive Relations for Multi-
uid Mixing and Combustion . . . . 18
2.2 Simulation Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.1 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.2 Time Marching Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.4 In
ow turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.4.1 3D In
ow Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.4.2 2D In
ow Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.5 Sponge layer and out
ow boundary condition . . . . . . . . . . . . . . . . . 32
2.6 Simulation cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.6.1 Passive scalar mixing in STI . . . . . . . . . . . . . . . . . . . . . . . 33
iv
2.6.2 Reactive 2D STI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
Chapter 3: Verication 36
3.1 Grid convergence for 3D STI . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2 Verication for 2D HIT simulations . . . . . . . . . . . . . . . . . . . . . . . 39
3.3 Verication for 2D STI simulations . . . . . . . . . . . . . . . . . . . . . . . 42
3.4 Verication for Combustion Capability . . . . . . . . . . . . . . . . . . . . . 46
3.4.1 Ignition Delay Time . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.4.2 1D Shock Induced Ignition . . . . . . . . . . . . . . . . . . . . . . . . 48
3.4.3 Shock Hydrogen Bubble Interaction . . . . . . . . . . . . . . . . . . . 49
Chapter 4: Results 57
4.1 Passive Scalar mixing in 3D STI . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.1.1 Shock-normal statistics . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.1.1.1 Reynolds stress, vorticity variance anisotropy, and relevant
streamwise locations . . . . . . . . . . . . . . . . . . . . . . 58
4.1.1.2 Scalar variance . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.1.1.3 Scalar variance budgets . . . . . . . . . . . . . . . . . . . . 61
4.1.1.4 Eect of Schmidt number . . . . . . . . . . . . . . . . . . . 65
4.1.1.5 Scalar Taylor-like microscale . . . . . . . . . . . . . . . . . . 66
4.1.1.6 Scalar dissipation rate budgets . . . . . . . . . . . . . . . . 67
4.1.1.7 Spectra and PDFs . . . . . . . . . . . . . . . . . . . . . . . 68
4.1.2 Flow topology and conditioned dissipation rate . . . . . . . . . . . . 71
4.1.2.1 Characterization of local
ow topology . . . . . . . . . . . . 71
4.1.2.2 One dimensional PDFs and conditioned scalar dissipation
distributions . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.1.2.3 Joint qr PDFs and conditioned scalar dissipation distributions 76
4.1.2.4 Proportion of topologies and their contribution to the scalar
dissipation . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.1.3 Alignment of passive scalar gradient, vorticity, shock-normal and strain
directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.1.3.1 Alignments conditioned on
ow topology across the shock . 83
4.1.3.2 Barycentric map representation of alignments . . . . . . . . 86
4.2 Autoignition in 2D HIT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.3 Reactive 2D Shock Turbulence Interaction . . . . . . . . . . . . . . . . . . . 91
Chapter 5: Conclusion 109
Bibliography 113
Appendices 135
Appendix A: Transport equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
Appendix A.1: Passive scalar variance . . . . . . . . . . . . . . . . . . . . . 135
Appendix A.2: Passive scalar variance dissipation rate . . . . . . . . . . . . 137
Appendix B: Regime diagram for reactive 2D STI simulations . . . . . . . . . . . 139
v
List Of Tables
2.1 3D passive scalar mixing simulation cases, line styles and symbols (for post-
shock state) used in gures. . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.2 2D Reactive STI simulation cases, line styles used in gures. Please note that
the resolution and domain size are for the rened part before the sponge layer. 34
vi
List Of Figures
1.1 Illustration of a scramjet engine highlighting the dierent physics at play,
including shock waves, turbulent boundary layers, fuel injection and mixing
with the oxidizer, recirculating cavity
ameholders, and nozzle
ow to provide
thrust. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1 Problem setup of STI and blended HIT precursor simulations. . . . . . . . . 23
2.2 Schematic representation of the variation of energy and enstrophy with time,
following Davidson (2015). . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.1 (a) Streamwise grid spacing and (b) k
max
proles for all the cases in table 2.1 37
3.2 Grid convergence test for caseF , with (M;M
t
;Re
) = (2:0; 0:3; 40) on increas-
ingly ner grids: 426 128
2
(dotted), 852 256
2
(dash-dotted), 1280 384
2
(dashed), 1700512
2
(solid). (a)k
max
, (b) transverse vorticity variance!
0
2
!
0
2
.
(c) scalar variance,
f
002
(d) scalar dissipation rate, ^ . Quantities in (b)-(d) are
normalized with their corresponding value immediately upstream of the shock. 38
3.3 Grid convergence test for caseK, with (M;M
t
;Re
) = (5:0; 0:3; 72) on increas-
ingly ner grids: 600 256
2
(dotted), 1200 512
2
(dash-dotted), 1800 768
2
(dashed), 2400 1024
2
(solid). (a) k
max
, (b) transverse vorticity variance
!
0
2
!
0
2
. (c) scalar variance, (d) scalar dissipation rate, ^ . Normalizations as in
gure 3.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.4 Streamwise proles comparing left- and right-hand side terms of (a) scalar
variance transport equation (3.2) for case J and (b) scalar dissipation trans-
port equation 3.1) for case A (at a higher resolution of 1700 512 512). . . 41
3.5 Proles of (a) time evolution of energy and enstrophy, compared with the
enstrophy decay in Fox & Davidson (2010) and (b) enstrophy spectra. . . . . 42
3.6 Streamwise proles of Reynolds averaged density (a), streamwise velocity (b),
pressure (c) and temperature (d). Each curve has been normalized by the value
immediately upstream of the unsteady shock region. Grey shaded regions and
horizontal line segments mark the extent of the unsteady shock region for each
case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
vii
3.7 Streamwise proles of streamwise Reynolds Stress (a), transverse Reynolds
stress (b), M
t
(c) and Re
(d), vorticity variance (e), Kolmogorov scale (f).
Each curve in (a); (b); (e); (f) has been normalized by the value immediately
upstream of the unsteady shock region. Grey shaded regions and horizontal
line segments mark the extent of the unsteady shock region for each case. . . 44
3.8 Streamwise proles of scalar variance (a), SDR (b). Each curve has been
normalized by the value immediately upstream of the unsteady shock region.
Grey shaded regions and horizontal line segments mark the extent of the
unsteady shock region for each case. . . . . . . . . . . . . . . . . . . . . . . . 45
3.9 Ignition delay time for 2H
2
: O
2
: 7Ar mixtures at dierent initial temperatures
and ambient pressure of 1 atm, comparing results from 0D simulations using
Cantera with 1D and 2D ignition-spot simulations using HybridX. . . . . . . 47
3.10 Illustration of the simulation set up of 1D shock induced ignition . . . . . . . 48
3.11 Instantaneous streamwise proles of
ow quantities at time 180s of 1D shock-
induced ignition simulations, comparing reference data by Zikoski (2011) with
present simulation results using the newly proposed time-marching scheme,
and Strang splitting. Each quantity is normalized by the initial value behind
the incident shock. The arrow indicates the direction of detonation front
propagation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.12 Instantaneous streamwise proles of
ow quantities at time 230s of 1D shock-
induced ignition simulations, comparing reference data by Zikoski (2011) with
present simulation results using the newly proposed time-marching scheme,
and Strang splitting. Each quantity is normalized by the initial value behind
the incident shock. The red arrow indicates the direction of the detonation
front propagation. The three letters mark A) detonation front, B) contact
discontinuity, C) rarefaction wave, respectively. . . . . . . . . . . . . . . . . 51
3.13 Illustration of the simulation set up of 2D shock hydrogen bubble interaction 51
3.14 Comparison between results of simulations of a shock/hydrogen-bubble in-
teraction obtained with HybridX in the present work (left) and the reference
data by Billet et al. (2008) (right) for increasing simulation time (1.5s, 2.0s,
2.5s, and 3.0s from top to bottom). Colormap shows the concentration of
H
2
with isolines of pressure superimposed. . . . . . . . . . . . . . . . . . . . 54
3.15 Comparison between results of simulations of a shock/hydrogen-bubble in-
teraction obtained with HybridX in the present work (left) and the reference
data by Billet et al. (2008) (right) for increasing simulation time (3.5s, 4.0s,
6.0s, and 8.8s from top to bottom). Colormap shows the concentration of
H
2
with isolines of pressure superimposed. . . . . . . . . . . . . . . . . . . . 55
viii
3.16 Comparison between results of simulations of a shock/hydrogen-bubble inter-
action obtained withHybridX in the present work (left) and the reference data
by Billet et al. (2008) (right) for increasing simulation time (13.6s, 16.6s,
25.6s, and 41.6s from top to bottom). Colormap shows the concentration
of H
2
with isolines of pressure superimposed. . . . . . . . . . . . . . . . . . . 56
4.1 Streamwise proles of Favre-averaged streamwise Reynolds stresses (a;b) and
vorticity variance anisotropy (c;d) for dierent simulation cases. Each curve
in (a;b) has been normalized by the value immediately upstream of the un-
steady shock region. Grey shaded regions and horizontal line segments in
(a;b) mark the extent of the unsteady shock region for each case, using the
corresponding line style. All unsteady shock regions are removed in (c;d) for
a clearer comparison among cases. Vertical segments starting from the hori-
zontal axis indicate the location of x
^ !y
. M, M
t
, and Re
values for each case
are given in parentheses in the legend following table 2.1. . . . . . . . . . . 59
4.2 Streamwise proles of Favre-averaged scalar variance for (a) cases with the
same M
t
and (b) cases with the same M. Each curve is normalized by the
value immediately upstream of the unsteady shock region, which is greyed out
and marked, for each case, by an horizontal segment with the same line style
as the corresponding curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3 Streamwise proles of non-negligible terms in the transport equation of scalar
variance for cases with the same M
t
(a;c;e) and cases with the same M
(b;d;f). Each curve is normalized by the value immediately upstream of the
unsteady shock region (removed for clarity). The preshock location is oset
to be at x = 0 for all cases, and the postshock location, x
post
, is marked
with symbols, (a;b) dilatational term; (c;d) turbulent diusion term; (e;f)
Dissipation rate of scalar variance. . . . . . . . . . . . . . . . . . . . . . . . 63
4.4 Streamwise proles of scaled scalar dissipation rate, ^ , for (a) cases with
the same M
t
and (b) cases with the same M. Each curve is normalized by
the value immediately upstream of the unsteady shock region (removed for
clarity). Symbols mark the postshock location, x
post
. Legend as in gure 4.3. 65
4.5 Streamwise proles of (a) Favre-averaged scalar variance and (b) its scaled rate
of dissipation for scalars with dierent Sc obtained from case H (M = 5:0,
M
t
= 0:3, Re
= 37). Each curve is normalized by the value immediately
upstream of the unsteady shock region, shaded in greyed. . . . . . . . . . . . 66
4.6 Streamwise proles of scalar Taylor-like microscale,
=
00
00
=^
1=2
, for
(a) cases with the same M
t
and (b) cases with the same M. Each curve is
normalized by the value immediately upstream of the unsteady shock region,
collapsed for better comparison of the eective jumps for all cases. . . . . . . 67
ix
4.7 Streamwise proles of right-hand-side terms in the Reynolds-averaged trans-
port equation for the scaled dissipation rate of scalar variance for cases with
the sameM
t
(a,c,e) and cases with the sameM (b,d,f). (a,b) Correlation be-
tween the
uctuating velocity and the gradient of dissipation rate,G; (c,d) In-
teraction between velocity gradient tensor and scalar gradient,H; (e,f) Molec-
ular diusivityI. Each curve is normalized byI
u
, the absolute value of the
molecular term immediately upstream of the unsteady shock, removed for
better comparison of the eective jumps for dierent cases. . . . . . . . . . . 69
4.8 Spectra of (a) passive scalar and (b) scalar variance for dierent simulation
cases at pre- and postshock locations. Each spectrum is normalized by its
integral over all wavenumbers. Spectra for cases with higher Re
( 70) have
been shifted up two decades. Postshock spectra are shifted by a factor of
3 relative to preshock counterparts. The black straight segment with5=3
slope marks the extent of the inertial range of scales for the higherRe
cases,
obtained from spectra of the kinetic energy (not shown). . . . . . . . . . . . 70
4.9 Comparison of time-averaged spectra of scalar dissipation on transverse planes
for dierent simulation cases in the (a) preshock and (b) postshock states.
Each spectrum is normalized to unitary integral. Spectra for cases with higher
Re
( 70) are shifted up one decade for clarity. (c) Comparison of spectra
of scalar dissipation at dierent streamwise locations for case K, where each
spectrum in (c) is normalized by the value at the same wave number in the
preshock spectrum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.10 PDFs of (a) passive scalar, (b) scalar variance, and (c) scalar dissipation for
case K at x
pre
(dotted), x
post
(solid), x
^ !y
(dashed), x
iso
(dash-dotted). Dash-
dot-dotted line in (a) corresponds to a Gaussian distribution with zero mean
and the standard deviation of the PDF of passive scalar at x
iso
. . . . . . . . 72
4.11 Flow topologies inpqr space mapped onto three representative planes with (a)
p< 0 (extension), (b)p = 0 (volume-preserving), and (c)p> 0 (contraction),
partitioned by the intersections of each plane with the dividing surfaces S1a
(dashed), S1b (dash-dotted), S2 (solid), and the r = 0 plane (dotted): SFS,
stable focus stretching; UFS unstable focus stretching; SFC, stable focus com-
pressing; UFC, unstable focus compressing; SNSS, stable-node/saddle/saddle;
UNSS, unstable-node/saddle/saddle; SN
3
, stable node/stable-node/stable-node;
UN
3
, unstable-node/unstable-node/unstable-node; SSN
3
, stable-star-node/stable-
star-node/stable-star-node; USN
3
, unstable-star-node/unstable-star-node/unstable-
star-node. The last two topologies correspond to the intersection of S1a and
S1b for p< 0 and p> 0, respectively. . . . . . . . . . . . . . . . . . . . . . . 73
x
4.12 PDFs ofp (a,b,c),q (d,e,f), andr (g,h,i), along with normalized distributions
of SDR conditioned on the respective invariant, extracted atx
pre
(a,d,g),x
post
(b,e,h), andx
iso
(c,f,i) for cases with dierent (M;M
t
;Re
). Each plot shows
the PDFs of the invariant on the bottom curves (left vertical axis), and the
conditioned SDR distributions on top (right vertical axis). . . . . . . . . . . 74
4.13 Joint qr-PDFs (black-white contour lines at 15% and 60% of the peak value
per plot) and conditioned SDR distributions (jet-scale contours masked below
1% of the peak of the qr-PDF) for case K at dierent streamwise locations
(left to right: x
pre
,x
post
,x
^ !y
,x
iso
) and values ofp (top to bottom: mean plus
1.5, 0 and -1.5 times the standard deviation per location). Short-dashed black
lines: S1a, S1b, S2 curves and q axis (see Suman & Girimaji, 2010). In the
color contour, red corresponds to the highest conditioned SDR in each gure,
and blue corresponds to zero. . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.14 Jointqr-PDFs (dashed black-white contour lines at 15% and 60% of the peak
value per plot) and conditioned SDR distributions (jet-scale contours masked
below 1% of the peak of the qr-PDF) for case I (left), J (middle), K (right)
at x
post
for p = 0. Color map represents the conditioned scalar dissipation
(blue to green to yellow to red, from zero to highest value in each plot). . . . 78
4.15 Streamwise evolution of the time-averaged proportion of (a)
ow topologies,
hT
p
i, and of (b) the scalar dissipation in each topology,h
pjT
i, for dierent sim-
ulations. Topologies ordered by largest relative contribution (top to bottom).
Unsteady shock regions are removed and streamwise axes are normalized by
the dissipation length scale immediately upstream of the shock. Symbols mark
the postshock location, x
post
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.16 PDFs of the cosine of the angle () between the passive scalar gradient,r,
and (a) the most extensional,, and (b) the most compressive,
, eigenvec-
tors of the velocity gradient tensor, conditioned on topology for dierent sim-
ulation cases, comparing preshock and postshock states. Green dotted lines
correspond to the preshock state for case E, representative of the preshock
state for all other simulation cases (not shown for clarity, as they nearly overlap). 85
4.17 PDFs of the cosine of the angle between the scalar gradient,r, and (a) the
intermediate eigenvector, , of the velocity gradient tensor, (b) the vorticity,
!, and (c) the streamwise axis, x, conditioned on the SFS topology (similar
PDFs of these quantities are observed for all other topologies). . . . . . . . . 86
xi
4.18 Barycentric map representation of the PDFs of alignment between the eigen-
vectors (most extensive,; intermediate,; and most compressive,
) of the
the strain-rate tensor and (a) passive scalar gradient (r), (b) vorticity (!),
and (c) the streamwise direction (x). Symbols mark the peak of the PDF,
whereas lines represent the isocontour at 80% of the peak value of each dis-
tribution. Segments normal to each triangle side mark the 1=4, 1=2, and 3=4
partition points. Insets zoom into the region of interest of each barycentric
map. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.19 Barycentric map representation of distributions of scalar dissipation rate con-
ditioned on the alignment between the,,
eigenvectors of the strain-rate
tensor and (a) scalar gradient (r) at x
pre
, (b) scalar gradient (r) at x
post
,
and (c) the streamwise direction (x) at x
post
. . . . . . . . . . . . . . . . . . . 90
4.20 The mean (top) and normalized standard deviation (bottom) for 2D reactive
HIT simulations of ignition delay time at dierent initial T , Re
and M
t
compared with 0D simulations for H
2
/O
2
/Ar mixture. . . . . . . . . . . . . . 92
4.21 Mean ignition delay time at dierent initial T and M
t
compared with 0D
simulation for JP7 and air mixture. . . . . . . . . . . . . . . . . . . . . . . . 93
4.22 Time evolution of normalized thermodynamic quantities in caseP , with y axis
normalized by preshock value, and x axis normalized by the domain length
before the sponge layer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.23 Comparison of dierent thermodynamic quantities between 1D laminar reac-
tive STI simulation and 2D turbulent reactive STI simulations at t = 230s,
with line style the same as table 2.2. . . . . . . . . . . . . . . . . . . . . . . 96
4.24 Comparison of dierent thermodynamic quantities between 1D laminar reac-
tive STI simulation and 2D turbulent reactive STI simulations at t = 370s,
with line style the same as table 2.2. . . . . . . . . . . . . . . . . . . . . . . 97
4.25 Time evolution of temperature contour of case L(2; 1250; 0:3) from t = 1:5
10
4
s to t = 3:5 10
4
s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.26 Normalized reactant consumption rate, _ m
R
= _ m
R;1
and
ame front velocity,
s=u
1
, for 1D laminar reactive STI before the sponge layer. Negative s implies
that the
ame propagates towards the inlet. . . . . . . . . . . . . . . . . . . 101
4.27 Normalized (a) reactant consumption rate and (b)
ame front speed from
caseL to caseP before the sponge layer, with linestyle the same as table 2.2.
The three vertical dotted lines correspond to t
1
, t
2
and t
3
, and the markers
correspond to t
peak
for dierent cases. . . . . . . . . . . . . . . . . . . . . . . 102
xii
4.28 Flame regime diagram presented for caseO(2; 600; 0:1) (top) and caseN(2; 600; 0:3)
(bottom) showing the joint distributions in theDa
t
-Re
t
plane of domain with
heat release larger than 10
8
J s
1
m
3
. The plot is colored by cell number for
snapshots t
1
, t
2
, t
peak
, t
3
from left to right. Regions with cell numbers below
50 are masked. The three contour lines represent 80%, 40% and 20% of the
peak cell number at the corresponding time instance. . . . . . . . . . . . . . 105
4.29 Flame regime diagram presented for case O(2; 600; 0:1) showing the joint dis-
tributions in the Da
t
-Re
t
plane of domain with heat release larger than 1e8
J/s/m
3
. The plot is colored by TFI (top) and HRR (bottom) for snapshotst
1
,
t
2
,t
peak
,t
3
from left to right. Regions with cell numbers below 50 are masked.
The three contour lines represent 80%, 40% and 20% of the peak cell number
at the corresponding time instance. . . . . . . . . . . . . . . . . . . . . . . . 107
4.30 Flame regime diagram for all the turbulent cases showing the joint distribu-
tions in theDa
t
-Re
t
plane of domain with heat release larger than 10
8
J s
1
m
3
at t
peak
. The plot is colored by TFI (top) and HRR (bottom). Regions with
cell numbers below 50 are masked. The three contour lines represent 80%,
40% and 20% of the peak cell number at the corresponding time instance. . . 107
5.1 Flame regime diagram presented for case L showing the joint distributions in
the Da
t
-Re
t
plane of domain with heat release larger than 1e8 J/s/m
3
. The
scatter plot is colored by cell number (top, in log scale), TKI (middle) and
HRR (bottom, in log scale) for snapshots t
1
, t
2
, t
peak
, t
3
from left to right.
Regions with cell numbers below 50 are masked. The three contour lines
represent 80%, 40% and 20% of the peak cell number at the corresponding
time instance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
5.2 Flame regime diagram presented for caseM showing the joint distributions in
the Da
t
-Re
t
plane of domain with heat release larger than 1e8 J/s/m
3
. The
scatter plot is colored by cell number (top, in log scale), TKI (middle) and
HRR (bottom, in log scale) for snapshots t
1
, t
2
, t
peak
, t
3
from left to right.
Regions with cell numbers below 50 are masked. The three contour lines
represent 80%, 40% and 20% of the peak cell number at the corresponding
time instance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
5.3 Flame regime diagram presented for caseN showing the joint distributions in
the Da
t
-Re
t
plane of domain with heat release larger than 1e8 J/s/m
3
. The
scatter plot is colored by cell number (top, in log scale), TKI (middle) and
HRR (bottom, in log scale) for snapshots t
1
, t
2
, t
peak
, t
3
from left to right.
Regions with cell numbers below 50 are masked. The three contour lines
represent 80%, 40% and 20% of the peak cell number at the corresponding
time instance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
xiii
5.4 Flame regime diagram presented for caseO showing the joint distributions in
the Da
t
-Re
t
plane of domain with heat release larger than 1e8 J/s/m
3
. The
scatter plot is colored by cell number (top, in log scale), TKI (middle) and
HRR (bottom, in log scale) for snapshots t
1
, t
2
, t
peak
, t
3
from left to right.
Regions with cell numbers below 50 are masked. The three contour lines
represent 80%, 40% and 20% of the peak cell number at the corresponding
time instance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
xiv
Abstract
Shock induced mixing and ignition are studied in the canonical shock-turbulence interac-
tion (STI) conguration via shock-capturing Direct Numerical Simulations (DNS). First, a
parametric study of passive scalar mixing in the STI conguration is performed, varying the
shock Mach number (M = 1:28 to 5), turbulence Mach number (M
t
= 0:1 to 0:4), Taylor
microscale Reynolds number (Re
40; 70), and Schmidt number (Sc = 0:5; 1; 2). Stream-
wise budgets across the shock of the transport of scalar variance and scalar dissipation rate
(SDR) are examined. The dominant quantity of scalar variance transport is SDR. There
are two dominant terms in SDR transport equation. The rst one represents the interaction
between scalar gradient and velocity gradient. The change of this term across the shock
increases with higher M, lower M
t
and higher Re
. The second term represents molecular
diusion. Both this term and the SDR show a larger change across the shock with higher
M, lower M
t
and lower Re
.
The signicance of dierent
ow topologies enhancing mixing across the shock is studied
from changes in probability density functions of the invariants of the velocity gradient tensor
and the distributions of scalar dissipation conditioned on these invariants. The stable fo-
cus stretching topology is dominant, while the unstable-node/saddle/saddle topology is the
most dissipative throughout the domain, despite variations across the shock. Pre- and post-
shock distributions of the alignment between strain-rate tensor eigenvectors and the scalar
gradient, the vorticity, and the mean streamwise vector conditioned on
ow topology are
studied. A novel barycentric map representation is introduced for more direct visualization
of alignments and conditioned scalar dissipation distributions. Interaction with the shock
xv
increases alignment of the scalar gradient with the most extensive eigenvector, decreasing
it with the most compressive, which is still dominant. Across the shock, the most prob-
able alignment between the passive scalar gradient and the eigenvectors of the strain-rate
tensor converges towards the alignment with the largest scalar dissipation, enhancing scalar
dissipation immediately downstream of the shock.
Shock induced ignition in the STI setting is studied with 2D DNS using nite-rate detailed
chemistry. To verify the implementation of the combustion capability, comparison with
several benchmark cases is performed. The ability of 2D STI simulation to reproduce relevant
physics present in the 3D STI conguration is assessed. The, 2D reactive homogeneous
isotropic turbulence (HIT) simulations are performed to explore the in
uence ofRe
andM
t
on the ignition delay time, nding that larger M
t
triggers early ignition.
Four reactive STI simulations are conducted with initial Re
= (1250; 600) and M
t
=
(0:1; 0:3) at M = 2 with H
2
, O
2
and Ar mixture, and compared with laminar simulation at
the same Mach number. Compressibility of the upstream turbulence leads to earlier ignition
compared with laminar simulation, and the peak values of thermodynamic quantities at the
ame front are found to be lower for the turbulent cases under consideration, owing to the
partially premixed nature of the mixture.An analysis of the temporal evolution of the
ame
regime reveals that the reaction happens mostly in the thin reaction zone regime, which is
characterizes by a broadened preheat zone. Lower M
t
brings a slightly higher probability
that the combustion happens in regime of corrugated
amelets. The time evolution of the
Takeno Flame Index (TFI) indicates that, as the
ame propagates upstream, the combus-
tion becomes increasingly non-premixed because of the larger variance of the species mass
fractions.
xvi
Chapter 1
Introduction
1.1 Motivation
In a sunny afternoon in 1967, a rocket-powered X-15 crossed the sky of southern California.
This X-15 hit a record in human history, the fastest manned, powered, winged aircraft at
that time, with a speed of Mach 6.7 (Thompson, 2003; Urzay, 2018).The pursuit of faster
aircraft did not stop after that. The record hitting X-15 was an engine
ight test of NASA's
Hypersonic Research Engine (HRE) project, which started in 1964 (Hallion et al., 1987).
HRE's aimed
ight engine test was canceled in 1968 and was replaced with ground tests of
full-scale engine models (Andrews, 1994). At the end of HRE in 1975, a structural model and
combustion/propulsion model was proposed based on wind tunnel tests up to Mach number 7
(Andrews, 1994). The National Aero-Space Plane (NASP; 1986{1993) attempted the design
of the Rockwell X-30 passenger spaceliner, capable of
ight at Mach 25 (Barthelemy, 1989).
This project aimed at the construction of a vehicle that can routinely
y from Earth to space
and back, from conventional airelds, in aordable ways (Barthelemy, 1991). However, the
project was canceled in 1993 after several engine ground tests, because of a series of technical
diculties, e.g., the incapability to build materials to survive severe thermal and aerody-
namic environment, and the lack of detailed CFD models and computational resources, which
are still major technology challenges today (Piland, 1987; Barthelemy, 1991; Leyva, 2017).
While NASP turned out to be a failure, it still provided insights into the issues to be solved
1
for further advancement of hypersonic vehicles. Right after NASP was canceled, NASA's
Hyper-X project was initiated to explore the potential and viability of scramjet propulsion in
hypersonic vehicles (Freeman et al., 1997). Through the Hyper-X project, NASA made his-
tory with the rst and second successful
ights of a scramjet-powered airplane at hypersonic
speeds (Peebles, 2009), and from then on, the study of scramjets caught more attention.
Scramjet is short for supersonic combustion ramjet. Compared with rocket propulsion,
scramjet powered vehicles do not need to carry oxidizer, as they retrieve (\breathe") it from
the air in the atmosphere. Hence, a larger fraction of the take-o weight can be used for
payload. Figure 1.1 shows an illustration of a scramjet engine. The surrounding air, which
moves relative to the vehicle at speeds exceeding ve times the sound speed, rst undergoes
a compression produced by the vehicle's bow shock, followed by further compression through
a series of oblique shocks forming a shock train in the intake system. Then fuel is injected
to mix with the air before it reaches the combustor. Then the supersonic partially premixed
mixture of fuel and air are burnt in the combustor. The exhaust gas is discharged through
the nozzle/diuser to produce thrust.
In recent years, the number of
ight and ground tests performed on scramjet engines
has increased. One of the most studied tests is the second
ight of the HIFiRE (Hypersonic
International Flight Research Experimentation) program. This
ight focused on the transi-
tion from dual-mode (ramjet) to supersonic combustion. More detailed information of the
conguration and experiment can be found in Jackson et al. (2015). Wall-modeled Large
Eddy Simulations (LES) were also conducted by Bermejo-Moreno et al. (2013b); Saghaan
et al. (2015) and compared with ground tests of Hass et al. (2011). Reasonable agreement of
pressure traces with experiments showed that CFD can be used as a predictive tool to help
with the design of scramjets, but many open questions on modeling the complex
ow physics
that take place inside a scramjet remain. Besides HIFiRE, there have been other successful
scramjet research programs, e.g. HyShot Program (Smart et al., 2006) and SCRAMSPACE
2
Tirtey et al. (2014) at the University of Queensland (Australia), as well as ground exper-
imentation on scramjet combustors such as work by Gamba & Mungal (2015) at Stanford
University, just to name a few. A more complete account of scramjet projects can be found
in a recent review by Urzay (2018) and references therein. Also, there is also a comprehensive
review of scramjet experimental data base (Seleznev et al., 2019).
As mentioned above, the speed of the mixture in the combustor of a scramjet is super-
sonic, which translates into a characteristic residence time of the mixture in the order of
milliseconds (Urzay, 2018). Without optimization and additional design, eective combus-
tion and mixing between the fuel and air cannot occur inside the combustor within this short
residence time. One approach to increase this residence time is the transverse injection of
fuel relative to the incoming air stream. This can create large adverse pressure gradients and
a near-wall recirculation zone (Gamba & Mungal, 2015; Urzay, 2018). The adverse pressure
gradient can be strengthened with the impingement of an oblique shock. Residence time
can also be increased by adding a cavity-based
ameholder in the combustor, as shown in
gure 1.1. The recirculation and subsonic region in the cavity can eectively increase the
residence time (Gruber et al., 2004; Micka & Driscoll, 2009; Barnes & Segal, 2015). Besides
increasing residence time, premixed, shock-induced combustion (PSIC) is also very attrac-
tive in boosting performance of scramjets (Rubins & Bauer, 1994; McGuire, 2007; Axdahl,
2013; Urzay, 2018). In PSIC, fuel is injected to the air upstream of the inlet. The series of
oblique shocks in the inlet nozzle can enhance mixing of air and fuel mixture. At the same
time, the intake mixture experiences a large increase in pressure and temperature, which
may induce autoignition. Hence, PSIC can help improve scramjet performance in two ways.
The rst one is shock-induced mixing, the second one is autoignition at elevated temperature
and pressure.
In this work, we split our investigation of shock-induced ignition into two steps. The
rst step is to study passive scalar mixing in shock-turbulence interaction (STI), where only
the coupled eects of the shock and the turbulence on mixing are considered, with all other
3
Figure 1.1: Illustration of a scramjet engine highlighting the dierent physics at play, in-
cluding shock waves, turbulent boundary layers, fuel injection and mixing with the oxidizer,
recirculating cavity
ameholders, and nozzle
ow to provide thrust.
eects eliminated, e.g. boundary layer and shock bifurcation in shock tube experiments (Lip-
kowicz et al., 2018). The second step is to add nite-rate chemistry to study shock-induced
ignition. Even though Bermejo-Moreno et al. (2013b); Saghaan et al. (2015) achieved mod-
erately good agreement with experimental data, they adopted wall-modelled LES (WMLES)
and
amelet-based combustion models, furthermore considering chemical reactions to occur
innitely fast. The aim of the present work is to explore the physics of shock-induced mix-
ing and ignition with Direct Numerical Simulations (DNS) and nite-rate chemistry in the
STI conguration. A second objective of this work is to provide comprehensive numerical
datasets that can contribute to the development of better turbulent combustion models and
LES model for more computationally aordable simulation approaches of scramjet engines,
such as large-eddy simulation (LES).
1.2 Background
1.2.1 Shock Induced Mixing
The enhanced mixing that results from the amplication of turbulence across a shock wave
can be critical in both natural and man-made applications. This enhanced mixing can be
4
both benecial and detrimental. Shock-induced mixing enhancement has been used in scram-
jets for decades (Marble et al., 1987; Yang et al., 1993). As discussed before, shock-induced
mixing is benecial in hypersonic propulsion, but it is detrimental in Inertial Connement
Fusion (ICF), which uses high-energy laser beams to create shocks that compress and heat
fuel in order to induce thermonuclear reactions. However, uneven compression may lead to
the formation of Rayleigh-Taylor and Richtmyer-Meshkov instabilities, which can in turn
result in premature mixing and reduced heating ecacy. The study of shock-induced mixing
is also crucial for supernova explosions (Hughes et al., 1999; De Avillez & Mac Low, 2002).
The canonical STI conguration has been extensively studied with theories, experiments, and
through numerical simulations. However, very few studies have concentrated on quantifying
the eects of STI on passive scalar mixing, in particular as the relevant physical parameters
and regime of the interaction are changed, which is a primary focus of the present work.
1.2.1.1 Theoretical work
Theoretical studies of the interaction between shock waves and turbulence were pioneered
by Ribner (1953, 1954) and Moore (1953), who proposed linear interaction analyses (LIA)
to study the postshock state of turbulence, neglecting nonlinear and viscous eects. In LIA,
the preshock turbulence is assumed to consist of a superposition of dierent modes of small-
amplitude disturbances (Kov asznay, 1953) that are amplied separately across the shock.
Kerrebrock (1956) showed that any of the three modes of
uctuations (acoustic, entropy
and vorticity waves), can lead to all three modes after the shock. LIA, proposed more than
sixty years ago, is still widely used in the community with a number of recent extensions.
Griond (2005) extended LIA to a mixture of two perfect gases, and later applied it to study
the Richtmyer-Meshkov Instability (RMI) (Griond, 2006). Wouchuk et al. (2009) adopted
LIA to study the changes of an isotropic turbulent vorticity eld across the shock. Huete
et al. (2011) used LIA to investigate three-dimensional random isotropic density eld across
shockwaves. LIA has also been extended to account for the interaction between turbulence
5
and detonation/reactive shocks: Jackson et al. (1990, 1993) studied the interaction between
detonation and vorticity waves with constant density using linear theory and the assumption
that the thickness of the reaction zone is much smaller than turbulence length scales Huete
et al. (2013, 2014, 2017) pushed the analysis of detonations to a higher Mach number, where
the eects of entropy waves become more dominant when both the reaction-zone thickness
is smaller and larger than turbulence length scales, respectively. Quadros et al. (2016a,b)
investigated turbulent energy
ux and streamwise velocity-temperature correlations in STI
comparing LIA predictions and datasets from numerical simulations. Using LIA, Livescu &
Ryu (2016) predicted the postshock turbulence at high Taylor microscale Reynolds numbers,
unattainable with current simulations. Recently, Hejranfar & Rahmani (2019) assessed the
eects of real gas on the STI with linear analysis.
The other widely used theoretical tool to study STI is Rapid Distortion Theory (RDT).
Lele (1992) used the rapid distortion assumption to solve the closure problem in the exact
jump relations across a shock in a turbulent mean
ow. Cambon et al. (1993) adopted RDT
and observed increased kinetic energy growth rate under the rapid distortion assumption.
Jacquin et al. (1993) used homogeneous RDT to predict the amplication of turbulent ki-
netic energy across the shock and found large dierences from LIA due to the inhomogeneous
nature of shock problems. Kitamura et al. (2016) applied RDT to obtain analytical expres-
sions for the amplication of turbulence kinetic energy (TKE), enstrophy and the changes
of turbulence lengthscales of homogeneous (but possibly anisotropic) turbulence across the
shock, nding a dependence on the upstream anisotropy.
Besides LIA and RDT, new theoretical developments using dimensional and similarity
analyses have been proposed in the last decade (Donzis, 2012a,b; Chen, 2018). In these
studies a quasi-equilibrium assumption of instantaneous adjustment of the shock to the local
ow conditions is invoked to treat the shock surface as a collection of innitesimal laminar
shocks at dierent conditions. This allows predictions of the jumps of thermodynamic quan-
tities across the shock, with a wider range of applicability than LIA, for stronger upstream
6
turbulence. Also, in contrast with LIA, the results not only depend on the shock Mach
number but also on the turbulence Mach number, which is shown to aect the jump of mean
quantities across the shock.
1.2.1.2 Experimental work
Experiments have also contributed to the study of STI (Andreopoulos et al., 2000) and
shock-driven mixing. Hesselink & Sturtevant (1988) studied the eect of shock propagation
on the properties of a random medium, observing enhanced mixing in the shock-processed,
strongly compressed
uid. Marble et al. (1987) considered the interaction of weak oblique
shocks and cylindrical jets of hydrogen in air to produce misalignment between pressure
gradient and density gradient, which results in faster mixing by baroclinic vorticity. Menon
(1989) studied the interaction between shock waves and shear layers, observing an accelerated
spreading of the shear layer downstream of the shock. Honkan & Andreopoulos (1990, 1992)
studied the passage of a normal shock wave through grid-generated turbulence in a shock
tube. Barre et al. (1996) adopted a multinozzle concept to produce quasihomogeneous
isotropic turbulence and studied its interaction with normal shocks at Mach number 3,
nding close agreement with LIA in the amplication of longitudinal velocity
uctuation.
Agui et al. (2005) found that the interaction with shock waves also brought an enhancement
of
uctuations of the smaller turbulence scales. Andreopoulos (2008) studied the change of
velocity and vorticity alignments across the shock, and concluded that the
ows downstream
of the shock become more vortex-dominated than their upstream counterparts.
1.2.1.3 Numerical work
Numerical simulations play a decisive role in the study of STI and shock-driven mixing. The
geometric simplicity of the canonical STI conguration renders Direct Numerical Simulation
(DNS) as an ideal computational tool to study the fundamental physical mechanisms at play.
For most combinations of Mach and Reynolds numbers, the large gap between the smallest
7
scales of the turbulence and those of the shock makes shock-capturing DNS, in which the
turbulence scales are resolved but the shockwaves are numerically captured, an attractive
methodology that allows to signicantly increase the attainable Reynolds number for the
same computational power, compared with shock-resolving DNS. Shock-tting methods have
also been used to simulate STI (Sesterhenn et al., 2005; Rawat & Zhong, 2010), although the
tracking of the shock makes these algorithms not optimal for the commonly used structured
grids (Pepe et al., 2014), such as the ones used in the present work.
Early DNS studies by Lee et al. (1993, 1994) observed amplication of TKE and trans-
verse vorticity variance, and decreased turbulent length scales across the shock, in agreement
with LIA predictions. However, they also found a rapid evolution of turbulence kinetic en-
ergy after the shock, which cannot be predicted with LIA. In the DNS of Hannappel &
Friedrich (1995) an increased compressibility of the in
ow turbulence reduced the ampli-
cation of turbulence kinetic energy and thermodynamic
uctuations across the shock, but
increased the amplication of vorticity
uctuations.
Shock-capturing DNS by Lee et al. (1997) pushed their earlier simulations to a higher
Mach number, nding that the amplication of turbulence kinetic energy saturates for Mach
numbers larger than 3, which could have implications on the achievable shock-induced mixing
enhancement as the shock strength increases. In the DNS of Mahesh et al. (1997), the role of
upstream entropy
uctuations in the amplication of turbulent kinetic energy was studied.
It was observed that negative correlation between velocity and temperature
uctuations can
enhance the amplication, while positive correlation suppresses it. Jamme et al. (2002)
looked into the eects of dierent upstream
uctuations on STI, conrming that upstream
entropy
uctuations enhance the amplication of turbulent kinetic energy and transverse
vorticity variances across the shock. Additionally, the vorticity-entropy
uctuations led to
an increased reduction of the transverse Taylor microscale. They observed less amplication
of the kinetic energy, less reduction of the transverse microscale and more amplication
8
of the transverse vorticity variance across the shock when upstream acoustic and vorticity
uctuations were dominant.
Larsson & Lele (2009); Larsson et al. (2013) explored the parameter space of STI with
high-order schemes and relatively higher Mach number and Taylor microscale Reynolds num-
ber, covering both the `wrinkled' and `broken' shock regimes. They proposed a criterion to
determine the interaction regime, based on the scaling of the PDF of the density shock
jump under grid renement. This interaction regime criterion is aligned with that of Donzis
(2012a,b), which was obtained instead from similarity arguments as the ratio of the lami-
nar shock thickness and the Kolmogorov scale of the incoming turbulence. When this ratio
is small, implying a large separation of scales between the shock and the turbulence, and
for small enough turbulence Mach numbers, Ryu & Livescu (2014) noted that viscous and
nonlinear eects across the shock become small, and found LIA to reliably predict the am-
plication of all Reynolds stress and vorticity components immediately downstream of the
shock obtained from a large set of shock-resolving DNS. The criterion in (Donzis, 2012a)
also predicts a third, new regime of interaction that corresponds to highly intense incoming
turbulence and a relatively weak shock. This newly termed `vanished regime', in which tur-
bulence quantities do not present amplication across the shock but a monotonic decay, was
observed in the large set of shock-resolving DNS performed by Chen (2018).
Besides DNS, currently restricted to small Reynolds numbers due to its computational
cost, several numerical studies have focused on lower delity simulation approaches, such as
large-eddy simulations (LES) (Bermejo-Moreno et al., 2010; Hickel et al., 2014; Braun et al.,
2019) and Reynolds (or Favre) Averaged Navier Stokes (RANS) simulations (Sinha et al.,
2003; Quadros & Sinha, 2016; Schwarzkopf et al., 2016). Insight from and comparisons with
available DNS datasets are often used in these studies to inform modications to subgrid-
scale (SGS) and RANS models for STI. Simulations targeting Reynolds numbers much higher
than existing DNS resort to LIA for validation.
9
1.2.1.4 Eects of Compressibility on Turbulent Mixing
Studies of scalar mixing in the canonical STI conguration are scarce, in comparison with
other compressible
ow congurations, such as compressible homogeneous isotropic turbu-
lence (Pan & Scannapieco, 2010, 2011; Ni et al., 2015; Ni, 2015, 2016; Danish et al., 2016) and
the shock-driven Richtmyer-Meshkov instability (Brouillette, 2002; Lombardini et al., 2012;
Orlicz et al., 2015; Reese et al., 2018; Desjardins et al., 2018). Tian et al. (2017) analyzed
the eects of density variations on STI and mixing by comparing single- and two-
uid cases
using shock-capturing DNS. Faster increase of turbulence stretching and an increased scalar
dissipation rate across the shock were identied as the causes for better mixing enhance-
ment for multi-
uid mixing. Boukharfane et al. (2018) compared the streamwise evolution
of scalar variance and its dissipation rate between the shocked and unshocked cases, and
observed an increase of scalar dissipation rate (SDR) across the shock. They also related the
increase of SDR with the change of alignments across the shock.
These pioneering studies of scalar mixing in STI were however limited in the Reynolds
number as well as the shock and turbulence Mach numbers of the interactions being con-
sidered, which were all in the wrinkled-shock regime. Besides the limited number of stud-
ies on scalar mixing under STI, other works have investigated compressibility eects on
mixing (Pantano & Sarkar, 2002; Vaghe & Madnia, 2015; Karimi & Girimaji, 2016; Ja-
hanbakhshi & Madnia, 2016; Arun et al., 2019; Yu & Lu, 2019). Although not focused
on passive scalars, these investigations have identied
ow patterns that are favored and
dominate dissipation in locally compressed regions, relevant to STI.
1.2.2 Shock-Induced Ignition
As discussed above, the sudden increase of temperature and pressure experienced across the
shock can help ignite the mixture of fuel and air. For that reason, shock-induced ignition
has been adopted to boost performance of scramjets (Rubins & Bauer, 1994; McGuire, 2007;
Axdahl, 2013; Urzay, 2018). However, shock-induced ignition is not desirable in applications
10
such as the storage of H
2
(Wolanski, 1973; Golub et al., 2007; Mogi et al., 2008; Kim et al.,
2013). In a high-pressureH
2
storage reservoir, if an accidental release ofH
2
happens, the high
pressure dierence between the storage reservoir and surrounding air will create a shock wave.
The shock wave will propagate and heat up the surrounding air. Then diusion happens
between air and H
2
, which induces ignition and explosion at the contact interface. Shock-
induced ignition has been extensively studied in shock-tube experiments and simulations.
However, the literature about shock-induced ignition in canonical STI congurations (which
remove the eect of the walls) is still sparse.
1.2.2.1 Experimental Work
Shock tube experiments have traditionally been an important tool to study shock-induced
ignition. Some of the most notable experiments to measure ignition delay time (IDT) with
shock tubes were carried out by Meyer & Oppenheim (1971); Oppenheim et al. (1975); Cheng
& Oppenheim (1984). They adopted pressure transducers and laser-schlieren measurements
to detect the passage of the shock wave and the
ame front. These studies categorized ig-
nition into two regimes, weak ignition and strong ignition. Weak ignition is characterized
by the formation of isolated
ame kernels in the early stage of ignition. Then these kernels
develop and become connected with each other, followed by the formation of planar com-
bustion. Strong ignition is characterized by the instantaneous formation of the planar
ame
front.
Pressure transducers and laser-schlieren measurements are still powerful tools to measure
the IDT in shock tubes nowadays (Penyazkov et al., 2005; Huang et al., 2006; Pang et al.,
2009; Shen et al., 2009), but these techniques only provide line-of-sight data. McMillin et al.
(1991) innovatively adopted Planar Laser-Induced Fluorescence (PLIF) to high-speed reac-
tive
ow experiments and obtained instantaneous 2D images of shock-induced ignition. PLIF
measurements were also performed by Morris et al. (1996) on reactive
ows over a wedge test-
body in an expansion tube for both hydrogen-based and ethylene-based reactive mixtures. In
11
these experiments, shock-induced de
agration was observed after oblique shocks. Penyazkov
et al. (2005) conducted experiments at dierent temperature, pressure and equivalence ratios
on propane-air mixture, inferring empirical correlations of IDT. Shen et al. (2009) performed
a shock-tube study of the IDT of n-heptane, n-decane, n-dodecane, and, for the rst time,
n-tetradecane. More recent experimental studies of shock-induced ignition include Cancino
et al. (2010); Yamashita et al. (2012); Mathieu et al. (2015); Figueroa-Labastida & Farooq
(2020); Merkel & Ciccarelli (2020), just to name a few.
Recent developments of experimental techniques, such as 3D time-resolved Particle Image
Velocimetry (PIV) and planar laser-induced
uorescence (PLIF), have enabled measurements
beyond those provided by single-point probes. Even though experiments can provide infor-
mation at high Reynolds numbers, which is not feasible for DNS, they are still limited by
spatial resolution and derivative quantities (Westerweel et al., 2013). In this regard, DNS can
be a good complement to experiments and provide more comprehensive access to volumetric
ow elds.
1.2.2.2 Numerical Work
Numerical simulations of shock-induced ignition range from 1D simulation with one-step
chemistry to 3D DNS with detailed chemical kinetics. Melguizo-Gavilanes et al. (2011)
adopted 1D simulation with one-step Arrhenius reaction to track the evolution of the shock-
initiated ignition process, from slow formation and rapid growth of the hot spot to the birth
of a secondary shock and its transition into a detonation wave. Sharpe (2002); Sharpe &
Short (2003) adopted two-step chain-branching reactions to compare with one-step chem-
istry, nding dierences in the initial ignition spot acceleration and on the secondary shock
formation. There are also several 2D and 3D simulation studies of ignition behind re
ected
shocks. Ihme et al. (2013) studied the interaction between the re
ected shock and reactive
boundary layer for Argon-diluted hydrogen/oxygen mixtures. In this study, the re
ected
shock creates a bifurcation, which ingests cold
uid from the upstream boundary layer into
12
the stagnant region. This thermally stratied
ow leads to inhomogeneities at the near-wall
region, and then acts as an ignition kernel. Khokhlov et al. (2015) adopted 3D DNS to study
the development of hot spots and ignition behind the re
ected shock. They observed that
entropy
uctuations can lead to preferential locations for hot spot generation Lipkowicz et al.
(2018) performed 3D LES simulations with detailed chemistry on re
ected shock-induced ig-
nition in a shock tube. They observed early occurrence of mild ignition away from the end
wall. With the help of Lagrangian tracking of temperature history, it was concluded that
early ignition is brought by the acceleration of the re
ected shock, which leads to bifurcation
growth and high temperature away from the end wall. Additionally, there are also other
works about shock-induced ignition considering boundary layer eects in shock tubes, to
name a few, Chaos & Dryer (2010); Dziemi nska & Hayashi (2013); Grogan & Ihme (2015);
Kiverin & Yakovenko (2019).
Autoignition induced by elevated pressure and hot ignition spots has also been extensively
studied. Im et al. (1998) used 2D DNS with detailed chemistry to study the autoignition
of hydrogen/air mixing layer with large temperature gradient. It was observed that the
dependence of the ignition delay on turbulence intensity is non-monotonic. Echekki & Chen
(2002, 2003) also performed 2D DNS with detailed chemistry to study autoignition, but
with a partially-mixed hydrogen/air mixture. In this study, combustion is observed to be
initiated in lean premixed
ame fronts that subsequently propagate into richer mixtures.
Also, ignition is induced by an increase in HO
2
, which leads to the production of OH radicals,
and nally results in chain-branching reactions. Yoo et al. (2011) conducted 2D DNS with
detailed chemistry to study the autoignition of methane/air mixture. They revealed that
larger temperature
uctuation induces a larger temporal spreading of mean heat release, and
ignition happens more often with lower temperature
uctuation
The rst 3D DNS of autoignition with detailed chemistry is found in Yu & Bai (2013).
This work compared the dierences between 3D DNS and 2D DNS, and it concluded that
heat transfer is enhanced in 3D simulations because of a higher strain rate. Yao et al.
13
(2018) studied autoignition at dierent elevated pressures with 2D DNS and with the same
simulation conguration as Im et al. (1998). They also identied HO
2
as an indicator for
ignition. Chi et al. (2018) conducted 1D, 2D and 3D DNS simulations with detailed chemistry
to study misre and ignition delay time with hydrogen/air mixtures. Besides DNS, there
are also some LES simulations about autoignition. Jones et al. (2007) performed LES on a
hydrogen jet issuing into a turbulent co-
owing air stream with detailed chemistry. They
compared their simulation results with experiments and conrmed the accuracy of their
sub-grid species concentration and temperature
uctuations obtained from a sub-grid joint
PDF for the reactive scalars. Stankovi c et al. (2011) combined LES simulations with a rst-
order Conditional Moment Closure (CMC) approach to simulate a nitrogen-diluted hydrogen
jet, igniting in a turbulent co-
owing hot air stream. They found that stronger turbulence
promotes ignition.
According to the previous literature review, it can be concluded that 3D DNS simulation
with detailed chemistry in canonical STI congurations is not explored yet. However, 3D
DNS simulation of shock-induced ignition in the STI conguration with detailed chemistry
is still computationally prohibitive for most conditions of practical relevance. In Lipkowicz
et al. (2018), a 3D LES simulation with 1.4 billion cells took 3.8 million core hours. If
meshes with three times of the cell number are needed to conduct DNS at the same physical
condition, a computation cost of the order of 100 million core hours will be required. Instead,
2D DNS simulations will be performed in this work. As it will be shown in section 3.3, 2D
STI simulations can adequately capture relevant physical phenomena to be studied compared
with 3D STI. For that reason, 2D STI appears as a good compromise to gain further physical
insight of to reactive STI problem.
14
1.3 Outline of present work
The governing equations and key parameters, along with the computational approach, are
described inx2, for both the simulations of passive scalar mixing and combustion appli-
cations. Validation and verication are done inx3, including grid convergence study for
3D STI simulations for passive scalar mixing inx3.1, verication for 2D HIT simulation in
x3.2, similarity between 2D STI and 3D STI inx3.3, 1D and 2D ignition delay time inx3.4.1,
ignition induced by re
ected shock inx3.4.2, 2D shock hydrogen bubble interaction inx3.4.3.
Simulation results are presented inx4. scalar variance and SDR budgets are analyzed
inx4.1.1 for dierent shock and turbulent Mach numbers, Taylor microscale Reynolds num-
bers, and Schmidt numbers. PDFs of the invariants of the velocity gradient tensor and the
corresponding conditioned SDR are examined inx4.1.2.x4.1.3 presents a parametric study
of the changes across the shock of the alignments between the scalar gradient, strain-rate
eigenvectors and vorticity for dierent topologies, introducing a novel barycentric map rep-
resentation of alignments. The eects of M
t
and Re
on the ignition delay time of reactive
HIT are examined inx4.2. Major ndings about reactive 2D STI simulations are included in
x4.3 .The main conclusions are summarized inx5. Derivation of budget equations, gures for
verication of 2D shock bubble interaction, and combustion regime diagram for 2D reactive
STI simulations can be found inxAppendices.
15
Chapter 2
Methodology
2.1 Governing Equations
Shock-capturing DNS are used to solve the conservative form of the equations of mass,
momentum, and total energy conservation, complemented with advection-diusion-reaction
equations.
@
@t
+
@u
j
@x
j
= 0; (2.1)
@u
i
@t
+
@u
i
u
j
@x
j
=
@
ij
@x
j
; (2.2)
@e
T
@t
+
@e
T
u
j
@x
j
=
@u
i
ij
@x
j
@q
j
@x
j
; (2.3)
@Y
@t
+
@Y
u
j
@x
j
=
@J
j
@x
j
+W
_ !
; (2.4)
where Einstein summation convention is used except for Greek subscripts, is the density,
u
i
is the i-th component of the velocity (i = 1; 2; 3), and Y
is the mass fraction of the -th
species for multi-
uid mixing and combustion applications, or the value for the -th passive
scalar for passive scalar mixing.
The stress tensor is dened as:
ij
=p
ij
+ 2(S
ij
S
kk
ij
=3); (2.5)
16
wherep is the thermodynamic pressure,S
ij
=
1
2
(@u
i
=@x
j
+@u
j
=@x
i
) is the strain-rate tensor,
and is the dynamic viscosity. The ideal gas law, p = RT , is used to relate static ther-
modynamic variables in terms of gas constant, R. J
j
represents the diusive
ux of species
or passive scalar in the j-th direction, whose formula will be introduced later for both
passive scalar mixing and combustion applications. W
and _ !
are molecular weight and
molar production rate of species , respectively. e
T
=e +
1
2
u
i
u
i
is the total energy per unit
mass, wheree is the internal energy per unit mass. q
j
is the heat
ux inj-th direction. The
formulation ofe andq
j
will be introduced separately later for passive scalar and combustion
applications.
Reynolds and Favre averages of any
ow quantity, f, are denoted by
f and
e
f = f=,
respectively, with corresponding
uctuations dened as f
0
= ff and f
00
= f
e
f. For
simulations in canonical STI, averaging is done spatially over the statistically homogeneous
transverse directions (y =x
2
andz =x
3
for 3D simulations,y =x
2
for 2D simulations) and
temporally for non-reacting simulations, since the
ow becomes statistically stationary after
an initial transient, by adequately setting out
ow back pressure (seex2.5).
The relevant dimensionless parameters of STI are the mean
ow (or shockwave) Mach
number, M =e u
1;u
=e c, the turbulence Mach number, M
t
=u
0
rms
=e c, and the Taylor microscale
Reynolds number,Re
= u
0
rms
= , where = [u
0
2
u
0
2
=(@u
2
=@x
2
)
2
]
1=2
is the Taylor microscale,
u
0
rms
= (u
0
i
u
0
i
=3)
1=2
is the root-mean square (rms) velocity
uctuation in 3D simulations,
u
0
rms
= (u
0
i
u
0
i
=2)
1=2
is the rms velocity
uctuation in 2D simulations, and c is the local
speed of sound. A second Reynolds number can be considered, Re
L
= u
0
rms
L
=, based
on the dissipation-based length scale characterizing the large eddies, L
= K
3=2
=, where
K =u
00
i
u
00
i
=2 is the turbulence kinetic energy (TKE), and = 2S
ij
S
ij
is the dissipation rate
of TKE. The diusion of each passive scalar or species introduces an additional dimensionless
parameter, given by the Schmidt number, Sc
=D
=.
17
2.1.1 Constitutive Relations for Passive Scalar Mixing
For the study of passive scalar mixing, we assume that the gas is calorically perfect, im-
plying that the internal energy per unit mass can be expressed as e = c
v
T , where T is
the temperature and c
v
is the specic heat capacity at constant volume, independent of T .
Fourier's law of heat conduction is used to express the heat
ux as q
j
=@T=@x
j
, where
is the thermal conductivity, chosen such that the Prandtl number Pr =c
p
= equals 0:7,
withc
p
the specic heat capacity at constant pressure, also independent of T . The dynamic
viscosity follows the power law =
ref
(T=T
ref
)
3=4
, where the subscript `ref' indicates ref-
erence values. A zero volumetric viscosity has been assumed in equation (2.5). For passive
scalar simulations, the Schmidt number for each scalar is assumed to be constant. In later
sections, instead of Y
,
is used to denote -th passive scalar. Furthermore, the diusive
ux of passive scalar follows Fick's law, J
;j
=D
@
@x
j
. The Schmidt number is dened as
Sc
= =D
, which is constant for each scalar eld for passive scalar mixing simulations.
Additionally, there is no chemical reaction in passive scalar mixing, so the reaction term in
equation (2.4) is removed for passive scalars. Hence, the transport equation (2.4) for passive
scalars can be expressed as below:
@
@t
+
@
u
j
@x
j
=
@
@x
j
D
@
@x
j
: (2.6)
2.1.2 Constitutive Relations for Multi-
uid Mixing and Combustion
For multi-
uid mixing and combustion,e =hp= is mixture internal energy per unit mass,
and h is the mixture enthalpy:
h =
Ns
X
=1
Y
h
; h
=h
0
+
Z
T
T
0
c
p;
dT; (2.7)
c
p
=
Ns
X
=1
Y
c
p;
; c
v
=
Ns
X
=1
Y
c
v;
; (2.8)
18
where h
is the enthalpy of species , h
0
is the enthalpy of formation of species at
temperature T
0
. c
p
and c
v
are the isobaric and isochoric heat capacities, respectively, which
are now functions ofT . N
s
is the total number of species in the gas mixture. While no longer
assumed as calorically perfect, the mixture here still follows the ideal gas law, p =R
mix
T ,
where R
mix
is the mixture gas constant:
R
mix
=R
u
=W
mix
; (2.9)
whereR
u
is the universal gas constant and W
mix
is the mixture molecular weight, given by:
W
mix
=
Ns
X
=1
Y
=W
!
1
=
Ns
X
=1
X
W
: (2.10)
Species mass fraction Y
and mole fraction X
are related by the following expression:
Y
X
=
W
W
mix
: (2.11)
In this work, only N
s
1 species transport equations are solved, instead of N
s
. The mass
fraction of the last species is obtained by:
Y
Ns
= 1
Ns1
X
i=1
Y
: (2.12)
The heat
ux vector can be written as:
q
j
=
@T
@x
j
| {z }
1
+
Ns
X
=1
J
j
h
| {z }
2
+
Ns
X
=1
d
j
p
Y
D
T
| {z }
3
; (2.13)
whereJ
j
andd
;j
represent the diusive
ux and diusion driving force of species in the
j-th direction, respectively. D
T
is the thermodiusion coecient. In equation (2.13), term 1
represents Fourier's conduction law, term 2 represents the energy
ux brought by species
19
ux, and term 3 is the heat
ux induced by species concentration gradient, also known as
thermophoresis or Soret eect.
J
j
and d
;j
are related through the equations below:
J
i
=Y
V
i
; (2.14)
V
i
=
1
X
Ns
X
=1
Y
X
D
d
i
D
T
Y
@lnT
@x
i
| {z }
4
; (2.15)
d
i
=
@X
@x
i
+ (X
Y
)
@lnp
@x
i
| {z }
5
: (2.16)
In equation (2.14),V
;i
is diusive velocity. In equation (2.15), term 4 represents Dufour
eect, which characterizes heat transfer brought by species concentration gradient. Term 5
in equation (2.16) represents baroclinic diusion. Neglecting Dufour and Soret eects (Billet
et al., 2008), V
i
can be approximated as (Kee et al., 2005):
V
i
=
D
mix,X
X
d
i
=
D
mix,X
X
@X
@x
i
D
mix,X
X
(X
Y
)
@lnp
@x
i
; (2.17)
where D
mix,X
is the mixture-averaged diusion coecient, which can be expressed as (Kee
et al., 2005):
D
mix,X
=
1X
P
6=
X
=D
: (2.18)
The above expression is called Hirschfelder and Curtiss approximation. With this approx-
imation, additional steps should be taken to ensure mass conservation. The rst approach
is to follow equation (2.12) and make the last species to be dilutant. Then the dilutant can
absorb all the inconsistencies brought by the Hirschfelder and Curtiss approximation. The
second approach is to add a correction velocity to the diusive velocityV
c
in equation (2.17):
V
i
=
D
mix,X
X
d
i
+V
c;i
: (2.19)
20
According to Kee et al. (2005), the formulation of V
c;i
has a formulation as below:
V
c;i
=
Ns
X
=1
D
mix;X
W
W
mix
@X
@x
i
: (2.20)
In equation (2.17), the gradient of mole fraction needs to be calculated. In our compu-
tational implementation of these equations, mass fractions are stored in the memory instead
of mole fractions, which makes it easier to calculate gradients of mass fractions compared
with those of mole fraction. Hence, (2.17) can be rewritten in terms of the gradient of the
mass fraction as:
V
i
=
D
mix,Y
Y
@Y
@x
i
D
mix,X
X
(X
Y
)
@lnp
@x
i
; (2.21)
where D
mix,Y
is the mixture-averaged mass diusion coecient, which can be expressed as
(Kee et al., 2005):
1
D
mix,Y
=
X
6=
X
D
+
X
1Y
X
6=
Y
D
: (2.22)
The thermal conductivity of the mixture
mix
is computed from the following mixture
rule:
mix
=
1
2
Ns
X
=1
X
+
1
P
Ns
=1
X
=
!
: (2.23)
The mixture diusivity
mix
is calculated following Wilke (1950):
mix
=
Ns
X
=1
X
P
Ns
=1
;
X
; (2.24)
;
=
"
1 +
r
q
W
W
#
2
2
p
2
p
1 +W
=W
: (2.25)
21
2.2 Simulation Setup
The main physical conguration under investigation in the present study is that of the
canonical shock turbulence interaction (STI), in which homogeneous isotropic turbulence
is processed by a nominally planar shock wave that is in turn corrugated by the action of
the turbulence. The reference frame is chosen such that the shock becomes statistically
stationary in the direction of the incoming turbulence and statistically homogeneous in the
transverse directions (by imposing periodicity to approximate an innitely long domain in
those directions). Precursor simulations of homogeneous isotropic turbulence (HIT) are
conducted to provide the necessary turbulent in
ow to the STI simulations.
The computational domain for all 3D STI simulations is a rectangular box of size 4
2 2 (see gure 2.1). The governing equations are discretized on rectilinear, Cartesian
grids, uniformly spaced in the transverse directions, and stretched in the mean streamwise
direction (x = x
1
) to increase the resolution near the average shock location. For all simu-
lation cases, x
2;3
=x
1;s
= 2:4, where x
2;3
is the grid spacing in the transverse (y and z)
directions, and x
1;s
is the grid spacing in the streamwise direction immediately after the
shock. Proles of the streamwise grid spacing stretching used for dierent simulations are
shown in gure 3.1(a).
In 2D STI simulations, the domain size in the transverse direction is the same as the size
for 2D turbulence in
ow, which is determined with the approach presented in section 2.4.2.
The domain length in the streamwise direction is discussed in section 2.6.2. Dierent from
3D STI simulations, the third dimension (z) has only one grid cell. To make sure the velocity
in the third dimension is always zero, a symmetry boundary condition is adopted. In the
transverse direction (y), the boundary condition is periodic. The boundary condition for the
streamwise direction (x) is the same as 3D STI simulations, and grid spacing is uniform in
both transverse and streamwise directions.
22
shock
wave
x
y
z
blended
region
blended
region
HIT
1
HIT
2
HIT
3
2
2
2 2 2
HIT in
ow database
in
ow
2
4, STI domain
sponge layer
Figure 2.1: Problem setup of STI and blended HIT precursor simulations.
2.3 Numerical Methods
2.3.1 Spatial Discretization
The governing equations are solved on structured Cartesian grids 3D for passive scalar mix-
ing simulations and 2D for reacting simulations. A solution-adaptive numerical method is
employed (Larsson & Gustafsson, 2008; Larsson & Lele, 2009) that combines a fth-order
weighted essentially non-oscillatory (WENO) scheme with HLL
ux splitting near shocks,
and a sixth-order central dierence scheme in the split form of Ducros et al. (2000) elsewhere.
Numerical stability and conservation across the interfaces introduced between WENO and
central dierencing schemes by the solution-adaptive approach are ensured following Piroz-
zoli (2002). Discontinuities in thermodynamic quantities are detected with the criterion
> 1:2!
k
!
k
, based on the relative importance of the rate of compression (quantied by
negative dilatation, = S
kk
= @u
k
=@x
k
) compared with the rate of rotation (quantied
by enstrophy, !
k
!
k
, and !
i
=
ijk
@u
k
=@x
j
is the vorticity). Besides discontinuities in ther-
modynamic quantities, we also need to detect discontinuities in species mass fraction and
passive scalar. These discontinuities are detected with a WENO-based sensor, dened with
smoothness
r;d
calculated following Shu (1998), where the subscript r represents the r-th
23
stencil, andd represents thed-th direction (d =x;y;z). 5th order WENO is adopted in this
work, so there are three stencils (r = 1; 2; 3). According to this WENO-based sensor, the
WENO scheme is also adopted when
S
WENO
= max
d
max
r
r;d
min
r
r;d
> 1:01: (2.26)
In
ow and out
ow boundaries along the streamwise direction use summation by parts (SBP)
operators with simultaneous approximation terms (SAT) (Fern andez, Hicken & Zingg, 2014;
Sv ard & Nordstr om, 2014). Periodic boundary conditions are imposed in the transverse
direction(s).
2.3.2 Time Marching Scheme
For passive scalar mixing, the equations are advanced in time using a fourth-order, four-stage
explicit Runge-Kutta scheme. One major challenge in combustion applications is how to deal
with the sti chemistry term in the time integration. One common approach to address this
challenge is to use Strang splitting (Strikwerda, 2004), which separates the integration of
the sti (reaction) and non-sti transport) parts. However, Lu et al. (2017) proved that
this separation will lead to early near-limit phenomena (ignition and extinction), as will
be shown in section 3.4.2. They suggested that the sti and non-sti parts should not be
separated completely, and proposed a semi-implicit method, where they treated the non-sti
part explicitly, but the sti part implicitly. Inspired by Lu et al. (2017), we propose a fully
explicit method for time integration, which ts well into the four-stage, 4th-order Runge-
Kutta method (RK4) in use. To illustrate more details of our new time integration method,
we rst introduce the RK4 algorithm. The governing equations can be rewritten in a more
compact form as:
@F
@t
=T (F ) +R(F ); (2.27)
24
F =
2
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
4
u
v
w
h
Y
1
.
.
.
Y
Ns1
3
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
5
;R =
2
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
4
0
0
0
0
0
W
1
!
1
.
.
.
W
Ns1
!
Ns1
3
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
5
; (2.28)
whereR represents the reaction (sti) part of the governing equations, which is nonzero only
for the species transport equations,T represents the transport (non-sti) part of the govern-
ing equations, which includes convective and diusive/viscous terms. The RK4 algorithm is
illustrated in Algorithm 1 for the governing equations with no reactions:
Algorithm 1 RK4
F
state
n
=F
n
b =f
1
2
;
1
2
; 1; 1g, g =f
1
6
;
1
3
;
1
3
;
1
6
g
for i = 0 to 3 do
f
f
= tT (F
state
n
)
F
state
n
=F
n
+f
f
b(i)
F
n+1
=F
n+1
+g(i)f
f
Update primitive variables based on F
state
n
end for
Update primitive variables based on F
n+1
Before introducing RK4 with chemical reactions (both with Strang splitting and the newly
proposed method), F is split into two parts, the
uid mechanics part, F
F
, and species part,
F
S
. Similarly,T (F ) is split intoT (F
F
) andT (F
S
), whereasR(F ) is split intoR(F
F
) and
R(F
S
).R(F
F
) is a zero vector, so it will be dropped from the discussion below. Algorithm 2
illustrates how Strang splitting works. For each time step, the sti chemistry integrator rst
integrates
@F
S
@t
=R(F
S
) over half of the time step. This zero-dimensional (0D) integration
conserves mass and makes sure each cell is independent from each other, so the density will
not be changed within each cell during the integration. However, the pressure, temperature
25
and mass fractions will change. Then the second step is to integrate
@F
@t
= T (F ) with
RK4. At the beginning of each stage, all the primitive variables are updated based on
an intermediate state (F
state
F
and F
state
S
) calculated in each stage. The last step of this
algorithm is to integrate
@F
S
@t
=R(F
S
) for the other half time step. As discussed in Lu et al.
(2017), in this algorithm, the integration of the reaction term and the transport term is
totally separated. This is why the name of the integrator in Algorithm 2 is IntegratorR,
which means this is the integrator only for the reaction term. As discussed in (Strikwerda,
2004; Lu et al., 2017), Strang splitting is only second order accurate in time, so RK4 is not
necessary for Strang splitting. However, to compare with the semi-coupled explicit method,
which is intertwined with RK4, we still combine Strang splitting and RK4 in Algorithm 2.
Dierent from Strang splitting in Algorithm 2, the name of the sti chemistry integrator
Algorithm 2 Strang Splitting
Reaction Step 1
f
s
=
1
2
tR(F
S;n
) = IntegratorR(
1
2
t;F
S;n
)
F
S;n
=F
S;n
+f
s
update mass fractions, pressure and temperature
Transport Step
b =f
1
2
;
1
2
; 1; 1g, g =f
1
6
;
1
3
;
1
3
;
1
6
g
F
state
F;n
=F
F;n
, F
state
S;n
=F
S;n
F
S;n+1
=F
S;n
, F
F;n+1
=F
F;n
for i = 0 to 3 do
f
f
= tT (F
state
F;n
)
F
state
F;n
=F
F;n
+f
f
b(i)
F
F;n+1
=F
F;n+1
+g(i)f
f
f
s
= tT (F
state
S;n
)
F
state
S;n
=F
S;n
+f
f
b(i)
F
S;n+1
=F
S;n+1
+g(i)f
s
Update primitive variables based on F
state
F;n
and F
state
S;n
end for
Update primitive variables based on F
F;n+1
and F
S;n+1
Reaction Step 2
f
s
=
1
2
tR(F
S;n+1
) = IntegratorR(
1
2
t;F
S;n+1
)
F
S;n+1
=F
S;n+1
+f
s
Update mass fractions, pressure and temperature
26
is IntegratorRT, which means both the reaction and transport terms are integrated at the
same time. This approach to coupling the time integration of the reaction and transport
terms will be shown later to solve the early ignition and distinction brought by Strang
splitting. We name this new time integration method \explicit semi-coupled" method, whose
algorithm is shown in Algorithm 3. This algorithm is fully coupled with RK4. The only
dierence is that
@F
S
@t
=T (F
S
) +R(F
S
) is solved with a sti chemistry ODE integrator in
each stage of RK4. In the sti chemistry integrator, the time step determined according to
uid mechanics is divided into many substeps. In each substep,R(F
S
) is updated, butT (F
S
)
is not updated in order to keep the zero-dimensional feature of the chemistry integrator. This
is also why this method is called semi-coupled.
Algorithm 3 Explicit Semi-coupled Method
F
S;n+1
=F
S;n
, F
F;n+1
=F
F;n
b =f
1
2
;
1
2
; 1; 1g, g =f
1
6
;
1
3
;
1
3
;
1
6
g
for i = 0 to 3 do
f
f
= tT (F
state
F;n
)
F
state
F;n
=F
F;n
+f
f
b(i)
F
F;n+1
=F
F;n+1
+g(i)f
f
f
s
= t(T (F
state
S;n
) +R(F
state
S;n
)) = IntegratorRT(t;F
state
S;n
)
F
state
S;n
=F
S;n
+f
f
b(i)
F
S;n+1
=F
S;n+1
+g(i)f
s
Update primitive variables based on F
state
F;n
and F
state
S;n
end for
Update primitive variables based on F
F;n+1
and F
S;n+1
2.4 In
ow turbulence
2.4.1 3D In
ow Turbulence
Precursor simulations of decaying homogeneous isotropic turbulence (DHIT) with transport
of passive scalars are performed in a periodic box of size (2)
3
to provide turbulence databases
that will be superimposed at the in
ow plane of the STI simulations to the mean advection
velocity by means of Taylor's hypothesis. The uniform, isotropic grid spacing in each DHIT
27
simulation is dictated by the transverse resolution of the corresponding STI simulation that
will use the in
ow dataset.
DHIT simulations are initialized following the method of Ristorcelli & Blaisdell (1997),
extended to include passive scalars. The initial spectra for both velocity and passive scalar
elds are set to:
E(k) =
16
k
5
0
r
2
k
4
e
2k
2
=k
2
0
; (2.29)
where k
0
= 4 is the wavenumber at which the energy peaks. Phases for Fourier coecients
of passive scalars are randomized, and a unitary rms is ensured before inverse-Fourier trans-
forming to physical space. Hence, the resulting initial passive scalar elds have zero mean
and unitary variance.
DHIT simulations are run until reaching the targeted M
t
and Re
, which account for
the additional streamwise decay that takes place in the STI simulation from the in
ow
boundary plane to the reference state immediately upstream of the unsteady shock re-
gion. The turbulence is assumed fully developed once the velocity derivative skewness
(S = (@u
2
=@x
2
)
3
=
(@u
2
=@x
2
)
2
3=2
) stabilizes around0:5 and the dissipation-based length
scale, L
, increases with time.
To obtain an in
ow database that is suciently long to allow for the computation of
statistically decorrelated time averages in the STI simulation, snapshots from multiple re-
alizations of DHIT, with the same physical parameters (M
t
, Re
, and Sc
) but dierent
initial conditions, are blended following Larsson (2009). Each pair of passive scalar elds
with zero mean, resulting from two dierent DHIT realizations (A and B), are blended as:
=
A
cos +
B
sin, with =
4
(1 + sin (x
b
=2l
b
)), where l
b
is the blending length and
jx
b
j < l
b
is the blending region, centered at x
b
= 0. In addition, from each DHIT snap-
shot (box) two additional boxes, obtained by a series of rotations, are subsequently blended.
From the original box, A, with axes denoted as (x;y;z) and velocity components (u;v;w),
the second box, B, is obtained by two consecutive =2 rotations around the y and z axes,
respectively, such that (x
B
;y
B
;z
B
) = (y;z;x) and (u
B
;v
B
;w
B
) = (v;w;u). Starting from
28
Figure 2.2: Schematic representation of the variation of energy and enstrophy with time,
following Davidson (2015).
the second box, the same rotations are applied again to obtain a third box of DHIT,C, such
that (x
C
;y
C
;z
C
) = (z;x;y) and (u
C
;v
C
;w
C
) = (w;u;v). Then these three DHIT boxes are
blended as the input turbulence database for STI simulations.
2.4.2 2D In
ow Turbulence
2D DHIT datasets are required as in
ow for the 2D STI simulations. 2D turbulence presents
signicant dierences with respect to 3D turbulence. One distinctive characteristic of 2D
turbulence is that vortex stretching is absent. As a consequence, in 2D turbulence, certain
types of vortical structures retain their general shape and intensity for time periods much
longer than the characteristic eddy turnover time. These vortices are called coherent vor-
tices (McWilliams, 1990; Fox & Davidson, 2010). After the initial transient period, eddies of
integral length scale increases in time, which leads to an increasedRe
with time (Davidson,
2015). These dierences make the criteria used to assess that the 2D turbulence is fully devel-
oped are dierent from the 3D case. As shown in Fox & Davidson (2010); Davidson (2015),
after initial transient period, there will be an accelerated decay of enstrophy (see gure 2.2).
For the 2D turbulence simulations presented in this work, fully-developed turbulence with
targeted M
t
, Re
and accelerated decay of enstrophy is fed into 2D STI simulations as in-
ow. The initialization of 2D turbulence follows the same procedure described in 2.4.1, but
29
with no third (transverse) dimension. The initial spectra for both velocity and passive scalar
elds are set to:
E(k) =k
3
exp(3k
2
=2k
2
0
); (2.30)
wherek
0
= 4 is the wavenumber at which the energy peaks. Then the rms velocity and pas-
sive scalar or species concentration variance are normalized to specied values. Additionally,
the species concentration generated with this approach has a mean of zero. For multi-
uid
mixing and combustion applications, the mean needs to be oset to be the value specied
for each species. In the present computational implementation, the mean and variance for
species concentration are specied for onlyN
s
1 species. The species concentration for the
last species is calculated with equation (2.12).
The approach above is similar to the initialization for passive scalar, so the thermody-
namic quantities generated are dimensionless, and the domain size is 2 2. The thermo-
dynamic quantities and domain size are then dimensionalized keeping M
t
andRe
constant
are picked to be constant during the process. In the transformation below, the subscript
`new' means dimensional quantities after the transformation, and subscript `old' means di-
mensionless quantities before the process. Before the transformation, T
old
= 1,
old
= 1,
and the gas constant is R
old
. According to the ideal gas law, p
old
=
old
R
old
T
old
= R
old
.
The transformation is dened by setting the target values T
new
and p
new
. Once the target
T
new
andp
new
are set, and preserving the original mean species concentrations, the resulting
mean gas constant,R
new
, and density
new
, of the mixture can be determined. AsM
t
is kept
constant throughout the transformation, we have:
M
2
t
=
hu
2
old
i=2
R
old
T
old
=
hu
2
new
i=2
R
new
T
new
; (2.31)
u
new
=u
old
s
R
new
T
new
R
old
T
old
: (2.32)
30
In this process, the specic heat capacity ratio,
, is the same before and after the trans-
formation, and it is obtained from T
new
, p
new
and the specied mean species concentration.
Angular bracketshi indicate spatial averaging. over the whole HIT domain.
The following equation is used to determine the domain size after the transformation:
L
old
old
=m
L
new
new
; (2.33)
where is Taylor micro scale, and L
old
= 2. m is a scaling factor that changes with k
0
.
In all the reactive
ow simulations here considered, m is set to be k
0
=4. In 2D reactive STI
simulations, the transverse direction averaging suers from a worse convergence compared
with 3D simulations because of the missing third dimension. One remedy is to increase
the amount of data in the transverse direction, keeping the other physical parameters of
the simulation unchanged. This can be done by increasing k
0
and transverse length at the
same time. For example, k
0
= 4 and k
0
= 8 indicate that there are four and eight large
eddies in the domain, respectively. If the transverse length of the simulation with k
0
= 8
is twice the transverse length of the simulation with k
0
= 4, then both simulations have
large eddies of the same size. By doubling k
0
, the transverse length and the number of grid
elements in the transverse direction, the amount of data available for spatial averaging in
the transverse direction can be eectively increased. The ratio between
old
and
new
are
obtained by assuming that Re
is not changed by the transformation:
Re
=
old
u
rms;old
old
old
=
new
u
rms;new
new
new
; (2.34)
new
old
=
u
rms;old
u
rms;new
rms;new
rms;old
rms;old
rms;new
(2.35)
where
new
is also determined from the speciedT
new
,p
new
and the unchanged mean species
concentrations of the mixture. The blending described in section 2.4.1 to create the in
ow
datasets for the STI simulations is also adopted for 2D simulations.
31
2.5 Sponge layer and out
ow boundary condition
A sponge layer is used in the STI simulations to eectively damp acoustic re
ections gener-
ated at the out
ow boundary, preventing upstream propagation into the subsonic region of
the domain. Each governing equation (2.1))-((2.6) adds a source term of the form:
k
0
u
1;u
L
y
xx
sp
x
max
x
sp
2
(fhfi
trans
(x)): (2.36)
For simulations of passive scalar mixing, L
y
= 2, x2 [x
sp
;x
max
] = [3; 4]. For multi-
uid
mixing, L is the domain length in the y direction, and the length for sponge layer is the
last quarter of the domain. In equation (2.36), u
1;u
is the mean advection velocity in the
streamwise direction upstream of the shock, andhfi
trans
represents instantaneous spatial
average in each transverse yz-plane for 3D simulations, and in each transverse y-line for 2D
STI, with f =f;u
i
;E;
(Y
)g.
To obtain a statistically stationary shock, the back pressure at the out
ow plane, p
b
,
is specied following Larsson & Lele (2009). Starting from the laminar Rankine-Hugoniot
pressure jump condition, successive adjustments to the back pressure are applied to minimize
the mean shock-drift velocity,U
s
, according to: p
b
4
u
e u
1;u
U
s
=(
+1), where the subscript
u denotes values upstream of the shock.
= c
p
=c
v
is the ratio of specic heat capacities
at constant pressure and volume, respectively, taken as 1:4 in this work for passive scalar
mixing. For the reactive 2D STI simulations, the value of
in cells within the sponge layer is
determined on-the-
y by spatially averaging the local
at each streamwise location within
that sponge layer.
32
Case Re
Re
L
M M
t
L
= Resolution Regime Line style Symbol
A 41 213 1.28 0.31 65 1280 384
2
broken
B 37 177 1.5 0.09 56 1700 512
2
wrinkled
C 38 191 1.5 0.20 60 1280 384
2
wrinkled
D 39 193 1.5 0.31 60 1280 384
2
wrinkled
E 40 201 1.5 0.42 62 1280 384
2
broken
F 38 177 2.0 0.30 57 1700 512
2
wrinkled
G 37 191 3.0 0.30 60 1280 384
2
wrinkled
H 37 204 5.0 0.31 63 1280 384
2
wrinkled
I 74 679 1.5 0.42 154 2400 1024
2
broken
J 66 548 3.0 0.29 132 2400 1024
2
wrinkled
K 72 655 5.0 0.31 151 2400 1024
2
wrinkled
Table 2.1: 3D passive scalar mixing simulation cases, line styles and symbols (for postshock
state) used in gures.
2.6 Simulation cases
2.6.1 Passive scalar mixing in STI
Table 2.1 presents a summary of the 3D simulations of passive scalar mixing under STI
performed in this study, varyingRe
between 40 (casesA toH) and 70 (I toK),M from
1:28 to 5:0, andM
t
from 0:09 to 0:42, which covers a wider parameter space than previously
investigated (Tian et al., 2017; Boukharfane et al., 2018). Wrinkled and broken-shock regimes
are thus considered in these simulations, according to the criterion that M
t
=(M 1)' 0:6
is required for the broken-shock regime (see Donzis, 2012b; Larsson et al., 2013). Simulation
casesA-H include passive scalars withSc = 0:5, 1, and 2, whereas casesI-K considerSc = 1
only.
The initial transient period needed to obtain a statistically stationary state (after which
the shock drift is negligible) is t
trans
3L
input
=e u
1;u
, or, in dimensionless form, k
0
e u
1;u
t
trans
=
72, whereL
input
is the length of the in
ow turbulence database obtained by blending snap-
shots from periodic (2)
3
DHIT boxes along the streamwise direction, ande u
1;u
is the mean
streamwise velocity upstream of the shock. After the initial transient, statistics are collected
over a time period L
input
=e u
1;u
, or k
0
e u
1;u
t
stats
= 24 in dimensionless form.
33
Case Re
;0
M M
t;0
Resolution Domain Size T
1
(K) p
1
(kPa) Line Style
L 1250 2 0.3 5120 1024 12 cm2.26 cm 560 60
M 1250 2 0.1 5120 1024 12 cm6.81 cm 560 60
N 600 2 0.3 5120 1024 12 cm2.12 cm 560 60
O 600 2 0.1 5120 1024 12 cm6.38 cm 560 60
P 0 2 0 51200 1 12 cm 560 30
Table 2.2: 2D Reactive STI simulation cases, line styles used in gures. Please note that the
resolution and domain size are for the rened part before the sponge layer.
2.6.2 Reactive 2D STI
Table 2.2 shows the ve cases to be presented in this section. From caseL to caseO, the mean
mole composition in the in
ow database is 0:2H
2
: O
2
: 7Ar, and the chemical mechanism is
the H
2
=O
2
reaction subset from GRI 3.0 (Smith et al., 1999). The species mass fractions are
initialized with equation (2.30) and choosing the rms
uctuation as 10% of the mean mass
fraction for each species. These four simulations cover two initial Re
= 1250; 600 and two
initialM
t
= 0:1; 0:3 at the inlet boundary. The domain length in the streamwise direction is
50 times larger than in the transverse length for each case. Uniform grid spacing is adopted
in the region before the sponge layer. In the sponge layer , the grid is stretched with a factor
of 1.1 in the streamwise direction. Please note that the resolution and domain size in table
2.2 do not include the sponge layer.
The STI domain is initialized without fuel (i.e., with only O
2
and dilutant). The rst 10
points in streamwise direction are blended between the STI initial condition and the input
HIT database to avoid discontinuities. The initial condition in the STI domain is generated
with M = 2 and an initial composition of O
2
: 7Ar. The initial composition is dierent
between the STI domain and the in
ow database, with the purpose to mimic the process of
the Ar-diluted fuel/oxidizer mixture as it is processed by the shock. In all the simulations,
the streamwise length is 1.2m. The shock is placed 3 cm behind the inlet boundary, and the
sponge layer starts 12 cm behind the inlet boundary. In the initial STI domain, the velocity
of the mixture before the shock is 875 m=s, and the velocity of the mixture after the shock is
367 m=s. The induction length, dened as the distance between the shock and the location
34
where ignition initiates, is approximately 3 cm when the Ar-diluted fuel/oxidizer mixture
is ignited after the shock, and is obtained by multiplying the velocity after the shock by
the ignition delay time obtained from 0D simulation under postshock conditions (1103 K,
283 kPa). Three ensemble runs are performed for each turbulent case. The computational
cost for caseL is 36864 core hours with Intel SKylake processors. If an additional transverse
direction with the same number of cells (1024) were considered, the resulting 3D counterpart
simulations of caseL would take approximately 38 million core hours (1024 36864), which
is quite expensive for the present study.
35
Chapter 3
Verication
The numerical simulations of the present study are performed with the HybridX research
code, which builds on the Hybrid code originally developed by Johan Larsson and exten-
sively veried and validated for STI studies in prior work (Larsson & Lele, 2009; Johnsen
et al., 2010; Larsson et al., 2013; Bermejo-Moreno et al., 2013a). The code is written in
C++, and parallelized through the Message Passing Interface (MPI), emphasizing the use
of nonblocking communication with concurrent computation. In this work, capabilities to
simulate mixing and combustion applications are added. Thermodynamic properties of ther-
mally perfect gases and reaction rates in the combustion are obtained by interfacing with
Cantera (David et al., 2018). CVODE (Cohen et al., 1996) is used as the stiif ODE solver for
the reaction terms.
This chapter focuses on the verication and validation of the implemented simulation
methodology used in this work. There are four sections in this chapter. Section 3.1 shows
the grid convergence for dierent quantities in dierent 3D STI simulations. Section 3.2
addresses the capability of the code to simulate 2D HIT. Good agreement with Fox &
Davidson (2010) is achieved. The similarity between 2D STI and 3D STI is assessed in
section 3.3, which conrms that 2D STI can re
ect certain relevant aspects of the physics
found in 3D STI, and provides a justication for using 2D STI simulations as a rst step in
the study of reactive STI conguration. Section 3.4 assesses the validity of the combustion
36
0:0 2:5 5:0 7:5 10:0
x
0:000
0:005
0:010
0:015
0:020
0:025
0:030
x
A
D
B
C
D
E
I
F
G
J
H
K
(a)
1 0 1 2 3
x=L
;u
1:0
1:5
2:0
2:5
3:0
3:5
4:0
4:5
k
max
(b)
Figure 3.1: (a) Streamwise grid spacing and (b) k
max
proles for all the cases in table 2.1
.
capability comparing results of our implemented methodology with relevant prior numerical
studies found in the literature.
3.1 Grid convergence for 3D STI
The grids used in the present simulations are stretched such that the resolution is highest
near the shock, where the Kolmogorov scale is the smallest (Larsson & Lele, 2009). Figure 3.1
shows proles of the grid spacing in the mean shock-normal direction for dierent simulations
along with proles of k
max
, exceeding unity at all locations. Grid convergence of statistics
is assessed by comparing results from series of simulations run at increasing grid resolution
for dierent combinations of dimensionless parameters. Figure 3.2 shows Kolmogorov scale,
transverse vorticity variance, passive scalar variance (for Sc = 1) and scalar dissipation
for case (M;M
t
;Re
) = (2:0; 0:3; 38), comparing results obtained from four independent
simulations performed at increasing resolution: 426128
2
, 852256
2
, 1280384
2
, 1700512
2
(caseF on table 1). Dierences between the 1280 384
2
and 1700 512
2
grids are less than
1% for all statistics, including background turbulent
ow quantities and those of the passive
scalar. Thus, the 1280 384
2
simulations are deemed grid-converged. A similar analysis is
37
2 0 2 4 6 8
x=L
;u
0
1
2
3
4
5
k
max
(a)
2 0 2 4 6 8
x=L
;u
0
1
2
3
4
5
6
!
02
y
=!
02
y;u
F (1700 512
2
)
F (1280 384
2
)
F (852 256
2
)
F (426 128
2
)
(b)
2 0 2 4 6 8
x=L
;u
0:00
0:25
0:50
0:75
1:00
1:25
1:50
f
002
=
f
002
u
(c)
2 0 2 4 6 8
x=L
;u
0:0
0:5
1:0
1:5
2:0
2:5
^ = ^
u
(d)
Figure 3.2: Grid convergence test for case F , with (M;M
t
;Re
) = (2:0; 0:3; 40) on increas-
ingly ner grids: 426 128
2
(dotted), 852 256
2
(dash-dotted), 1280 384
2
(dashed),
1700 512
2
(solid). (a)k
max
, (b) transverse vorticity variance!
0
2
!
0
2
. (c) scalar variance,
f
002
(d) scalar dissipation rate, ^ . Quantities in (b)-(d) are normalized with their corresponding
value immediately upstream of the shock.
38
presented for simulations at the highest Reynolds number considered in this study. Figure 3.3
shows results for simulations with (M;M
t
;Re
) = (5:0; 0:3; 72) at resolutions 600 256
2
,
1200 512
2
, 1800 768
2
, and 2400 1024
2
in which the convergence of statistics can also be
observed. The eects of using a nite computational box size and varying the wavenumber
at which the energy peaks, k
0
, were discussed in Larsson, Bermejo-Moreno & Lele (2013).
In the present simulations we choose the same box size, 4 (2)
2
, and peak-energy mode
k
0
= 4, which were shown to have a minimal in
uence on the simulation results. Please note
that in reactive STI simulations, the domain sizes are dierent, which are converted to be
dimensional values based on the method proposed in section 2.4.2.
Another relevant metric to examine the convergence of averaged simulation results is the
dierence between the sum of all terms on the left and right hand sides of each transport
equation presented inx4.1. Figure 3.4(a) shows the streamwise proles of non-negligible
terms in the scalar variance transport equation as well as the total contribution from the
left- and right-hand sides, for a representative case at (M;M
t
;Re
) = (1:28; 0:31; 41). Fig-
ure 3.4(b) shows equivalent results for the scalar dissipation transport equation. In both
cases, the dierence between each side of the equation is negligible. However, convergence of
the left and right hand sides of the scalar dissipation transport equation generally requires an
additional level of renement compared with the other statistics, as it involves higher-order
spatial derivatives (notice the 1700 512 512 grid resolution of the simulation used for
that gure).
3.2 Verication for 2D HIT simulations
Two-dimensional homogeneous isotropic turbulence has been studied for the past few decades
(Brachet et al., 1988; Fox & Davidson, 2010; Davidson, 2015; Zhou & Xu, 2020). In this
section, a 2D HIT simulation is conducted with the same nondimensional initial condition
as Fox & Davidson (2010) to verify the correctness of the code. Whereas the simulation in
39
2 0 2 4 6 8
x=L
;u
0
1
2
3
k
max
(a)
2 0 2 4 6 8
x=L
;u
0
2
4
6
8
10
12
!
02
y
=!
02
y;u
2400 1024
2
1800 768
2
1200 512
2
600 256
2
(b)
2 0 2 4 6 8
x=L
;u
0:00
0:25
0:50
0:75
1:00
1:25
f
002
=
f
002
u
(c)
2 0 2 4 6 8
x=L
;u
0
1
2
3
4
^ = ^
u
(d)
Figure 3.3: Grid convergence test for case K, with (M;M
t
;Re
) = (5:0; 0:3; 72) on increas-
ingly ner grids: 600 256
2
(dotted), 1200 512
2
(dash-dotted), 1800 768
2
(dashed),
2400 1024
2
(solid). (a) k
max
, (b) transverse vorticity variance !
0
2
!
0
2
. (c) scalar variance,
(d) scalar dissipation rate, ^ . Normalizations as in gure 3.2.
40
2 0 2 4 6 8
x=L
;u
1:00
0:75
0:50
0:25
0:00
Scalar variance budgets
L1
E
C
A
ABC +DEF
(a)
2 0 2 4 6 8
x=L
;u
500
250
0
250
500
Scalar dissipation budgets
I
H
G
L2
GH +I
(b)
Figure 3.4: Streamwise proles comparing left- and right-hand side terms of (a) scalar vari-
ance transport equation (3.2) for case J and (b) scalar dissipation transport equation 3.1)
for case A (at a higher resolution of 1700 512 512).
Fox & Davidson (2010) is in the incompressible
ow regime, all the simulations conducted
in this work account for compressibility.
In the present 2D HIT simulation, the initial spectra is given by equation (2.30) with the
peak energy at k
0
= 8. The dimensionless initial rms velocity is unity. The (2)
2
domain
has a resolution of 8192
2
grid points. As in Fox & Davidson (2010), an inverse time scale
!
0
=h!
2
i
1=2
is dened based on the initial enstrophy. The time scale in this simulation
is nondimensionalized as = t!
0
. In Fox & Davidson (2010), the kinematic viscosity
is determined with the initial Re = !
0
=(k
2
0
) = 29000. In the present compressible
ow
simulation, the reference dynamic viscosity
ref
is set to be
ref
=hi, wherehi = 1 is
the initial dimensionless mean density. Figure 3.5(a) shows the time evolution of energy and
enstrophy, together with the enstrophy decay with time in Fox & Davidson (2010), with
both quantities averaged over the whole domain. An initial transient period of 30 is
observed, followed by an exponential decay, as shown in Lowe & Davidson (2005), Fox &
Davidson (2010), and Davidson (2015). A region of E
!
(k) k
1
for enstrophy spectrum
has been found in many previous studies (Kraichnan, 1967; Batchelor, 1969; Dritschel et al.,
2007; Fox & Davidson, 2010). In gure 3.5(b), E
!
(k) is multiplied by k, so the k
1
region
41
10
2
10
0
10
2
10
0:4
10
0:2
10
0
(a)
Enstrophy
Energy
Fox et al. (2010)
10
0
10
1
10
2
10
3
k
0
5
10
15
kE
!
(b)
= 30
= 50
= 70
Figure 3.5: Proles of (a) time evolution of energy and enstrophy, compared with the en-
strophy decay in Fox & Davidson (2010) and (b) enstrophy spectra.
appears as a plateau. As the
ow develops and decays, the plateau becomes less obvious,
because a relatively large Re is needed to observe a wide k
1
region. A clear dissipation
range to the right of the plateau is also observed in gure 3.5(b).
3.3 Verication for 2D STI simulations
To show that the 2D STI setting can re
ect important elements of the physics of 3D STI,
the 2D version of simulation D (see table 2.1) is performed (and denoted as case L). Two-
and three-dimensional simulations (D and L) are then compared. The in
ow database of
the 2D simulation is blended from 100 independent 2D HIT planes with similar M
t
and
Re
. Time-averaged statistics are collected for the time it takes the mean upstream velocity
to go over the STI domain 100 times. The domain size is 24 2, and the grid size is
4200 512. The mesh is stretched in the streamwise direction similarly to how it is done
in the 3D STI simulations to ensure enough resolution near the shock region. The mean
shock location is placed 4 behind the inlet boundary, which is the distance for the 2D
turbulence to decay approximately to the targeted M
t
and Re
. Figure 3.6 compares the
change of thermodynamic quantities across the shock. The streamwise evolution of these
42
Figure 3.6: Streamwise proles of Reynolds averaged density (a), streamwise velocity (b),
pressure (c) and temperature (d). Each curve has been normalized by the value immediately
upstream of the unsteady shock region. Grey shaded regions and horizontal line segments
mark the extent of the unsteady shock region for each case.
43
Figure 3.7: Streamwise proles of streamwise Reynolds Stress (a), transverse Reynolds
stress (b), M
t
(c) and Re
(d), vorticity variance (e), Kolmogorov scale (f). Each curve
in (a); (b); (e); (f) has been normalized by the value immediately upstream of the unsteady
shock region. Grey shaded regions and horizontal line segments mark the extent of the
unsteady shock region for each case.
44
Figure 3.8: Streamwise proles of scalar variance (a), SDR (b). Each curve has been normal-
ized by the value immediately upstream of the unsteady shock region. Grey shaded regions
and horizontal line segments mark the extent of the unsteady shock region for each case.
quantities across the shock is very similar between 2D and 3D simulations, but the plots for
2D simulation are closer to the laminar shock jump relation.
Figure 3.7(a;b) shows the streamwise evolution of streamwise and transverse Reynolds
stresses. 3D and 2D results show similar trends. However, the amplication of transverse
Reynolds stress in 2D simulation across the shock is notably larger than that of 3D simu-
lation. The reason is not clear to the author. Further analysis of budgets of the transport
equation of Reynolds stress can shed light on the reason behind this observation, which is
not the focus of the current work. The change ofM
t
and Kolmogorov scale across the shock
are presented in gure 3.7(c;f), respectively. Both of them show similar trends between
2D and 3D simulations. The evolution of Re
is shown in gure 3.7(d), which re
ects a
dierence between 3D turbulence and 2D turbulence. In 2D turbulence, there is a criti-
cal Reynolds number, above which the Reynolds number increases with time, as proved by
numerical simulations by Chasnov (1997). Figure 3.7(e) shows the streamwise evolution of
vorticity variance. The dierence between 2D and 3D simulation in vorticity variance is
brought by the missing vortex stretching production in 2D simulations, which may explain
why the amplication of vorticity variance is smaller for 2D simulations.
45
The evolution of scalar variance and SDR across the shock are presented in gure 3.8. 2D
and 3D simulations show similar trends for both quantities. The increase of SDR is larger
for 3D simulation, consistent with the larger increase of vorticity variance across the shock.
Larger increase in vorticity variance means more energy is injected to small scales, so the
more energy in small scales in 3D STI can bring a better enhancement for turbulent mixing
across the shock compared with 2D STI. In summary, from the presented comparison between
2D and 3D STI simulations, we can conclude that the 2D STI simulations can qualitatively
re
ect relevant aspects of the physics of 3D STI simulations.
3.4 Verication for Combustion Capability
This section introduces several test cases intended to verify and assess the validity of our
implementation of the combustion solver. The rst case includes 0D, 1D and 2D simulations
to calculate ignition delay time of dierent reactive mixtures under varying conditions. The
0D simulation veries that the interface between Cantera and HybridX is correct. 1D and
2D simulations of autoignition at high pressure and high temperature are performed to assess
the accuracy of the newly propose time-integration method (introduced in section 2.3.2) to
capture ignition delay time. The second test case considers ignition induced by a re
ected
shock in a 1D shock tube. The third one is a 2D simulation of shock/hydrogen bubble
interaction. The results from our simulations for all test cases are compared with existing
literature (Fedkiw et al., 1997; Zikoski, 2011; Billet et al., 2008; Houim & Kuo, 2011).
3.4.1 Ignition Delay Time
In this part, 0D Cantera simulations are conducted rst. These simulations cover a range of
initial temperature (from 900 K to 1400 K) under a constant pressure of 1 atm. The chemical
mechanism adopted here is a subset of GRI Mech 3.0 (Smith et al., 1999), which involves
H
2
=O
2
mixtures diluted with Ar. Afterwards, 1D and 2D simulations are conducted with
46
Figure 3.9: Ignition delay time for 2H
2
: O
2
: 7Ar mixtures at dierent initial temperatures
and ambient pressure of 1 atm, comparing results from 0D simulations using Cantera with
1D and 2D ignition-spot simulations using HybridX.
ignition spot at dierent temperatures. In 1D simulations, the size of the domain is 10 cm.
The size of the ignition spot is 6 cm, and it is placed n the middle of the domain. In 2D
simulations, the size of the domain is 10 cm by 10 cm. A circular ignition spot with a radius
of 3 cm is centered in the domain. In both 1D and 2D simulations, the initial pressure of
the whole domain is 1 atm. The temperature outside of the ignition spot is 300 K, and the
temperature inside of the ignition spot varies from simulation to simulation. The initial mole
fractions for all the simulations are given by the ratios 2H
2
: O
2
: 7Ar. In both 1D and 2D
simulations, the resolution is 512 points in each dimension, found to provide grid-converged
results, and periodic boundary conditions are considered in all coordinate direction. For
each simulation, the ignition delay time is dened by the time instant with the fastest mean
temperature increase (T=t) in the entire computational domain between two consecutive
time steps. Comparisons between 0D simulations using Cantera and 1D, 2D simulations
47
Figure 3.10: Illustration of the simulation set up of 1D shock induced ignition
using HybridX are shown in gure 3.9. It can be observed that the ignition delay times for
1D and 2D simulations match really well with 0D simulations in the entire range of initial
temperature considered.
3.4.2 1D Shock Induced Ignition
As the incident shock produced in a shock tube containing a reactive mixture re
ects from
the end wall, the high temperature and high pressure created by the re
ected shock leads
to the formation of a reaction front, which eventually catches up and merges with the shock
wave, forming three waves. From the re
ective (end wall) boundary to the outlet boundary
considered at the other end (from right to left), the three waves formed are a rarefaction, a
contact discontinuity, and a detonation front. An illustration of the computational set up of
this simulation is shown in gure 3.10.
In this case, the 1D multi-species Navier-Stokes equations with chemical reactions are
solved. The computational domain is 12 cm long and is discretized into a grid containing
32000 points distributed evenly. The
uid is a stoichiometric mixture of hydrogen and
oxygen, diluted by 70% argon, i.e., 2H
2
: O
2
: 7Ar. The chemical mechanism from Westbrook
(1982) and Zikoski (2011) is used here. This mechanism includes 37 elementary reactions
involving 9 species: H; O; Ar; H
2
; O
2
; OH; H
2
O; HO
2
and H
2
O
2
. The shock is initially placed
at x = 6 cm, and moves to the right. Behind the shock (x 6 cm), the state is set to be
1
= 0:1875 kg=m
3
, u
1
= 487:34 m=s, p
1
= 35594 Pa, values that are also set at the in
ow
48
boundary condition on the left end of the domain. The right boundary is a re
ective wall.
Ahead of the shock (x> 6 cm), the
uid state is set to be
2
= 0:072 kg=m
3
, p
2
= 7173 Pa,
and u
2
= 0 m=s. Similar simulations are conducted in Fedkiw et al. (1997) and Zikoski
(2011). Data extracted from Zikoski (2011) is used here for comparison with the present
simulations.
Figures 3.11 and 3.12 show scaled plots of the density, velocity, pressure, and temperature
at simulation times of 180s and 230s, respectively. Each quantity is scaled by its post-
incident-shock state. The x axis is normalized by the total domain length. The results
match well with Zikoski (2011). The three waves (rarefaction, contact discontinuity, and
detonation front) can be observed in gure 3.12. The detonation front generated in HybridX
is slightly ahead of the reference. Two main factors may be responsible for this discrepancy.
In Zikoski (2011), the sti term and the non-sti term are updated together, so they are
fully coupled. In HybridX, when integrating the species transport equations with substeps,
the non-sti part of the equation is xed. Additionally, Zikoski (2011) considered thermal
diusion, which will diuse the peak of OH radical and delay ignition. Results obtained
from Strang splitting are also shown at t = 180s. It can be observed that the detonation
front catches up with the re
ected shock already, which indicates an early ignition.
3.4.3 Shock Hydrogen Bubble Interaction
This test case aims to assess the ability of the code to simulate the interaction between
shock wave and chemical reactions in 2D settings, by comparing with benchmark reference
data from prior simulations (Zikoski, 2011; Houim & Kuo, 2011; Billet et al., 2008). An
illustration of the simulation setup is present in gure 3.13. In this test case, a hydrogen
bubble placed in air passes through a shock wave, and the elevated temperature and pressure
after the shock will lead to ignition of the hydrogen bubble in air. There are two setups of
the simulation found in previous work. In Zikoski (2011), the bubble location is xed, and
the shock moves and passes through the hydrogen bubble. In Billet et al. (2008) and Houim
49
Figure 3.11: Instantaneous streamwise proles of
ow quantities at time 180s of 1D shock-
induced ignition simulations, comparing reference data by Zikoski (2011) with present sim-
ulation results using the newly proposed time-marching scheme, and Strang splitting. Each
quantity is normalized by the initial value behind the incident shock. The arrow indicates
the direction of detonation front propagation.
50
Figure 3.12: Instantaneous streamwise proles of
ow quantities at time 230s of 1D shock-
induced ignition simulations, comparing reference data by Zikoski (2011) with present sim-
ulation results using the newly proposed time-marching scheme, and Strang splitting. Each
quantity is normalized by the initial value behind the incident shock. The red arrow indi-
cates the direction of the detonation front propagation. The three letters mark A) detonation
front, B) contact discontinuity, C) rarefaction wave, respectively.
Figure 3.13: Illustration of the simulation set up of 2D shock hydrogen bubble interaction
51
& Kuo (2011), the shock location is xed, and the hydrogen bubble passes through the
shock. In the present simulation we replicate the latter setup. Because of the symmetry of
this problem, only half of the bubble is simulated, and both the top and bottom boundaries
apply symmetry boundary conditions. The size of the domain is 2 m 0:75 cm. The shock
is placed at x = 1:1 cm. Initially, the center of the hydrogen bubble is at y = 0 cm and
x = 0:8 cm. The radius of the bubble is 0:28 cm. The sponge layer starts at x = 40 cm. The
ow goes from left (in
ow) to right (out
ow) in the computational domain. Before the shock
(unshocked state),T
1
= 1000 K,p
1
= 1 atm,u
1
= 1240 m=s, which is also the state set at the
in
ow boundary condition. Behind the shock (shocked state), T
2
= 1565 K, p
2
= 4:45 atm,
u
2
= 435 m=s, which is also the state set at the out
ow boundary condition. There is no
velocity in the y direction initially. The grid spacing in the y direction is uniform, which is
10
3
cm. The grid is rened in the x direction from x = 0:5 cm to x = 1:8 cm, where the
grid spacing is also 10
3
cm. Outside of the rened domain, the grid spacing increases with
a ratio of 1:1 towards the in
ow and the out
ow boundary. In this simulation, the time step
is xed to be 10
9
s.
The chemical mechanism used here is the same as in Billet et al. (2008) and Zikoski
(2011), with 9 species (H; O; N
2
; H
2
; O
2
; OH; H
2
O; HO
2
and H
2
O
2
) and 19 reversible reactions.
Outside of the hydrogen bubble, the composition of air is 768N
2
: 232O
2
in mole fraction,
and inside of the bubble, the composition is pure H
2
, with the bubble interface dened by:
Y
H
2
= 0:5 + 0:5tanh
r
b
r
C
; (3.1)
Y
O
2
= 0:232(1Y
H
2
); (3.2)
Y
N
2
= 0:768(1Y
H
2
); (3.3)
wherer
b
is the radius of the bubble, which is 0:28 cm in this case. C represents the thickness
of the interface, which is 0:03 mm. r is the distance between a point in the computational
domain and the center of the bubble.
52
Contour maps of H
2
concentration with superimposed pressure contour lines are shown
in gures 3.14, 3.15, 3.16. Our results are compared with both Houim & Kuo (2011) and
Billet et al. (2008). Good agreements are obtained with both papers. However, Billet et al.
(2008) conducted simulation with a longer evolution time, so only comparison with Billet
et al. (2008) are also included. In each gure, there are 100 pressure contour lines, whose
maximum and minimum values follow Billet et al. (2008). Hydrogen concentration below
5% is masked. The reference data shown in gures 3.14-3.16 has been digitized from the
contour plots provided in Billet et al. (2008) for comparison with our present simulations.
Figure 3.14 shows the initial evolution of the hydrogen bubble from t = 1:5s to
t = 3:0s. At t = 1:5s, the bubble already interacts with the shock. The cold bubble
moves towards the hot air on the other side of the shock and generates a re
ected shock,
which travels left inside of the bubble. The shock is also deformed by the interaction with
the bubble. Betweent = 2:0s andt = 2:5s, some waves are re
ected from the left bound-
ary of the bubble and travels to the right. In gure 3.15, from t = 3:0s to t = 8:8s,
we can observe the generation of one major vortex and a secondary vortex resulting from
the misalignment between pressure gradient and density gradients. All these phenomena
discussed above are captured by all of the gures. Figure 3.16 shows the later evolution of
the hydrogen bubble. During this period, the main vortex keeps straining the bubble until
it is wrapped. After the comparison with Billet et al. (2008) at various time steps, we can
conclude that there is an overall good agreement between these two works.
53
Figure 3.14: Comparison between results of simulations of a shock/hydrogen-bubble inter-
action obtained with HybridX in the present work (left) and the reference data by Billet
et al. (2008) (right) for increasing simulation time (1.5s, 2.0s, 2.5s, and 3.0s from top
to bottom). Colormap shows the concentration ofH
2
with isolines of pressure superimposed.
54
Figure 3.15: Comparison between results of simulations of a shock/hydrogen-bubble inter-
action obtained with HybridX in the present work (left) and the reference data by Billet
et al. (2008) (right) for increasing simulation time (3.5s, 4.0s, 6.0s, and 8.8s from top
to bottom). Colormap shows the concentration ofH
2
with isolines of pressure superimposed.
55
Figure 3.16: Comparison between results of simulations of a shock/hydrogen-bubble interac-
tion obtained with HybridX in the present work (left) and the reference data by Billet et al.
(2008) (right) for increasing simulation time (13.6s, 16.6s, 25.6s, and 41.6s from top to
bottom). Colormap shows the concentration of H
2
with isolines of pressure superimposed.
56
Chapter 4
Results
This chapter presents the major ndings of the current work. Section 4.1 performs parametric
study to assess the eects of Sc, M
t
and Re
on passive scalar mixing in STI. The eects
of M
t
on ignition delay time of reactive 2D HIT simulations are investigated in section 4.2.
Results of 2D reactive STI are presented in section 4.3.
4.1 Passive Scalar mixing in 3D STI
This section focuses on passive scalar mixing in 3D STI. First, shock-normal statistics are
studied in section 4.1.1, including streamwise/transverse Reynolds stress/vorticity variance,
scalar variance and its budgets, scalar dissipation rate and its budgets. In section 4.1.2, the
change of
ow topology across the shock, and its eects on mixing are investigated. Align-
ments of passive scalar gradient, vorticity, shock-normal and strain directions are studied in
4.1.3, with eorts to explain why turbulent mixing is enhanced across the shock with the
change of alignments through the interaction.
4.1.1 Shock-normal statistics
In this section we focus on the streamwise evolution of scalar mixing statistics averaged on
transverse planes and in time. Statistics of the turbulent background
ow quantities (e.g.,
Reynolds stresses and vorticity variances) have been extensively documented for comparable
57
values of physical parameters in previous studies (e.g., Larsson et al., 2013) and are only
minimally discussed here, mainly to introduce relevant streamwise locations.
The region of shock unsteadiness diers among simulations, as the shock corrugation and
strength of the interaction vary with the change of dimensionless parameters (M, M
t
, and
Re
). To compare results from dierent simulations in the gures to follow, for each case:
1) the streamwise origin is shifted to the start of the unsteady shock region; 2) the streamwise
coordinate is rescaled by the dissipation-based length scale taken immediately upstream of
the unsteady shock region,L
;u
; 3) the unsteady shock region is either highlighted by a grey
shade and a horizontal line segment marking its width, or removed for better comparison of
the eective jumps of
ow quantities across the averaged unsteady shock regions; 4) when
the unsteady shock regions are removed, symbols mark the postshock location for each case,
following Table 2.1.
4.1.1.1 Reynolds stress, vorticity variance anisotropy, and relevant streamwise
locations
Figure 4.1(a,b) shows averaged proles of streamwise Reynolds stresses, R
xx
= u
002
, where
u
00
is the Favre streamwise velocity
uctuation, averaged in spanwise planes and time. For
each simulation case, R
xx
rst peaks in the region of shock unsteadiness, generally followed
by a second peak farther downstream (outside the unsteady shock region) that results from
an acoustic-to-vortical energy transfer (Lee et al., 1993; Larsson & Lele, 2009). For each
simulation, we dene the averaged preshock,x
pre
, and postshock,x
post
, locations by the rst
two local minima found in the prole of R
xx
. The unsteady shock region lies in between
these preshock and postshock locations, as marked by the greyed areas and the horizontal line
segments in gure 4.1(a,b). The width of the unsteady region is seen to decrease withM, and
increase withM
t
andRe
. The peak ofR
xx
downstream of the shock increases in magnitude
with M, tending toward saturation for M > 3, and decreases with Re
and M
t
, which also
bring it closer to the shock. Note that for case I, with (M;M
t
;Re
) = (1:5; 0:4; 74), the
58
1 0 1 2 3
0:50
0:75
1:00
1:25
1:50
1:75
f
u
002
=
f
u
002
u
(a)
0 2 4
0:4
0:6
0:8
1:0
1:2
1:4
f
u
002
=
f
u
002
u
(b)
0 1 2 3
x=L
;u
10
0
10
1
f
!
002
y
=
f
!
002
x
K(5:0; 0:3; 72)
H(5:0; 0:3; 37)
J(3:0; 0:3; 66)
G(3:0; 0:3; 37)
F (2:0; 0:3; 38)
D(1:5; 0:3; 39)
A(1:28; 0:3; 41)
(c)
0 2 4
x=L
;u
1:0
1:5
2:0
2:5
3:0
f
!
002
y
=
f
!
002
x
I(1:5; 0:4; 74)
E(1:5; 0:4; 40)
D(1:5; 0:3; 39)
C(1:5; 0:2; 38)
B(1:5; 0:1; 37)
(d)
Figure 4.1: Streamwise proles of Favre-averaged streamwise Reynolds stresses (a;b) and
vorticity variance anisotropy (c;d) for dierent simulation cases. Each curve in (a;b) has
been normalized by the value immediately upstream of the unsteady shock region. Grey
shaded regions and horizontal line segments in (a;b) mark the extent of the unsteady shock
region for each case, using the corresponding line style. All unsteady shock regions are
removed in (c;d) for a clearer comparison among cases. Vertical segments starting from the
horizontal axis indicate the location of x
^ !y
. M, M
t
, and Re
values for each case are given
in parentheses in the legend following table 2.1.
59
unsteady shock region is the widest and encloses the region of acoustic-to-vortical energy
transfer, such that the downstream evolution of R
xx
is monotonic, without a second peak.
Thus, caseI does not present a local minimum ofR
xx
immediately donwstream of the shock,
and the postshock location for this case is dened instead by the discontinuous change in
slope of R
xx
, as observed in gure 4.1(b).
The anisotropy of the small scales, represented by the ratio of vorticity variance in trans-
verse and streamwise directions, !
002
y
=!
002
x
, is plotted in gure 4.1(c,d). The unsteady shock
region has been removed from these plots for better comparison of the eective jump across
the shock for dierent cases. As in earlier studies (Larsson et al., 2013), the anisotropy
is highest immediately downstream of the shock, increasing for larger M, and smaller M
t
and Re
. In contrast with R
xx
, the anisotropy does not saturate for the range of M here
considered, and increasing Re
for cases in the wrinkled shock regime steepens the nega-
tive slope immediately downstream of the shock, toward a faster recovery of isotropy of the
small scales. In this study, we dene the location of recovery of small-scale isotropy, x
iso
,
at the streamwise station downstream of the shock where the vorticity variance anisotropy,
!
002
y
=!
002
x
, reduces to a value of 1.01.
4.1.1.2 Scalar variance
Figure 4.2 shows the eect ofM,M
t
, andRe
on the scalar variance,
00
00
, for all simulation
cases performed, considering Sc = 1 only rst. Whereas there is no local jump of passive
scalar variance across the shock, the averaging results into an eective jump across the
unsteady shock region, owing to its nite width. M has a nearly negligible in
uence on this
eective jump of scalar variance, due to the combined eect of a reduced unsteady shock
width and a steepening of the decay rate across the shock. This faster decay of scalar variance
downstream of the shock, evidenced by the increased negative slope of the proles, is the
result of the combined eect of an increased scalar dissipation rate (discussed inx4.1.1.3)
and a decreased mean velocity downstream of the shock.
60
1 0 1 2 3
x=L
;u
0:25
0:50
0:75
1:00
1:25
1:50
f
002
=
f
002
u
K(5:0; 0:3; 72)
H(5:0; 0:3; 37)
J(3:0; 0:3; 66)
G(3:0; 0:3; 37)
F (2:0; 0:3; 38)
D(1:5; 0:3; 39)
A(1:28; 0:3; 41)
(a)
1 0 1 2 3
x=L
;u
I(1:5; 0:4; 74)
E(1:5; 0:4; 40)
D(1:5; 0:3; 39)
C(1:5; 0:2; 38)
B(1:5; 0:1; 37)
(b)
Figure 4.2: Streamwise proles of Favre-averaged scalar variance for (a) cases with the same
M
t
and (b) cases with the same M. Each curve is normalized by the value immediately
upstream of the unsteady shock region, which is greyed out and marked, for each case, by
an horizontal segment with the same line style as the corresponding curve.
Increasing M
t
widens the unsteady shock region, due to the increased shock corruga-
tion, and translates into a larger eective jump of scalar variance across the shock. Larger
Re
slows down the decay rate of normalized scalar variance downstream of the shock, a
consequence of the reduced jump of scalar dissipation rate to be shown inx4.1.1.3).
4.1.1.3 Scalar variance budgets
The transport equation for the Favre-averaged passive scalar variance,
00
00
is (seexAppendix
A.1):
@
00
00
@t
+e u
i
@
00
00
@x
i
=
00
00
@e u
i
@x
i
| {z }
A
2u
00
i
00
@
e
@x
i
| {z }
B
@
@x
i
u
00
i
00
00
| {z }
C
+ 2
@
@x
i
D
00
@
@x
i
| {z }
D
2D
@
00
@x
i
@
00
@x
i
| {z }
E
2D
@
00
@x
i
@
e
@x
i
| {z }
F
;
(4.1)
61
where D is the diusivity of scalar . The left-hand side includes the local time derivative,
which vanishes once statistical stationarity is reached, and the convection of scalar variance
by the Favre-averaged mean velocity. The right-hand side includes six terms: 1) the di-
latational term,A, proportional to the dilatation of the Favre-averaged mean velocity; 2)
the production term,B, proportional to the Favre-averaged mean scalar gradient; 3-4) the
turbulent and molecular diusion terms,C andD, respectively; 5-6) the (turbulent) scalar
variance dissipation rate terms,E andF. For all simulation cases considered, termsB,D,
andF have a negligible contribution outside the unsteady shock region, and will not be
shown.
Figure 4.3 shows streamwise proles of the dilatational termA, the turbulent diusion
termC, and the scalar dissipation rate termE. Each prole has been normalized by the scalar
dissipation rate immediately upstream of the shock,E
u
, of each simulation. It is observed that
the scalar dissipation rate is the dominant term, in agreement with Tian et al. (2017), while
all three terms increase across the shock. Tian et al. (2017) did not perform a parametric
study of dierent terms in equation 4.1, while Boukharfane et al. (2018) only studied the
eects ofM (for 1.7, 2.0, and 2.3) on the scalar dissipation term. Dilatational and turbulent
diusion terms reach about ten percent of the scalar dissipation rate term downstream of
the shock, in contrast to the results of Tian et al. (2017), who found these two terms to
be negligible for the mixing of passive scalars, but not for multi-
uid mixing. Our present
simulations predict the dilatational and turbulent diusion terms to be also signicant in the
mixing of passive scalars, more so for larger M and smaller Re
. Boukharfane et al. (2018)
also observed a non-negligible turbulent diusion term. The dilatational term reaches its
peak shortly downstream of the the unsteady shock region, closer to the postshock location
for larger M. The normalized peak increases with M for M 3 but appears to saturate
or slightly decrease beyond that shock Mach number. Note that the eect of an increased
density jump across the shock for increasing M is included inA=E
u
.
62
0:0
0:1
0:2
0:3
A=E
u
K(5:0; 0:3; 72)
H(5:0; 0:3; 37)
J(3:0; 0:3; 66)
G(3:0; 0:3; 37)
F (2:0; 0:3; 38)
D(1:5; 0:3; 39)
A(1:28; 0:3; 41)
(a)
0:1
0:0
0:1
0:2
A=E
u
I(1:5; 0:4; 74)
E(1:5; 0:4; 40)
D(1:5; 0:3; 39)
C(1:5; 0:2; 38)
B(1:5; 0:1; 37)
(b)
0:1
0:0
0:1
0:2
0:3
C=E
u
(c)
0:10
0:05
0:00
0:05
0:10
C=E
u
(d)
0 1 2 3
x=L
;u
0:0
2:5
5:0
7:5
10:0
12:5
15:0
E=E
u
(e)
0 1 2 3
x=L
;u
0:0
0:5
1:0
1:5
2:0
E=E
u
(f)
Figure 4.3: Streamwise proles of non-negligible terms in the transport equation of scalar
variance for cases with the sameM
t
(a;c;e) and cases with the sameM (b;d;f). Each curve
is normalized by the value immediately upstream of the unsteady shock region (removed
for clarity). The preshock location is oset to be at x = 0 for all cases, and the postshock
location, x
post
, is marked with symbols, (a;b) dilatational term; (c;d) turbulent diusion
term; (e;f) Dissipation rate of scalar variance.
63
Downstream of the shock, a rapid decay of the dilatational and turbulent diusion terms
follows along a streamwise distance shorter thanL
;u
, and becomes slightly negative and then
recovering to a negligible value about 2L
;u
behind the shock.. Increasing Re
lowers the
reached peak of both dilatational and turbulent diusion terms downstream of the unsteady
shock region, and results in a faster decay and crossover to negative values. The latter is also
observed for increasing M
t
, although it can be considered an eect of the wider unsteady
shock region. We note that case I, which corresponds to the broken shock regime (for the
largest M
t
and Re
considered) does not show an increase of these two terms (dilatational
and turbulent diusion) immediately downstream of the shock, where both terms are already
negative instead, possibly due to a wider shock region that required a dierent denition of
the postshock streamwise location than all other cases. Despite the
uctuations observed in
the proles of the turbulent diusion term, which stem from the diculty to converge its
statistics in time, general trends can be observed.
The eective jump of scalar dissipation rateE across the shock monotonically increases
with M (in agreement with Boukharfane et al., 2018), and decreases with M
t
and Re
,
showing no sign of saturation in the range of M = 5 considered in the present study. The
lower jump for increasingRe
is responsible for the slower decay rate of scalar variance seen
in gure 4.2. To decouple the eect of increased density and viscosity across the shock for
cases with dierent M, we introduce a scaled scalar dissipation rate (SDR), ^ =E=2D =
D(@
00
=@x
i
)
2
=D. This quantity, plotted in gure 4.4 shows a nearly linear jump with
M across the unsteady shock region, up to M = 3 with signs of saturation beyond that
Mach number. A faster downstream decay is also observed as M increases. A larger Re
mainly lowers the postshock value. The wider shock region brought in by increasedM
t
leads
to eective jumps of ^ =^
u
below unity for the broken shock regime. It is expected that
as M
t
is further increased, possibly entering the `vanished (shock) regime', amplication of
scalar-related quantities will be progressively reduced, eventually leading to a decay across
the shock, as seen for turbulence quantities in Chen & Donzis (2019).
64
0 1 2 3
x=L
;u
0
1
2
3
4
^ = ^
u
(a)
0 1 2 3
x=L
;u
0:0
0:5
1:0
1:5
2:0
^ = ^
u
(b)
Figure 4.4: Streamwise proles of scaled scalar dissipation rate, ^ , for (a) cases with the
sameM
t
and (b) cases with the sameM. Each curve is normalized by the value immediately
upstream of the unsteady shock region (removed for clarity). Symbols mark the postshock
location, x
post
. Legend as in gure 4.3.
4.1.1.4 Eect of Schmidt number
The eect of the Schmidt number, Sc, is examined in gure 4.5, which shows streamwise
proles of scalar variance and SDR for three Schmidt numbers (Sc = 0:5, 1, and 2) extracted
from simulation case F . It is observed that Sc in this range has a negligible eect on the
scalar variance normalized by its corresponding value immediately upstream of the shock.
A small dierence is seen in the jump of normalized SDR across the unsteady shock region,
slightly larger in magnitude for increasing Sc, that persists downstream in the otherwise
similar evolution of SDR for all Sc. Similarly, Sc was found to have a negligible eect on
other scalar quantities (not shown) including other terms in the transport equations of scalar
variance and SDR, on the dissipation conditioned by
ow topology, and on the alignment of
the scalar gradient with vorticity and strain-rate tensor eigenvectors, to be discussed later.
In what follows, only results for Sc = 1 will be considered, and the eects of the other
dimensionless numbers (M, M
t
, and Re
) will be studied in detail.
65
0 1 2 3
x=L
;u
0:00
0:25
0:50
0:75
1:00
f
002
=
f
002
u
Sc = 2:0
Sc = 1:0
Sc = 0:5
(a)
0 1 2 3
x=L
;u
0
2
4
=
u
Sc = 2:0
Sc = 1:0
Sc = 0:5
(b)
Figure 4.5: Streamwise proles of (a) Favre-averaged scalar variance and (b) its scaled rate
of dissipation for scalars with dierent Sc obtained from case H (M = 5:0, M
t
= 0:3,
Re
= 37). Each curve is normalized by the value immediately upstream of the unsteady
shock region, shaded in greyed.
4.1.1.5 Scalar Taylor-like microscale
Figure 4.6 shows the in
uence of M, M
t
and Re
on the scalar Taylor scale, dened as
= (
00
00
=^ )
1=2
. It is observed that
decreases across the shock, more so for increasing
M and decreasing Re
. Downstream of the unsteady shock region,
grows at a rate that
increases withM andRe
. Cases with relatively lowM (1.28, 1.5) show a monotonic growth
rate, with a nearly uniform
for downstream distances below approximately L
;u
. Cases
with largerM (2-5) show a non-monotonic growth rate, with a fast growth very close to the
unsteady shock region x < L
;u
=3, followed by a slowdown and a subsequent acceleration.
Increasing M
t
shows a change in behavior from cases in the wrinkled shock regime, which
share a similar eective jump and downstream recovery of the normalized
, to cases in
the broken shock regime, for which the widened unsteady shock region results in smaller
eective jump and faster downstream recovery.
66
0 1 2 3
x=L
;u
0:4
0:6
0:8
1:0
=
;u
(a)
0 1 2 3
x=L
;u
0:8
1:0
1:2
1:4
=
;u
(b)
Figure 4.6: Streamwise proles of scalar Taylor-like microscale,
=
00
00
=^
1=2
, for (a)
cases with the same M
t
and (b) cases with the same M. Each curve is normalized by the
value immediately upstream of the unsteady shock region, collapsed for better comparison
of the eective jumps for all cases.
4.1.1.6 Scalar dissipation rate budgets
The Reynolds-averaged transport equation for the (scaled) scalar dissipation rate, =
(@=@x
i
)
2
, can be expressed as (seexAppendix A.2):
@
@t
+u
k
@
@x
k
=u
0
k
@
@x
k
| {z }
G
2
@u
k
@x
i
@
@x
i
@
@x
k
| {z }
H
+ 2
@
@x
i
@
@x
i
1
@
@x
k
D
@
@x
k
| {z }
I
; (4.2)
Budgets of scalar dissipation were previously analyzed in compressible DHIT (Danish et al.,
2016, e.g.), but the shock-normal evolution of these terms under STI has not been reported
to date. The local time derivative on the left-hand side vanishes in the statistically stationary
regime reached after the initial transient. The convection of scalar variance by the Reynolds-
averaged mean velocity is then balanced by the terms on the right-hand side. The rst term,
G, represents the correlation between Reynolds velocity
uctuations and the gradient of the
scalar dissipation rate. The second term,H, results from the interaction of spatial gradients
of the velocity and passive scalar elds. The alignment between the scalar gradient and the
eigenvectors of the velocity gradient tensor thus plays a key role in the resulting amplication
67
or reduction of scalar dissipation and in the overall mixing. The eects of such alignment,
as well as the scalar dissipation rate conditioned on
ow topologies (dened from invariants
of the velocity gradient tensor) will be studied in detail inx4.1.3 and x4.1.2. The last term,
I, accounts for molecular diusion processes.
Figure 4.7 compares the contribution ofG,H, andI for all simulated cases. Each prole
is normalized with the preshock magnitude of the molecular diusion term,jE
u
j. Note that
the velocity-scalar interaction term has a net positive contribution (i.e., amplication of
dissipation rate hence leading to better mixing) whereas the molecular diusion term has a
net negative contribution to the right-hand side equation 4.2.
The contribution from the
uctuating-velocity/dissipation-gradient correlation term is
predominantly positive, but much smaller (below 1%) than the other two terms. Whereas
the velocity-scalar-gradient and molecular diusion terms are of the same order, the latter
is always the largest in magnitude (for all cases and streamwise locations). The molecular
diusion term shows an eective jump in magnitude across the shock that is monotonically
increasing with largerM, and smallerRe
andM
t
. The streamwise proles show a monotonic
shallowing of the slope downstream of the shock. In contrast, the velocity-scalar-gradient
interaction term experiences a jump of magnitude across the shock that appears to saturate
for M' 3 and increases with larger Re
and smaller M
t
, in the wrinkled shock regime.
4.1.1.7 Spectra and PDFs
Time-averaged spectra of passive scalar calculated on transverse planes,
E
, at pre- and
postshock locations are compared in gure 4.8(a) for dierent simulation cases. Within
the limited inertial range of lengthscales of the present simulations, the spectra of kinetic
energy (not shown) upstream of the shock exhibit an approximate
E(k)/
2=3
k
5=3
scaling,
consistent with the prediction of Kolmogorov (1941) for incompressible HIT. However, the
passive scalar spectra in the inertial-convective range depart from the
E
(k)/
1=3
k
5=3
prediction of M. (1949) and Corrsin (1951) for passive scalar mixing in incompressible HIT.
68
0:03
0:02
0:01
0:00
0:01
0:02
0:03
G=jI
u
j
(a)
0:000
0:005
0:010
0:015
0:020
G=jI
u
j
I(1:5; 0:4; 74)
E(1:5; 0:4; 40)
D(1:5; 0:3; 39)
C(1:5; 0:2; 38)
B(1:5; 0:1; 37)
(b)
0
1
2
3
4
H=jI
u
j
K(5:0; 0:3; 72)
H(5:0; 0:3; 37)
J(3:0; 0:3; 66)
G(3:0; 0:3; 37)
F (2:0; 0:3; 38)
D(1:5; 0:3; 39)
A(1:28; 0:3; 41)
(c)
0:0
0:5
1:0
1:5
2:0
H=jI
u
j
(d)
0 1 2 3
x=L
;u
0:0
2:5
5:0
7:5
10:0
12:5
I=jI
u
j
(e)
0 1 2 3
x=L
;u
0:0
0:5
1:0
1:5
2:0
2:5
I=jI
u
j
(f)
Figure 4.7: Streamwise proles of right-hand-side terms in the Reynolds-averaged transport
equation for the scaled dissipation rate of scalar variance for cases with the same M
t
(a,c,e)
and cases with the same M (b,d,f). (a,b) Correlation between the
uctuating velocity and
the gradient of dissipation rate,G; (c,d) Interaction between velocity gradient tensor and
scalar gradient,H; (e,f) Molecular diusivityI. Each curve is normalized byI
u
, the absolute
value of the molecular term immediately upstream of the unsteady shock, removed for better
comparison of the eective jumps for dierent cases.
69
10
0
10
1
10
2
k=k0
10
9
10
7
10
5
10
3
10
1
10
1
5=3
(a) (a)
5=3
(a) (a)
5=3
(a) (a)
5=3
(a) (a)
5=3
(a) (a)
5=3
(a) (a)
E
(k)
E(1.5,0.4,40)
G(3,0.3,37)
H(5,0.3,37)
I(1.5,0.4,74)
J(3,0.3,63)
K(5,0.3,72)
10
0
10
1
10
2
k=k0
10
8
10
6
10
4
10
2
10
0
x
pre
x
post
x
pre
x
post
Re
40
Re
70
(b) (b)
x
pre
x
post
x
pre
x
post
Re
40
Re
70
(b) (b)
x
pre
x
post
x
pre
x
post
Re
40
Re
70
(b) (b)
x
pre
x
post
x
pre
x
post
Re
40
Re
70
(b) (b)
x
pre
x
post
x
pre
x
post
Re
40
Re
70
(b) (b)
x
pre
x
post
x
pre
x
post
Re
40
Re
70
(b) (b)
E
0
0(k)
Figure 4.8: Spectra of (a) passive scalar and (b) scalar variance for dierent simulation
cases at pre- and postshock locations. Each spectrum is normalized by its integral over all
wavenumbers. Spectra for cases with higher Re
( 70) have been shifted up two decades.
Postshock spectra are shifted by a factor of 3 relative to preshock counterparts. The black
straight segment with5=3 slope marks the extent of the inertial range of scales for the
higher Re
cases, obtained from spectra of the kinetic energy (not shown).
A slope shallower than5=3 is observed in the inertial-convective range upstream of the
shock for all cases, which may be attributed to the moderate Reynolds numbers of our
simulations, following Danaila & Antonia (2009). For all cases, the spectral content of passive
scalar and its variance increases across the shock for small scales (gure 4.8a,b), plausibly a
consequence of the immediate amplication of transverse vorticity. In contrast, spectra of
scalar dissipation,
E
, display a shift of spectral content across the shock from smaller to
larger scales for increasingM and lowerRe
(see gure 4.9a,b) Taking caseK as an example
(gure 4.9c) amplication of all wavenumbers is rst observed in the postshock state (x
post
),
predominantly at very small scales, followed by large and intermediate scales. Betweenx
post
andx
^ !y
, the readjustment of transverse into streamwise vorticity counteracts the postshock
viscous decay at small scales, such that the spectral content of scalar dissipation is reduced
much more signicantly at large and intermediate scales than at small scales. Downstream
of x
^ !y
, viscous decay aects all scales similarly, preserving the shape of the spectra.
Time-averaged, transverse-plane PDFs of passive scalar, its variance and its dissipation
for case K at dierent streamwise locations are shown in gure 4.10.
70
10
0
10
1
10
2
k=k0
10
5
10
4
10
3
10
2
10
1
(a)
Re
40
Re
70
(a)
Re
40
Re
70
(a)
Re
40
Re
70
(a)
Re
40
Re
70
(a)
Re
40
Re
70
(a)
Re
40
Re
70
E
(k;x
pre
)
E(1.5,0.4,40)
G(3,0.3,37)
H(5,0.3,37)
I(1.5,0.4,74)
J(3,0.3,63)
K(5,0.3,72)
10
0
10
1
10
2
k=k0
10
5
10
4
10
3
10
2
10
1
(b)
Re
40
Re
70
(b)
Re
40
Re
70
(b)
Re
40
Re
70
(b)
Re
40
Re
70
(b)
Re
40
Re
70
(b)
Re
40
Re
70
E
(k;x
post
)
E(1.5,0.4,40)
G(3,0.3,37)
H(5,0.3,37)
I(1.5,0.4,74)
J(3,0.3,63)
K(5,0.3,72)
10
0
10
1
10
2
k=k
0
10
0
10
1
(c)
E
(k;x)=E
(k;x
pre
)
x
pre
x
post
x
^ !y
x
iso
Figure 4.9: Comparison of time-averaged spectra of scalar dissipation on transverse planes
for dierent simulation cases in the (a) preshock and (b) postshock states. Each spectrum is
normalized to unitary integral. Spectra for cases with higher Re
( 70) are shifted up one
decade for clarity. (c) Comparison of spectra of scalar dissipation at dierent streamwise
locations for caseK, where each spectrum in (c) is normalized by the value at the same wave
number in the preshock spectrum.
Super-Gaussian passive scalar PDFs, indicative of intermittency, show a monotonic nar-
rowing downstream as the scalar variance decreases, as previously observed in Tian et al.
(2017) and Boukharfane et al. (2018). The streamwise evolution of PDFs of scalar variance
and scalar dissipation (gure 4.10b,c) was not reported in those previous studies. PDFs of
scalar dissipation shift signicantly across the shock toward larger values, followed down-
stream by a progressive recovery to preshock shapes, completed between the streamwise
locations of maximum enstrophy (x
^ !y
) and small-scale isotropy recovery (x
iso
). Similar ob-
servations hold for other cases (not shown), with higherRe
resulting in wider tails of scalar
dissipation, consistent with the study of passive scalar mixing in compressible HIT
4.1.2 Flow topology and conditioned dissipation rate
4.1.2.1 Characterization of local
ow topology
The deformation, rotation, and stretching (hence, mixing) of material elements can be char-
acterized by the study of local
ow topologies introduced by Chong et al. (1990), which
71
2 0 2
10
7
10
5
10
3
10
1
(a) (a) (a) (a)
0:00 0:25 0:50 0:75 1:00
log(1 +
0
0
)
10
7
10
5
10
3
10
1
(b) (b) (b) (b)
0 2 4 6
log(1 +)
10
7
10
5
10
3
(c) (c) (c) (c)
x
pre
x
post
x
^ !
y
x
iso
Figure 4.10: PDFs of (a) passive scalar, (b) scalar variance, and (c) scalar dissipation for
caseK atx
pre
(dotted), x
post
(solid),x
^ !y
(dashed),x
iso
(dash-dotted). Dash-dot-dotted line
in (a) corresponds to a Gaussian distribution with zero mean and the standard deviation of
the PDF of passive scalar at x
iso
.
relates the underlying pattern of streamlines in a local frame moving with the
uid to the
invariants of the velocity gradient tensor using critical point theory. Compared with the
incompressible
ow case, compressibility introduces a richer set of
ow topologies (Suman
& Girimaji, 2010; Pirozzoli & Grasso, 2004) along with several possible formulations of the
topology analysis. In the present study, following Suman & Girimaji (2010), we normalize the
velocity gradient tensorA
ij
=@u
j
=@x
i
asa
ij
=A
ij
=
p
A
ij
A
ij
, from which normalized strain-
rate and rotation-rate tensors can be obtained as s
ij
= (a
ij
+a
ji
)=2 and w
ij
= (a
ij
a
ji
)=2,
respectively. The three invariants of a
ij
are then expressed as
p =tr(a
ij
) =a
ii
; (4.3)
q =
1
2
(p
2
s
ij
s
ji
w
ij
w
ji
); (4.4)
r = det(a
ij
) =
1
3
(p
3
+ 3pqs
ij
s
jk
s
ki
3w
ij
w
jk
s
ki
); (4.5)
72
0:01 0:00 0:01
r
0:1
0:0
0:1
q
SFS
UFS
SNSS
UN
3
USN
3
UFC
UNSS
(a)
0:01 0:00 0:01
r
SNSS UNSS
SFS UFC
(b)
0:01 0:00 0:01
r
SNSS
SFS
UNSS
SN
3
SSN
3
SFC
UFC
(c)
Figure 4.11: Flow topologies in pqr space mapped onto three representative planes with
(a) p < 0 (extension), (b) p = 0 (volume-preserving), and (c) p > 0 (contraction),
partitioned by the intersections of each plane with the dividing surfaces S1a (dashed),
S1b (dash-dotted), S2 (solid), and the r = 0 plane (dotted): SFS, stable focus stretch-
ing; UFS unstable focus stretching; SFC, stable focus compressing; UFC, unstable fo-
cus compressing; SNSS, stable-node/saddle/saddle; UNSS, unstable-node/saddle/saddle;
SN
3
, stable node/stable-node/stable-node; UN
3
, unstable-node/unstable-node/unstable-
node; SSN
3
, stable-star-node/stable-star-node/stable-star-node; USN
3
, unstable-star-
node/unstable-star-node/unstable-star-node. The last two topologies correspond to the in-
tersection of S1a and S1b for p< 0 and p> 0, respectively.
and determine the character (real or complex) of the three eigenvalues of A
ij
, thus dening
the local streamline pattern and
ow topology. Thepqr space is partitioned by three dividing
surfaces:
S1a;b :
1
3
p
q
2
9
p
2
2
27
3q +p
2
1=2
r = 0; S2 :pqr = 0 (4.6)
and the r = 0 plane into ten dierent
ow topologies (see gure 4.11). S1a;b separate
regions with three real eigenvalues (non-focal topologies) from regions with one real and two
complex conjugate eigenvalues (focal topologies). S2 contains purely imaginary eigenvalues.
Stable (unstable) topologies indicate that the local streamlines are directed toward (away
from) the critical point.
73
0:5 0:0 0:5
p
0
10
PDF(p)
^
E(jp)
(a)
PDF(p)
^
E(jp)
(a)
PDF(p)
^
E(jp)
(a)
PDF(p)
^
E(jp)
(a)
PDF(p)
^
E(jp)
(a)
PDF(p)
^
E(jp)
(a)
x
pre
0:6
0:8
1:0
1:2
0:5 0:0 0:5
p
0
5
PDF(p)
^
E(jp)
(b)
PDF(p)
^
E(jp)
(b)
PDF(p)
^
E(jp)
(b)
PDF(p)
^
E(jp)
(b)
PDF(p)
^
E(jp)
(b)
PDF(p)
^
E(jp)
(b)
x
post
0:6
0:8
1:0
1:2
0:5 0:0 0:5
p
0
5
PDF(p)
^
E(jp)
(c)
PDF(p)
^
E(jp)
(c)
PDF(p)
^
E(jp)
(c)
PDF(p)
^
E(jp)
(c)
PDF(p)
^
E(jp)
(c)
PDF(p)
^
E(jp)
(c)
x
iso
0:6
0:8
1:0
1:2
0:5 0:0 0:5
q
0
1
PDF(q)
^
E(jq)
(d)
PDF(q)
^
E(jq)
(d)
PDF(q)
^
E(jq)
(d)
PDF(q)
^
E(jq)
(d)
PDF(q)
^
E(jq)
(d)
PDF(q)
^
E(jq)
(d)
0:5
1:0
1:5
0:5 0:0 0:5
q
0
1
2
PDF(q)
^
E(jq)
(e)
PDF(q)
^
E(jq)
(e)
PDF(q)
^
E(jq)
(e)
PDF(q)
^
E(jq)
(e)
PDF(q)
^
E(jq)
(e)
PDF(q)
^
E(jq)
(e)
0:5
1:0
1:5
0:5 0:0 0:5
q
0
1
PDF(q)
^
E(jq)
(f)
PDF(q)
^
E(jq)
(f)
PDF(q)
^
E(jq)
(f)
PDF(q)
^
E(jq)
(f)
PDF(q)
^
E(jq)
(f)
PDF(q)
^
E(jq)
(f)
0:5
1:0
1:5
0:18 0:00 0:18
r
0
5
10
PDF(r)
^
E(jr)
(g)
PDF(r)
^
E(jr)
(g)
PDF(r)
^
E(jr)
(g)
PDF(r)
^
E(jr)
(g)
PDF(r)
^
E(jr)
(g)
PDF(r)
^
E(jr)
(g)
0:6
0:8
1:0
1:2
0:18 0:00 0:18
r
0
20
PDF(r)
^
E(jr)
(h)
PDF(r)
^
E(jr)
(h)
PDF(r)
^
E(jr)
(h)
PDF(r)
^
E(jr)
(h)
PDF(r)
^
E(jr)
(h)
PDF(r)
^
E(jr)
(h)
0:6
0:8
1:0
1:2
0:18 0:00 0:18
r
0
5
10
PDF(r)
^
E(jr)
(i)
PDF(r)
^
E(jr)
(i)
PDF(r)
^
E(jr)
(i)
PDF(r)
^
E(jr)
(i)
PDF(r)
^
E(jr)
(i)
PDF(r)
^
E(jr)
(i)
0:6
0:8
1:0
1:2
D(1:5; 0:3; 39)
E(1:5; 0:4; 40)
G(3:0; 0:3; 37)
J(3:0; 0:3; 66)
H(5:0; 0:3; 37)
K(5:0; 0:3; 72)
Figure 4.12: PDFs ofp (a,b,c),q (d,e,f), andr (g,h,i), along with normalized distributions of
SDR conditioned on the respective invariant, extracted at x
pre
(a,d,g),x
post
(b,e,h), andx
iso
(c,f,i) for cases with dierent (M;M
t
;Re
). Each plot shows the PDFs of the invariant on
the bottom curves (left vertical axis), and the conditioned SDR distributions on top (right
vertical axis).
74
4.1.2.2 One dimensional PDFs and conditioned scalar dissipation distributions
For each case and streamwise location, x, we denote by E(j f) the distribution of scalar
dissipation rate, , conditioned on the invariant f 2 (p;q;r). This distribution is subse-
quently time-averaged and normalized by the Reynolds-averaged scalar dissipation rate, ,
at the same streamwise location,x, to obtain
^
E(jf) =hE(jf)i=, wherehi is the time-
averaging operator. Figure 4.12 shows time-averaged PDFs of the p,q andr invariants, and
the corresponding normalized conditional distributions of scalar dissipation rate,
^
E(j p),
^
E(j q), and
^
E(j r), extracted at the three streamwise locations introduced inx4.1.1.1:
1) preshock,x
pre
, where the streamwise Reynolds stress reaches the minimum before the in-
crease in the shock region; 2) postshock,x
post
, where the streamwise Reynolds stress reaches
the minimum before it increases after the shock region; and 3) return to vorticity-variance
(i.e., small-scale) isotropy,x
iso
. Results from only a subset of the simulation cases are plotted
for clarity, to highlight the eect of M, M
t
, and Re
.
PDFs of p (gure 4.12a,b,c) are nearly symmetric around p = 0 (i.e., incompressible
events), where they peak. The PDFs are wider for larger M
t
(owing to its larger velocity
uctuations relative to the sound speed) and narrower for increasing M and Re
. Imme-
diately downstream of the unsteady shock region, the peak of the p-PDF is reduced to
about half its preshock magnitude. The normalized distributions of SDR conditioned on
p,
^
E(j p), are nearly symmetric around p = 0, where they also peak, for all cases in
the wrinkled shock regime. For case E, the only case shown in gure 4.12 in the broken
shock regime, the preshock SDR distribution conditioned on p is asymmetric and skewed
toward p < 0. At x
post
the SDR distributions signicantly
atten, as the p-PDFs widen,
which indicates a more even distribution of the dissipation across compressible
ow topolo-
gies compared to the preshock state. The broken-shock regime simulation still shows an
asymmetric SDR-distribution conditioned on p at x
post
. At x
iso
all distributions become
predominantly symmetric about p = 0, including the broken-shock regime case, and wider
than in the preshock state.
75
Preshock q-PDFs (gure 4.12d,e,f) are negatively skewed (i.e., strain-rate dominated),
more so for increasing Re
. Immediately downstream of the unsteady shock region, the
shape of theq-PDFs is altered, especially for largerM, increasing the magnitude of the peak
(found at q 0), which narrows the core of the distributions. At x
iso
, the q-PDFs recover
the preshock shape, nearly overlapping for all simulation cases considered. The normalized
distributions of SDR conditioned on q,
^
E(j q), are skewed toward q < 0, indicating that
scalar dissipation occurs predominantly in strain-dominated topologies. At x
pre
andx
iso
the
distributions nearly collapse for all simulations. However, at x
post
, largerM tends to reduce
the skewness, and a local maximum atq 0 develops. IncreasingRe
reduces the previously
described eects of a larger M on both the PDF(q) and
^
E(jq) immediately downstream
of the unsteady shock region.
PDFs ofr (gure 4.12g,h,i) atx
pre
andx
iso
locations are nearly insensitive to the change
of M, M
t
, and Re
, and appear skewed toward positive values. At x
post
, higher M and
smaller Re
tend to narrow the PDFs, which attain larger peaks at r = 0. A tendency
toward symmetrization of the PDF of this third invariant was observed using LIA and DNS
by Ryu & Livescu (2014). The distribution of SDR conditioned on r,
^
E(j r), is largely
skewed toward r> 0 at all streamwise locations, more so for larger Re
, peaking at r 0:1
for x
pre
and x
iso
, while
attening and bringing its peak closer to r = 0 at x
post
.
4.1.2.3 Joint qr PDFs and conditioned scalar dissipation distributions
Distinct regions in the qr planes for dierent values of p dene a variety of
ow topologies
schematically represented in gure 4.11. The predominance of dierent
ow topologies and
their relative role on the dissipation of passive scalar at dierent streamwise locations is
explored next. Figure 4.13 shows, for the case K, with (M;M
t
;Re
= (5; 0:3; 72), joint qr-
PDFs (contour lines) and conditional SDR distributions (colored contourmap) atx
pre
,x
post
,
andx
iso
locations, for three dierent values ofp: negative (extension), zero (incompressible),
and positive (contraction). The teardrop shape commonly found in incompressible
ows is
76
Figure 4.13: Jointqr-PDFs (black-white contour lines at 15% and 60% of the peak value per
plot) and conditioned SDR distributions (jet-scale contours masked below 1% of the peak of
the qr-PDF) for case K at dierent streamwise locations (left to right: x
pre
, x
post
, x
^ !y
, x
iso
)
and values of p (top to bottom: mean plus 1.5, 0 and -1.5 times the standard deviation per
location). Short-dashed black lines: S1a, S1b, S2 curves and q axis (see Suman & Girimaji,
2010). In the color contour, red corresponds to the highest conditioned SDR in each gure,
and blue corresponds to zero.
recovered for p = 0 at the three streamwise locations considered. SDR concentrates on the
right lower tail (q< 0 andr> 0) of the teardrop, below and near theS1b line, which conrms
the UNSS (unstable node saddle saddle) topology as the most dissipative at those streamwise
locations, consistent with previous ndings for compressible and incompressible turbulence
(Danish et al., 2016; Brethouwer et al., 2003; O'Neill & Soria, 2005). Both positive and
negativep values modify the shape of theqr-PDFs, tending to a symmetrization of its outer
edge, whereas the distribution of SDR remains largely concentrated in the UNSS topology.
Figure 4.14 compares the jointqr-PDF and the corresponding conditioned SDR distributions
for p = 0 among the three cases with the higher Re
( 70). As M increases, the joint qr-
PDF preserves the teardrop shape but considerably narrows, concentrating near the origin.
The conditional SDR distribution still appears predominantly concentrated in the fourth
quadrant for the higher M cases, spreading between the r = 0 axis and S1b. The peak of
77
Figure 4.14: Joint qr-PDFs (dashed black-white contour lines at 15% and 60% of the peak
value per plot) and conditioned SDR distributions (jet-scale contours masked below 1% of
the peak of the qr-PDF) for case I (left), J (middle), K (right) at x
post
for p = 0. Color
map represents the conditioned scalar dissipation (blue to green to yellow to red, from zero
to highest value in each plot).
the SDR distribution shows a progression toward the origin as M increases, as indicated
previously in the parametric study 1D distributions of SDR conditioned on q and r.
4.1.2.4 Proportion of topologies and their contribution to the scalar dissipation
Figure 4.15(a) shows a comparison, for dierent simulation cases, of the streamwise evolution
of the proportion of grid cells with a given
ow topologyT , averaged in time, denoted
hT
p
(x)i, where brackets represent the time-averaging operator. The corresponding time-
averaged proportion of scalar dissipation rate carried by each topology, denotedh
pjT
i, is
shown in Figure 4.15(b). The instantaneous proportions are calculated as:
T
p
(x) =
0
@
X
c2C
T
(x)
V
c
1
A
,
0
@
X
c2C(x)
V
c
1
A
;
pjT
(x) =
0
@
X
c2C
T
(x)
V
c
1
A
,
0
@
X
c2C(x)
V
c
1
A
(4.7)
whereC(x) is the set of grid cells intersected by the transverse plane with streamwise coordi-
natex,C
T
(x) is the subset of cells withinC(x) with a topologyT obtained from the invariants
of the normalized velocity gradient tensor (p, q, and r) calculated at that cell, and V
c
is the
volume of a cell c. For all simulation cases considered, SFS is the most dominant topology,
accounting for over one third of the volume, followed by UNSS ( 25%), UFC ( 20%),
SNSS ( 10%), SFC ( 8%), and UFS ( 1%). Upstream of the shock, these proportions
78
1 0 1 2 3 4
x=L
;u
0:01
0:04
1 0 1 2 3 4
x=L
;u
0:01
0:04
0:07
0:08
0:04
0:06
0:08
0:12
0:09
0:12
0:18
0:19
0:15
0:17
0:25
0:29
0:30
0:34
0:31
0:37
hT
p
i
0:29
0:41
SFS UNSS
UNSS SFS
UFC
UFC
SNSS SNSS
SFC SFC
UFS UFS
h
pjT
i
(a) (b)
B(1:5; 0:1; 37)
D(1:5; 0:3; 39)
G(3:0; 0:3; 37)
J(3:0; 0:3; 66)
H(5:0; 0:3; 37)
K(5:0; 0:3; 72)
Figure 4.15: Streamwise evolution of the time-averaged proportion of (a)
ow topologies,
hT
p
i, and of (b) the scalar dissipation in each topology,h
pjT
i, for dierent simulations.
Topologies ordered by largest relative contribution (top to bottom). Unsteady shock regions
are removed and streamwise axes are normalized by the dissipation length scale immediately
upstream of the shock. Symbols mark the postshock location, x
post
.
79
of topologies remain nearly uniform with x for each simulation case. Downstream of the
unsteady shock region, the proportion of SFS and UNSS decreases, whereas that of UFS
and SNSS increases by comparable amounts, respectively. These changes are amplied with
larger M. Increasing Re
mainly shifts the proportions throughout the domain, enhancing
those of non-focal topologies (UNSS and SNSS), while reducing all other. The dominance
of SFS and UNSS topologies has been observed in previous numerical studies of compress-
ible HIT (e.g., Suman & Girimaji, 2010; Wang et al., 2012) and STI (Ryu & Livescu, 2014;
Boukharfane et al., 2018), the latter for lowerM,M
t
, andRe
than the ones here considered.
A characterization of the downstream evolution of these topologies has not been reported to
date.
An initial recovery of the proportion of each topology toward its preshock state follows
shortly downstream of the shock, within the region of acoustic-to-vortical energy transfer (see
gure 4.1a,b). Increasing Re
shortens this initial recovery distance. Farther downstream,
cases with lowerM recover the upstream proportion for each topology within 3L
;u
, whereas
cases with largerM show a sustained downstream reduction of SFS and UNSS topologies at
the expense of an increased proportion of SFC and UFC topologies. The order of dominant
topologies remains practically unchanged within the computational domain for all simula-
tions. UNSS and SFS topologies are the most sensitive to changes of dilatation (Parashar
et al., 2019), which is consistent with the larger changes experienced by these topologies
across the shock (gure 4.15a). From gure 4.15(b), UNSS presents the largest proportion
of scalar dissipation ( 40%), followed by SFS ( 33%), swapping order with respect to
the dominant topology proportions,T
p
, discussed earlier, while all other topologies maintain
the same order. Results (not shown) for the ratioh
pjT
i=T
p
, indicate that UNSS, followed
by SNSS and then SFS, have the largest relative mean scalar dissipation when conditioned
on the topology. Thus, nonfocal topologies are more dissipative than focal ones, as also
found for the kinetic energy dissipation in compressible mixing layers (Vaghe & Madnia,
2015). Multiscale analyses of compressible HIT have shown that UNSS is favored at small
80
scales (Danish & Meneveau, 2018). The larger amplication of small scales seen across the
shock, combined with the dominance of UNSS as the most dissipative topology can help
explain the postshock increase of scalar dissipation.
Among focal topologies, those that are fully compressing (SFC) or fully stretching (UFS)
are less dissipative than those with simultaneous compressing and stretching (SFS and UFC).
These results remain valid downstream of the unsteady shock region and are consistent with
ndings in compressible HIT (Danish et al., 2016). Besides, UNSS and SFS are known to
present the maximum vortex stretching rate on average from studies of mixing layers (Vaghe
& Madnia, 2015). Although associated with weaker enstrophy content, UNSS regions present
high vorticity stretching rates, explaining their predominant dissipation of TKE and scalar
variance. Despite the contraction of vorticity that dominates regions of compression, SFS
(characterized by radial contraction and axial stretching), presents a net positive stretching
rate of vorticity and thus contributes to the dissipation of TKE and scalar variance. The
tendency of vortex tubes (associated mainly with SFS topologies) to orient parallel to the
shock downstream could then be linked to the increase of dissipation (of TKE and scalar
variance). Studies tracking vortex tubes across the shock would help elucidate this issue.
However, immediately downstream of the unsteady shock region, UNSS and SFS topolo-
gies reduce their contribution to the scalar dissipation, whereas all other topologies increase
theirs. Larger M amplies these eects for all topologies except for SFS. At all streamwise
locations, largerRe
shifts the contribution of UNSS and SNSS towards a larger proportion
of the scalar dissipation, reducing the relative contribution of all other topologies, consistent
with the eects onT
p
seen earlier. Likewise,h
pjT
i proles also show an initial recovery
shortly downstream of the shock, followed by a slower asymptotic evolution toward preshock
values, except for a sustained decrease of scalar dissipation proportion seen for SFS at the
expense of an increase for UFS.
To investigate dominant paths of transition between dierent topologies, Lagrangian
tracking of 10 000
uid particles seeded upstream of the shock is conducted for case D
81
(M;M
t
;Re
1:5; 0:3; 40), following their trajectories across and downstream of the shock.
The time origin of each trajectory is shifted to match its shock-crossing instant. From the
ensemble of trajectories, the probability of transition from each
ow topology into every
possible topology is calculated, as a function of non-dimensional time. From this study,
it is concluded that before and after the unsteady shock region, each topology has a high
probability of not changing its topology. SFS and UNSS are the least prone to change,
particularly upstream of the shock, whereas UFS, SNSS and SFC are most likely to transition,
although still with small probability. Across the shock, most topologies (SFS, UNSS, UFC,
SFC) exhibit a high probability of transitioning into SNSS. This is most noticeable for SFS
and UNSS, which are the most dominant topologies in the
ow, and also with the largest
contributions to the passive scalar dissipation. Although SNSS shows the largest probability
of not changing across the shock, transitions into SFS and UNSS are frequent, and remain
signicant farther downstream.
In summary, preshock
ow topologies and their contribution to the scalar dissipation are
altered across the shock. The widening of the nearly symmetric PDF of dilatation across
the shock enhances the role of compressible
ow topologies (compression-only, SFC, and
expansion-only, UFS) on the scalar dissipation, at the expense of a decreased contribution
from those topologies that are dominant in incompressible
ows (SFS and UNSS). Despite the
symmetry of the PDFs of dilatation (gure 4.12), preshock fractions of SFC and UFS topolo-
gies are asymmetric, favoring the compression-only SFC pattern. Immediately downstream
of the shock, the widened (still symmetric) PDFs of dilatation enhance the importance of the
expansion-only UFS topology. UFC (involving radial stretching and axial contraction with
outward spiraling streamlines) and SNSS (bi-axial contraction) are also favored and increase
their contribution to scalar dissipation across the shock. The observed changes become more
signicant for stronger shocks (larger M).
82
4.1.3 Alignment of passive scalar gradient, vorticity, shock-normal
and strain directions
Preferential alignment of vorticity and passive scalar gradient with the eigenvectors of the
strain-rate tensor is known to play a central role in the dissipation of turbulence kinetic
energy and passive scalar variance (Kerr, 1985; Ashurst et al., 1987; Jimenez, 1992), and
forms the basis of some structure-based models of turbulence and passive scalar ne scales
(Misra & Pullin, 1997; Pullin, 2000). We denote the eigenvectors of the strain-rate tensor,
S
ij
= (A
ij
+A
ji
)=2, as , and
, with their respective eigenvalues ordered from most
extensional to most compressive,
(note that their sum, + +
, need not be
zero for compressible
ows).
The amplication of the passive scalar dissipation rate brought by the net positiveness
of the termH =rAr in the transport equation for (discussed inx4.1.1.6) has
been related to the alignments of the principal strain directions and the passive scalar gra-
dient. In their numerical study of compressible decaying HIT, Danish et al. (2016) found
that this nonlinear amplication of the scalar dissipation rate is dominated by two non-focal
topologies, UNSS and SNSS, resulting from an enhanced alignment between the scalar gra-
dient and the most compressive eigenvector of the strain-rate tensor,
. In this section, we
present parametric studies, for dierent simulation cases, of the change of alignments across
the unsteady shock region, rst conditioned on
ow topology, and then followed by a global
assessment of the relation between these alignments and their contribution to the dissipation
of passive scalar
uctuations. The eects of variable M and Re
will be highlighted, as M
t
is seen to play a lesser role.
4.1.3.1 Alignments conditioned on
ow topology across the shock
For dierent simulation cases, gure 4.16 shows the change across the unsteady shock region
of alignments between the gradient of passive scalar,r, and both the most extensive,,
and the most compressive,
, principal strain directions, by means of time-averaged PDFs of
83
the cosine of the angles between those directions at pre- and postshock streamwise locations.
Preferential alignment with
is observed for all topologies, ordered left-to-right from most to
least aligned with
in that gure (UNSS>SNSS>SFS>UFC>UFS>SFC). For clarity, only
case E, with (M;M
t
;Re
) = (1:5; 0:4; 40), is shown at the preshock location, x
pre
, since all
other simulation cases share nearly the same PDFs at that location (i.e., the preshock state
is practically insensitive to M, M
t
and Re
). UNSS presents the most preferred alignment
ofr and
, followed by SNSS and SFS, whereas SFC shows the least alignment (although
still dominant) between those two directions, in favor of an increasing alignment ofr with
. Note that two of the topologies with better alignment ofr and
at x
pre
, UNSS and
SFS, were found inx4.1.2.4 to have the largest contribution to the scalar dissipation rate.
These preshock results are consistent with the ndings for compressible decaying HIT in
Danish et al. (2016), who attributed variations in the amplication (production) of scalar
dissipation to the ability of a topology to enhance the alignment between the scalar gradient
and the most compressive strain-rate eigendirection. However, immediately downstream of
the shock, atx
post
, the alignments show a strong dependence onM andRe
, but not onM
t
:
alignment with
becomes less dominant for all topologies, favoring an increased alignment
with . This trend becomes more obvious with larger M and smaller Re
. The change
of alignments for simulation cases in the broken shock regime is signicantly smaller than
for the wrinkled shock regime, likely resulting from the wider unsteady shock region in the
broken shock cases (the same trend was observed for case A, in the broken-shock regime,
not shown in gure 4.16).
Figure 4.17 shows the change, across the unsteady shock region, of the alignments between
r and the intermediate strain-rate eigenvector, , the vorticity, !, and the streamwise
direction, x. All topologies present nearly identical preshock and postshock PDFs, so only
SFS, the most dominant topology found in 4.1.2, will be discussed in the parametric study,
with the same conclusions applying to all other topologies. Alignment betweenr and
at x
pre
is small, and becomes even smaller immediately downstream of the unsteady
84
0:5
1:0
1:5
PDF[j cos(;r)j]
(a)
UNSS SNSS SFS UFC UFS SFC
0 0.5 1
j cosj
1
2
3
4
PDF[j cos(
;r)j]
(b)
0 0.5 1
j cosj
0 0.5 1
j cosj
0 0.5 1
j cosj
0 0.5 1
j cosj
0 0.5 1
j cosj
E(1:5; 0:4; 40)
x
pre
E(1:5; 0:4; 40)
x
post
I(1:5; 0:4; 74)
x
post
G(3; 0:3; 37)
x
post
J(3; 0:3; 63)
x
post
H(5; 0:3; 37)
x
post
K(5; 0:3; 72)
x
post
Figure 4.16: PDFs of the cosine of the angle () between the passive scalar gradient,r, and
(a) the most extensional, , and (b) the most compressive,
, eigenvectors of the velocity
gradient tensor, conditioned on topology for dierent simulation cases, comparing preshock
and postshock states. Green dotted lines correspond to the preshock state for case E,
representative of the preshock state for all other simulation cases (not shown for clarity, as
they nearly overlap).
85
0 0.5 1
j cos(;r)j
0:5
1:0
1:5
2:0
2:5
PDF
(a)
0 0.5 1
j cos(!;r)j
0
1
2
3
4
5
(b)
0 0.5 1
j cos(x;r)j
0
2
4
6
8
10
(c)
E(1:5; 0:4; 40)
x
pre
E(1:5; 0:4; 40)
x
post
I(1:5; 0:4; 74)
x
post
G(3; 0:3; 37)
x
post
J(3; 0:3; 63)
x
post
H(5; 0:3; 37)
x
post
K(5; 0:3; 72)
x
post
Figure 4.17: PDFs of the cosine of the angle between the scalar gradient,r, and (a) the
intermediate eigenvector, , of the velocity gradient tensor, (b) the vorticity, !, and (c) the
streamwise axis, x, conditioned on the SFS topology (similar PDFs of these quantities are
observed for all other topologies).
shock region, in favor of better alignment with. Larger M and smaller Re
signicantly
reduce the alignment ofr with at x
post
, but almost no eect is observed for varying
M
t
. Upstream of the shock, the passive scalar gradient is predominantly orthogonal to
the vorticity, consistent with previous studies in incompressible and compressible HIT and
STI (Kerr, 1985; Ashurst et al., 1987; Boukharfane et al., 2018), whereas, as expected,
no preferential alignment exists betweenr and x. At x
post
, the scalar gradient becomes
slightly more orthogonal to! and signicantly aligned withx, more so for increasingM and
decreasingRe
. A similar trend withM was also observed by Boukharfane et al. (2018), who,
however, did not explore Re
eects nor condition the alignments on topologies, presently
analyzed.
4.1.3.2 Barycentric map representation of alignments
We propose a combined representation of the alignments between a vector of interest (e.g.,
vorticity, scalar gradient, unit vector inx) and the three eigenvectors of the strain-rate tensor,
, and
by the barycentric coordinates of a point inside an equilateral triangle. The
86
barycentric coordinates are calculated from the cosines of the angles formed by the vector of
interest,v, and the strain-rate eigenvectors as (a;b;c) = (cos(;v); cos(;v); cos(
;v))=,
where = cos(;v) + cos(;v) + cos(
;v), such that a +b +c = 1. The corresponding
Cartesian coordinates on the plane of the triangle are (;) = (b=2;
p
3b=2 +c). The
vertices of the triangle, with Cartesian coordinates given by (0; 0), (1=2;
p
3=2), and (1; 0),
represent perfect alignment with,, and
, respectively.
At each streamwise location and for eachv of interest, we calculate the -joint PDF of
the alignments ofv with,, and
combined, along with the distribution of the passive
scalar dissipation rate, , conditioned on those alignments. Joint PDFs and conditioned
distributions are then time-averaged, and can be represented by contours on a barycentric
map.
Figure 4.18 shows the barycentric map representation of the alignments between the
scalar gradient, the vorticity, and the unit x vector with the principal strain directions,
through time-averaged joint PDFs, for dierent simulation cases. As in guresThesis (4.16,4.17),
the preshock state is shown only for case E, with (M;M
t
;Re
) = (1:5; 0:4; 40), since it is
representative of all other simulation cases at x
pre
. All PDFs have a common, single-peaked
shape, monotonically decreasing away from the peak, and thus are suitable for representa-
tion by a reduced signature: the peak of each PDF is marked by a symbol, whereas the
isoline of the PDF at 80% of its peak value illustrates how concentrated the PDF is about
its peak value (i.e., the higher concentration the closer the 80% isoline is to the peak). The
peak location together with the shape and spread of the 80% isoline are used for comparison
among dierent simulation cases. In the preshock state,r shows (gure 4.18a) preferential
alignment with
, followed by, and a much smaller alignment with. Atx
post
, the peak of
the PDFs moves toward the center of the
side of the triangle, implying a relative increased
alignment with at the expense of a reduced alignment with
and, increasing the or-
thogonality with the latter. LargerM drastically enhances this trend, not only bringing the
peak of the postshock PDF closer to the
side, but also signicantly narrowing the PDF.
87
(a)
(b)
(c)
E(1:5; 0:4; 40)
x
post
I(1:5; 0:4; 74)
x
post
G(3; 0:3; 37)
x
post
J(3; 0:3; 63)
x
post
H(5; 0:3; 37)
x
post
K(5; 0:3; 72)
x
post
E(1:5; 0:4; 40)
x
pre
Figure 4.18: Barycentric map representation of the PDFs of alignment between the eigen-
vectors (most extensive,; intermediate,; and most compressive,
) of the the strain-rate
tensor and (a) passive scalar gradient (r), (b) vorticity (!), and (c) the streamwise direc-
tion (x). Symbols mark the peak of the PDF, whereas lines represent the isocontour at 80%
of the peak value of each distribution. Segments normal to each triangle side mark the 1=4,
1=2, and 3=4 partition points. Insets zoom into the region of interest of each barycentric
map.
An increasing Re
counteracts the eect of larger M, widening the PDF, and returning its
peak toward its preshock state. These ndings are in agreement with the individual trends
described inx4.1.3.1 for the alignments of the scalar gradient conditioned on dierent
ow
topologies.
Preferential alignment of vorticity with the intermediate strain eigenvector, is clearly
inferred from gure 4.18(b). The shape of the joint PDF along the side of the triangle
is indicative of the predominant orthogonality of ! with the most compressive principal
direction of strain,
. Larger M brings the postshock PDF closer to the vertex, while
increasing Re
opposes this trend. Prior studies in HIT (Danish & Meneveau, 2018) and
turbulent shear
ows (Fiscaletti et al., 2016) found better alignment of ! and at small
scales. Moreover, Gonzalez (2012) observed more ecient scalar mixing in HIT when! and
are better aligned. The larger amplication of small-scale turbulence that occurs across
the shock may thus explain the overall increased alignment shown by these postshock PDFs,
and, in turn, the enhanced scalar mixing.
88
The PDF of alignment between the scalar gradient and x (gure 4.18c) is centered and
symmetric upstream of the shock, as expected for isotropic turbulence. Downstream of
the shock, the PDF shifts toward the
side of the triangle, implying an almost equally
probable alignment ofx with the most compressive and extensional strain eigenvectors, while
the intermediate eigenvector becomes more orthogonal to the average shock-normal direction.
Larger M and lower Re
strengthen this eect.
Joint distributions of scalar dissipation rate conditioned on the
-alignments withr
are given in gure 4.19. The distribution concentrates near the
vertex of the triangle,
expected as preferential alignment with this eigenvector is responsible for an amplication
of the scalar gradient (note that the scalar dissipation rate is proportional to the square of
the scalar gradient magnitude). At x
pre
, the shape of the conditioned joint distribution is
rather insensitive to M, Re
, and M
t
, although the peak location shifts slightly toward
for larger M and smaller Re
. At x
post
, each conditioned distribution closely follows the
joint PDF of scalar gradient seen in gure 4.18 concentrating around its peak, which moves
closer to. The postshock distribution of scalar dissipation conditioned on x shows peaks
very near the
vertex, whereas the 80% isolevels spreads farther away from that vertex with
smallerM and higherRe
, nearly equipartitioned between and, with a slight preference
for the former. Comparing gure 4.18(a) and gure 4.18(b), it can be gured out that the
most preferred scalar gradient alignment becomes better matched with the most dissipative
scalar gradient alignment in the barycentric map across the shock, which explains why there
is an increase of SDR across the shock. Summarizing, in the preshock state, the scalar
gradient is mostly aligned with the most compressive strain-rate eigendirection, more so for
ow topologies responsible for most of the scalar dissipation. Interaction with the shock
better aligns the scalar gradient with the shock-normal direction, balancing the alignment
with both the most extensional and the most compressive eigendirections, while increasing
the orthogonality with the intermediate eigendirection (and in turn, with the vorticity). The
convergence between the most preferred alignment and the most dissipative alignment of
89
(a)
(b)
(c)
E(1:5; 0:4; 40)
I(1:5; 0:4; 74)
G(3; 0:3; 37)
J(3; 0:3; 63)
H(5; 0:3; 37)
K(5; 0:3; 72)
Figure 4.19: Barycentric map representation of distributions of scalar dissipation rate con-
ditioned on the alignment between the , ,
eigenvectors of the strain-rate tensor and
(a) scalar gradient (r) at x
pre
, (b) scalar gradient (r) at x
post
, and (c) the streamwise
direction (x) at x
post
.
passive scalar gradient and strain-rate principal directions is seen to be responsible for the
enhancement of mixing across the shock. These eects increase with larger M and smaller
Re
.
4.2 Autoignition in 2D HIT
To evaluate the eect of turbulence and compressibility on autoignition, numerical experi-
ments are conducted for mixtures with the same initial mole fraction 2H
2
: O
2
: 7Ar. These
mixtures have dierent initial mean temperatures (T = 1000 K, 1100 K, 1200 K, 1300 K,
1400 K), dierent initial M
t
(0:1; 0:2; 0:3; 0:4) and dierent Re
(800; 1600; 3200). There are
60 cases in total, and the ignition delay time is obtained for each case by ensemble aver-
aging of 10 simulations. At the moment of ignition, the nondimensional k
max
for all the
simulations is larger than 1:5, ensuring that the turbulence is adequately resolved. Both
initial scalar and velocity elds are initialized with the spectra given in equation (2.30), with
initial rms of mass fractions of H
2
and O
2
set as 10% of its mean, while the mass fraction
of Ar is obtained considering that the mass fractions of all species must add up to unity, by
denition. The initial elds are obtained following the procedure described in section 2.4.2.
90
The mean ignition delay timet
ig
for dierent cases is compared with 0D simulation results
from Cantera with the same initial temperatures (T = 1000 K, 1100 K, 1200 K, 1300 K,
1400 K). The introduction of turbulent
uctuations of dierent degree of compressibility
leads to local temperature that can be higher than the mean temperature, which will shorten
the t
ig
, as observed in gure 4.20. The gure also shows that larger M
t
can bring shorter
t
ig
because of higher local temperature. The same observation is made by Stankovi c et al.
(2011) from numerical simulations of nitrogen-diluted hydrogen jets in turbulent co-
ows
of heated air. The trend of change of t
ig
with Re
is not conclusive, so the comparison at
dierentRe
is not presented. The ratio between standard deviationt
ig;
and mean ignition
delay timet
ig
is also presented in gure 4.20. It can be observed that high standard deviation
happens when temperature and M
t
are both high, and higher Re
also makes the standard
deviation higher. Similar simulations are also conducted on a stoichiometric mixture of a
partially-cracked JP7 surrogate (consisting of a binary mixture of 36% CH
4
and 64% C
2
H
4
by volume) and air, with results shown in gure 4.21. A reduced mechanism with 32 species
and 206 reactions from Luo et al. (2012) is adopted. The comparison between t
ig
obtained
from 0D simulations and 2D turbulent simulations shows that a higherM
t
will make ignition
delay time shorter also for JP7. Compared with gure 4.20, the same M
t
shows a larger
reduction in IDT when JP7 is used as fuel.
4.3 Reactive 2D Shock Turbulence Interaction
To investigate the eects of turbulence on shock-induced ignition, a laminar 1D simulation
(case P ) is conducted rst as a baseline. This case corresponds to the laminar 1D counter-
part of casesL toO, with the sameT
1
,p
1
and length before the sponge layer. However, this
1D laminar simulation has a ner resolution (51200 points before the sponge layer) in the
streamwise direction. This 1D laminar simulation is also used to decide the adequate resolu-
tion for the reactive STI simulations . 1D simulations with 5120 points and 51200 points are
91
Figure 4.20: The mean (top) and normalized standard deviation (bottom) for 2D reactive
HIT simulations of ignition delay time at dierent initial T ,Re
andM
t
compared with 0D
simulations for H
2
/O
2
/Ar mixture.
92
Figure 4.21: Mean ignition delay time at dierent initial T and M
t
compared with 0D
simulation for JP7 and air mixture.
conducted. The comparison is not shown, because the line plots show no dierence, which
indicates that the grid is already converged at 5120 points. Results for the simulation with
51200 points will be presented as the 1D laminar case.
The time evolution of dierent
ow quantities is shown in gure 4.22. In the time
evolution, three stages of the simulation can be observed. First, ignition starts at t = 180s
at a localized spot of high temperature emerging downstream of the shock. Then the ignition
and
ame evolves, and propagates upstream towards the shock. At t = 290s, a clear
ame
front is formed, and is about to catch up with the shock. The formation of a secondary
shock is also present in Sharpe (2002) and Melguizo-Gavilanes et al. (2011). When the
secondary shock catches up with the shock, a rarefaction wave is re
ected towards the outlet
boundary. After the
ame catches up with the shock, the peaks of pressure, velocity and
density catch up with the
ame front, and a detonation front is formed, which is also observed
in Sharpe (2002). The
ame keeps propagating towards the inlet boundary. At t = 370s,
we can observe a
ame front with peak pressure, density and velocity. All the reacting
93
Figure 4.22: Time evolution of normalized thermodynamic quantities in case P , with y axis
normalized by preshock value, and x axis normalized by the domain length before the sponge
layer.
94
STI simulations (from case L to case P ) end when the
ame reaches the inlet boundary.
The phenomena that the
ame propagates back towards the inlet boundary is called \
ame
ashback". This has been observed in scramjets in dierent experimental and numerical
works (O'Byrne et al., 2005; Micka, 2010; Sun et al., 2015; Zhao et al., 2019). In O'Byrne
et al. (2005), the
ashback is concluded to be initiated by the transition from diusive
combustion to premixed combustion. In Sunami et al. (2006); Sunami & Kodera (2012);
Wang et al. (2012), it is observed that the
ashback started from the well-premixed region in
a fashion similar to DDT (De
agration-to-Detonation Transition). This might be explained
by the nding in Ma et al. (2005) that the acoustic wave generated by unsteady heat release
can propagate upstream in the subsonic region between the
ame and the shock wave. These
acoustic waves cannot penetrate the shock into the supersonic region, so they are trapped
between the
ame and the shock. The pressure in this region is then increased with a
peak, and the
ame front is formed in a DDT fashion. Sharpe (2002) also attributes the
propagation of
ame towards the leading shock to the acoustic waves between the leading
shock and the ignition spot. Figure 4.23 shows dierent thermodynamic quantities at
t = 230s, a moment after a clear ignition spot is emerged and the
ame front has not
caught up with the shock yet for all the cases. The lines for the turbulent cases are ensemble
averaged over three runs with the same parameters. The quantities for cases with M
t
= 0:3
are ltered with a Wiener lter with a lter width of 0:47 cm, and quantities for cases with
M
t
= 0:1 are ltered with a Wiener lter with a lter width of 0:07 cm. Wiener lter has
already been adopted successfully to denoise data from
uid mechanics experiments (V etel
et al., 2011; Oxlade et al., 2012; Brindise & Vlachos, 2017). It is obvious that the cases with
higher M
t
have
ame fronts closer to the shock, and the cases with lower M
t
have
ame
fronts further away from the shock. This nding is consistent with what was observed in the
autoignition in 2D HIT. A higher M
t
will lead to higher local temperature, and this higher
local temperature will lead to earlier ignition. In 2D STI simulations, this earlier ignition
will be re
ected as a shorter induction length, which is the major reason why the simulations
95
Figure 4.23: Comparison of dierent thermodynamic quantities between 1D laminar reactive
STI simulation and 2D turbulent reactive STI simulations at t = 230s, with line style the
same as table 2.2.
96
Figure 4.24: Comparison of dierent thermodynamic quantities between 1D laminar reactive
STI simulation and 2D turbulent reactive STI simulations at t = 370s, with line style the
same as table 2.2.
97
with higher M
t
have
ame fronts closer to the shock. Even though the spatially averaged
ame front of cases with M
t
= 0:3 does not develop a clear peak as observed in gure 4.22,
instantaneous contour plots show higher local peaks compared with the M
t
= 0:1 and 1D
laminar cases.
It is also observed that the cases with M
t
= 0:1 do not show
ame fronts closer to the
shock compared with the 1D laminar case. This is attributed to the partially premixed
combustion that is present in the turbulent cases, in contrast to the premixed combustion
of the 1D laminar case. In premixed combustion, the combustion rate is limited by chemical
kinetics, while in non-premixed combustion, the combustion rate is limited by diusion. The
time scale of chemical kinetics is much smaller than the diusion time scale, so the reaction
rate for a non-premixed
ame is typically lower than a premixed
ame, and that of a partially
premixed
ame is in between. The 1D laminar simulation represents premixed combustion,
and the turbulent simulations represents partially-premixed combustion. For simulations
withM
t
= 0:1, ignition still occurs before the 1D laminar case, which will be shown later in
the discussion of global consumption rate for reactants, but the partially premixed nature
of the
ame makes the reaction rate and accumulation of thermal energy slower than in
the 1D laminar case. Hence, at t = 230s, even though the two cases with M
t
= 0:1 were
ignited earlier than the 1D laminar case, the slower reaction rate still makes its spatially-
averaged peak velocity, density and pressure lower than the 1D laminar case. The fact that
the averaged peak of cases with M
t
= 0:3 is closer to the shock means that the cases with
higher M
t
are ignited much earlier than the lower M
t
and 1D laminar cases, as discussed
in section 4.2. The same trend is observed at t = 370s in gure 4.24, in which lines for
turbulent cases are also ensemble-averaged over three runs with the same parameters. At
that time instant, the
ame fronts of all the cases have caught up with the shock. Cases with
M
t
= 0:3 have
ame fronts closer to the inlet boundary, then followed by the 1D laminar
case, and then the cases with lower M
t
.
98
According to Oppenheim et al. (1975), shock-induced ignitions can be divided into two
regimes, weak ignition and strong ignition. Weak ignition is characterized by the formation
of isolated
ame kernels in the early stage of ignition. Then these kernels develop and
become connected with each other, followed by the formation of planar combustion. Strong
ignition is characterized by the instantaneous formation of the planar
ame front. A very fast
formation of the
ame indicates that current simulations fall into the strong ignition regime.
Meyer & Oppenheim (1971) also showed where the strong ignition regime of H
2
/O
2
mixture
is located in the Arrhenius plot (t
ig
[O
2
] vs 1000=T , where [O
2
] is molar concentration of
O
2
). The current simulations are located at point (0.91K
1
, 0:2 mols l
1
) in the Arrhenius
plot, which falls into the strong ignition regime. However, in some of our simulation cases,
it appears that
ame fronts are initiated from dispersed ignition kernels, which indicates
weak ignition. The time evolution of temperature contours for case L is shown in gure
4.25. At t = 1:5 10
4
s, ignition regions start to emerge. Then at t = 2:0 10
4
s, the
ignition regions start to expand both upstream and downstream. At t = 3:0 10
4
s, the
ame region reaches the shock, but the detonation front is not formed yet. Att = 3:510
4
s,
after the
ame region passes through the shock, a detonation front is nally formed. This
process indicates a weak ignition. The H
2
/O
2
mixture in Meyer & Oppenheim (1971) is
stoichiometric and undiluted, which diers from the present work. More experiments and
simulations are needed to dene criteria that can determine whether ignition is strong or
weak for non-stoichiometric and diluted mixtures. Figure 4.26 shows the time evolution of
reactant (H
2
and O
2
) consumption rate, _ m
R
, and
ame front velocity, s, for case P within
the domain before the sponge layer. The consumption rate is normalized by the in
ow rate
of reactant _ m
R;1
, the
ame front velocity is normalized by the
uid
ow velocity upstream
of the shock, and the time is normalized by the ignition delay time. The time origin, t = 0,
is oset to the moment when ignition happens. The ignition delay time is dened as the
moment when the normalized reactant consumption rate exceeds 0:02. To automate the
calculation of the
ame front velocity, the streamwise location of the
ame front is dened
99
Figure 4.25: Time evolution of temperature contour of caseL(2; 1250; 0:3) fromt = 1:510
4
s
to t = 3:5 10
4
s
100
Figure 4.26: Normalized reactant consumption rate, _ m
R
= _ m
R;1
and
ame front velocity,s=u
1
,
for 1D laminar reactive STI before the sponge layer. Negative s implies that the
ame
propagates towards the inlet.
where the minimum
ow velocity occurs. Before ignition happens, it is hard to detect the
lowest velocity downstream of the shock, because the velocity prole is nearly constant.
Hence, the velocity of the
ame front is only measured after ignition (i.e., for t> 0).
The time evolution of the normalized reactant consumption rate shows ve stages. The
rst stage is the pre-ignition stage, where the consumption rate is nearly zero. Then the
consumption rate starts to increase, and at the same time, the
ame front starts to move
towards the outlet boundary. The third stage begins when the
ame front starts to move
towards the shock. During this stage, the
ame front moves through the partially reacted
mixture, but progressively less so, which explains why the reactant consumption rate in-
creases during this stage. The initial downstream movement of the ignition spot, and the
subsequent movement toward the shock is also conrmed in Melguizo-Gavilanes et al. (2011).
The fourth stage takes place as the
ame front catches up with the shock, which coincides
with rapid increases toward peaks of the consumption rate and the
ame front speed. After-
wards, the fth and nal stage is characterized by a rapid decrease of consumption rate and
ame speed due to the low pressure and temperature encountered by the
ame upstream of
the shock. In previous studies of DDT, the evolution of the
ame front velocity starts with
an initial exponential acceleration and is followed by an almost linear increase to supersonic
101
Figure 4.27: Normalized (a) reactant consumption rate and (b)
ame front speed from case
L to case P before the sponge layer, with linestyle the same as table 2.2. The three vertical
dotted lines correspond to t
1
, t
2
and t
3
, and the markers correspond to t
peak
for dierent
cases.
values (Valiev et al., 2009, 2010; Han et al., 2017). The rst stage described for the 1D
laminar case ts the tail of the exponential acceleration stage for DDT, and the second stage
for the 1D laminar case also shows a linear acceleration.
Figure 4.27 compares ensemble-averaged normalized consumption rate and
ame front
velocity from caseL to caseP . The time axis is normalized by the ignition delay time (t
ig;P
)
of caseP , the 1D laminar case, with the time origin (t = 0) oset to the moment of ignition
for the 1D laminar case. The streamwise location of the
ame front for the turbulent cases
is dened where the minimum transverse-averaged
ow velocity occurs. The trend for
ame
front speed is the same for all the cases, with a peak reached at the moment the
ame
passes through the shock. Four characteristic snapshots are marked in gure 4.27. Three
of them are marked by vertical dashed lines. From left to right, the rst one corresponds
to t=t
ig;P
= 0:05, when the ignition just happened, and the
ame front moves towards the
out
ow boundary. The second one corresponds to t=t
ig;P
= 0:075, when the
ame front
moves towards the shock. The third one corresponds tot=t
ig;P
= 1:5, when the
ame front is
about to exit the domain. These three snapshots are denoted as t
1
, t
2
and t
3
. There is one
more characteristic snapshot for these simulations, which is the time when the consumption
102
rate reaches its peak. This peak is marked by symbols of the same color as the linestyle for
each case in gure 4.27, and denoted as t
peak
. These four snapshots for each case will be
used to study the temporal evolution of the Takeno Flame Index (TFI) (Yamashita et al.,
1996), the heat release rate (HRR), and the combustion regime later.
Regarding the global consumption rate, turbulent cases start to increase from zero earlier
than the laminar 1D case, and higher M
t
will bring earlier departure from zero, which
coincides with our ndings that turbulence and, in particular, increased compressibility of
turbulence
uctuations lead to earlier ignition (section 3.2). The time sequence when these
cases reach peak consumption rates is the same as the sequence when their
ame fronts
catch up with the shock, which is within our expectation. Even though cases with M
t
= 0:1
are ignited earlier compared with the 1D laminar case, the reactant consumption rate of
the 1D laminar case increases much faster than the turbulent cases. Thus, the 1D laminar
case catches up with turbulent cases with M
t
= 0:1 shortly after ignition, and exceeds the
consumption rate of cases with M
t
= 0:3 at t=t
ig;P
= 0:4. This can also explain why the 1D
laminar case has the highest peak velocity, pressure and density of all the cases, together
with the facts that in the present simulations, combustion in the 1D laminar case occurs in
the fully-premixed regime, whereas partially-premixed combustion is present in the turbulent
cases under investigation.
Flame regime diagrams have been widely used to study combustion problems (Borghi,
1985; Poinsot, 1996; Quinlan et al., 2014; Nordin-Bates et al., 2017). The gures inxAppendix
B show
ame regime for the four characteristic snapshots dened above of each case in
a Da
t
vs Re
t
diagram, colored by the number of cells with heat release above threshold
(10
8
J s
1
m
3
), conditioned heat release rate (HRR) and conditioned Takeno Flame Index
(TFI). The heat release rate is dened as the rate of change of the enthalpy in the mixture:
HRR =
Ns
X
=1
W
_ !
h
: (4.8)
103
The Da
t
vs Re
t
diagrams shown include only regions of the domain with a HRR larger
than 10
8
J s
1
m
3
, which marks the region with more active reactions. TFI measures the
alignment between the gradient of mass fraction of fuel and the gradient of mass fraction of
oxygen:
TFI =
rY
H
2
rY
O
2
jrY
H
2
rY
O
2
j
: (4.9)
According to this denition, positive and negative unitary values of the TFI correspond to
fully premixed and non-premixed combustion, respectively.
Da
t
denotes turbulent Damk ohler number, which represents the ratio between the tur-
bulent time scale,
t
, and the chemical time scale,
c
:
Da
t
=
t
=
c
= (L
=u
rms
)=(l
c
=s
l
): (4.10)
L
=K
3=2
= is the turbulent length scale, where K is the turbulence kinetic energy (TKE),
and is the dissipation rate of TKE. u
rms
is the rms velocity. l
c
= =s
l
is the combustion
length scale, where = = is the kinematic viscosity, and s
l
is the laminar
ame speed.
In the present calculations, the laminar
ame speed, s
l
, is dened by s
2
l
= 2 _ !
H
2
=
1;H
2
,
according to a two-zone model for a laminar premixed
ame (Grogan et al., 2015; Turns,
2000), where is the thermal diusivity, _ !
H
2
is the local reaction rate of H
2
, and
1;H
2
is
the mean partial density of H
2
at the in
ow boundary. The turbulent Reynolds number is
dened as Re
t
=u
rms
l
t
= = (u
rms
=s
l
)(L
=l
c
).
Figure 4.28 compares the
ame regime diagram for casesO(2; 600; 0:1) andN(2; 600; 0:3).
The diagrams are divided into six regions: wrinkled
amelets (WF), corrugated
amelets
(CF), thin reaction zone (TRZ), broken reaction zone (BRZ), distributed reaction zone
(DRZ), and well stirred reactor (WSR). In WF, the laminar
ame speed is larger than
the turbulent velocity. This indicates that the turbulence is not strong enough to compete
with laminar speed of
ame propagation, so
ame propagation dominates over the
ame
corrugation brought by turbulence. In the regime of CF, the Kolmogorov scale, is smaller
104
Figure 4.28: Flame regime diagram presented for case O(2; 600; 0:1) (top) and case
N(2; 600; 0:3) (bottom) showing the joint distributions in theDa
t
-Re
t
plane of domain with
heat release larger than 10
8
J s
1
m
3
. The plot is colored by cell number for snapshots t
1
,
t
2
, t
peak
, t
3
from left to right. Regions with cell numbers below 50 are masked. The three
contour lines represent 80%, 40% and 20% of the peak cell number at the corresponding
time instance.
than l
c
. In this regime, the
ame is usually considered to be innitely thin. In TRZ, the
Kolmogorov length scale is smaller thanl
c
, but larger than the thickness of the reaction zone.
One characteristic of the TRZ regime is that eddies can penetrate into the preheat zone of
the
ame and enhance the mixing between the material in the preheat zone and unburnt
gas, which means that the
ame structure will go through instantaneous changes, and the
preheat zone will be broadened. One distinctive feature of BRZ is that the turbulence length
scale is smaller than the reaction zone. The eddies may penetrate into the reaction zone, and
break down chemical kinetics by transporting heat into the preheat zone. DRZ and WSR
are characterized by the large diusion time scale compared with the combustion time scale.
Dierent from WSR, the reactions in DRZ have turbulence scales smaller than the reaction
zone thickness, which leads to distributed and volumetric reaction zone. A comprehensive
overview of the dierent
ame regimes and their corresponding physical meaning can be
found in Peters (2000). Figure 4.28 shows that most reactions reside in TRZ. At t
1
, case O
has fewer number of points above the HRR threshold. This coincides with the fact that the
reactant consumption rate is low for caseO att
1
in gure 4.27. Att
2
andt
peak
, the case with
105
lower M
t
shows non-negligible region in the CF regime, while the diagram for the case with
higher M
t
does not reside in that CF regime, which means that lower M
t
shifts the regime
in the direction of higherDa
t
. We hypothesize that, since the ignition spot appears closer to
the shock for simulations with higher M
t
, these cases do not have enough time to develop a
high temperature when the
ame front catches up with the shock as the 1D laminar case and
2D turbulent cases with lower M
t
. Hence, this lower temperature when the
ame catches
up with the shock makes the combustion time scale larger, which leads to a smaller Da
t
. In
the TRZ regime, the preheat zone is broadened. In gure 4.23, a larger preheat zone can
be observed for simulations with larger M
t
, which is an indicator that the reactions in these
simulations reside more in the TRZ regime compared with simulations with lowerM
t
. Att
3
,
there is a drop inDa
t
because of the less vigorous combustion brought by the low temperature
and pressure in the unshocked mixture. In all the diagrams colored by cell numbers, there
is a single peak. The average of the peak is around the order of (Da
t
;Re
t
) O(10
2
; 10
5
),
which is of the same order as previous studies about the hydrogen/air HyShot II scramjet
engine (Nordin-Bates et al., 2017). The regime diagrams of simulations with higher Re
behave similarly in terms of the eects of M
t
and time evolution. There are no obvious
eects of Re
in the combustion regime diagram, because the dierence between the two
Re
considered is only a factor of approximately 2. Since case L(2; 1250; 0:3) and case
M(2; 1250; 0:1) do not provide more physical insights, the
ame regime diagrams of these
two cases are provided inxAppendix B. Figure 4.29 presents the time evolution of TFI
and HRR for case O(2; 600; 0:1). Similar observations can be made for all simulations, so
only one case is shown here. Again, the time evolution of TFI and HRR of other cases
for the dierent stages of ignition and
ame propagation can be found inxAppendix B. In
gure 4.29, HRR increases from t
1
to t
peak
, but decreases from t
peak
to t
4
, which follows
the same trend as the reactant consumption rate in gure 4.27. The same order of heat
release is also presented in many previous studies about combustion in scramjets, including
Moule et al. (2014), Wang et al. (2016), Zhang et al. (2016), and Nordin-Bates et al. (2017),
106
Figure 4.29: Flame regime diagram presented for case O(2; 600; 0:1) showing the joint dis-
tributions in the Da
t
-Re
t
plane of domain with heat release larger than 1e8 J/s/m
3
. The
plot is colored by TFI (top) and HRR (bottom) for snapshots t
1
, t
2
, t
peak
, t
3
from left to
right. Regions with cell numbers below 50 are masked. The three contour lines represent
80%, 40% and 20% of the peak cell number at the corresponding time instance.
Figure 4.30: Flame regime diagram for all the turbulent cases showing the joint distributions
in the Da
t
-Re
t
plane of domain with heat release larger than 10
8
J s
1
m
3
at t
peak
. The
plot is colored by TFI (top) and HRR (bottom). Regions with cell numbers below 50 are
masked. The three contour lines represent 80%, 40% and 20% of the peak cell number at
the corresponding time instance.
107
just to name a few. From t
1
to t
3
, TFI also indicates that the reaction becomes more and
more non-premixed. This is because, as the
ame propagates upstream, the species mass
fraction variances increase, which implies a more non-premixed nature of the mixture. At
t
3
, where non-premixed combustion and premixed combustion coexist, the premixed region
shows a higher HRR. This explains why the 1D laminar case has a faster increase in reactant
consumption rate, and why it has a larger peak of thermodynamic quantities. The contour
of HRR shows an increase fromt
1
tot
peak
, and a decrease fromt
peak
tot
3
, which is consistent
with the time evolution of reactant consumption rate in gure 4.27. Figure 4.30 compares
TFI and HRR of all the cases att
peak
. It can be observed that cases with lowerM
t
will lead
to a more non-premixed combustion and a higher HRR when the
ame front catches up with
the shock. The more non-premixed combustion in cases with lowerM
t
att
peak
indicates that
the decay of species mass concentration variance is slower with lower M
t
. The higher HRR
of cases with lower M
t
coincides with higher reactant consumption rate seen in gure 4.27.
In summary, we found out that higher M
t
will lead to earlier ignition behind the shock.
Flame fronts in simulations with higher M
t
catch up with the shock earlier than the 1D
laminar case, but cases with lowerM
t
catch up with the shock later than the 1D laminar case.
In the presence of turbulence, the transverse-averaged peak values of pressure, temperature
and velocity are lower than those in the 1D laminar case. Also, the reactant consumption
rate of the 1D laminar case increases faster than the turbulent cases. The
ame regime
diagram shows that most reactions happen in the TRZ regime. A lowerM
t
will also increase
the chance for the reaction to happen in the CF
ame regime slightly, which will lead to a
narrower preheat zone. The eect of Re
is not obvious for the range considered in these
simulations.
108
Chapter 5
Conclusion
In this study, DNS of passive scalar mixing in shock-turbulence interaction is performed
under dierent parameters. The DNS data includes Mach number (M) from 1.28 to 5,
turbulent Mach number (M
t
) from 0.1 to 0.4 and Taylor microscale Reynolds number (Re
)
from 40 to 70. These simulations cover relatively high M, M
t
,, Re
and also interactions in
the broken-shock regime, which were not addressed in previous literature in the context of
mixing in STI.
The analysis of the budget equation of scalar variance reveals that the dissipation term
dominates, dilatational term and turbulent diusion terms are not negligible, while all other
terms are negligible. The eective jump of scalar dissipation rate across the shock mono-
tonically increases with M, and decreases with M
t
and Re
, showing no sign of saturation
in the range of M = 5 considered in the present study, contrary with the amplication of
TKE, which is known to saturate forM 3. A scalar Taylor-like microscale, representative
of mixing at an intermediate length scale between dissipation and energy-containing eddies
is also investigated for dierent cases. There is a decrease of this scale across the shock, and
this decrease is larger for increasing M and decreasing Re
. The budget equation of scalar
dissipation rate is also explored. It is observed that velocity-scalar-gradient and molecu-
lar diusion terms dominate, both terms being amplied in magnitude across the shock, a
trend that increases monotonically with larger M, smaller M
t
and Re
. It is interesting to
109
note that the jump of the velocity-scalar-gradient interaction term seems to saturate around
M = 3, coinciding with the saturation of TKE.
The PDFs of invariants (p;q;r) of the velocity gradient tensor and the corresponding
conditioned SDR distribution are studied. At all streamwise locations, the teardrop shape
of the PDF in the qr plane, prevalent in many turbulent
ow types, is observed for p = 0.
The shape is narrowed across the shock. Conditioned SDR concentrates on the right lower
tail (q < 0 and r > 0) of the teardrop, which conrms that UNSS is the most dissipative
topology. At p > 0 (contraction) and p < 0 (extension), the shape of the joint qr-PDFs
becomes more symmetrized. The proportion of topologies and their contribution to the
scalar dissipation are also investigated. It is found that SFS is the most dominant topology
both before and after the shock. Downstream of the unsteady shock region, the proportion
of SFS and UNSS decreases, whereas that of UFS and SNSS increases. However, the relative
order of these topologies does not change across the shock. UNSS contributes the largest
proportion of scalar dissipation. Fully compressing (SFC) or fully stretching (UFS) topologies
are less dissipative than those with simultaneous compressing and stretching (SFS and UFC).
Immediately after the shock, UNSS and SFS topologies reduce their contribution to the scalar
dissipation, whereas all other topologies increase theirs.
Alignments between scalar gradient and the three eigenvectors of strain rate, vorticity and
unit streamwise vector conditioned on dierent
ow topologies at preshock and postshock
locations are compared. In the scalar dissipation rate budget equation, the velocity-scalar-
gradient term is related to the alignments of the principal strain directions and the passive
scalar gradient, which is of importance for structure-based SGS models of passive scalar
transport. The alignments between scalar gradient and the most extensive () and the most
compressive (
) strain-rate eigendirections shows a large variation conditioned on dierent
topologies, while the alignments between scalar gradient and the intermediate strain-rate
eigendirection (), vorticity and unit streamwise vector do not show an obvious dependence
on topology. Among all the topologies, UNSS is the most aligned with
, and least aligned
110
with , and SFC is the least aligned with
, and most aligned with . The changes of
the alignments across the shock show a negligible variation with M
t
, but larger changes are
observed with higherM and lowerRe
. In general, the scalar gradient becomes more aligned
with and less aligned with
across the shock, which is also conrmed by the barycentric
map representation of these alignments. Also, across the shock, scalar gradient becomes more
aligned with the unit streamwise vector, and more orthogonal with and vorticity. Also, the
barycentric map representation shows that vorticity becomes more aligned with across the
shock. All these changes are more obvious with higherM and lowerRe
. Another important
observation is that the most dissipative alignment and the most preferred alignment of scalar
gradient become better matched across the shock, which explains why there is an increase
of SDR and better mixing across the shock.
Additionally, it is also found that the Schmidt number, which ranges from 0.5 to 2 in
these simulations, has little eects on all the quantities discussed in this work.
Four reactive 2D STI simulations are conducted with initialRe
= (1250; 600) andM
t
=
(0:1; 0:3) and M = 2 with a mixture of H
2
, O
2
, and Ar, and compared with a laminar 1D
simulation at the same Mach number. Three stages are found in these simulations. In the
rst stage, ignition spot is formed downstream of the shock, and is transported towards
the out
ow boundary. In the second stage, combustion becomes more and more vigorous,
and the
ame front starts to move upstream towards the shock, and catches up with the
shock in the end. In the last stage, after the
ame front passes the shock, it keeps moving
towards the inlet boundary. The study of time evolution of streamwise proles of dierent
thermodynamic quantities reveals that higher M
t
will lead to earlier ignition behind the
shock. Flame fronts in simulations with higher M
t
catch up with the shock earlier than the
1D laminar case, but cases with lowerM
t
catch up with the shock later than the 1D laminar
case. In the presence of turbulence, the peak values of pressure, temperature and velocity
are lower compared with the laminar case. The analysis of the reactant consumption rate
and the speed of the
ame front shows that the consumption rate of the 1D laminar case
111
increases faster than in the turbulent cases. The study of
ame regime of the 2D reactive
STI simulations indicates that that a higherM
t
will make reactions more likely to happen in
the TRZ regime, which will lead to a wider preheat zone. The eects of Re
is not obvious
in these simulations.
In summary, this work extends the DNS of passive scalar mixing in STI setting to a wider
parametric space, with the eects of both M
t
and Re
on mixing enhancement across the
shock. The mechanism behind shock-induced mixing is explained with alignments and
ow
topologies. Also, DNS of shock-induced ignition in the STI setting is carried out. This work
reveals the eects of turbulence and compressibility on shock-induced ignition, compared
with previous 1D laminar simulations. The understanding of turbulence eects on shock-
induced ignition and the database generated in this work can be used to develop better
reduced-order models for future scramjet simulations.
112
Bibliography
Agui, J. H., Briassulis, G. & Andreopoulos, Y. 2005 Studies of interactions of a
propagating shock wave with decaying grid turbulence: velocity and vorticity elds. J.
Fluid Mech. 524, 143{195.
Andreopoulos, Y. 2008 Vorticity and velocity alignment in compressible
ows: An exper-
imental study of helicity density in turbulence and vortices. Russ. J. Electrochem. 44 (4),
390.
Andreopoulos, Y., Agui, J. H. & Briassulis, G. 2000 Shock wave-turbulence inter-
actions. Annu. Rev. Fluid Mech. 32 (1), 309{345.
Andrews, E. H. 1994 NASA's Hypersonic Research Engine Project: A Review. Tech. Rep.
TM-107759. National Aeronautics and Space Administration, Langley Research Center.
Arun, S., Sameen, A., Srinivasan, B. & Girimaji, S. S. 2019 Topology-based char-
acterization of compressibility eects in mixing layers. J. Fluid Mech. 874, 38{75.
Ashurst, W. T., Kerstein, A. R, Kerr, R. M. & Gibson, C. H. 1987 Alignment
of vorticity and scalar gradient with strain rate in simulated Navier{Stokes turbulence.
Phys. Fluids 30 (8), 2343{2353.
Axdahl, E. L. 2013 A study of premixed, shock-induced combustion with application to
hypervelocity
ight. PhD thesis, Georgia Institute of Technology.
Barnes, F. W. & Segal, C. 2015 Cavity-based
ameholding for chemically-reacting
supersonic
ows. Prog. Aerosp. Sci. 76, 24{41.
113
Barre, S., Alem, D. & Bonnet, J. P. 1996 Experimental study of a normal
shock/homogeneous turbulence interaction. AIAA J. 34 (5).
Barthelemy, R. R. 1989 Recent progress in the National Aero-Space Plane Program.
IEEE Aero. El. Sys. Mag. 4 (5), 3{12.
Barthelemy, R. R. 1991 The National Aero-Space Plane Program-a revolutionary con-
cept. Tech. Rep.. SAE Technical Paper.
Batchelor, G. K. 1969 Computation of the energy spectrum in homogeneous two-
dimensional turbulence. Phys. Fluids 12 (12), II{233.
Bermejo-Moreno, I., Bodart, J., Larsson, J., Barney, B. M., Nichols, J. W. &
Jones, S. 2013a Solving the compressible Navier-Stokes equations on up to 1.97 million
cores and 4.1 trillion grid points. In Proc. Intl. Conf. on High Performance Computing,
Networking, Storage and Analysis, p. 62:1{62:10. ACM.
Bermejo-Moreno, I., Larsson, J., Bodart, J. & Vicquelin, R. 2013b Wall-modeled
large-eddy simulations of the HIFiRE-2 scramjet. Center for Turbulence Research, Annual
Research Briefs pp. 3{19.
Bermejo-Moreno, I., Larsson, J. & Lele, S. K. 2010 LES of canonical shock-
turbulence interaction. Center for Turbulence Research, Annual Research Briefs pp. 209{
222.
Billet, G., Giovangigli, V. & De Gassowski, G. 2008 Impact of volume viscosity on
a shock hydrogen-bubble interaction. Combust. Theory Model. 12 (2), 221{248.
Borghi, R. 1985 On the structure and morphology of turbulent premixed
ames. In Recent
advances in the Aerospace Sciences, pp. 117{138. Springer.
Boukharfane, R., Bouali, Z. & Mura, A. 2018 Evolution of scalar and velocity dy-
namics in planar shock-turbulence interaction. Shock Waves pp. 1{25.
114
Brachet, M. E., Meneguzzi, M., Politano, H. & Sulem, P. L. 1988 The dynamics
of freely decaying two-dimensional turbulence. J. Fluid Mech. 194, 333{349.
Braun, N. O., Pullin, D. I. & Meiron, D. I. 2019 Large eddy simulation investigation
of the canonical shock{turbulence interaction. J. Fluid Mech. 858, 500{535.
Brethouwer, G., Hunt, J. C. R. & Nieuwstadt, F. T. M. 2003 Micro-structure and
Lagrangian statistics of the scalar eld with a mean gradient in isotropic turbulence. J.
Fluid Mech. 474, 193{225.
Brindise, M.C. & Vlachos, P.P. 2017 Proper orthogonal decomposition truncation
method for data denoising and order reduction. Exp. Fluids 58 (4), 28.
Brouillette, M. 2002 The Richtmyer-Meshkov instability. Annu. Rev. Fluid Mech. 34 (1),
445{468.
Cambon, C., Coleman, G. N. & Mansour, N. N. 1993 Rapid distortion analysis and
direct simulation of compressible homogeneous turbulence at nite Mach number. J. Fluid
Mech. 257, 641{665.
Cancino, L. R., Fikri, M., Oliveira, A. A. M. & Schulz, C. 2010 Measurement
and chemical kinetics modeling of shock-induced ignition of ethanol-air mixtures. Energy
& Fuels 24 (5), 2830{2840.
Chaos, M. & Dryer, F. L. 2010 Chemical-kinetic modeling of ignition delay: Consider-
ations in interpreting shock tube data. Int. J. Chem. Kinet. 42 (3), 143{150.
Chasnov, J. R. 1997 On the decay of two-dimensional homogeneous turbulence. Phys.
Fluids 9 (1), 171{180.
Chen, C. 2018 Shock-turbulence interactions at high turbulence intensities: Theory and
direct numerical simulations. PhD thesis, Texas A&M University.
115
Chen, C. H. & Donzis, D. A. 2019 Shock{turbulence interactions at high turbulence
intensities. J. Fluid Mech. 870, 813{847.
Cheng, R. K. & Oppenheim, A. K. 1984 Autoignition in methane-hydrogen mixtures.
Combust.
ame 58 (2), 125{139.
Chi, C., Abdelsamie, A. & Th evenin, D. 2018 Direct numerical simulations of hotspot-
induced ignition in homogeneous hydrogen-air pre-mixtures and ignition spot tracking.
Flow Turbul. Combus. 101 (1), 103{121.
Chong, M. S., Perry, A. E. & Cantwell, B. J. 1990 A general classication of
three-dimensional
ow elds. Phys. Fluids 2 (5), 765{777.
Cohen, S. D., Hindmarsh, A. C. & Dubois, P. F. 1996 CVODE, a sti/nonsti ODE
solver in C. Comput. Phys. 10 (2), 138{143.
Corrsin, S. 1951 On the spectrum of isotropic temperature
uctuations in an isotropic
turbulence. J. Appl. Phys. 22 (4), 469{473.
Danaila, L. & Antonia, R. A. 2009 Spectrum of a passive scalar in moderate Reynolds
number homogeneous isotropic turbulence. Phys. Fluids 21 (11), 111702.
Danish, M. & Meneveau, C. 2018 Multiscale analysis of the invariants of the velocity
gradient tensor in isotropic turbulence. Phys. Rev. Fluids 3 (4), 044604.
Danish, M., Suman, S. & Girimaji, S. S. 2016 In
uence of
ow topology and dilatation
on scalar mixing in compressible turbulence. J. Fluid Mech. 793, 633{655.
David, G. G., Raymond, L. S., Harry, K. M. & Bryan, W. W. 2018 Cantera:
An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport
processes. https://www.cantera.org, version 2.4.0.
Davidson, P. A. 2015 Turbulence: An Introduction for Scientists and Engineers. Oxford
University Press.
116
De Avillez, M. A. & Mac Low, M. 2002 Mixing timescales in a supernova-driven
interstellar medium. Astrophys. J. 581 (2), 1047.
Desjardins, T., Merritt, E. C., Flippo, K. A., Kline, J. L., Di Stefano, C.,
DeVolder, B. G., Doss, F. W., Fierro, F., Randolph, R. B., Schmidt, D. W.
& others 2018 Mshock, a thin layer Richtmyer-Meshkov instability experiment. Tech.
Rep.. Los Alamos National Lab.(LANL), Los Alamos, NM (United States).
Donzis, D. A. 2012a Amplication factors in shock-turbulence interactions: Eect of shock
thickness. Phys. Fluids 24 (1), 011705.
Donzis, D. A. 2012b Shock structure in shock-turbulence interactions. Phys. Fluids 24 (12),
126101.
Dritschel, D. G., Tran, C. V. & Scott, R. K. 2007 Revisiting Batchelor's theory of
two-dimensional turbulence. J. Fluid Mech. 591, 379{391.
Ducros, F., Laporte, F., Souleres, T., Guinot, V., Moinat, P. & Caruelle, B.
2000 High-order
uxes for conservative skew-symmetric-like schemes in structured meshes:
application to compressible
ows. J. Comput. Phys. 161 (1), 114{139.
Dziemi nska, E. & Hayashi, A. K. 2013 Auto-ignition and DDT driven by shock wave{
boundary layer interaction in oxyhydrogen mixture. Int. J. Hydrog. Energy 38 (10), 4185{
4193.
Echekki, T. & Chen, J. H. 2002 High-temperature combustion in autoigniting non-
homogeneous hydrogen/air mixtures. Proc. Combust. Inst. 29 (2), 2061{2068.
Echekki, T. & Chen, J. H. 2003 Direct numerical simulation of autoignition in non-
homogeneous hydrogen-air mixtures. Combust. Flame 134 (3), 169{191.
Fedkiw, R. P., Merriman, B. & Osher, S. 1997 High accuracy numerical methods for
thermally perfect gas
ows with chemistry. J. Comput. Phys. 132 (2), 175{190.
117
Fern andez, D. C. D. R., Hicken, J. E. & Zingg, D. W. 2014 Review of summation-
by-parts operators with simultaneous approximation terms for the numerical solution of
partial dierential equations. Comput. Fluids 95, 171{196.
Figueroa-Labastida, M. & Farooq, A. 2020 Simultaneous lateral and endwall high-
speed visualization of ignition in a circular shock tube. Combust. Flame 214, 263{265.
Fiscaletti, D., Elsinga, G. E., Attili, A., Bisetti, F. & Buxton, O. R. H. 2016
Scale dependence of the alignment between strain rate and rotation in turbulent shear
ow. Phys. Rev. Fluids 1 (6), 064405.
Fox, S. & Davidson, P. A. 2010 Freely decaying two-dimensional turbulence. J. Fluid
Mech. 659, 351{364.
Freeman, Jr. D. C., McClinton, C. R., Rausch, V. L. & Crawford, J. L. 1997
The NASA Hyper-X Program. Tech. Rep.. NASA Langley Technical Report Server.
Gamba, M. & Mungal, M. G. 2015 Ignition,
ame structure and near-wall burning in
transverse hydrogen jets in supersonic cross
ow. J. Fluid Mech. 780, 226{273.
Golub, V. V., Baklanov, D. I., Bazhenova, T. V., Bragin, M. V., Golovastov,
S. V., Ivanov, M. F. & Volodin, V. V. 2007 Shock-induced ignition of hydrogen
gas during accidental or technical opening of high-pressure tanks. J. Loss Prevent Proc.
20 (4-6), 439{446.
Gonzalez, M. 2012 Analysis of scalar dissipation in terms of vorticity geometry in isotropic
turbulence. J. Turbul. 13 (1), N41.
Griffond, J. 2005 Linear interaction analysis applied to a mixture of two perfect gases.
Phys. Fluids 17 (8), 086101.
Griffond, J. 2006 Linear interaction analysis for Richtmyer-Meshkov instability at low
Atwood numbers. Phys. Fluids 18 (5), 054106.
118
Grogan, K. P., Goldsborough, S. S. & Ihme, M. 2015 Ignition regimes in rapid
compression machines. Combust. Flame 162 (8), 3071{3080.
Grogan, K. P. & Ihme, M. 2015 Weak and strong ignition of hydrogen/oxygen mixtures
in shock-tube systems. Proc. Combust. Inst. 35 (2), 2181{2189.
Gruber, M. R., Donbar, J. M., Carter, C. D. & Hsu, K. Y. 2004 Mixing and
combustion studies using cavity-based
ameholders in a supersonic
ow. Combust. Flame
20 (5), 769{778.
Hallion, R. P., Becker, John., Vitalli, J. & Young, J. 1987 The hypersonic rev-
olution, Vol. 2: From scramjet to the National Aero-Space Plane. Tech. Rep.. Wright
Patterson Airf. Base, OH: Aeronaut. Syst. Div.
Han, W., Gao, Y. & Law, C. K. 2017 Flame acceleration and de
agration-to-detonation
transition in micro-and macro-channels: An integrated mechanistic study. Combust. Flame
176, 285{298.
Hannappel, R. & Friedrich, R. 1995 Direct numerical simulation of a Mach 2 shock
interacting with isotropic turbulence. Appl. Sci. Res. 54 (3), 205{221.
Hass, Neal, Cabell, Karen, Storch, Andrea & Gruber, Mark 2011 HIFiRE
direct-connect rig (HDCR) phase I scramjet test results from the NASA Langley Arc-
Heated Scramjet Test Facility. In 17th AIAA International Space Planes and Hypersonic
Systems and Technologies Conference, p. 2248.
Hejranfar, K. & Rahmani, S. 2019 An assessment of shock-disturbances interaction
considering real gas eects. J. Fluids Eng. 141 (1), 011201.
Hesselink, L. & Sturtevant, B. 1988 Propagation of weak shocks through a random
medium. J. Fluid Mech. 196, 513{553.
119
Hickel, S., Egerer, C. P. & Larsson, J. 2014 Subgrid-scale modeling for implicit
large eddy simulation of compressible
ows and shock-turbulence interaction. Phys. Fluids
26 (10), 106101.
Honkan, A. & Andreopoulos, J. 1990 Experiments in a shock wave/homogeneous
turbulence interaction. In 21st Fluid Dynamics, Plasma Dynamics and Lasers Conference,
p. 1647. AIAA.
Honkan, A. & Andreopoulos, J. 1992 Rapid compression of grid-generated turbulence
by a moving shock wave. Phys. Fluids 4 (11), 2562{2572.
Houim, R. W. & Kuo, K. K. 2011 A low-dissipation and time-accurate method for
compressible multi-component
ow with variable specic heat ratios. J. Comput. Phys.
230 (23), 8527{8553.
Huang, J., Bushe, W. K., Hill, P. G. & Munshi, S. R. 2006 Experimental and
kinetic study of shock initiated ignition in homogeneous methane{hydrogen{air mixtures
at engine-relevant conditions. Int. J. Chem. Kinet. 38 (4), 221{233.
Huete, C esar, Jin, Tai, Mart nez-Ruiz, Daniel & Luo, Kun 2017 Interaction of a
planar reacting shock wave with an isotropic turbulent vorticity eld. Phys. Rev. E 96 (5),
053104.
Huete, C., S anchez, A. L. & Williams, F. A. 2013 Theory of interactions of thin
strong detonations with turbulent gases. Phys. Fluids 25 (7), 076105.
Huete, C., S anchez, A. L. & Williams, F. A. 2014 Linear theory for the interaction
of small-scale turbulence with overdriven detonations. Phys. Fluids 26 (11), 116101.
Huete, C., Velikovich, A. L. & Wouchuk, J. G. 2011 Analytical linear theory for
the interaction of a planar shock wave with a two-or three-dimensional random isotropic
density eld. Phys. Rev. E 83 (5), 056320.
120
Hughes, J. P., Rakowski, C. E., Burrows, D. N. & Slane, P. O. 1999 Nucleosyn-
thesis and mixing in Cassiopeia A. Astrophys. J. 528 (2), L109.
Ihme, M., Sun, Y. & Deiterding, R. 2013 Detailed simulations of shock-bifurcation
and ignition of an argon-diluted hydrogen/oxygen mixture in a shock tube. In 51st AIAA
Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition,
p. 538.
Im, H. G., Chen, J. H. & Law, C. K. 1998 Ignition of hydrogen-air mixing layer in
turbulent
ows. Proc. Combust. Inst. 27 (1), 1047{1056.
Jackson, K. R., Gruber, M. R., & Buccellato, S. 2015 Mach 6{8+ hydrocarbon-
fueled scramjet
ight experiment: the HIFiRE Flight 2 project. J. Propul. Power 31 (1),
36{53.
Jackson, T. L., Hussaini, M. Y. & Ribner, H. S. 1993 Interaction of turbulence with
a detonation wave. Phys. Fluids 5 (3), 745{749.
Jackson, T. L., Kapila, A. K. & Hussaini, M. Y. 1990 Convection of a pattern of
vorticity through a reacting shock wave. Phys. Fluids 2 (7), 1260{1268.
Jacquin, L., Cambon, C. & Blin, E. 1993 Turbulence amplication by a shock wave
and rapid distortion theory. Phys. Fluids 5 (10), 2539{2550.
Jahanbakhshi, R. & Madnia, C. K. 2016 Entrainment in a compressible turbulent shear
layer. J. Fluid Mech. 797, 564{603.
Jamme, S., Cazalbou, J., Torres, F. & Chassaing, P. 2002 Direct numerical simu-
lation of the interaction between a shock wave and various types of isotropic turbulence.
Flow Turbul. Combust. 68 (3), 227{268.
Jimenez, J. 1992 Kinematic alignment eects in turbulent
ows. Phys. Fluids 4 (4), 652{
654.
121
Johnsen, Eric, Larsson, Johan, Bhagatwala, Ankit V, Cabot, William H,
Moin, Parviz, Olson, Britton J, Rawat, Pradeep S, Shankar, Santhosh K,
Sj ogreen, Bj orn, Yee, Helen C & others 2010 Assessment of high-resolution meth-
ods for numerical simulations of compressible turbulence with shock waves. J. Comput.
Phys. 229 (4), 1213{1237.
Jones, W. P., Navarro-Martinez, S. & R ohl, O. 2007 Large eddy simulation of
hydrogen auto-ignition with a probability density function method. Proc. Combust. Inst.
31 (2), 1765{1771.
Karimi, M. & Girimaji, S. S. 2016 Suppression mechanism of Kelvin-Helmholtz insta-
bility in compressible
uid
ows. Phys. Rev. E 93 (4), 041102.
Kee, R. J., Coltrin, M. E. & Glarborg, P. 2005 Chemically reacting
ow: theory
and practice. John Wiley & Sons.
Kerr, R. M. 1985 Higher-order derivative correlations and the alignment of small-scale
structures in isotropic numerical turbulence. J. Fluid Mech. 153, 31{58.
Kerrebrock, J. L. 1956 The interaction of
ow discontinuities with small disturbances in
a compressible
uid. PhD thesis, California Institute of Technology.
Khokhlov, A., Austin, J. & Knisely, A. 2015 Development of hot spots and ignition
behind re
ected shocks in 2H
2
+O
2
. In Proceedings of the 25th International Colloquium
on the Dynamics of Explosions and Reactive Systems, p. 20.
Kim, Y. R., Lee, H. J., Kim, S. & Jeung, I. 2013 A
ow visualization study on self-
ignition of high pressure hydrogen gas released into a tube. Proc. Combust. Inst. 34 (2),
2057{2064.
Kitamura, T., Nagata, K., Sakai, Y., Sasoh, A. & Ito, Y. 2016 Rapid distortion
theory analysis on the interaction between homogeneous turbulence and a planar shock
wave. J. Fluid Mech. 802, 108{146.
122
Kiverin, A. D. & Yakovenko, I. S. 2019 Ignition and detonation onset behind incident
shock wave in the shock tube. Combust. Flame 204, 227{236.
Kolmogorov, A. N. 1941 The local structure of turbulence in incompressible viscous
uid
for very large Reynolds numbers. Cr Acad. Sci. URSS 30, 301{305.
Kov asznay, L. S. G. 1953 Turbulence in supersonic
ows. J. Aeronaut. Sci. 20, 657{674.
Kraichnan, R. H. 1967 Inertial ranges in two-dimensional turbulence. Phys. Fluids 10 (7),
1417{1423.
Larsson, J. 2009 Blending technique for compressible in
ow turbulence: algorithm local-
ization and accuracy assessment. J. Comput. Phys. 228 (4), 933{937.
Larsson, J., Bermejo-Moreno, I. & Lele, S. K. 2013 Reynolds- and Mach-number
eects in canonical shock{turbulence interaction. J. Fluid Mech. 717, 293{321.
Larsson, Johan & Gustafsson, Bertil 2008 Stability criteria for hybrid dierence
methods. J. Comput. Phys. 227 (5), 2886{2898.
Larsson, J. & Lele, S. K. 2009 Direct numerical simulation of canonical shock/turbulence
interaction. Phys. Fluids 21 (12), 126101.
Lee, S., Lele, S. K. & Moin, P. 1993 Direct numerical simulation of isotropic turbulence
interacting with a weak shock wave. J. Fluid Mech. 251, 533{562.
Lee, S., Lele, S. K. & Moin, P. 1994 Corrigendum: Direct numerical simulation of
isotropic turbulence interacting with a weak shock wave. J. Fluid Mech. 264, 373{374.
Lee, S., Lele, S. K. & Moin, P. 1997 Interaction of isotropic turbulence with shock
waves: eect of shock strength. J. Fluid Mech. 340, 225{247.
Lele, S. K. 1992 Shock-jump relations in a turbulent
ow. Phys. Fluids 4 (12), 2900{2905.
Leyva, I. 2017 The relentless pursuit of hypersonic
ight. Physics Today 70 (11), 30{36.
123
Lipkowicz, J. T., Wlokas, I. & Kempf, A. M. 2018 Analysis of mild ignition in
a shock tube using a highly resolved 3D-LES and high-order shock-capturing schemes.
Shock Waves pp. 1{11.
Livescu, D. & Ryu, J. 2016 Vorticity dynamics after the shock{turbulence interaction.
Shock Waves 26 (3), 241{251.
Lombardini, M., Pullin, D. I. & Meiron, D. I. 2012 Transition to turbulence in
shock-driven mixing: a Mach number study. J. Fluid Mech. 690, 203{226.
Lowe, A. J. & Davidson, P. A. 2005 The evolution of freely-decaying, isotropic, two-
dimensional turbulence. Eur. J. Mech. B Fluids 24 (3), 314{327.
Lu, Z., Zhou, H., Li, S., Ren, Z., Lu, T. & Law, C. K. 2017 Analysis of operator
splitting errors for near-limit
ame simulations. J. Comput. Phys. 335, 578{591.
Luo, Z., Yoo, C. S., Richardson, E. S., Chen, J. H., Law, C. K. & Lu, T. 2012
Chemical explosive mode analysis for a turbulent lifted ethylene jet
ame in highly-heated
co
ow. Combust. Flame 159 (1), 265{274.
M., Obukhov A. 1949 Structure of temperature eld in turbulent
ow. Ser. Geogr. i Gooz.
XIII (1), 58{69.
Ma, F., Li, J., Yang, V., Lin, K. & Jackson, T. 2005 Thermoacoustic
ow instability
in a scramjet combustor. In 41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference
& Exhibit, p. 3824.
Mahesh, K., Lele, S. K. & Moin, P. 1997 The in
uence of entropy
uctuations on the
interaction of turbulence with a shock wave. J. Fluid Mech. 334, 353{379.
Marble, F. E., Hendricks, G. J. & Zukoski, E. E. 1987 Progress toward shock
enhancement of supersonic combustion processes. In 23rd AIAA, SAE, ASME, and ASEE,
Joint Propulsion Conference, pp. 1{8. San Diego, CA.
124
Mathieu, O., Pemelton, J. M., Bourque, G. & Petersen, E. L. 2015 Shock-induced
ignition of methane sensitized by NO
2
and N
2
O. Combust. Flame 162 (8), 3053{3070.
McGuire, J. R. 2007 Ignition enhancement for scramjet combustion. PhD thesis, Univer-
sity of New South Wales.
McMillin, B.K., Lee, M.P., Paul, P.H. & Hanson, R.K. 1991 Planar laser-induced
uorescence imaging of shock-induced ignition. Symposium (International) on Combustion
23 (1), 1909{1914.
McWilliams, J. C 1990 The vortices of two-dimensional turbulence. J. Fluid Mech. 219,
361{385.
Melguizo-Gavilanes, J., Rezaeyan, N., Tian, M. & Bauwens, L. 2011 Shock-
induced ignition with single step Arrhenius kinetics. Int. J. Hydrog. Energy 36 (3), 2374{
2380.
Menon, S. 1989 Shock-wave-induced mixing enhancement in scramjet combustors. AIAA
Paper 89-0104 .
Merkel, A. C. & Ciccarelli, G. 2020 Visualization of lean methane-air ignition behind
a re
ected shock wave. Fuel 271, 117617.
Meyer, J. W. & Oppenheim, A. K. 1971 On the shock-induced ignition of explosive
gases. Proc. Combust. Inst. 13 (1), 1153{1164.
Micka, D. J. 2010 Combustion stabilization, structure, and spreading in a laboratory
dual-mode scramjet combustor. PhD thesis, University of Michigan.
Micka, D. J. & Driscoll, J. F. 2009 Combustion characteristics of a dual-mode scramjet
combustor with cavity
ameholder. Proc. Combust. Inst. 32 (2), 2397{2404.
Misra, A. & Pullin, D. I. 1997 A vortex-based subgrid stress model for large-eddy
simulation. Phys. Fluids 9 (8), 2443{2454.
125
Mogi, T., Kim, D., Shiina, H. & Horiguchi, S. 2008 Self-ignition and explosion during
discharge of high-pressure hydrogen. J. Loss Prevent Proc. 21 (2), 199{204.
Moore, F. K. 1953 Unsteady oblique interaction of a shock wave with a plane disturbances.
Tech. Rep.. NACA Tech. Rep. 1165.
Morris, C., Kamel, M., Stouklov, I. & Hanson, R. 1996 PLIF imaging of the super-
sonic reactive
ows around projectiles in an expansion tube. In 34th Aerospace Sciences
Meeting and Exhibit, p. 855.
Moule, Y., Sabel'nikov, V., Mura, A. & Smart, M. 2014 Computational
uid dy-
namics investigation of a mach 12 scramjet engine. J. Propul. Power 30 (2), 461{473.
Ni, Q. 2015 Compressible turbulent mixing: Eects of Schmidt number. Phys. Rev. E 91 (5),
053020.
Ni, Q. 2016 Compressible turbulent mixing: Eects of compressibility. Phys. Rev. E 93 (4),
043116.
Ni, Q., Shi, Y. & Chen, S. 2015 A numerical investigation on active and passive scalars
in isotropic compressible turbulence. arXiv:1505.02685 .
Nordin-Bates, K., Fureby, C., Karl, S. & Hannemann, K. 2017 Understanding
scramjet combustion using les of the HyShot II combustor. Proc. Combust. Inst. 36 (2),
2893{2900.
O'Byrne, S., Stotz, I., Neely, A., Boyce, R., Mudford, N. & Houwing, F.
2005 OH PLIF imaging of supersonic combustion using cavity injection. In AIAA/CIRA
13th International Space Planes and Hypersonics Systems and Technologies Conference,
p. 3357.
126
O'Neill, P. & Soria, J. 2005 The relationship between the topological structures in
turbulent
ow and the distribution of a passive scalar with an imposed mean gradient.
Fluid Dyn. Res. 36 (3), 107.
Oppenheim, A. K., Cohen, L. M., Short, J. M., Cheng, R. K. & Hom, K. 1975
Dynamics of the exothermic process in combustion. Proc. Combust. Inst. 15 (1), 1503{
1513.
Orlicz, G. C., Balasubramanian, S., Vorobieff, P. & Prestridge, K. P. 2015
Mixing transition in a shocked variable-density
ow. Phys. Fluids 27 (11), 114102.
Oxlade, A.R., Valente, P.C., Ganapathisubramani, B. & Morrison, J.F. 2012
Denoising of time-resolved PIV for accurate measurement of turbulence spectra and re-
duced error in derivatives. Exp. Fluids 53 (5), 1561{1575.
Pan, L. & Scannapieco, E. 2010 Mixing in supersonic turbulence. Astrophys. J. 721 (2),
1765.
Pan, L. & Scannapieco, E. 2011 Passive scalar structures in supersonic turbulence. Phys.
Rev. E 83 (4), 045302.
Pang, G. A., Davidson, D. F. & Hanson, R. K. 2009 Experimental study and modeling
of shock tube ignition delay times for hydrogen{oxygen{argon mixtures at low tempera-
tures. Proc. Combust. Inst. 32 (1), 181{188.
Pantano, C. & Sarkar, S. 2002 A study of compressibility eects in the high-speed
turbulent shear layer using direct simulation. J. Fluid Mech. 451, 329{371.
Parashar, N., Sinha, S. S. & Srinivasan, B. 2019 Lagrangian investigations of velocity
gradients in compressible turbulence: lifetime of
ow-eld topologies. J. Fluid Mech. 872,
492{514.
127
Peebles, C. 2009 Learning from experience: case studies of the Hyper-X project. In 47th
AIAA Aerospace Sciences Meeting including The New Horizons Forum and Aerospace
Exposition, p. 1523.
Penyazkov, O. G., Ragotner, K. A., Dean, A. J. & Varatharajan, B. 2005
Autoignition of propane{air mixtures behind re
ected shock waves. Proc. Combust. Inst.
30 (2), 1941{1947.
Pepe, R., Bonfiglioli, A., D
0
Angola, A., Colonna, G. & Paciorri, R. 2014 Shock-
tting versus shock-capturing modeling of strong shocks in nonequilibrium plasmas. IEEE
T. Plasma. Sci. 42 (10), 2526{2527.
Peters, Norbert 2000 Turbulent Combustion. Cambridge University Press.
Piland, W. M. 1987 Technology challenges for the National Aero-Space Plane. Tech. Rep..
NASA National Aero-Space Plane Program Management Oce.
Pirozzoli, S. 2002 Conservative hybrid compact-WENO schemes for shock-turbulence in-
teraction. J. Comput. Phys. 178 (1), 81{117.
Pirozzoli, S. & Grasso, F. 2004 DNS of isotropic compressible turbulence: in
uence of
compressibility on dynamics and structures. Phys. Fluids 16 (12), 4386{4407.
Poinsot, Thierry 1996 Using direct numerical simulations to understand premixed tur-
bulent combustion. Symposium (International) on Combustion 26 (1), 219{232.
Pullin, D. I. 2000 A vortex-based model for the subgrid
ux of a passive scalar. Phys.
Fluids 12 (9), 2311{2319.
Quadros, R. & Sinha, K. 2016 Modelling of turbulent energy
ux in canonical shock-
turbulence interaction. Int. J. Heat Fluid Fl. 61, 626{635.
128
Quadros, R., Sinha, K. & Larsson, J. 2016a Kovasznay mode decomposition of
velocity-temperature correlation in canonical shock-turbulence interaction. Flow Turbul.
Combust. 97 (3), 787{810.
Quadros, R., Sinha, K. & Larsson, J. 2016b Turbulent energy
ux generated by
shock/homogeneous-turbulence interaction. J. Fluid Mech. 796, 113{157.
Quinlan, J., McDaniel, J. C., Drozda, T. G., Lacaze, G. & Oefelein, J. C. 2014
A priori analysis of
amelet-based modeling for a dual-mode scramjet combustor. In 50th
AIAA/ASME/SAE/ASEE Joint Propulsion Conference, p. 3743.
Rawat, P. S. & Zhong, X. 2010 On high-order shock-tting and front-tracking schemes
for numerical simulation of shock{disturbance interactions. J. Comp. Phys. 229 (19),
6744{6780.
Reese, D. T., Ames, A. M., Noble, C. D., Oakley, J. G., Rothamer, D. A. &
Bonazza, R. 2018 Simultaneous direct measurements of concentration and velocity in
the Richtmyer{Meshkov instability. J. Fluid Mech. 849, 541{575.
Ribner, H. S. 1953 Convection of a pattern of vorticity through a shock wave. Tech. Rep..
NACA Tech. Rep. TN-2864.
Ribner, H. S. 1954 Shock-turbulence interaction and the generation of noise. Tech. Rep..
NACA Tech. Rep. 1233.
Ristorcelli, J. R. & Blaisdell, G. A. 1997 Consistent initial conditions for the DNS
of compressible turbulence. Phys. Fluids 9 (1), 4{6.
Rubins, P. M. & Bauer, R. C. 1994 Review of shock-induced supersonic combustion
research and hypersonic applications. J. Propuls. Power 10 (5), 593{601.
Ryu, J. & Livescu, D. 2014 Turbulence structure behind the shock in canonical shock{
vortical turbulence interaction. J. Fluid Mech. 756.
129
Saghafian, A., Shunn, L., Philips, D. A. & Ham, F. 2015 Large eddy simulations of
the HIFiRE scramjet using a compressible
amelet/progress variable approach. P. Com-
bust. Inst. 35 (2), 2163{2172.
Schwarzkopf, J. D., Livescu, D., Baltzer, J. R., Gore, R. A. & Ristorcelli,
J. R. 2016 A two-length scale turbulence model for single-phase multi-
uid mixing. Flow
Turbul. Combust. 96 (1), 1{43.
Seleznev, R.K., Surzhikov, S.T. & Shang, J.S. 2019 A review of the scramjet exper-
imental data base. Prog. Aerosp. Sci. 106, 43{70.
Sesterhenn, J., Dohogne, J. & Friedrich, R. 2005 Direct numerical simulation of
the interaction of isotropic turbulence with a shock wave using shock-tting. C. R. M ec.
333 (1), 87{94.
Sharpe, G. J. 2002 Shock-induced ignition for a two-step chain-branching kinetics model.
Phys. Fluids 14 (12), 4372{4388.
Sharpe, G. J. & Short, M. 2003 Detonation ignition from a temperature gradient for a
two-step chain-branching kinetics model. J. Fluid Mech. 476, 267{292.
Shen, H. S., Steinberg, J., Vanderover, J. & Oehlschlaeger, M. A. 2009 A
shock tube study of the ignition of n-heptane, n-decane, n-dodecane, and n-tetradecane
at elevated pressures. Energy & Fuels 23 (5), 2482{2489.
Shu, C. W. 1998 Essentially non-oscillatory and weighted essentially non-oscillatory
schemes for hyperbolic conservation laws. In Advanced numerical approximation of non-
linear hyperbolic equations, pp. 325{432. Springer.
Sinha, K., Mahesh, K. & Candler, G. V. 2003 Modeling shock unsteadiness in
shock/turbulence interaction. Phys. Fluids 15 (8), 2290{2297.
130
Smart, M. K., Hass, N. E. & Paull, A. 2006 Flight data analysis of the HyShot 2
scramjet
ight experiment. AIAA J. 44 (10), 2366{2375.
Smith, G. P., Golden, D. M., Frenklach, M., Moriarty, N. W., Eite-
neer, B., Goldenberg, M., Bowman, C. T., Hanson, R. K., Song, S.,
Gardiner, W. C., Lissianski, V. V. & Qin, Z. 1999 GRI 3.0 Mechanism.
http://combustion.berkeley.edu/gri-mech/version30/text30.html.
Stankovi c, I., Triantafyllidis, A., Mastorakos, E., Lacor, C. & Merci, B.
2011 Simulation of hydrogen auto-ignition in a turbulent co-
ow of heated air with LES
and CMC approach. Flow Turbul. and Combust. 86 (3-4), 689{710.
Strikwerda, J. C. 2004 Finite dierence schemes and partial dierential equations. SIAM.
Suman, S. & Girimaji, S. S. 2010 Velocity gradient invariants and local
ow-eld topology
in compressible turbulence. J. Turb. (11), N2.
Sun, M., Cui, X., Wang, H. & Bychkov, V. 2015 Flame
ashback in a supersonic
combustor fueled by ethylene with cavity
ameholder. J. Propul. Power 31 (3), 976{981.
Sunami, T., Itoh, K., Satoh, K. & Komuro, T. 2006 Mach 8 ground tests of the
hypermixer scramjet for HyShot-IV
ight experiment. In 14th AIAA/AHI Space Planes
and Hypersonic Systems and Technologies Conference, p. 8062.
Sunami, T. & Kodera, M. 2012 Numerical investigation of a detonation wave system
in a scramjet combustor. In 18th AIAA/3AF International Space Planes and Hypersonic
Systems and Technologies Conference, p. 5861.
Sv ard, M. & Nordstr om, J. 2014 Review of summation-by-parts schemes for initial{
boundary-value problems. J. Comput. Phys. 268, 17{38.
Thompson, M. 2003 At the Edge of Space: The X-15 Flight Program. Smithsonian.
131
Tian, Y., Jaberi, F. A., Li, Z. & Livescu, D. 2017 Numerical study of variable density
turbulence interaction with a normal shock wave. J. Fluid Mech. 829, 551{588.
Tirtey, Sandy, Boyce, Russell, Brown, Melrose, Capra, Bianca, Creagh,
Michael, Dedman, Amy, Dimitrijevic, Igor, Pudsey, Adrian, Sharp,
Bradley, Van Standen, Paul & others 2014 Scramspace: Radical-farming scramjet
for access to space. In STO-EN-AVT-234 Hypersonic Flight Testing. Von Karman Institute
for Fluid Dynamics (VKI).
Turns, S.R. 2000 An Introduction to Combustion: Concepts and Applications. McGraw-
hill.
Urzay, J. 2018 Supersonic combustion in air-breathing propulsion systems for hypersonic
ight. Annu. Rev. Fluid Mech. 50, 593{627.
Vaghefi, N. S. & Madnia, C. K. 2015 Local
ow topology and velocity gradient invariants
in compressible turbulent mixing layer. J. Fluid Mech. 774, 67{94.
Valiev, D., Bychkov, V., Akkerman, V., Law, C. K. & Eriksson, L. 2010 Flame
acceleration in channels with obstacles in the de
agration-to-detonation transition. Com-
bust. Flame 157 (5), 1012{1021.
Valiev, D. M., Bychkov, V., Akkerman, V. & Eriksson, L. 2009 Dierent stages
of
ame acceleration from slow burning to chapman-jouguet de
agration. Phys. Rev. E
80 (3), 036317.
V etel, J., Garon, A. & Pelletier, D. 2011 Denoising methods for time-resolved PIV
measurements. Exp. Fluids 51 (4), 893{916.
Wang, J., Shi, Y., Wang, L., Xiao, Z., He, X. & Chen, S. 2012 Eect of compress-
ibility on the small-scale structures in isotropic turbulence. J. Fluid Mech. 713, 588{631.
132
Wang, Z., Cai, Z., Sun, M., Wang, H. & Zhang, Y. 2016 Large eddy simulation of the
ame stabilization process in a scramjet combustor with rearwall-expansion cavity. Int. J.
Hydrog. Energy 41 (42), 19278{19288.
Westbrook, C. K. 1982 Chemical kinetics of hydrocarbon oxidation in gaseous detona-
tions. Combust. Flame 46, 191{210.
Westerweel, J., Elsinga, G. E. & Adrian, R. J. 2013 Particle image velocimetry for
complex and turbulent
ows. Annu. Rev. Fluid Mech. 45, 409{436.
Wilke, C. R. 1950 A viscosity equation for gas mixtures. J. Chem. Phys. 18 (4), 517{519.
Wolanski, P. 1973 Investigation into the mechanism of the diusion ignition of a com-
bustible gas
owing into an oxidizing atmosphere. In Fourteenth Symposium (Interna-
tional) on Combustion, 1973 .
Wouchuk, J. G., Huete, C. & Velikovich, A. L. 2009 Analytical linear theory for the
interaction of a planar shock wave with an isotropic turbulent vorticity eld. Phys. Rev.
E 79 (6), 066315.
Yamashita, H., Kasahara, J., Sugiyama, Y. & Matsuo, A. 2012 Visualization study
of ignition modes behind bifurcated-re
ected shock waves. Combust. Flame 159 (9), 2954{
2966.
Yamashita, H., Shimada, M. & Takeno, T. 1996 A numerical study on
ame stability at
the transition point of jet diusion
ames. In Symposium (International) on Combustion,
, vol. 26, pp. 27{34. Elsevier.
Yang, J., Kubota, T. & Zukoski, E. E. 1993 Applications of shock-induced mixing to
supersonic combustion. AIAA J. 31 (5), 854{862.
133
Yao, T., Yang, W. H. & Luo, K. H. 2018 Direct numerical simulation study of hydro-
gen/air auto-ignition in turbulent mixing layer at elevated pressures. Computers & Fluids
173, 59{72.
Yoo, C. S., Lu, T., Chen, J. H. & Law, C. K. 2011 Direct numerical simulations of
ignition of a lean n-heptane/air mixture with temperature inhomogeneities at constant
volume: Parametric study. Combust. Flame 158 (9), 1727{1741.
Yu, J. & Lu, X. 2019 Topological evolution near the turbulent/non-turbulent interface in
turbulent mixing layer. J. Turbul. 20 (5), 300{321.
Yu, R. & Bai, X. 2013 Direct numerical simulation of lean hydrogen/air auto-ignition in
a constant volume enclosure. Combust. Flame 160 (9), 1706{1716.
Zhang, S., Li, J., Qin, F., Huang, Z. & Xue, R. 2016 Numerical investigation of
combustion eld of hypervelocity scramjet engine. Acta Astronaut. 129, 357{366.
Zhao, G., Sun, M., Wu, J., Cui, X. & Wang, H. 2019 Investigation of
ame
ash-
back phenomenon in a supersonic cross
ow with ethylene injection upstream of cavity
ameholder. Aerosp. Sci. and Tech. 87, 190{206.
Zhou, Z., Fang L. Ouellette N. T. & Xu, H. 2020 Vorticity gradient stretching in the
direct enstrophy transfer process of two-dimensional turbulence. Phys. Rev. Fluids 5 (5),
054602.
Zikoski, Z. J. 2011 A parallel adaptive wavelet method for multidimensional simulations
of hypersonic propulsion. PhD thesis, University of Notre Dame.
134
Appendices
Appendix A: Transport equations
Appendix A.1: Passive scalar variance
The advection-diusion equation for a passive scalar can be written in conservation (or
ux) form as:
@
@t
+
@
@x
i
(u
i
) =
@M
i
@x
i
(5.1)
where M
i
=D@=@x
i
, and in advection form (using the continuity equation) as:
@
@t
+u
i
@
@x
i
=
@M
i
@x
i
(5.2)
Multiplying equation (5.1) by , it results, after rearranging:
@
2
@t
+
@
@x
i
(u
i
2
) =
@
@t
+u
i
@
@x
i
+
@M
i
@x
i
= 2
@M
i
@x
i
=
@M
i
@x
i
M
i
@
@x
i
(5.3)
where equation (5.2) has been used to obtain the second equality.
Reynolds-averaging equation (5.1) and using a =e a and ac
=e a
e
b +a
00
b
00
:
@
@t
(
e
) +
@
@x
i
(e u
i
e
) +
@
@x
i
(u
00
i
00
) =
@M
i
@x
i
(5.4)
135
Multiplying the previous equation by the Favre-averaged scalar,
e
, after rearranging, it
results:
@
2
@t
+
@
@x
i
(e u
i
e
2
)
e
@
@t
+e u
i
@
e
@x
i
!
+
e
@
@x
i
(u
00
i
00
) =
e
@M
i
@x
i
Reynolds-averaging equation (5.3):
@
@t
2
+
@
@x
i
u
i
2
=
@
@t
2
+
@
@t
00
00
+
@
@x
i
h
e u
i
e
2
+e u
i
00
00
+ 2
e
u
00
i
00
+u
00
i
00
00
i
= 2
@M
i
@x
i
Subtracting the last two equations, we obtain:
@
@t
00
00
+
@
@x
i
e u
i
00
00
+
@
@x
i
2
e
u
00
i
00
+u
00
i
00
00
+
e
@
e
@t
+e u
i
@
e
@x
i
!
e
@
@x
i
u
00
i
00
= 2
@M
i
@x
i
e
@M
i
@x
i
(5.5)
The parenthesis in the forth term on the left-hand side can be expressed as:
@
e
@t
+e u
i
@
e
@x
i
=
@
@t
e
+
@
@x
i
e u
i
e
e
@
@t
+
@e u
i
@x
i
| {z }
0
=
@
@x
i
u
00
i
00
+
@M
i
@x
i
where the Reynolds-averaged continuity equation has been used to zero the last term on the
left hand side and equation (5.4) has been used in the last equality.
Therefore, equation (5.5) results:
@
@t
00
00
+
@
@x
i
e u
i
00
00
+
@
@x
i
2
e
u
00
i
00
+u
00
i
00
00
2
e
@
@x
i
u
00
i
00
= 2
@M
i
@x
i
2
e
@M
i
@x
i
= 2
00
@M
i
@x
i
(5.6)
136
where the last equality is obtained noting thatab =e ab +a
00
b. Rewriting the following terms:
@
@x
i
e u
i
00
00
=e u
i
@
@x
i
00
00
+
00
00
@e u
i
@x
i
;
e
@
@x
i
u
00
i
00
=
@
@x
i
e
u
00
i
00
u
00
i
00
@
e
@x
i
;
00
@M
i
@x
i
=
@
@x
i
D
00
@
@x
i
"
D
@
e
@x
i
@
00
@x
i
+D
@
00
@x
i
@
00
@x
i
#
;
and rearranging, the transport equation for the Favre-averaged scalar variance, results:
@
00
00
@t
+e u
i
@
00
00
@x
i
=
00
00
@e u
i
@x
i
2u
00
i
00
@
e
@x
i
@
@x
i
u
00
i
00
00
+ 2
@
@x
i
D
00
@
@x
i
2D
@
00
@x
i
@
00
@x
i
2D
@
00
@x
i
@
e
@x
i
;
(5.7)
Appendix A.2: Passive scalar variance dissipation rate
From the transport equation for a passive scalar :
@
@t
+
@
@x
k
(u
k
) =
@
@x
k
D
@
@x
k
The LHS can be expanded as:
@
@t
+
@
@t
+u
k
@
@x
k
+
@
@x
k
(u
k
) =
@
@t
+u
k
@
@x
k
=
@
@x
k
D
@
@x
k
;
where continuity has been used in the rst equality. Dividing by the density, it results:
@
@t
+u
k
@
@x
k
=
1
@
@x
k
D
@
@x
k
(5.8)
Taking @(5:8)=@x
i
:
@
@t
@
@x
i
+
@u
k
@x
i
@
@x
k
+u
k
@
@x
k
@
@x
i
=
@
@x
i
1
@
@x
k
D
@
@x
k
(5.9)
137
Multiplying (5:9)@=@x
i
:
@
@x
i
@
@t
@
@x
i
+
@u
k
@x
i
@
@x
k
@
@x
i
+u
k
@
@x
i
@
@x
k
@
@x
i
=
@
@x
i
@
@x
i
1
@
@x
k
D
@
@x
k
(5.10)
Developing the rst and third terms of the LHS of the previous equation:
@
@t
@
@x
i
@
@x
i
@
@x
i
@
@t
@
@x
i
+
@u
k
@x
i
@
@x
k
@
@x
i
+ (5.11)
u
k
@
@x
k
@
@x
i
@
@x
i
u
k
@
@x
i
@
@x
k
@
@x
i
= (5.12)
@
@x
i
@
@x
i
1
@
@x
k
D
@
@x
k
(5.13)
Combining equations (5.10) and (5.13) and reorganizing terms:
@
@t
@
@x
i
@
@x
i
+u
k
@
@x
k
@
@x
i
@
@x
i
= (5.14)
2
@u
k
@x
i
@
@x
k
@
@x
i
+ 2
@
@x
i
@
@x
i
1
@
@x
k
D
@
@x
k
(5.15)
Reynolds-averaging the previous equation, the sought transport equation for the (scaled)
rate of dissipation of the scalar variance is obtained:
@
@t
@
@x
i
@
@x
i
+u
k
@
@x
k
@
@x
i
@
@x
i
= (5.16)
u
0
k
@
@x
k
@
@x
i
@
@x
i
2
@u
k
@x
i
@
@x
k
@
@x
i
+ 2
@
@x
i
@
@x
i
1
@
@x
k
D
@
@x
k
(5.17)
138
Appendix B: Regime diagram for reactive 2D STI simulations
Figure 5.1: Flame regime diagram presented for caseL showing the joint distributions in theDa
t
-Re
t
plane of domain with heat
release larger than 1e8 J/s/m
3
. The scatter plot is colored by cell number (top, in log scale), TKI (middle) and HRR (bottom,
in log scale) for snapshots t
1
,t
2
,t
peak
,t
3
from left to right. Regions with cell numbers below 50 are masked. The three contour
lines represent 80%, 40% and 20% of the peak cell number at the corresponding time instance.
139
Figure 5.2: Flame regime diagram presented for case M showing the joint distributions in the Da
t
-Re
t
plane of domain with
heat release larger than 1e8 J/s/m
3
. The scatter plot is colored by cell number (top, in log scale), TKI (middle) and HRR
(bottom, in log scale) for snapshots t
1
, t
2
, t
peak
, t
3
from left to right. Regions with cell numbers below 50 are masked. The
three contour lines represent 80%, 40% and 20% of the peak cell number at the corresponding time instance.
140
Figure 5.3: Flame regime diagram presented for case N showing the joint distributions in the Da
t
-Re
t
plane of domain with
heat release larger than 1e8 J/s/m
3
. The scatter plot is colored by cell number (top, in log scale), TKI (middle) and HRR
(bottom, in log scale) for snapshots t
1
, t
2
, t
peak
, t
3
from left to right. Regions with cell numbers below 50 are masked. The
three contour lines represent 80%, 40% and 20% of the peak cell number at the corresponding time instance.
141
Figure 5.4: Flame regime diagram presented for case O showing the joint distributions in the Da
t
-Re
t
plane of domain with
heat release larger than 1e8 J/s/m
3
. The scatter plot is colored by cell number (top, in log scale), TKI (middle) and HRR
(bottom, in log scale) for snapshots t
1
, t
2
, t
peak
, t
3
from left to right. Regions with cell numbers below 50 are masked. The
three contour lines represent 80%, 40% and 20% of the peak cell number at the corresponding time instance.
142
Abstract (if available)
Abstract
Shock induced mixing and ignition are studied in the canonical shock-turbulence interaction (STI) configuration via shock-capturing Direct Numerical Simulations (DNS). First, a parametric study of passive scalar mixing in the STI configuration is performed, varying the shock Mach number (M=1.28 to 5), turbulence Mach number (Mₜ=0.1 to 0.4), Taylor microscale Reynolds number (Reλ ≈ 40, 70), and Schmidt number (Sc=0.5, 1, 2). Streamwise budgets across the shock of the transport of scalar variance and scalar dissipation rate (SDR) are examined. The dominant quantity of scalar variance transport is SDR. There are two dominant terms in SDR transport equation. The first one represents the interaction between scalar gradient and velocity gradient. The change of this term across the shock increases with higher M, lower Mₜ and higher Taylor microscale Reynolds number. The second term represents molecular diffusion. Both this term and the SDR show a larger change across the shock with higher M, lower Mₜ and lower Taylor microscale Reynolds number. ❧ The significance of different flow topologies enhancing mixing across the shock is studied from changes in probability density functions of the invariants of the velocity gradient tensor and the distributions of scalar dissipation conditioned on these invariants. The stable focus stretching topology is dominant, while the unstable-node/saddle/saddle topology is the most dissipative throughout the domain, despite variations across the shock. Pre- and post-shock distributions of the alignment between strain-rate tensor eigenvectors and the scalar gradient, the vorticity, and the mean streamwise vector conditioned on flow topology are studied. A novel barycentric map representation is introduced for more direct visualization of alignments and conditioned scalar dissipation distributions. Interaction with the shock increases alignment of the scalar gradient with the most extensive eigenvector, decreasing it with the most compressive, which is still dominant. Across the shock, the most probable alignment between the passive scalar gradient and the eigenvectors of the strain-rate tensor converges towards the alignment with the largest scalar dissipation, enhancing scalar dissipation immediately downstream of the shock. ❧ Shock induced ignition in the STI setting is studied with 2D DNS using finite-rate detailed chemistry. To verify the implementation of the combustion capability, comparison with several benchmark cases is performed. The ability of 2D STI simulation to reproduce relevant physics present in the 3D STI configuration is assessed. The, 2D reactive homogeneous isotropic turbulence (HIT) simulations are performed to explore the influence of Taylor microscale Reynolds number and Mₜ on the ignition delay time, finding that larger Mₜ triggers early ignition. ❧ Four reactive STI simulations are conducted with initial Taylor microscale Reynolds number to be 1250, 600 and Mₜ = (0.1, 0.3) at M = 2 with H₂, O₂, and Ar mixture, and compared with laminar simulation at the same Mach number. Compressibility of the upstream turbulence leads to earlier ignition compared with laminar simulation, and the peak values of thermodynamic quantities at the flame front are found to be lower for the turbulent cases under consideration, owing to the partially premixed nature of the mixture. An analysis of the temporal evolution of the flame regime reveals that the reaction happens mostly in the thin reaction zone regime, which is characterizes by a broadened preheat zone. Lower Mₜ brings a slightly higher probability that the combustion happens in regime of corrugated flamelets. The time evolution of the Takeno Flame Index (TFI) indicates that, as the flame propagates upstream, the combustion becomes increasingly non-premixed because of the larger variance of the species mass fractions.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Tracking and evolution of compressible turbulent flow structures
PDF
Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations
PDF
Large eddy simulations of turbulent flows without use of the eddy viscosity concept
PDF
Accuracy and feasibility of combustion studies under engine relevant conditions
PDF
Supersonic flow study of a capsule geometry using large-eddy simulation
PDF
Re-assessing local structures of turbulent flames via vortex-flame interaction
PDF
Investigations of fuel and hydrodynamic effects in highly turbulent premixed jet flames
PDF
On the simulation of stratified turbulent flows
PDF
Numerical study of focusing effects generated by the coalescence of multiple shock waves
PDF
CFD design of jet-stirred chambers for turbulent flame and chemical kinetics experiments
PDF
Numerical simulations of linearly stratified flow around a sphere
PDF
Hydrodynamic stability of two fluid columns of different densities and viscosities subject to gravity
PDF
Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles
PDF
The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems
PDF
Model based design of porous and patterned surfaces for passive turbulence control
Asset Metadata
Creator
Gao, Xiangyu
(author)
Core Title
Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Mechanical Engineering
Publication Date
07/14/2020
Defense Date
06/19/2020
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
combustion,direct numerical simulation,mixing,OAI-PMH Harvest,shock,turbulence
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Bermejo-Moreno, Ivan (
committee chair
), Domaradzki, Julian (
committee member
), Egolfopoulos, Fokion (
committee member
), Nakano, Aiichiro (
committee member
), Pantano-Rubino, Carlos (
committee member
)
Creator Email
xiangyug@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c89-330182
Unique identifier
UC11664340
Identifier
etd-GaoXiangyu-8674.pdf (filename),usctheses-c89-330182 (legacy record id)
Legacy Identifier
etd-GaoXiangyu-8674.pdf
Dmrecord
330182
Document Type
Dissertation
Rights
Gao, Xiangyu
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
combustion
direct numerical simulation
mixing
turbulence