NUMERICAL SIMULATIONS AND ANALYSES OF
SHOCK WAVE-BOUNDARY LAYER INTERACTIONS
By
Avinash Jammalamadaka

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Mechanical Engineering – Doctor of Philosophy
2013

ABSTRACT
NUMERICAL SIMULATIONS AND ANALYSES OF SHOCK
WAVE-BOUNDARY LAYER INTERACTIONS
By
Avinash Jammalamadaka
Shock-boundary layer interaction (SBLI) is becoming one of the benchmark problems
in the high-speed flow modeling and simulation community. The interaction of shock wave
with the boundary layer is a very complex phenomenon that requires high-fidelity numerical
methods like direct numerical simulation (DNS) and large-eddy simulation (LES) to capture
the flow physics. In this study, SBLI is examined for various flow conditions using DNS and
LES.
In the first part of this study, DNSs are conducted for a Mach 2.75 turbulent boundary layer interacting with an impinging shock at three different shock incidence angles.
Instantaneous flow visualizations show the effect of shock on turbulence structure in the
shock-boundary layer interaction zone and also in the flow downstream of the interaction
region. The separation bubbles exhibit highly unsteady and three-dimensional behavior and
are larger for stronger shocks but the maximum probability of flow separation is found to
be independent of the shock strength. A detailed examination of the terms in the turbulent
kinetic equation shows a strong coupling exists between the mean and turbulent fields in
the interaction region with energy being continuously exchanged from one field to another.
However, the compressibility-related terms in the transport equations for turbulent kinetic
energy and enstrophy are found to be small for the simulated flows.
In the second part of this study, data generated by DNS for a Mach 2.75 zero-pressure
gradient turbulent boundary layer interacting with shocks of different intensities are used for

a priori analysis of subgrid-scale (SGS) turbulence and various terms in the compressible
filtered Navier-Stokes equations. The behavior of SGS stresses and their components, namely,
Leonard, Cross and Reynolds components, is examined in various regions of the flow for
different shock intensities and filter widths. The backscatter in various regions of the flow is
found to be significant only instantaneously while the ensemble-averaged statistics indicate
no significant backscatter. The budgets for the SGS kinetic energy equation are examined for
a better understanding of shock-tubulence interactions at the subgrid level and also with the
aim of providing useful information for one-equation LES models. A term-by-term analysis
of the SGS terms in the filtered total energy equation indicate that while each term in this
equation is significant by itself, the net contribution by all of them is relatively small. This
observation is consistent with our a posteriori analysis.
In the third part of the study, assessments of several existing SGS LES models along with
a new model are made for SBLI through systematic a priori and a posteriori analyses. In
the problem chosen for this study, the incident shock is strong enough to generate a marginal
separation of the boundary layer near the interaction region, hence providing the SGS models
with a non-trivial challenge. The effect of SGS stress term on the resolved velocity field is
found to be significant and shock-dependent. The subgrid-scale models tested include the
mixed-time-scale model, the dynamic Smagorinsky model, the dynamic mixed model and a
new dynamic model, termed the compressible serial decomposition model. A priori analysis
indicate that the new dynamic model is more accurate than other SGS closures. A posteriori
tests also indicate better predictions of DNS results by the LES employing the compressible
serial decomposition model.

ACKNOWLEDGMENTS

Firstly, I would like to sincerely thank Dr. Farhad Jaberi, the driving force behind this study,
for giving me an opportunity to work with him. This work would not have been successful
if it weren’t for his patience and constant encouragement. It has been a truly wonderful
learning experience to work for Dr. Jaberi.
I am grateful to Dr. Ahmed Naguib, Dr. Giles Brereton and Dr. Guowei Wei, for
agreeing to be a part of my PhD committee. I thank them for their time and effort put into
reading this dissertation and for their helpful suggestions and comments. The discussions I
had with Dr. Naguib were very helpful.
This work was primarily supported by the Air Force Research Laboratory (AFRL) under
the agreement number FA8650-06-2-3625 and I wish to sincerely thank Dr. John Benek of
AFRL for his guidance.
I would like to sincerely thank Dr. Zhaorui Li for all his help and suggestions. The
numerous discussions I had with him were very fruitful. I would also like to gratefully
acknowledge the help received from Dr. Sergio Pirozzoli of University of Rome, La Sapienza,
Italy, for his help in implementing the hybrid code used in this study.
I express my gratitude to the ICER/HPCC staff for their support in providing the computational facilities for this study.
Finally, I am ever-indebted to my parents and my sister who have always supported me.

iv

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

Chapter 1 Parametric Investigations of Shock Wave/Boundary Layer Interactions Using Direct Numerical Simulation . . . . . . . . . . .
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Computational Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2 Numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.3 Inflow turbulence generation . . . . . . . . . . . . . . . . . . . . . . .
1.2.4 Computational setup for DNS of SBLI . . . . . . . . . . . . . . . . .
1.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Accuracy of DNS . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2 Flow Visualizations . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.3 Reynolds- and Favre-Averaging . . . . . . . . . . . . . . . . . . . . .
1.3.4 Flow Reversal Comparisons . . . . . . . . . . . . . . . . . . . . . . .
1.3.5 Mean and Fluctuating Wall-Pressure Comparisons . . . . . . . . . . .
1.3.6 Thermodynamic Variables . . . . . . . . . . . . . . . . . . . . . . . .
1.3.7 Turbulent Kinetic Energy Budget . . . . . . . . . . . . . . . . . . . .
1.3.8 Enstrophy budget . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
4
4
5
9
10
12
13
16
21
24
31
34
39
51
57

Chapter 2 A Priori Analysis of Subgrid-Scale Turbulence for Shock-Boundary
Layer Interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
2.2 DNS data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
2.3.1 SGS terms in the filtered momentum equation . . . . . . . . . . . . . 66
2.3.2 SGS terms in the filtered energy equation . . . . . . . . . . . . . . . 75
2.3.3 Resolved and residual kinetic energies . . . . . . . . . . . . . . . . . . 83
2.3.4 Backscatter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
2.3.5 Directional filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
2.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
Chapter 3 SGS Models and their Applications to LES of Shock-Boundary
Layer Interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

v

3.1
3.2
3.3
3.4
3.5

3.6
3.7

Introduction . . . . . . . . . . . . . . . . . .
Governing Equations . . . . . . . . . . . . .
Subgrid-Scale Models . . . . . . . . . . . . .
Numerical Method . . . . . . . . . . . . . .
Results . . . . . . . . . . . . . . . . . . . . .
3.5.1 Computational Setup and Parameters
3.5.2 A Priori Assessment . . . . . . . . .
3.5.3 A Posteriori Assessment . . . . . . .
3D effects of lateral walls on SBLI . . . . . .
Conclusions . . . . . . . . . . . . . . . . . .

BIBLIOGRAPHY

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

vi

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

109
112
115
119
120
120
125
139
148
153

. . . . . . . . . . . . . . . . . 155

LIST OF TABLES

Table 1.1

Simulation Parameters for DNS of SBLI. . . . . . . . . . . . . . . .

12

Table 1.2

Boundary layer parameters for the incoming boundary layer. . . . .

13

Table 1.3

Separation and Interaction lengths. . . . . . . . . . . . . . . . . . .

31

Table 1.4

Legend and colors for turbulent kinetic energy budget plots. . . . . .

42

Table 1.5

Legend and colors for enstrophy budget plots.

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

52

Table 2.1

Color and legend for filtered total energy plots. . . . . . . . . . . . .

79

Table 2.2

Color and legend for ksgs and normal SGS stress budget plots. . . .

91

Table 3.1

Simulation Parameters for DNS of Mach 2 SBLI . . . . . . . . . . . 120

Table 3.2

Simulation Parameters for LESof Mach 2 SBLI . . . . . . . . . . . . 121

Table 3.3

Separation region sizes . . . . . . . . . . . . . . . . . . . . . . . . . 144

vii

LIST OF FIGURES

Figure 1.1

Instantaneous density contours overlaid with modified Ducros shock
sensor values. The critical regions identified by the shock sensor are
shown in red. For interpretation of the references to color in this and
all other figures, the reader is referred to the electronic version of this
dissertation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

Figure 1.2

The wall pressure spectra at x = 23δref for the “auxiliary” DNS. . .

11

Figure 1.3

(a) Van Driest transformed mean streamwise velocity profile and,
(b) density-scaled Reynolds stresses of the incoming boundary layer.
Symbols represent DNS data of Pirozzoli and Bernardini [1] . . . . .

15

Instantaneous streamwise densitygradient contours for Case III. Contour levels are evenly spaced between -1 and 1. . . . . . . . . . . . .

17

Instantaneous numerical schlieren visualizations at an x−y symmetry
plane. (a) Case I, (b) Case II, (c) Case III. . . . . . . . . . . . . . .

18

(a) Instantaneous velocity fluctuations and, (b) instantantaneous temperature fluctuations along the mean sonic plane for Case II. . . . .

20

Figure 1.4

Figure 1.5

Figure 1.6

Figure 1.7

Figure 1.8

Figure 1.9

Figure 1.10

Figure 1.11

Differences between Reynolds- and Favre-averaged quantities at X ∗ =
−1.1. (a) velocity, (b) temperature. . . . . . . . . . . . . . . . . . .

22

Differences between Reynolds- and Favre-averaged quantities at X ∗ =
−0.5. (a) velocity, (b) temperature. . . . . . . . . . . . . . . . . . .

23

(a) Comparison of Cf computed using Reynolds- and Favre-averaged
quantities. (b) Percentage difference in Cf computed by Reynoldsand Favre-averaged quantities. The results shown are for Case II. . .
Comparisons for ρ ui uj . (a) for different shock intensities. X ∗ =
−1.1, (b) X ∗ = 0.25. . . . . . . . . . . . . . . . . . . . . . . . . . .
Instantaneous streamwise skin-friction coefficient contours. The contour levels are evenly spaced between −6 × 10−3 and 0. (a) Case I,
(b) Case II, (c) Case III. . . . . . . . . . . . . . . . . . . . . . . . .
viii

25

26

28

Figure 1.12

Figure 1.13

Figure 1.14

Figure 1.15

Figure 1.16

Figure 1.17

Figure 1.18

Figure 1.19

Instantaneous spanwise skin-friction coefficient contours. The contour levels are evenly spaced between −5 × 10−3 and 5 × 10−3 . (a)
Case I, (b) Case II, (c) Case III. . . . . . . . . . . . . . . . . . . . .

29

Comparisons for the probability of flow reversal for different shock
intensities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

Comparisons for mean skin-friction coefficient for different shock intensities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

Mean wall pressure distribution for different shock intensities plotted
with different scalings. . . . . . . . . . . . . . . . . . . . . . . . . . .

33

(a) Instantaneous wall pressure fluctuations for Case II. (b) Streamwise evolution of r.m.s. of wall pressure fluctuations for different
shock intensities. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

R.M.S. of thermodynamic fluctuations at different streamwise locations for Case II. (a) pressure, (b) temperature. . . . . . . . . . . . .

37

Contours of turbulent fluctuations of thermodynamic variables for
Case II. (a) r.m.s. of pressure fluctuations, (b) r.m.s. of density
fluctuations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

Turbulent and fluctuating Mach numbers at different streamwise locations for Case II. (a) Mt , (b) M . . . . . . . . . . . . . . . . . .

40

Figure 1.20

Turbulent kinetic energy budget for Case II at X ∗ = −1.1. . . . . .

Figure 1.21

Turbulent kinetic energy budgets at X ∗ = −0.5. (a) Case II, (b)
Case III. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

Turbulent kinetic energy budgets at X ∗ = 0.25. (a) Case II, (b) Case
III. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

Turbulent kinetic energy budgets at X ∗ = 1.5. (a) Case II, (b) Case
III. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

(a) Streamwise evolution of turbulent kinetic energy along the mean
sonic line. (b) Streamwise evolution of turbulent kinetic energy budget along the sonic line. . . . . . . . . . . . . . . . . . . . . . . . . .

48

Figure 1.22

Figure 1.23

Figure 1.24

ix

42

Figure 1.25

Evolution of dissipation term and its components along the mean
sonic line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

Figure 1.26

Enstrophy budget for Case II at X ∗ = −1.1. . . . . . . . . . . . . .

53

Figure 1.27

Enstrophy budgets at X ∗ = −0.5. (a) Case II, (b) Case III. . . . . .

54

Figure 1.28

Enstrophy budgets at X ∗ = 0.25. (a) Case II, (b) Case III. . . . . .

55

Figure 1.29

Enstrophy budgets at X ∗ = 1.5. (a) Case II, (b) Case III. . . . . . .

56

Figure 2.1

SGS stresses at different streamwise locations as obtained from the
filtered DNS data for θ = 7.75o and ∆ = 8∆DN S . (a) X ∗ = −1.1,
(b) X ∗ = −0.5. Black solid line for τ11 , black dashed line for τ22 ,
black dash-dot line for τ33 , black dots for τ12 , red solid line for ksgs.

69

SGS stresses at different streamwise locations as obtained from the
filtered DNS data for θ = 7.75o and ∆ = 8∆DN S . (a) X ∗ = 0.25,
(b) X ∗ = 1.5. Black solid line for τ11 , black dashed line for τ22 , black
dash-dot line for τ33 , black dots for τ12 , red solid line for ksgs. . . .

70

PDFs of τ12 SGS stresse as obtained from the instantaneous filtered
DNS data for ∆ = 8∆DN S . . . . . . . . . . . . . . . . . . . . . . . .

71

SGS stress decomposition at different streamwise locations as obtained from the filtered DNS data for θ = 7.75o and ∆ = 4∆DN S .
(a) X ∗ = −1.1, (b) X ∗ = −0.5. Black solid line for τ11 , black dashed
line for L11 , black dash-dot line for C11 , black dots for R11 . . . . . .

73

SGS stress decomposition at different streamwise locations as obtained from the filtered DNS data for θ = 7.75o and ∆ = 8∆DN S .
(a) X ∗ = 0.25, (b) X ∗ = 1.5. Legend same as in Fig. 2.4. . . . . . .

74

Figure 2.2

Figure 2.3

Figure 2.4

Figure 2.5

Figure 2.6

Figure 2.7

SGS stress decomposition at different streamwise locations as obtained from the filtered DNS data for θ = 7.75o and ∆ = 16∆DN S .
(a) X ∗ = −0.5, (b) X ∗ = 0.25. Legend same as in Fig. 2.4. . . . . .
Various SGS terms in the filtered total energy equation at different
streamwise locations as obtained from the filtered DNS data for θ =
7.75o and ∆ = 8∆DN S . (a) X ∗ = −1.1, (b) X ∗ = −0.5. See
Table. 2.1 for legend. . . . . . . . . . . . . . . . . . . . . . . . . . .

x

76

80

Figure 2.8

Figure 2.9

Figure 2.10

Skin-friction coefficient predictions as obtained from a posteriori LES
calculations for θ = 7.75o and ∆ = 8∆DN S . Solid black line for
filtered DNS, red dashed line for LES1, blue dash-dot line for LES2.

82

Turbulent Reynolds stress at different streamwise locations as obtained from a posteriori LES calculations for θ = 7.75o and ∆ =
8∆DN S . (a) X ∗ = −1.1, (b) X ∗ = −0.5. Solid black line for filtered
DNS, red dashed line for LES1, blue dash-dot line for LES2. . . . .

84

Turbulent heat flux at different streamwise locations as obtained from
a posteriori LES calculations for θ = 7.75o and ∆ = 8∆DN S . (a)
X ∗ = −1.1, (b) X ∗ = −0.5. Legend same as in Fig. 2.9. . . . . . . .

85

Figure 2.11

The ratio of SGS kinetic energy to the total kinetic energy (ksgs /(ksgs +
Er )) at different streamwise locations as obtained from the filtered
DNS data for θ = 7.75o and ∆ = 4∆DN S . (a) X ∗ = −1.1, (b)
X ∗ = −0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

Figure 2.12

ksgs /(ksgs + Er ) at different streamwise locations as obtained from
the filtered DNS data for θ = 7.75o and ∆ = 8∆DN S . (a) X ∗ = −1.1,
(b) X ∗ = −0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.13

Figure 2.14

Figure 2.15

Figure 2.16

ksgs /(ksgs +Er ) at different streamwise locations as obtained from the
filtered DNS data for θ = 7.75o and ∆ = 16∆DN S . (a) X ∗ = −1.1,
(b) X ∗ = −0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ksgs /(ksgs + Er ) as obtained from the filtered DNS data for different
shock intensities with ∆ = 8∆DN S . (a) θ = 6o , (b) θ = 7.75o , (c)
θ = 9.25o . The contours levels are uniformly spaced between 0 and
0.5. The figures are shown not to scale for clarity. . . . . . . . . . .
SGS kinetic energy (ksgs ) budget at different streamwise locations as
obtained from the filtered DNS data for θ = 7.75o and ∆ = 8∆DN S .
(a) X ∗ = −1.1, (b) X ∗ = −0.5. See Table. 2.2 for legend. . . . . . .
SGS kinetic energy (ksgs ) budget at different streamwise locations as
obtained from the filtered DNS data for θ = 7.75o and ∆ = 8∆DN S .
(a) X ∗ = 0.25, (b) X ∗ = 1.5. See Table. 2.2 for legend. . . . . . . .

xi

88

89

90

93

94

Figure 2.17

Figure 2.18

Instantaneous contours of SGS production as obtained from the filtered DNS data for θ = 7.75o and ∆ = 8∆DN S . Only the backscatter
regions are shown. 50 uniformly spaced contour levels between 0.0
and 0.005. (a) X ∗ = −1.1, (b) X ∗ = −0.5, (c) X ∗ = 0.25, (d)
X ∗ = 1.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97

Probability of backscatter in the x − y symmetry plane as obtained
from filtered DNS data for θ = 7.75o and ∆ = 8∆DN S . . . . . . . .

98

Figure 2.19

Directionally filtered SGS stresses as obtained from the filtered DNS
data for θ = 7.75o and ∆ = 8∆DN S at X ∗ = −1.1. (a) τ11 , (b) τ22 .
solid black line for x-only filter, dashed black line for y-only filter, red
dash-dot line for z-only filter . . . . . . . . . . . . . . . . . . . . . . 101

Figure 2.20

Directionally filtered SGS stresses as obtained from the filtered DNS
data for θ = 7.75o and ∆ = 8∆DN S at X ∗ = −1.1. (a) τ33 , (b) τ12 .
solid black line for x-only filter, dashed black line for y-only filter, red
dash-dot line for z-only filter . . . . . . . . . . . . . . . . . . . . . . 102

Figure 2.21

Directionally filtered SGS stresses as obtained from the filtered DNS
data for θ = 7.75o and ∆ = 8∆DN S at X ∗ = −0.5. (a) τ11 , (b) τ22 .
solid black line for x-only filter, dashed black line for y-only filter, red
dash-dot line for z-only filter . . . . . . . . . . . . . . . . . . . . . . 103

Figure 2.22

Directionally filtered SGS stresses as obtained from the filtered DNS
data for θ = 7.75o and ∆ = 8∆DN S at X ∗ = −0.5. (a) τ33 , (b) τ12 .
solid black line for x-only filter, dashed black line for y-only filter, red
dash-dot line for z-only filter . . . . . . . . . . . . . . . . . . . . . . 104

Figure 2.23

Directionally filtered SGS stresses as obtained from the filtered DNS
data for θ = 7.75o and ∆ = 8∆DN S at X ∗ = 0.25. (a) τ11 , (b) τ22 .
solid black line for x-only filter, dashed black line for y-only filter, red
dash-dot line for z-only filter . . . . . . . . . . . . . . . . . . . . . . 105

Figure 2.24

Directionally filtered SGS stresses as obtained from the filtered DNS
data for θ = 7.75o and ∆ = 8∆DN S at X ∗ = 0.25. (a) τ33 , (b) τ12 .
solid black line for x-only filter, dashed black line for y-only filter, red
dash-dot line for z-only filter . . . . . . . . . . . . . . . . . . . . . . 106

Figure 3.1

An x − y slice of the sample DNS grid. . . . . . . . . . . . . . . . . 122

Figure 3.2

A schematic of the computational procedure for DNS and LES of SBLI.123

xii

Figure 3.3

Flow structure in a typical shock/boundary layer interaction configuration. Vortex structures identified using λ2 criterion [2] and colored
based on the local values of instantaneous velocity magnitude. . . . 124

Figure 3.4

Van Driest transformed mean streamwise velocity profile of the incoming boundary layer in inner scaling calculated from DNS data at
+
x = 6. Black dashed line represents UV D = y + and Eq. (1.14) is
represented by black dash-dot line. . . . . . . . . . . . . . . . . . . . 125

Figure 3.5

Density-scaled profiles of Reynolds stresses obtained by DNS at x = 6.126

Figure 3.6

Leonard, Cross and Reynolds parts of the spanwise-averaged τ11 SGS
stress component at x = 6 for ∆ = 4∆DN S . total stress - black solid
line; L11 - dashed black line; C11 - dash-dot red line; R11 - blue
dotted line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128

Figure 3.7

Leonard, Cross and Reynolds parts of the spanwise-averaged τ11 SGS
stress component at x = 6 for ∆ = 8∆DN S . total stress - black solid
line; L11 - dashed black line; C11 - dash-dot red line; R11 - blue
dotted line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129

Figure 3.8

Spanwise-averaged values of the trace-free τ11 component of SGS
stress tensor at (a) x = 6, (b) x = 10. Legend : open triangles
= True stresses, dotted black line = Modeled-MTSM, dash-dot green
line = Modeled-DSM, dashed blue line = Modeled-DMM, solid red
line = Modeled-CSDM. . . . . . . . . . . . . . . . . . . . . . . . . . 131

Figure 3.9

Spanwise-averaged isocontours of trace-free τ11 SGS stress component. (a) True stress, (b) Modeled-DSM, (c) Modeled-DMM, (d)
Modeled-CSDM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

Figure 3.10

Spanwise-averaged values of the τ12 SGS stress component at (a)
x = 10, (b) x = 15. See Fig. 3.8 for legend. . . . . . . . . . . . . . . 133

Figure 3.11

Spanwise-averaged values of q1 SGS heat flux component at (a) x =
10, (b) x = 15. See Fig. 3.8 for legend. . . . . . . . . . . . . . . . . . 135

Figure 3.12

Correlation coefficients between the true and the modeled values of
τ11 SGS stress for different values of density gradient parameter ψ.
Legend : open circles = DSM, open triangles = DMM , open squares
= CSDM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

xiii

Figure 3.13

Correlation coefficients between the true and the modeled values of q1
heat flux component for different values of density gradient parameter
ψ. Legend : open circles = DSM, open triangles = DMM , open
squares = CSDM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

Figure 3.14

Variations of backscatter with the density gradient parameter ψ. . . 139

Figure 3.15

Correlation coefficients for the SGS dissipation, εsgs , as a function of
the density gradient parameter ψ as predicted by CSDM. . . . . . . 140

Figure 3.16

Turbulent kinetic energy profiles at x = 6. Legend : open circles
= DNS, open triangles = Filtered DNS, crosses = LES-No-Model,
dotted black line = LES-MTSM, dash-dot green line = LES-DSM,
dashed blue line = LES-DMM, solid red line = LES-CSDM. . . . . . 141

Figure 3.17

One-dimensional energy spectra in spanwise direction at x = 6 and
y ≈ 0.07. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142

Figure 3.18

Evolution of mean skin-friction coefficient, Cf . See Fig. 3.16 for legend.143

Figure 3.19

Evolution of mean wall-pressure, Pw . See Fig. 3.16 for legend.

Figure 3.20

Van-Driest transformed mean streamwise velocity profiles in inner
scaling. (a) x = 10, (b) x = 15. Black dashed line represents law of
the wall. See Fig. 3.16 for legend. . . . . . . . . . . . . . . . . . . . 147

Figure 3.21

Comparisons of R11 Reynolds stress component obtained by DNS
and LES. (a) x = 6, (b) x = 10. See Fig. 3.16 for legend. . . . . . . 149

Figure 3.22

Comparisons of R22 Reynolds stress component obtained by DNS
and LES. (a) x = 10, (b) x = 15. See Fig. 3.16 for legend. . . . . . . 150

Figure 3.23

Instantaneous slices of temperature contours at different streamwise
locations for 3DSBLI. . . . . . . . . . . . . . . . . . . . . . . . . . . 152

xiv

. . . 145

Chapter 1
Parametric Investigations of Shock
Wave/Boundary Layer Interactions
Using Direct Numerical Simulation

1.1

Introduction

Shock-turbulent boundary layer interaction (SBLI) is becoming one of the benchmark problems for high-speed flow modeling and simulations. This problem has mainly been studied
in two flow configurations, a compression corner and a flat plate turbulent boundary layer
with an incident shock. In the former, the shock is generated as a result of the geometry of
the surface on which the boundary layer is developing. In the latter, the shock is generated
from an external source and impinges on the boundary layer. In either case, turbulence
is amplified across the shock and the boundary layer experiences a strong adverse pressure
gradient which may lead to its separation. An important feature of the separated boundary
layer is the low-frequency, large-scale unsteadiness of the shock that has been observed in
experiments and numerical studies, but is yet to be fully understood. The characteristic
frequency of the shock oscillation is generally an order or two magnitudes smaller than the
characteristic frequency of the incoming boundary layer.

1

The reflected SBLI problem has been receiving significant attention in recent years even
though one of the earliest experimental studies was performed more than four decades ago [3].
In this experimental study, the shock deflection angle (θ) varied between 2o and 10.5o while
the Mach number was kept constant at 2.5. Measurements were made using pitot probes
and wall-pressure taps and flow visualizations were made using surface oil and optical techniques like schlieren and shadowgraph techniques. The experimental results indicated that
the boundary layer remains attached for weaker incident shocks but separates for stronger
incident shocks. The interaction region for the separated boundary layer was shown to exhibit a complex pattern with the primary reflected shock followed by an expansion wave and
then a series of compression waves near the reattachment point.
Dupont et al.[4] studied the interaction of an impinging oblique shock with a flat-plate
turbulent boundary layer at M = 2.3 using PIV technique. The primary focus was on the
unsteadiness of the reflected shock. Measurements were made for different incidence shock
angles θ which varied between 7o and 9.5o . By analyzing the wall-pressure signals in the
interaction region, the authors found some low-frequency oscillations in the second part of
the recirculation bubble. The characteristic frequency of these oscillations was found to be
close to that of the reflected shock suggesting that the reflected shock unsteadiness could be
related to the unsteadiness in the recirculation bubble. The effect of interaction strength on
the unsteadiness of separation bubble was experimentally investigated by Souverein et al. [5].
From their experimental data, the authors inferred the following: when no separation occurs,
unsteadiness of the shock motion is correlated with the incoming boundary layer; when mean
flow separation occurs, the shock motion exhibits a strong correlation with the unsteadiness
of the separation bubble; when incipient separation occurs, the shock unsteadiness appeared
to be related to a superposition of both upstream and downstream flow conditions.
2

Touber and Sandham[6] studied the effect of shock strength on SBLI using large-eddy
simulation (LES). They replicated three different experimental configurations and found
a good agreement between LES and experiments. They showed that for low-to-moderate
shock intensities, the length of SBLI region, when scaled by the upstream boundary layer
thickness, is linearly proportional to the shock intensity scaled by the upstream wall shear
stress. Using LES, Morgan et al. [7] carried out extensive parametric studies of SBLI for
various wedge angles, Reynolds numbers, streamwise and spanwise domain sizes. They found
that the Reynolds number did not have any significant effect on the separation bubble size.
Increasing the wedge angle led to an increased size of the interaction region but did not
necessarily cause a higher probability of reversed flow.
A state-of-the-art DNS calculations were performed by Pirozzoli and Bernardini [8] for a
Mach 2.25 boundary layer interacting with an oblique shock. The numerical strategy relied
on a hybrid scheme that used a high-order non-dissipative central scheme [9] in the shockfree regions and a fifth-order WENO scheme in the shock regions. A very long streamwise
domain size was used to avoid possible adverse effects of the recycling-rescaling procedure
and also to capture the recovery of boundary layer to an equilibrium state downstream of the
interaction region. While the differences in Reynolds numbers prevented a definite one-toone comparison with experiments, the DNS was able to reproduce the overall flow structure
of SBLI observed in experiments.
The main objectives of the present study are to examine SBLI under various flow conditions and to improve the basic understanding of SBLI by a detailed analysis of DNS data.
The Mach number of the incoming boundary layer is kept fixed at 2.75 and the strength of
the interacting shock is varied by varying the flow deflection angle, θ. Three different flow
deflection angles, θ = 6o , 7.75o and 9.25o, are considered. The selected flow parameters lie
3

on the higher end of DNS studies performed to date on SBLIs. The DNS calculations are
performed using a hybrid numerical methodology that uses a minimally dissipative central
scheme [9] in the shock-free regions of the flow and a very robust monotonicity preserving
scheme [10, 11] in the shock regions. More details on the numerical method are described
in Sec. 1.2. The effect of interacting shock strength on various flow quantities, mean and
turbulent statistics is examined. The contributing terms to the turbulent kinetic energy and
enstrophy variations are examined with the aim of explaining the effect of shock on turbulence, for providing useful information for high-order turbulence models and also for testing
the accuracy of DNS calculations.

1.2
1.2.1

Computational Methodology
Governing equations

For the current DNS computations, we solve the compressible, three-dimensional NavierStokes equations in conservative form. The dimensionless equations read as:

∂
∂ρ
+
ρuj = 0
∂t ∂xj

(1.1)

∂σij
∂
∂p
∂
(ρui ) +
ρui uj = −
+
∂t
∂xj
∂xi
∂xj

(1.2)

∂Qj
∂E
∂
∂
+
(E + p) uj =
σij ui −
,
∂t
∂xj
∂xj
∂xj

(1.3)

where the viscous stress and heat flux terms are obtained from the following Newtonian
models:
σij =

2µ 1
Re 2

∂ui ∂uj
+
∂xj
∂xi
4

−

1 ∂uk
δ
3 ∂xk ij

(1.4)

Qj =

∂T
−µ
2 ∂x
(γ − 1)ReP rM∞ j

(1.5)

Here, the temperature dependence of µ is based on Sutherland’s law, µ(T ) = T 1.5 (1 +
C)/(T + C) with C = 0.82. In Eq. (3.3), E is the total energy defined as

E=

1
p
+ ρui ui
γ−1 2

(1.6)

Finally, the thermodynamic variables, ρ, p and T are linked by the ideal gas law,

p=

1.2.2

ρT
.
2
γM∞

(1.7)

Numerical scheme

DNS and LES of turbulent flows require low dissipation numerical schemes to avoid unphysical damping of “small” scales due to numerical dissipation. In the presence of shocks,
it becomes necessary to use shock-fitting schemes to prevent spurious oscillations near the
shock regions introduced due to the Gibbs phenomenon. While standard shock-capturing
schemes like WENO or monotonicity preserving scheme (MP) [10, 11] successfully capture
the shock without any spurious oscillations, they can be dissipative in the regions away from
the shock and may reduce the accuracy of the solution. It should be noted that DNS computations of SBLI have been performed using WENO [12, 13, 14, 15] and MP [16] schemes
with a good degree of accuracy. However, these schemes are not suitable for flows at higher
Reynolds numbers or flows with strong shock/turbulence interactions. An ideal solver for
numerical simulations of shock/turbulence interactions should have very low numerical dissipation in the regions away from the shock and should also successfully capture the shock

5

without introducing any spurious oscillations. One way to achieve this goal is to use a hybrid scheme which incorporates a very low dissipation central-difference scheme in shock-free
regions and a shock-capturing algorithm in regions near the shock. The switch between the
two algorithms can be made by using a shock sensor. This hybrid methodology has been
successfully used for various shock/turbulence interaction simulations [17, 8, 18].
In the present study, we use a similar hybrid methodology to discretize the convective
terms in Eqs. (3.1), (3.2) and (3.3). In the shock regions, we employ a fifth-order MP scheme.
(It should be noted that the order of accuracy degenerates to first-order near discontinuities).
This scheme has proved to be very robust for a wide variety of problems and has lesser
numerical dissipation and faster grid convergence when compared to a standard WENO
scheme [11]. The details of MP scheme can be found in Ref. [11]. In smooth regions of the
flow, the convective terms are cast in a split form [19, 9, 18]:
∂uj ϕ
∂ρuj
∂ρuj ϕ 1 ∂ρuj ϕ 1
∂ρϕ
=
+
uj
+ρ
+ϕ
∂xj
4 ∂xj
4
∂xj
∂xj
∂xj
∂ϕ
∂u
∂ρ
1
ρuj
+ ρϕ k + uj ϕ
+
4
∂xj
∂xj
∂xj

(1.8)

where ϕ is unity for the continuity equation, velocity components (uj ) for the momentum
equations and H = γ/(γ − 1)p/ρ + uk uk /2 for the total energy equation. The casting
of convective terms in spilt form ensures global conservation of total kinetic energy in the
limit of inviscid, incompressible flow and ensures strong numerical stability. Pirozzoli [9] has
shown that the above split form for convective terms can be cast in a locally conservative
form:
∂ρuj ϕ
1 ˆ
ˆ
=
f
− fj−1/2
∂xj x=x
∆x j+1/2
j

6

,

(1.9)

and be approximated using a central-difference operator of the form:
L

Dfj =
l=1

al
f
− fj−l
∆x j+l

.

(1.10)

This is achieved by defining the numerical flux as:
L

ˆ
fj+1/2 = 2

l−1

al
l=1

(ρ, u, ϕ)j−m,l

(1.11)

m=0

where
1
(f, g, h)j,l = (fl + fj+l )(gl + gj+l )(hl + hj+l ) .
8

(1.12)

Since the above formulation is based on central differencing, there is no additional numerical
dissipation induced as in upwind schemes. Also, since the split form ensures strong numerical
stability, there is no need for the use of spectral-type filters [20] which may also introduce
artificial dissipation. In the present study, we use a sixth-order accurate central scheme
(L = 3) for the formulation of convective terms in split form.
The switch between MP scheme and the central-difference scheme is made using a modified form of Ducros sensor [8]:

Θ=

(
(

u)2 + |

u)2
, 0
× u|2 + (u∞ /δ)2

Θ

1

(1.13)

In the present study, we use a threshold value of Θ = 0.075 to distinguish between shock and
shock-free regions, i.e., for regions where Θ

0.075, MP scheme is implemented and in all

the other regions, the central-difference scheme is used. The hybridization is achieved by first
computing the values of Θ in the computational domain and pre-marking the critical regions.

7

Figure 1.1: Instantaneous density contours overlaid with modified Ducros shock sensor values. The critical regions identified by the shock sensor are shown in red. For interpretation
of the references to color in this and all other figures, the reader is referred to the electronic
version of this dissertation.

Next, the fluxes are computed in the entire domain using the central-difference scheme. In
critical regions, the fluxes are then selectively replaced by using the MP scheme. For the
treatment of interfaces between the shock and shock-free regions, a padding of 3 points on
either side of the interface point is used to prevent the incursion of non-dissipative centraldifference scheme into the shock regions [8, 21]. In Fig. 1.1, we show the instantaneous density
contours overlaid with the shock sensor values in an x − y symmetry plane for θ = 7.75o .
The regions highlighted in red are identified as critical regions by the shock sensor and a
switch is made to MP scheme for computing the convective fluxes in these regions.
The viscous and diffusion terms are cast in a “Laplacian” form [22] for proper representation of diffusion effects at the highest resolved wavenumbers and also for improved numerical

8

stability. We found that when using a non-dissipative scheme for discretizing the convective
terms, as done in the present study, it is necessary to use the non-conservative form for discretizing the viscous and diffusion terms which otherwise could lead to numerical oscillations
in the solution, particularly near the wall. The viscous and diffusion terms are computed
using a sixth-order central compact scheme [20] for both first and second derivatives. The
time integration is performed using an explicit third-order Runge-Kutta scheme.

1.2.3

Inflow turbulence generation

The generation of a realistic inflow boundary condition for turbulent boundary layers is
very important. Several methods are described in the literature, each with their own merits
and demerits. In the present study, we adopt a two-step approach to prescribe the inflow
turbulence for SBLI simulations. This is similar to the method used in Refs. [16, 15, 23].
First, we first perform an auxiliary DNS for a canonical zero-pressure gradient turbulent
boundary layer. The domain size for the auxiliary simulation was 37.7δref × 15δref × 3.5δref
in the streamwise (x), wall-normal (y) and spanwise (z) directions, respectively. Here, δref
is the boundary layer thickness at the inlet of the auxiliary DNS. The useful part of the
domain in the streamwise direction was only until x = 24δref after which is the sponge
region [24] which is used to drive the flow to a uniform state towards the exit of the domain.
The corresponding grid resolution for the auxiliary simulation was 1464 × 202 × 193. The
grid was uniform in the streamwise direction until x = 24δref and stretched in the sponge
region. The grid was also uniform in the spanwise direction. In the wall-normal direction,
grid points were clustered near the wall using a hyperbolic sine transformation. In terms of
wall units, the streamwise and spanwise grid spacings were ∆x+ = ∆z + ≈ 4.5 (the reference
quantities are taken at x = 23δref ). In the wall-normal direction, 100 points were located
9

within y/δref = 1 and with the node next to the wall located at a distance of ∆y + ≈ 0.6.
For prescribing the inflow boundary condition for the auxiliary DNS, we used the recyclingrescaling method similar to the one described in Ref. [25]. The recycling station was located
at a distance of 23δref from the inlet of the domain. Once the flow reached a statistically
steady state, a time-history of all the primitive variables were stored in a y − z plane located
at a distance of 23δref from the inlet for a total length of approximately 360δref /U∞ . This
“inflow box” was then used as the inflow boundary condition for SBLI simulations. The
total length of the inflow box (360δref /U∞ ) was sufficient to carry out the SBLI simulations
without having to recycle the inflow box. To establish that the recycling procedure used
for auxiliary DNS does not introduce any spurious periodicity, we show the spectra of wall
pressure at the recycling station in Fig. 1.2. Evidently, we do not see any phase locking due
to the recycling procedure. As explained in Ref. [25], if a phase locking were to occur, a
peak in the spectra would appear at frec = Uc /xrec where Uc ≈ 0.8U∞ and xrec = 23δref .
This would correspond to frec = 0.0348U∞/δref (indicated by a vertical line in Fig. 1.2).
At the outflow of auxiliary DNS, a supersonic exit condition was used. At the top
boundary, a symmetry boundary condition was applied to all the primitive variables and a
periodic boundary condition was used in the spanwise direction.

1.2.4

Computational setup for DNS of SBLI

The primary motivation for using the inflow box method in our DNS computations was to
keep the incoming flow identical for all three selected flow deflection angles. Additionally,
this inflow box can also be used for other DNS and LES computations. The computational
parameters for the SBLI simulations conducted in this paper for various shock intensities
are given in Table 1.1. In this table, P2 /P1 is the pressure ratio across the incident shock.
10

−3

10

−4

E(f )

10

−5

10

−6

10

−7

10

−3

10

−2

10

−1

0

10
10
f δ ref /U∞

1

10

Figure 1.2: The wall pressure spectra at x = 23δref for the “auxiliary” DNS.

For all simulations, the grid is uniform in streamwise and spanwise directions. Similar to
the auxiliary DNS, the useful part of the domain in the streamwise direction extends only
up to x = 23.5δref after which is the sponge region. In the wall-normal direction, the
nodes are clustered near the wall using a hyperbolic sine function to resolve the boundary
layer. It should be noted that the inflow box generated from the auxiliary DNS had 202 and
193 grid points in the wall-normal and spanwise directions, respectively. However, for the
SBLI simulations, we increased the number points in both y and z to improve the resolution
(See Table 1.1). Hence, the inflow box was interpolated to the SBLI DNS grid in y and
z directions using a third-order polynomial. Since the time step for SBLI computations is
controlled by the CFL condition, a third-order polynomial interpolation was used to obtain
11

Table 1.1: Simulation Parameters for DNS of SBLI.
Parameter

Case I

Case II

Case III

M∞
θ
P2 /P1
Lx
Ly
Lz
Nx
Ny
Nz
∆x+
+
∆yw
∆z +

2.75
6o
1.51
23.5δref
15δref
3.48δref
1440
252
234
4.1
0.6
3.45

2.75
7.75o
1.69
23.5δref
13δref
3.48δref
1632
302
234
3.45
0.6
3.45

2.75
9.25o
1.86
23.5δref
13δref
3.48δref
1632
302
234
3.45
0.6
3.45

the inflow values at an arbitrary time t. As in the case for the auxiliary DNS, a sponge
layer is used at the outlet to damp out the fluctuations, thus enabling the use of a simple
extrapolation for all variables. At the top boundary, a symmetry boundary condition is
applied to all variables. The bottom boundary is an isothermal wall with wall temperature
set equal to the adiabatic wall temperature of the undisturbed boundary layer. A no-slip
condition is used for the three velocity components and the density at the wall is computed
using the unsteady characteristic boundary condition [26, 27].

1.3

Results

The DNS results for three different shock intensities or deflection angles, θ = 6o , 7.75o and
9.25o are discussed in this section. For convenience, we refer to simulations with θ = 6o , 7.75o
and 9.25o as Cases I, II and III, respectively. The ensemble-averaged statistics reported
in this paper are obtained by averaging in time and the homogeneous spanwise direction.
12

Table 1.2: Boundary layer parameters for the incoming boundary layer.
δ0

δ1

δ2

Reδ2

Cf

1.484δref

0.4418δref

0.102δref

2070

2.42 × 10−3

The time integration was for at least one time period of the low-frequency oscillation of
the reflected shock system. The time period of the shock system was estimated a priori
using the data presented in Fig. 7 of Ref. [4] and assuming a characteristic frequency
of StL = 0.03, which is the typical value observed in most experimental and numerical
studies. For calculating/presenting various statistics, we use either Reynolds decomposition
( φ = φ − φ ) or Favre decomposition ({φ} = φ − φ , {φ} = ρφ / ρ ) of the variables.
In the discussion of results for SBLIs, it is customary to use a non-dimensional streamwise
coordinate defined as X ∗ = (x−Xo )/L. Here Xo is the nominal shock impingement location
and L is the length of the interaction region, defined as the region between the nominal shock
impingement location (Xo ) and the apparent origin of the reflected shock. As done in other
studies, the apparent origin of the reflected shock is obtained by extrapolating the mean
reflected shock to the wall. In the present study, the nominal shock impingement location
was kept fixed at Xo = 15δref for all three DNS calculations. Using X ∗ as the streamwise
coordinate, the interaction region is located between −1 ≤ X ∗ ≤ 0.
The flow parameters for the incoming (undisturbed) boundary layer are shown in Table 1.2. For evaluating these parameters, data at the reference station x = 9δref was chosen.
For the purpose of this study, the boundary layer thickness (δ0 ) is defined as the location
where the mean spanwise vorticity reaches 0.5% of its maximum value.

13

1.3.1

Accuracy of DNS

To assess the accuracy of our DNS calculations, we first compare the mean and turbulence
quantities of the incoming boundary layer with the DNS results of Pirozzoli and Bernardini [1]
which were at Mach 2. For this purpose x = 9.0δref , which is located just upstream of the
interaction region for Case III, is chosen as the reference station. Figure 1.3(a) shows the
Van-Driest transformed mean streamwise velocity profile of the incoming boundary layer.
The universal linear relationship, U + = y + , (shown using dash-dot line in the figure) is
exhibited by the incoming boundary layer for y +
between 30

y+

5 and there is a logarithmic behavior

60. The logarithmic behavior is in a very good agreement with the

log-law profile (shown using dashed line in the figure) given by Eq. (1.14)

+
UV D =

1
ln(y + ) + C
K

(1.14)

where

K = 0.41,

C = 5.1

The differences in the wake region are due to the differences in Reynolds numbers computed using momentum thickness and wall properties (Reδ2w ). In the present DNS calculations, Reδ2w ≈ 1000 while in the DNS calculations of Pirozzoli and Bernardini [1],
Reδ2w ≈ 1300. The density-scaled Reynolds stresses are shown in Fig. 1.3(b) in outer scaling. Once again, a very good agreement is obtained between our current DNS results and
DNS results of Pirozzoli and Bernardini [1].
The accuracy of the incoming boundary layer statistics does not necessarily guarantee

14

25
20

+
UV D

15
10
5
0 0
10

10

1

2

3

10

10

y+
(a)

8
6

Rij

4

R11

R33
2
0

R22
R12

−2
0

0.2

0.4
0.6
y/δ 99

0.8

1

(b)
Figure 1.3: (a) Van Driest transformed mean streamwise velocity profile and, (b) densityscaled Reynolds stresses of the incoming boundary layer. Symbols represent DNS data of
Pirozzoli and Bernardini [1]
15

the accuracy of results in the interaction region and also downstream of the interaction
region. Since there are no experimental data available in the range of current simulations, it
is very important to establish the numerical accuracy of present DNS calculations in different
ways. The numerical accuracy of current DNS calculations, both in the upstream and the
downstream regions of the interaction, were assessed by checking the sum of all the terms in
the turbulent kinetic energy equation, which in theory should be equal to zero for a converged
solution. The net sum of all the terms in turbulent kinetic energy equation was found to be
negligible compared to the leading order terms. These results, discussed in detail in a later
section, and all the other results obtained for different grids indicate the numerical accuracy
of our DNS data through out the simulated domain.

1.3.2

Flow Visualizations

The complex features of SBLI are shown through the streamwise density gradient contours
for Case III in Fig. 1.4. The strong amplification of turbulence behind the reflected shock is
apparent from the density gradient contours in this figure. Also notice that the flow appears
more fine-grained downstream of the interaction region, indicating the effect of shock on
different turbulent length scales. The formation of expansion fan (represented by negative
values of the density gradient) is also evident behind the reflected shock. The expansion fan
region is followed by a series of compression waves which arise due to the reattachment of
boundary layer which mimics the effect of a compression corner.
The effect of shock strength on the shock-boundary layer interaction can be qualitatively examined by looking at the instantaneous numerical schlieren visualizations of the
flow obtained from DNS data (Fig. 1.5). We use the density gradient parameter defined
in Ref. [13] for flow visualizations in Fig. 1.5. The dark regions indicate regions of high
16

Figure 1.4: Instantaneous streamwise densitygradient contours for Case III. Contour levels
are evenly spaced between -1 and 1.

density gradients. One can clearly observe that the origin of reflected shock shifts upstream
with increasing shock intensity and the reflected shock is followed by a series of stronger
compression waves. Also notice that the strength of the reflected shock increases with an
increase in shock intensity.
The overall behavior of the boundary layer can be examined by looking at the flow in a
plane parallel to the wall, i.e., in an x − z plane. It is recommended to examine the streamwise evolution of flow structures in a fixed y/δ plane [1]. However, the presence of shock
makes the computation of boundary layer thickness (δ) in the interaction region difficult.
Hence, we chose the mean sonic plane, i.e., the plane wall-normal location at which the
mean Mach number is equal to 1, instead of a fixed y/δ plane. Additionally, the mean sonic
plane approximately represents the beginning of the supersonic regime in the boundary layer.

17

y

10

5

0
0

5

10

x

15

20

15

20

15

20

(a)

y

10

5

0
0

5

10

x

(b)

y

10

5

0
0

5

10

x

(c)
Figure 1.5: Instantaneous numerical schlieren visualizations at an x − y symmetry plane. (a)
Case I, (b) Case II, (c) Case III.

18

Figure 1.6 shows the instantaneous velocity and temperature fluctuations along the mean
sonic plane. For the incoming boundary layer, the velocity and temperature fields exhibit
a typical streaky structure with alternating low and high momentum regions observed in
a canonical zero-pressure-gradient boundary layer. These structures can sometimes extend
to lengths of around 40δ and are referred to as superstructures. Some SBLI studies [28]
have postulated that these superstructures may be responsible for the low-frequency oscillation of the shock system. In the region upstream of the interaction region, the velocity
and temperature fields are negatively correlated, i.e., high-momentum regions correspond to
low-temperature regions and vice-versa. In the separation region (indicated by the two vertical lines), the streaky flow structure is observed in the initial part of the region but it soon
becomes “chaotic” and “isotropic” with a finer structure compared to the incoming flow.
From a numerical perspective, this is very important because conventional shock capturing
or upwind schemes cannot accurately capture the turbulence downstream of the interaction
region due to their inherent numerical dissipation. The shock appears to enhance mixing in
the boundary layer too. The anti-correlation between the velocity and temperature fields
drops in the interaction region and also downstream of the interaction region. For the incoming boundary layer, the correlation between velocity and temperature field is approximately
−0.85 along the mean sonic line but drops to −0.55 in the interaction region and regions
downstream of the shock. This indicates the weakening of Reynolds analogy assumptions in
the interaction region and regions downstream of the shock. We also examined the behavior
of velocity and temperature fluctuations in the x − z plane along which the turbulent kinetic
energy was maximum and found similar results.

19

y

−0.4

−0.2

0

0.2

2
0
5

10

x

15

20

(a)

y

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

2
0
5

10

x

15

20

(b)
Figure 1.6: (a) Instantaneous velocity fluctuations and, (b) instantantaneous temperature
fluctuations along the mean sonic plane for Case II.

20

1.3.3

Reynolds- and Favre-Averaging

For a turbulent variable f , the Reynolds-averaged value f is related the Favre-averaged
variable {f } by the following equation:
f − {f } = f

.

(1.15)

Both Reynolds- and Favre-averaged forms of the flow statistics may be used for analysis of
compressible flows. Huang et al. [29] compared Reynolds- and Favre-averaged quantities for
a supersonic channel flow at Mach numbers of 1.5 and 3.0 and observed small variations
between them in regions close to the wall. Here we examine the two forms of averaging for
different shock intensities before discussing the flow statistics for SBLI.
Figure 1.7 shows the difference between the two kinds of averaging for the mean streamwise velocity and temperature at X ∗ = −1.1, which corresponds to the incoming boundary
layer not yet affected by the shock. For the streamwise velocity, the maximum difference is
about −2.5% w.r.t. the Reynolds-averaged velocity and it is observed to be near the wall. For
the static temperature, the maximum difference is about 1.5% w.r.t. the Reynolds-averaged
temperature but the difference persists to a greater distance along the wall-normal direction.
The difference between Reynolds- and Favre-averaged total temperature was found to be
less than 0.4% w.r.t. the Reynolds-averaged quantity and can be considered negligible. In
the interaction region (X ∗ = −0.5), the maximum difference between Reynolds-averaged
and Favre-averaged velocity increases to 5% for all three shock intensities (Fig. 1.8(a)). For
the temperature (Fig. 1.8(b)), there is only a marginal increase to 1.75%. In summary, differences between the two forms of averaging appear to be largely independent of the shock
intensities and relatively small for the selected cases.
21

−0.5
−1
−1.5

u

/ u × 100

0

θ = 6o o
θ = 7. 75 o
θ = 9. 25

−2
−2.5
0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

y /δ 0
(a)

1.5
1
0.5

T

/ T

× 100

2

0
−0.5
0

0.2

0.4
y /δ 0
(b)

Figure 1.7: Differences between Reynolds- and Favre-averaged quantities at X ∗ = −1.1. (a)
velocity, (b) temperature.

22

0
−5
−10

u

/ u × 100

5

θ = 6o o
θ = 7. 75 o
θ = 9. 25

−15
−20
0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

y /δ 0
(a)

1.5
1
0.5

T

/ T

× 100

2

0
−0.5
0

0.2

0.4
y /δ 0
(b)

Figure 1.8: Differences between Reynolds- and Favre-averaged quantities at X ∗ = −0.5. (a)
velocity, (b) temperature.

23

Since the differences in the two forms of averaging are more significant near the wall,
it can influence the prediction of skin-friction coefficient (Cf ). In Fig. 1.9, the Cf values
computed from Reynolds- and Favre-averaged velocities for Case II are compared. As evident
from the figure, the differences mainly occur in the interaction region and downstream of
the interaction region.
The relationship between the Reynolds- and Favre-averaged turbulent stresses is given
by
{ui uj } = ui uj − ui uj + ρ ui uj / ρ

(1.16)

Consistent with supersonic channel flow DNS results [29], our results show that the term
associated with the product of two mean Favre-averaged fluctuations, ui uj is approximately 1% of the Favre-averaged stress value in all the regions of flow. On the other hand,
the triple correlation term ρ ui uj is not negliglible, particularly in the outer regions of the
boundary layer, as can be seen in Fig. 1.10. Although, it should be noted that the magnitude
of stresses is small in the outer regions of the boundary layer. Similar results are observed
for turbulent heat flux terms.

1.3.4

Flow Reversal Comparisons

Figure 1.11 shows instantaneous streamwise skin-friction coefficient (Cf x ) contours for the
three cases considered in this paper. The streamwise skin-friction coefficient is defined as

Cf x =

2σxy
2
ρ∞ U∞

24

(1.17)

−3

4

x 10

3

Cf

2
1
0
−1
−1.5

Cf
{C f }
−1

−0.5

0
X∗

0.5

1

1.5

0
X∗

0.5

1

1.5

( C f − {C f }) / C f × 100

(a)

5

0

−5
−1.5

−1

−0.5

(b)
Figure 1.9: (a) Comparison of Cf computed using Reynolds- and Favre-averaged quantities.
(b) Percentage difference in Cf computed by Reynolds- and Favre-averaged quantities. The
results shown are for Case II.

25

ρ u u / ρ {u u } × 100

0.05
0

−0.05
−0.1

−0.15
0

θ = 6o o
θ = 7. 75 o
θ = 9. 25
0.2

0.4

0.6

0.8

1

0.6

0.8

1

y /δ 0

ρ u u / ρ {u u } × 100

(a)

0.05
0

−0.05
−0.1

−0.15
0

θ = 6o o
θ = 7. 75 o
θ = 9. 25
0.2

0.4
y /δ 0
(b)

Figure 1.10: Comparisons for ρ ui uj . (a) for different shock intensities. X ∗ = −1.1, (b)
X ∗ = 0.25.

26

where σxy is the x − y component of the shear stress at the wall. Only regions with Cf x < 0
are shown to illustrate the extent of boundary layer separation. The locations of mean
separation and mean reattachment points in this figure are indicated by the dashed vertical
lines. The effect of shock intensity on the extent of flow separation is evident from the
larger separation regions created by shocks with higher intensities. The instantaneous data
show much larger separation zones than the corresponding mean separation zones (denoted
by dashed vertical lines). Additionally, the separation regions are not continuous and are
distributed non-uniformly in the spanwise direction. This shows that the separation bubbles
are highly unsteady and three-dimensional.
To further illustrate the 3D nature of flow in the interaction region, instantaneous spanwise skin-friction coefficient (Cf z ) contours are shown in Fig. 1.12. Cf z is defined as

Cf z =

2σyz
2
ρ∞ U∞

(1.18)

The incoming boundary layer is statistically homogeneous in the spanwise direction and
more-or-less “uniform” but we observe several patches of high Cf z values in the separation
region as well as in the regions downstream of the reattachment point. This indicates that
the flow is three-dimensional, at least instantaneously, not only in the separation region but
also downstream of the reattachment point, where the boundary layer is recovering back to
an equilibrium state.
For an additional characterization of the unsteady and intermittent nature of flow separation, we computed the statistical probability of flow reversal (γ) at the wall. According to
the classification by Simpson [30], incipient detachment (ID) occurs when the instantaneous
flow reversal happens about 1% of the time, intermittent transitory detachment (ITD) oc-

27

y

2
0
8

10

12

14
x

16

18

20

14
x

16

18

20

14
x

16

18

20

y

(a)

2
0
8

10

12

y

(b)

2
0
8

10

12

(c)
Figure 1.11: Instantaneous streamwise skin-friction coefficient contours. The contour levels
are evenly spaced between −6 × 10−3 and 0. (a) Case I, (b) Case II, (c) Case III.

28

y

2
0
8

10

12

14
x

16

18

20

14
x

16

18

20

14
x

16

18

20

y

(a)

2
0
8

10

12

y

(b)

2
0
8

10

12

(c)
Figure 1.12: Instantaneous spanwise skin-friction coefficient contours. The contour levels
are evenly spaced between −5 × 10−3 and 5 × 10−3 . (a) Case I, (b) Case II, (c) Case III.

29

1

0.6

TD

γ

0.8

θ = 6o o
θ = 7. 75 o
θ = 9. 25

0.4
IT D
0.2
0
−1.5

−1

−0.5

0
X∗

0.5

1

1.5

Figure 1.13: Comparisons for the probability of flow reversal for different shock intensities.

curs when the instantaneous backflow happens 20% of the time and transitory detachment
(TD), which in practice coincides with the mean flow reversal, occurs when the flow reversal
happens 50% of the time. The streamwise variation of γ for the three cases with different
shock intensities is shown Fig. 1.13. All three cases show a maximum flow reversal probability between 65 − 70%. As evident from the figure, with an increase in shock strength, the
probability of flow reversal extends to regions further downstream of the interaction region.
However, by increasing the strength of the incident shock, the maximum probability of flow
reversal does not change very much. This behavior was also observed by Morgan et al. [7]
in their LES studies.
Henceforth, for convenience, the streamwise skin-friction coefficient (Cf x ) will be referred
30

Table 1.3: Separation and Interaction lengths.
Parameter

Case I

Case II

Case III

Lsep
L/δ0

0.915δref
2.04

1.76δref
2.49

3.3δref
3.44

to as just skin-friction coefficient and will be denoted by Cf . Figure 1.14 shows the streamwise distribution of mean skin-friction coefficient for the three cases simulated with different
shock strengths. The abscissa of the plot has been modified to coincide the mean separation
points for all three cases and all three cases exhibit appreciable mean separation bubbles.
The separation lengths (Lsep ) and the interaction lengths (L) for Cases I, II, and III are
given in Table 3.3. Dupont et al. [4] proposed a linear relationship (for weak-to-moderate
interactions) between the length of the interaction region (L) and the shock intensity. The
interaction lengths computed in the present simulations agree well with the data presented
in Fig.7 of Ref. [4]. Interestingly, Fig. 1.14 shows that the rate of decrease of Cf ahead of
the mean separation point is similar for all three cases, even though the recovery rates vary.

1.3.5

Mean and Fluctuating Wall-Pressure Comparisons

The streamwise variations of mean wall-pressure for Cases I, II and III are shown in Fig. 1.15.
In Fig. 1.15(a), the mean wall-pressure is again plotted against the streamwise coordinate
whose origin has been shifted to the mean separation point. The mean wall pressure for all
three cases begin to rise upstream of the mean separation point, i.e., before x − xsep = 0.
Also, similar to the trends shown for Cf , the rate of increase of wall-pressure is similar for all

31

−3

4

x 10

3

Cf

2
1
θ = 6o
θ = 7. 75 o
θ = 9. 25 o

0
−1
−5

0

x − xsep

5

10

Figure 1.14: Comparisons for mean skin-friction coefficient for different shock intensities.

three cases during the initial part of the interaction. When the non-dimensional mean wall
pressure, P ∗ = (Pw − p1 )/(p2 − p1 ) (p1 and p2 are the pressures upstream and downstream
of the incident shock, respectively) is plotted against X ∗ (Fig. 1.15(b)), the profiles show a
reasonable collapse across and downstream of the interaction region. However, the collapse
becomes poor as the shock intensity increases [4] since the non-dimensionalization is valid
only for the inviscid limit.
Figure 1.16(a) shows the instantaneous wall-pressure fluctuations for Case II. In the
incompressible separated boundary layer DNS study by Na and Moin [31], the wall pressure fluctuations were found to be lower in the separated region compared to the region
downstream of the reattachment point. However, in the present study, we find that the
32

0.3
0.25

Pw

0.2
0.15
θ = 6o o
θ = 7. 75 o
θ = 9. 25

0.1
0.05
−5

0

x − xsep

5

10

(a)

3
2.5

P∗

2
1.5
1
θ = 6o o
θ = 7. 75 o
θ = 9. 25

0.5
0
−1.5

−1

−0.5

0
X∗

0.5

1

1.5

(b)
Figure 1.15: Mean wall pressure distribution for different shock intensities plotted with
different scalings.

33

wall-pressure fluctuations are high even in the separation region (between the dashed vertical lines in the figure). The streamwise evolution of r.m.s. of wall pressure fluctuations,
normalized by the mean wall pressure and plotted against X ∗ , are shown in Fig. 1.16(b).
The data show a reasonable collapse, particularly for Case II and Case III, in the regions
downstream of the interaction region. Also note that the peak of the r.m.s. wall pressure
fluctuation moves upstream with increasing shock strength when plotted in this scaling.

1.3.6

Thermodynamic Variables

Figure 1.17 shows the variation of r.m.s. values of thermodynamic variables along the wallnormal direction at different streamwise locations for Case II. At all the streawise stations,
the pressure fluctuations (normalized by the mean pressure) are nearly constant in the nearwall regions, as shown in Fig. 1.17(a). For the incoming boundary layer (X ∗ = −1.1), the
maximum pressure fluctuations intensity is approximately 4% of the mean value. However,
in the interaction region (X ∗ = −0.5), the r.m.s. value of pressure increases to 12%. The
spike in the r.m.s. pressure profile seen at this station is due to the shock, and as shown in
the figure, there is a very sharp increase in pressure fluctuations across the shock. As we
move to downstream stations, the intensity of pressure fluctuations decreases, although the
peak values keep shifting to the outer part of the boundary layer. The variations in static
temperature fluctuations (Fig. 1.17(b)) due to the shock are not as significant as those of
the pressure. However, for the incoming boundary layer (X ∗ = −1.1), the r.m.s. values of
temperature fluctuations are quite significant with a peak value reaching to about 12% of
the mean temperature. In the interaction region (X ∗ = −0.5) there is a modest increase
in the peak value of r.m.s. temperature fluctuations. Even the jump in the temperature
plot across the shock (indicated by the spike) is not as strong as that in the pressure plot.
34

y

−0.05

0

0.05

0.1

2
0
5

10

x

15

20

(a)

0.12

θ = 6o o
θ = 7. 75 o
θ = 9. 25

p w /P w

0.1
0.08
0.06
0.04
0.02
−1.5

−1

−0.5

0
X∗

0.5

1

1.5

(b)
Figure 1.16: (a) Instantaneous wall pressure fluctuations for Case II. (b) Streamwise evolution of r.m.s. of wall pressure fluctuations for different shock intensities.

35

In the downstream stations, the r.m.s. temperature fluctuation values drop to the values
comparable to those found in the incoming boundary layer. Also, similar to the r.m.s. of
pressure fluctuations, the peak values of temperature r.m.s. keep shifting towards the outer
part of the boundary layer as we move to downstream stations. The fluctuations in the
total temperature (not shown) are not negligible, even for the incoming boundary layer.
Interestingly, for the selected streamwise stations, the maximum value of total temperature
fluctuations occur at X ∗ = 0.25, and not at X ∗ = −0.5, which is in the interaction region.
And, unlike the fluctuations for pressure and static temperature, the fluctuations in total
temperature are not so much affected by the shock, which is evident from the absence of
spikes in the total temperature plots.
To get a better understanding of the global behavior of various thermodynamic variables
and their turbulent fluctuations, we show the contours of their r.m.s. in Fig. 1.18 for Case
II. Figures. 1.18(a) and 1.18(b) show contours of pressure and density, respectively. The
highest fluctuations in the thermodynamic variables are observed in the interaction region.
The oscillations in pressure and density are shown to propagate along the reflected shock and
the strength of the fluctuations propagating along the reflected shock are shown to increase
with increasing incident shock strengths. Although the intensity of pressure and density
fluctuations decreases as we go further along the reflected shock towards the freestream,
they could be partly because of the grid coarsening.
A convenient and important parameter for predicting the effects of compressibility in
supersonic turbulent flows is the turbulent Mach number, Mt . It is defined as the ratio of
twice the turbulent kinetic energy and the mean speed of sound,

Mt =

{ui ui }
.
c
36

(1.19)

p / p

0.2
0.15

X∗
X∗
X∗
X∗

=
=
=
=

− 1. 1
− 0. 5
0. 25
1. 5

0.1
0.05
0

−2

−1

10

10

0

10

y /δ 0
(a)

0.2

0.1

T

/ T

0.15

0.05
0

−2

−1

10

10

0

10

y /δ 0
(b)
Figure 1.17: R.M.S. of thermodynamic fluctuations at different streamwise locations for Case
II. (a) pressure, (b) temperature.

37

y

4

2

2

1
−3

0
5

10

15

x

20

x 10
0

(a)

y

4
2
0
5

10

15

x

20

0.08
0.06
0.04
0.02
0

(b)
Figure 1.18: Contours of turbulent fluctuations of thermodynamic variables for Case II. (a)
r.m.s. of pressure fluctuations, (b) r.m.s. of density fluctuations.

Another useful quantity is the fluctuation Mach number ( M ) which is the r.m.s. of the
Mach number fluctuations. Compressibility effects in turbulence are predicted to become
important when Mt exceeds 0.3. In Fig. 1.19, we show the profiles for Mt and M

at various

streamwise locations for Case II. Both Mt and M

exhibiting

show similar trends with M

slightly larger values. The incoming boundary layer shows a peak Mt value of approximately
0.3 and it increases to 0.55 at X ∗ = −0.5. As we go to regions downstream of the interaction
region, the peak value of Mt reduces, but the location of peak moves towards the outer part
38

of the boundary layer, similar to pressure and temperature fluctuations. Unlike Mt , M
experiences a sharp jump across the shock, which can be seen in its profile at X ∗ = −0.5.
Similar trends are observed for the other two cases with the difference being the magnitude
of peak values which are higher for stronger shocks, as expected.

1.3.7

Turbulent Kinetic Energy Budget

Useful information for turbulence modeling can be obtained by looking at various terms in
the turbulent kinetic energy (TKE) transport equation. The transport equation for turbulent
kinetic energy (k = 1/2{ui ui }) is given by
∂ ρk
=C +P − +T +V +Π+D
∂t

(1.20)

where C, P, , T, V, Π, D are the contributions from mean convection/advection, production,
viscous dissipation, transport, viscous diffusion, pressure dilatation and mass diffusion terms,
respectively. The expressions for these terms are given below.

C=−

∂ ρ {uj }k
∂xj

∂{ui }
P = − ρ {ui uj }
∂xj
= σij
T = Tt + Tp = −

∂ui
∂xj

∂ 1
∂
{ui ui uj } −
∂xj 2
∂xj
V =

∂ σij ui
∂xj
39

p uj

0.8

Mt

0.6
0.4
0.2
0

−2

−1

10

10

0

10

y /δ 0
(a)

0.8

M

0.6

X∗
X∗
X∗
X∗

=
=
=
=

− 1. 1
− 0. 5
0. 25
1. 5

0.4
0.2
0

−2

−1

10

10

0

10

y /δ 0
(b)
Figure 1.19: Turbulent and fluctuating Mach numbers at different streamwise locations for
Case II. (a) Mt , (b) M .

40

Π= p

D = ui

∂ui
∂xi

∂ σij
∂ p
−
∂xj
∂xi

The production term (P ) generates the turbulent kinetic energy through the mean velocity
gradient. This term is responsible for the energy transfer between mean flow kinetic energy and the turbulent kinetic energy. The transport term consists of two parts, turbulent
transport term (Tt ) and pressure transport term (Tp ). The contributions from the pressure transport term (Tp ) and the pressure dilatation term (Π) together are known as the
pressure work. The pressure dilatation term is responsible for the energy transfer between
the turbulent kinetic energy and the mean internal energy and is significant only when the
compressibility effects become important. The pressure transport term is responsible for the
redistribution of energy to the inhomogeneous flow regions. As stated in Ref. [32], the mass
diffusion term (D) is usually associated with variable inertia effects for flows with low turbulent Mach numbers. But at higher turbulent Mach numbers, this term can also include the
effects of compressibility. Note that the mass diffusion term involves the Reynolds-averaged
value of the Favre fluctuation of the velocity which, as we have shown earlier, is small and
primarily confined to regions close to the wall.
Figure 1.20 shows the turbulent kinetic energy budget for the incoming boundary layer,
i.e., at X ∗ = −1.1. Since the profiles are similar for all the three cases at this location,
we show the results only for Case II. The legend for various terms in the turbulent kinetic
energy equation is given in Table 1.4. The results at X ∗ = −1.1 are representative of a typical
zero-pressure gradient boundary layer. Near the wall, viscous dissipation is balanced by the
viscous diffusion term and away form the wall, the production term largely balances the

41

T KE budget

0.01
0.005
0
−0.005
−0.01
−2

−1

10

10

0

10

y/δ 0
Figure 1.20: Turbulent kinetic energy budget for Case II at X ∗ = −1.1.
Table 1.4: Legend and colors for turbulent kinetic energy budget plots.
Term

Color Line style

convection (C)
production (P)
dissipation ( )
turbulent transport (Tt )
pressure transport (Tp )
viscous diffusion (V )
pressure dilatation (Π)

green
black
red
blue
blue
black
black

solid
solid
solid
solid
dashed
solid line with open circles
open squares

viscous dissipation term and the turbulent transport term. The advection term is negligible
for the incoming boundary layer due to the absence of any significant mean strain rates in
the flow direction. The sum of all the terms in Eq. 1.20 is shown by the dashed red line and
it indicates excellent numerical accuracy of our DNS calculations.
The profiles at X ∗ = −0.5 are shown in Fig. 1.21. At this location, the pressure transport
42

term becomes important near the wall (indicating the redistribution of TKE to the inhomogeneous directions) and cancels the other two important terms, the viscous diffusion and the
viscous dissipation terms. The dissipation near the wall increases with generally higher dissipation for higher shock intensities. Since dissipation is largely associated with small scales
of the flow, larger dissipation indicates more small scales are formed due to the interaction
of the boundary layer with the shock. As a result of the boundary layer interaction with
the shock, the peak value of production term increases by almost four times which results in
the amplification of turbulent kinetic energy at this location. Also, the peak of production
term shifts away from the wall into the mixing layer region developing over the separation
bubble. As shown earlier, the height of the separation bubble increases with an increase in
shock intensity and this explains the shifting of the peak away from the wall with increasing
shock intensity. Away from the wall, the convection term becomes significant, again due to
the formation of mixing layer over the separation bubble, and helps in balancing out the
production term in the outer part of the boundary layer. The turbulent transport term has
both positive and negative contributions in the boundary layer. Near the wall, it has a positive contribution and compensates for the draining of the energy due to viscous dissipation.
Near the regions where the production term attains its peak, the turbulent transport term
is negative to counter the production term but once again its contribution becomes positive
in the outermost region of the boundary layer to balance out the convective term. Once
again note the excellent balance between all the contributing terms which is indicated by
the dashed red line.
At X ∗ = 0.25 ( Fig. 1.22 ), again the pressure transport term, the viscous diffusion term
and the viscous dissipation term cancel each other near the wall. For Case I ( not shown ),
the production term exhibits two peaks, one in the inner layer and the other in the outer
43

T KE budget

0.02
0.01
0
−0.01
−0.02

−2

−1

10

10

0

10

y/δ 0
(a)

T KE budget

0.02
0.01
0
−0.01
−0.02

−2

−1

10

10

0

10

y/δ 0
(b)
Figure 1.21: Turbulent kinetic energy budgets at X ∗ = −0.5. (a) Case II, (b) Case III.

44

layer, with the inner peak having a larger value. For Case II ( Fig. 1.22(a)) and Case III (
Fig. 1.22(b)), the inner peaks are just beginning to form. This indicates the varying rates
of boundary layer recovery towards an equilibrium state for different shock intensities. The
peaks of production in the outer part of the boundary layer are balanced by the turbulent
transport and viscous dissipation terms. The convection term is positive in the inner part
of the boundary layer to negate the effect of viscous dissipation but becomes negative in the
outer part to balance the production and turbulent transport terms.
At X ∗ = 1.5 ( Fig. 1.23), a recovery towards an equlibrium state is still in process with
the peak of production term moving back towards the wall, although the peak in the outer
layer still exists. The convection term becomes small again at X ∗ = 1.5 compared to the
previous two stations and the pressure transport term continues to be important near the
wall. The significance of the pressure transport term in the interaction region and also
downstream of the interaction region indicates the effect of shock in redistributing the TKE
in the inhomogeneous directions. In other words, the shock tends to make the flow more
homogeneous near the wall. The compressibility related terms like Π and D are found to be
negligible at all the four selected streamwise locations.
The behavior of TKE and various terms in the turbulent kinetic energy budget across
the interaction region can be examined by looking at their evolution along the mean sonic
line (Fig. 1.24). We show the results only for Case II since similar observations are made
for the other two cases. Figure 1.24(a) shows the evolution of TKE along the mean sonic
line. During the first half of the interaction region, TKE increases significantly, it then drops
sharply before increasing again. This behavior can be explained by examining various terms
in the TKE transport equation.
Along the mean sonic line, the production term in the incoming boundary layer is bal45

T KE budget

0.01
0.005
0
−0.005
−0.01
−2

−1

10

10

0

10

y/δ 0
(a)

T KE budget

0.01
0.005
0
−0.005
−0.01
−2

−1

10

10

0

10

y/δ 0
(b)
Figure 1.22: Turbulent kinetic energy budgets at X ∗ = 0.25. (a) Case II, (b) Case III.

46

T KE budget

0.01
0.005
0
−0.005
−0.01
−2

−1

10

10

0

10

y/δ 0
(a)

T KE budget

0.01
0.005
0
−0.005
−0.01
−2

−1

10

10

0

10

y/δ 0
(b)
Figure 1.23: Turbulent kinetic energy budgets at X ∗ = 1.5. (a) Case II, (b) Case III.

47

t u r b u l e nt k i n e t i c e n e r gy

2.5
2
1.5
1
0.5
−1.5

−1

−0.5
X∗

0

0.5

0

0.5

(a)

T K E B u d ge t

0.03
0.02
0.01
0
−0.01
−0.02
−1.5

−1

−0.5
X∗
(b)

Figure 1.24: (a) Streamwise evolution of turbulent kinetic energy along the mean sonic line.
(b) Streamwise evolution of turbulent kinetic energy budget along the sonic line.

48

anced by the turbulent transport and the viscous dissipation terms. Due to the interaction
of the turbulence with the shock, the production term increases in the interaction region.
The production term couples the mean flow kinetic energy and the turbulent kinetic energy.
Hence, a positive increase in the production term means energy is transferred from the mean
flow to the turbulent field (which explains the increase in TKE in Fig. 1.24(a)). The loss of
energy from the mean flow is due to the bulk flow deceleration caused by the shock-generated
adverse pressure gradient. The adverse pressure gradient also causes boundary layer separation and, hence the lift-up of vortical structures which leads to an increased contribution
from the convection term. In the first half of the interaction region, the convection and
viscous dissipation terms balance the contributions by the production and turbulent transport terms. At approximately X ∗ = −0.4, the production term becomes negative and the
convection term becomes positive. This is due to the formation of the expansion region in
the boundary layer that causes a local flow acceleration and the energy transfer from the
turbulent field to the mean field (This explains the sharp drop in TKE in Fig. 1.24(a)). The
pressure dilatation term is both negative and positive in the interaction region. The negative
values of Πd (which indicates a loss in turbulent kinetic energy and gain in mean internal
energy in this region) occur near the foot of the incident shock while positive values occur
around the expansion fan region (which indicates a loss of energy from the mean internal
energy to the turbulent field). Hence, the compressibility effects are confined only to a small
part of the interaction region. Similar to the pressure dilatation term, the pressure transport
term is also significant in the expansion and compression regions of the interaction zone. As
we go downstream of the interaction zone, production is once again balanced by turbulent
transport and viscous dissipation.
The streamwise evolution of various terms in the TKE transport equation along the locus
49

of maximum TKE (which in the interaction region would approximately correspond to the
centerline of the mixing layer) exhibit behavior similar to that along the mean sonic line
except for the pressure-dilation term whose contribution was found to be negligible along
the locus of maximum TKE.
The viscous dissipation term in the TKE equation can be decomposed into solenoidal
( s ), dilatation ( d ) and inhomogeneous components ( I ).

= s+ d+ I

(1.21)

where

= µ ωi ωi
∂uk ∂ul
4
µ
d=
3
∂xk ∂xl
2 uu
∂uj
∂ i j
∂
−2
ui
=2 µ
I
∂xi xj
∂xi
∂xj
s

Figure 1.25 shows the streamwise variation of viscous dissipation and its components along
the mean sonic line. As evident from the figure, the solenoidal part is the most significant
component of the three. Similar to the pressure-dilatation term in the TKE equation, the
dilatational dissipation is only noticeable in the interaction region. From Fig. 1.25, we also
see that the solenoidal dissipation, and the total dissipation, increases across the interaction
region which suggests the amplification of vorticity variance or enstrophy across the interaction region. The effect of shock on various terms in the enstrophy equation are discussed in
the next section.

50

v i sc ou s d i ssi p at i on

−3

3

x 10

s
d
I

2
1
0
−1
−2
−3
−4
−1.5

−1

−0.5
X∗

0

0.5

Figure 1.25: Evolution of dissipation term and its components along the mean sonic line.

1.3.8

Enstrophy budget

The transport equation for enstrophy (Ω = 1/2 ωi ωi ) is given below.

uj

∂u
∂ ui
∂ uk
∂Ω
= ωi ωj
+ ωi ωj i − Ω
−
∂xj
∂xj
∂xj
∂xk
I

Ω

II

IV

III

∂uk
ωi ∂ρ ∂p
∂ 1 ∂σkq
+ εijk 2
+ εijk ωi
∂xk
∂xj ρ ∂xq
ρ ∂xj ∂xk
V

VI

−

(1.22)

VII

∂u
∂u
∂
∂ ωi
Ωuj − ωi uj
+ ωj ωi i − ωi ωi k
∂xj
∂xj
∂xj
∂xk
VIII

IX

X

XI

In this equation, term I represents the convection term, term II is the production/destruction
of enstrophy due to the mean velocity gradient, term III is the self-stretching due to the

51

Table 1.5: Legend and colors for enstrophy budget plots.
Term

Color

Line style

I
II
III
VII
X

black
black
black
black
black

dotted
dash-dot
solid
dashed
solid with open triangles

fluctuating vorticity field, term IV is the vorticity-mean dilatation, term V is the vorticityfluctuating dilatation, term VI represents the generation of vorticity due to baroclinic torques,
term VII is the viscous diffusion term, term VIII is the turbulent transport term, term IX
is the production of enstrophy due to mean vorticity field, term X is associated with the
interaction of the fluctuating stretching term (ωi ui ) with the mean vorticity and term XI
represents the interaction of the mean vorticity with the stretching term associated with the
dilatation term.
Not all terms are important in the enstrophy transport equation before, within and after
the interaction zone. In Fig. 1.26, we show the dominant terms in the enstrophy transport
equation at X ∗ = −1.1. The legend for the enstrophy budget plots is given in Table 1.5.
Near the wall, terms III and X have positive contributions to the enstrophy budget and they
are balanced by terms II and VII. As we go away from the wall, the stretching term (X)
becomes zero and the production term (II) becomes positive to compensate for the stretching
term (X).
In the interaction region at X ∗ = −0.5 (Fig. 1.27), the two most dominant terms are
the self-stretching term (III) and the viscous dissipation term (VII) and they balance each
others effects. The increase of enstrophy (and viscous dissipation in the TKE equation)
52

50

0

−50

−2

−1

10

10
y/δ 0

0

10

Figure 1.26: Enstrophy budget for Case II at X ∗ = −1.1.
in the interaction region is primarily due to the stretching term (III). There is a small
contribution from the convection term (I) and as the shock strength increases, the production
term (II) becomes less significant. At X ∗ = 0.25 (Fig. 1.28), the stretching term (X) and
the production term (II) again become significant near the wall with varying degrees of
contributions for different shock intensities. The convection term (I) becomes negligible at
this location. As we go further downstream to X ∗ = 1.5 (Fig. 1.29), where the boundary
layer is in the process of reaching an equilibrium state once again, the dominant terms in
the enstrophy transport equation show behavior similar to the one observed at X ∗ = −1.1.
Again, terms III and X have positive contributions to the enstrophy evolution near the wall
and are balanced by terms II and VII. The contribution of the baroclinic term (VI) was found
to be negligible at all the selected streamwise locations. Similar to the pressure-dilatation
term (Π) in the transport equation for the turbulent kinetic energy, the baroclinic term (VI)
is also associated with the effects of compressibility.

53

100
50
0
−50
−100

−2

−1

10

10
y/δ 0

0

10

(a)

100
50
0
−50
−100

−2

−1

10

10
y/δ 0

0

10

(b)
Figure 1.27: Enstrophy budgets at X ∗ = −0.5. (a) Case II, (b) Case III.

54

100
50
0
−50
−100

−2

−1

10

10
y/δ 0

0

10

(a)

100
50
0
−50
−100

−2

−1

10

10
y/δ 0

0

10

(b)
Figure 1.28: Enstrophy budgets at X ∗ = 0.25. (a) Case II, (b) Case III.

55

100
50
0
−50
−100

−2

−1

10

10
y/δ 0

0

10

(a)

100
50
0
−50
−100

−2

−1

10

10
y/δ 0

0

10

(b)
Figure 1.29: Enstrophy budgets at X ∗ = 1.5. (a) Case II, (b) Case III.

56

1.4

Summary and Conclusions

Direct numerical simulations are carried out for shock waves of different intensities interacting
with a Mach 2.75 flat-plate turbulent boundary layer. The numerical method used for discretizing the Euler terms is based on a hybrid strategy that uses a sixth-order, non-dissipative
central-differencing scheme in the shock-free regions and a robust fifth-order monotonicity
preserving scheme in the regions around the shock. The accuracy of DNS calculations are
verified by comparing the incoming boundary layer statistics with the published literature
and also by checking the balance of all the terms in the turbulent kinetic energy transport
equation.
Visualizations of instantaneous flow velocity and temperature in the planes parallel to the
wall clearly indicate the significant and complex effects of shock on the turbulence structure.
As expected, the flow downstream of the shock becomes more chaotic with a fine-grained
structure and the mixing is also enhanced. The change of length scales downstream of the
shock is very significant from a numerical perspective since it places additional constraints
on the downstream grid and also the dissipative properties of the numerical scheme. The
anti-correlation between velocity and temperature fields are found to be reduced across the
interaction region indicating a breakdown of the strong Reynolds analogy. The differences
between Reynolds- and Favre-averaged velocities are more significant than those for temperature but are only important near the wall and in the interaction region, as reflected in
the skin-friction coefficient computations. The increase in shock intensity does not appear
to increase the differences between Reynolds- and Favre-averaging.
As expected, the separation region grows with increasing shock intensity due to an increased adverse pressure gradient. Instantaneous flow reversal analysis indicate that the flow

57

reversals are highly unsteady and three dimensional. The effect of shock strength on the
boundary layer interaction is clearly visible in the flow reversal comparisons. However, the
maximum probability of flow reversal does not increase with increasing shock intensity, an
observation also made in the LES studies by Morgan et al. [7]. The length of interaction
regions computed from the present DNS calculations show very good agreement with experimental results (Fig. 7 in Ref. [4]). The wall pressure fluctuations measured by the r.m.s. of
pressure show a reasonable collapse downstream of the interaction region when normalized
by mean wall-pressure and plotted against the non-dimensional streamwise coordinate, X ∗ .
The fluctuations of thermodynamic properties are intensified in the interaction region with
pressure showing the largest amplification in the interaction region.
The turbulent kinetic energy (TKE) budgets are examined for a better understanding of
shock effects on turbulence. Various terms in the TKE equation showed similar behavior for
different shock intensities when measured at the same non-dimensional streamwise locations,
X ∗ . The streamwise evolution of TKE and terms contributing to its evolution are also
examined along the mean sonic line. The increase in TKE in the interaction region is
primarily due to an increase in the production term. The production term is found to be
negative in the expansion fan region after the impinging shock, indicating a loss of energy
from the turbulence to the mean flow. This effect is stronger for higher shock intensities.
The pressure dilatation term and the dilatational dissipation are found to be significant
only in the regions surrounding the shock. The pressure dilatation term is negative near
the foot of the reflected shock which indicates an energy transfer from turbulence to the
mean internal energy. In the expansion fan region, the pressure dilatation is positive and
responsible for a net energy transfer from mean internal energy to the turbulence. All these
indicate a strong and dynamic coupling exists between the turbulent and mean fields in
58

the interaction region with energy being continuously exchanged between the two fields.
Despite the relatively high values of incoming Mach number and interacting shocks, the
compressibility-related terms like the pressure dilatation term are found to be negligible
in the TKE equation in most of the flow regions although the turbulent Mach number
exceeded 0.3, a conventional threshold for the turbulent compressibility effects to become
important. Examination of the enstrophy budgets also revealed that the baroclinic term,
a term associated with compressibility effects, was negligible even for the highest shock
intensity. The vortex stretching term in the enstrophy transport equation is found to be
primarily responsible for the net increase in enstrophy in the shock-turbulence interaction
region.

59

Chapter 2
A Priori Analysis of Subgrid-Scale
Turbulence for Shock-Boundary Layer
Interactions

2.1

Introduction

Large-eddy simulations (LES) are becoming increasingly popular for high speed flow computations due to the current restrictions on performing direct numerical simulations (DNS)
and the short-comings of Reynolds-averaged Navier-Stokes simulations (RANS). In LES,
only the large scales in the flow are explicitly computed and the small scales are modeled
using a subgrid-scale (SGS) model. Unlike the turbulence models used in RANS, which are
mostly of the two-equation type, the commonly used subgrid-scale (SGS) models for LES
are the zero-equation eddy-viscosity type models. These models, despite their simplicity,
may still perform well because they model only the small (or unresolved) scales which are
generally more homogeneous than the large turbulent scales.
Shock-boundary layer interaction (SBLI) is becoming a benchmark problem in the highspeed flow modeling community due to its unique characteristics. Firstly, the setup for SBLI
flow is simple enough for DNS and LES but it is of great practical significance. Secondly,

60

the interaction of shock with the boundary layer is very complex and requires high-fidelity
numerical methods like DNS and LES to correctly capture the flow physics. With the obvious
limitations of DNS, the focus is on LES to provide a reasonably accurate description of the
flow. The development of models for LES of high speed (supersonic) flows, including SBLI,
is a continuing research. There have been several LES computations of SBLIs. Garnier et
al. [33] were the first to conduct LES of SBLI for a flow Mach number of 2.3 and a flow
deflection angle of θ = 8o . Their LES results were in good agreement with the available
experimental data. Using the mixed time scale SGS closure, Touber and Sandham [34] were
able to reproduce the experimentally-observed characteristic frequency of the reflected shock
system. In another separate LES study [6], they compared their LES results with those from
three different experimental configurations. Pirozzoli et al. [35] examined the capabilities of
LES and RANS models to predict unsteady flow characteristics of SBLI. Jammalamadaka
et al. [16] evaluated the performances of several SGS models for the problem of SBLI by
comparing DNS and LES results and found that dynamic similarity-type models gave good
predictions. LES has also been used to study SBLI in compression corners. [36, 37]
In LES of incompressible flows, the only terms requiring models are the SGS stresses
which appear in the filtered momentum equations. Several types of models are employed for
the SGS stresses in incompressible flows, the most common being the zero-equation eddyviscosity type models. These models have also been used for modeling the SGS stresses in
compressible flows with varying degrees of success. However, in LES of compressible flows,
in addition to SGS terms in the filtered momentum equations, there are several terms in
the filtered energy equation which also require modeling. The modeling of SGS terms in the
energy equation has been a bit of a grey area for LES of compressible flows where some terms
are modeled while the others are neglected. There have been some studies that examined the
61

behavior of SGS terms in the filtered energy equation but these studies have been limited
to relatively simple subsonic flows [38, 39]. In complex flows like SBLI that involve both
shock and turbulence, the behavior of SGS terms in the LES equations is largely unexplored.
The main objective of the present study is to examine the behavior of various SGS terms
appearing in the compressible filtered Navier-Stokes equations for SBLIs. This is primarily
done through a priroi analysis of DNS data for SBLI [40]. The description of the DNS
database is described in the next section. The contribution of various components of SGS
fluxes in various regions of the flow is examined. We also examine the budgets for SGS
kinetic energy to study the effect of shock on subgrid-level turbulence and with the aim of
providing useful information for one-equation SGS LES models. To the authors’ knowledge,
this is the first comprehensive attempt to examine the behavior of SGS correlations/terms
in supersonic flows of the type seen in SBLIs.

2.2

DNS data

This study is based on a priori analysis of DNS data we have recently generated for SBLI [40].
The setup for the DNS of SBLI was a typical impinging oblique shock interacting with
a zero-pressure gradient turbulent boundary layer developing on a flat plate. The DNS
computations were performed for three different shock intensities to examine the effect of
shock intensity on SBLI. The incoming freestream Mach number for the turbulent boundary
layer was 2.75 and the flow deflection angles for impinging oblique shock were θ = 6o ,
7.75o and 9.25o . For the selected flow deflection angles, all three cases exhibited mean
flow separation. It was also shown in the DNS study [40] that the shock wave significantly
affects the mean and turbulent quantities and these effects were amplified for higher shock

62

intensities. With these results in view, we can expect the shock to have a similar effect on
various SGS terms.
The computational algorithm used for generating the DNS data was based on a hybrid
methodology. Hybrid numerical methodologies are popular for shock/turbulence simulations. In regions away from the shock, we used a minimally dissipative central scheme [9] to
compute the Euler terms in the Navier-Stokes equations. For enhanced numerical stability
and accuracy, the Euler terms were cast in a split-form [41, 9]. Interested readers are referred
to Refs. [9, 41] for a detailed description of the method. In regions surrounding the shock, a
fifth-order monotonicity-preserving (MP) scheme [10, 11] was used. Taking advantage of the
hyperbolic nature of the Euler terms, and also for improved numerical stability, the fluxes
were computed in the characteristic space and then transformed back to the physical space.
The flux-splitting in characteristic space was done using Lax-Friedrich’s method. Details on
characteristic transformation can be found in Refs. [42, 43, 11]. The MP scheme has proven
to be very robust for a wide variety of problems involving discontinuities like shocks and also
has lesser numerical dissipation in comparison to the standard WENO scheme [44]. However,
for the present DNS setup, MP scheme alone was not sufficient to resolve all scales in the
flow. The switch between the central scheme and MP scheme was made using a dilatationbased shock sensor [8]. More details on the hybridization can be found in Refs. [21, 8]. The
diffusive terms of the Navier-Stokes equations were computed using sixth-order compact central schemes [20]. The diffusive terms were cast in a non-conservative form [22] for enhanced
numerical stability and also to ensure finite molecular dissipation at all wavelengths. In our
experience, it is necessary to use the non-conservative form for the diffusive terms when the
Euler terms are cast in a split form and computed using very low dissipation central schemes.
Otherwise, unphysical oscillations appear in the solution, particularly in regions close to the
63

wall.
While the details of the DNS can be found in Ref. [40], some computational details are
given here for a better understanding of the results presented in this paper. The grid spacings
in wall-units for DNS of SBLI were as follows: in the streamwise direction, ∆x+ = 4.1 for
θ = 6o and ∆x+ = 3.45 for θ = 7.75o and 9.25o . In the homogeneous spanwise direction,
∆z + = 3.45. The grid was clustered near the wall and the node immediately next to the wall
was located at a distance of ∆y + ≈ 0.6. For the incoming boundary layer, Reynolds number
based on the momentum thickness and freestream quantities was 2070. The accuracy of the
incoming boundary layer statistics was verified by comparing them against results published
in the literature [1]. However, the accuracy of incoming boundary layer characteristics does
not necessarily guarantee the accuracy of solution in the interaction region and regions
downstream. In the absence of any experimental data, the accuracy of DNS calculations in
these regions was verified by grid-sensitivity assessment and by checking the sum of all terms
in the turbulent kinetic energy (TKE) equation. An accurate and converged solution would
result in a near zero sum of all the contributing terms in the TKE equation. For the DNS
calculations, the sum was indeed found to be negligible compared to the leading order terms
in the TKE equation in all the flow regions. In our opinion, this is a better way to ensure
numerical accuracy for DNS of turbulent flows when compared to the standard practice
of conducting grid refinement studies. Our studies showed that carrying grid refinement
studies with MP scheme (which is based on an upwind scheme) alone gave converged results
for various second-order statistics at ∆x+ = ∆z + = 4.1, but the contributing terms in the
turbulent kinetic energy equation were not found to be in full balance, even after further
refinement of the grid, possibly because of the numerical dissipation of the scheme. Although
increasing the grid resolution would theoretically increase our ability to capture smallest
64

scales in the flow, the numerical dissipation of the scheme may prevent us from doing so.
The hybrid method employed here [40] generated fully accurate data in all flow regions with
the selected grids.

2.3

Results

The results presented in this section are obtained from either using a single realization of
instantaneous filtered DNS data or an ensemble average of the data. For ensemble averaging,
averages were taken in time and homogeneous spanwise direction. The duration for temporal
averaging was T ≈ 15δ/U∞ , with δ being the boundary layer thickness of the incoming
boundary layer and U∞ the free stream velocity. The primary motivation to do ensemble
averaging was to reduce the numerical noise in the results for a better interpretation. We
stress that we were not looking to obtain fully converged statistics at low frequencies, which
would need much longer integration.
One of the parameters studied in this paper is the filter width. Three different filter
widths are considered: ∆ = 4∆DN S , ∆ = 8∆DN S and ∆ = 16∆DN S , where ∆ is the
characteristic size of the filter and ∆DN S is the characteristic size of the DNS grid. In
this study, filtering operations were performed mostly in the streamwise (x) and spanwise
(z) directions where the grid was uniform. Filtering in the wall-normal (y) direction was
avoided because of the complexities associated with filtering on non-uniform grids. It should
be noted that this is a standard practice for wall-bounded flows. Filtering was performed
in the physical space using a top-hat filter with the filtered variable, φ, obtained using the
relation:
1 ∆/2
φ=
φdξ
∆ −∆/2
65

(2.1)

The integral was evaluated using the trapezoidal rule and the resulting 1-D discrete filter is
given by:
l=N/2

φi =

al φi+l

(2.2)

l=−N/2

where the coefficients al are the corresponding coefficients of the trapezoidal rule and N = 4,
8 and 16 for ∆ = 4∆DN S , 8∆DN S and 16∆DN S , respectively. Multi-dimensional filtering
was performed by applying the above 1-D filter consecutively in each direction. For LES
of compressible flows, it is customary to use a Favre-filtered variable, φ, which is given by
φ = ρφ/ρ.
Since the simulated SBLI has two inhomogeneous directions (x and y), various SGS
terms are compared along the wall-normal direction at selected streamwise locations. The
streamwise stations were located at X ∗ = −1.1, −0.5, 0.25 and 1.5. Here, X ∗ is the nondimensional streamwise coordinate defined as X ∗ = (x − XI )/L, where XI is the nominal
shock impingement location and L is the length of the interaction region, defined as the
distance between the nominal shock impingement location and the extrapolated wall-location
of the reflected shock. The interaction region is then defined as the region between −1 ≤
X ∗ ≤ 0. Hence, X ∗ = −1.1 would correspond to the the incoming (pre-shock) boundary
layer, X ∗ = −0.5 is located in the middle of the interaction region and X ∗ = 0.25 and 1.5
are located downstream of the interaction region.

2.3.1

SGS terms in the filtered momentum equation

There are no unclosed terms in the filtered continuity equation when Favre-filtering is used.

∂
∂ρ
+
ρuj = 0
∂t ∂xj
66

(2.3)

However, unclosed terms appear in the filtered momentum equation which require modeling. The filtered momentum equation is obtained by applying the filtering operator to the
compressible Navier-Stokes equations. The dimensionless form is given by:
∂σij
∂τij
∂
∂
∂
∂p
(ρui ) +
=
−
+
ρui uj +
σ − σij ,
∂t
∂xj
∂xi
∂xj
∂xj
∂xj ij

(2.4)

τij = ρ(ui uj − ui uj )
σij =

2µ(T ) 1
Re
2

∂uj
∂ui
+
∂xj
∂xi

−

1 ∂uk
δ
3 ∂xk ij

σ ij =

2µ(T ) 1
Re
2

∂uj
∂ui
+
∂xj
∂xi

−

1 ∂uk
δ
3 ∂xk ij

where µ is the dynamic viscosity of the fluid and Re is the reference Reynolds number. The
tilde denotes Favre-filtered value, φ = ρφ/ρ. The last term in Eq. (3.9) is due to the nonlinearities and differences in Favre and regular filtering operations. This term was found
to be negligible by Vreman et al. [38] and also in the present study. Like in low-speed
flows, the most important unclosed term in the momentum equation is the SGS stress, τij .
The profiles of τij at various streamwise locations are shown in Fig. 2.1 and Fig. 2.2. The
results are shown only for θ = 7.75o and ∆ = 8∆DN S since similar trends are observed for
other shock intensities and filter widths. In all figures, the wall-normal coordinate (y) is
normalized by the incoming boundary layer thickness (δ0 ). At X ∗ = −1.1, the profiles look
very similar to the turbulent Reynolds stresses with τ11 being the dominant stress. Similar
to the turbulent Reynolds stresses, peaks for various SGS stresses occur close to the wall
and are unaffected by the shock at this streamwise location. In the interaction region at
X ∗ = −0.5, however, the effect of shock is significant on the SGS stresses. All the SGS

67

stress components are amplified and their peaks move away from the wall into the mixing
layer developing over the separation bubble. This trend is once again similar to that seen for
the turbulent Reynolds stresses. At X ∗ = 0.25, all the three normal SGS stress components
exhibit very similar peak values. While the peak value for τ11 decreases in magnitude from
X ∗ = −0.5 to X ∗ = 0.25, those for τ22 and τ33 are increased between these two locations.
The shock appears to have made the subgrid-level turbulence more isotropic at this location.
The reason for this behavior is discussed in a later section where the transport equations for
normal SGS stress are examined. At X ∗ = 1.5, the peak value for τ11 shifts back to the near
wall region, but those for τ22 and τ33 are still located away from the wall. Hence, we see a
very complex behavior exhibited by various components of the SGS stress tensor.
The effect of shock strength on SGS stresses can be better understood by looking at the
PDFs of different stress components. The PDFs were obtained using a single realization of
the filtered DNS data in the regions between 5.0 ≤ x ≤ 20.0, 0.0 ≤ y ≤ 1.0 and 0.0 ≤ z ≤ 3.5.
Figure 2.3 shows the PDFs for τ12 . It is clear that as the shock strength increases, the PDFs
become wider, indicating the probability for large SGS stress values to occur in the flow is
higher. Hence, for stronger shocks more energy is present at the subgrid level. This was
also observed in the R.M.S. values of SGS stresses for different shock strengths. For the
same filter width, the peak R.M.S. of SGS stresses were higher for stronger shocks in the
interaction region and downstream of the shock.
The SGS stress, τij , can be decomposed into generalized Leonard stress (Lij ), Cross
stress (Cij ) and Reynolds stress (Rij ) parts [45, 46].

τij = Lij + Cij + Rij ,

68

(2.5)

−3

3

x 10

τij

2
1
0
−1

−2

−1

10

10

0

10

y/δ 0
(a)
−3

8

x 10

6
τij

4
2
0
−2

−2

−1

10

10

0

10

y/δ 0
(b)
Figure 2.1: SGS stresses at different streamwise locations as obtained from the filtered DNS
data for θ = 7.75o and ∆ = 8∆DN S . (a) X ∗ = −1.1, (b) X ∗ = −0.5. Black solid line for
τ11 , black dashed line for τ22 , black dash-dot line for τ33 , black dots for τ12 , red solid line
for ksgs.

69

−3

6

x 10

τij

4
2
0
−2

−2

−1

10

10

0

10

y/δ 0
(a)
−3

6

x 10

τij

4
2
0
−2

−2

−1

10

10

0

10

y/δ 0
(b)
Figure 2.2: SGS stresses at different streamwise locations as obtained from the filtered DNS
data for θ = 7.75o and ∆ = 8∆DN S . (a) X ∗ = 0.25, (b) X ∗ = 1.5. Black solid line for τ11 ,
black dashed line for τ22 , black dash-dot line for τ33 , black dots for τ12 , red solid line for
ksgs.

70

0

10

−2

p df

10

−4

10

−6

10
−0.05

0
τ 12

0.05

Figure 2.3: PDFs of τ12 SGS stresse as obtained from the instantaneous filtered DNS data
for ∆ = 8∆DN S .

where

Lij = ρ ui uj − ui uj
Cij = ρ ui uj + ui uj − ui uj − ui uj
Rij = ρ ui uj − ui uj
Lij is closed as it is represented by only the resolved variables. Rij is only a function of
unresolved components and Cij is a function of both resolved and unresolved components.
Hence, when it comes to modeling SGS stresses, the mixed-type models [46] are expected to

71

perform better because they include the Leonard stress part in the model formulation.
The decomposition of τ11 stress for ∆ = 4∆DN S are shown in Fig. 2.4. At all selected
streamwise locations, L11 has the maximum contribution towards the total SGS stress. The
contribution of R11 is the smallest at X ∗ = −1.1 but increases as we go to downstream
stations. Since R11 is represented by only the unresolved scales, an increase in its value indicates generation of scales smaller than the characteristic filter width by the shock. However,
the contribution of C11 remains largely unchanged through all the regions of the flow for
filter width of the size ∆ = 4∆DN S . The shapes of all the individual components of the
SGS stress are very similar to the total SGS stress. This suggests similar kind of models,
say similarity-type model, can be used to model individual stress parts. The compressible
serial decomposition model [16] is one such example. LES calculations with filter widths of
the size ∆ = 4∆DN S can be termed as “fine-scale” LES. The LES studies in Refs. [34, 6]
are some examples of“fine-scale” LES.
The Leonard, Cross and Reynolds components of τ11 stress for ∆ = 8∆DN S are shown
in Fig. 2.5. Once again, at X ∗ = −1.1 and X ∗ = −0.5 (not shown), the Leonard stress
component has the largest contribution. However, at X ∗ = 0.25, the Reynolds part of the
total SGS stress has nearly the same contribution as the Leonard part. At X ∗ = 1.5, the
contribution of R11 exceeds that of L11 . This once again highlights the effect shock has
on the small scales in the flow and subsequently, the SGS stresses. For the filter width
∆ = 8∆DN S , accurate modeling of Rij becomes important to correctly predict the behavior
of SGS stresses, and consequently the flow. LES calculations conducted with filter widths of
the size ∆ = 8∆DN S can be termed as “coarse-scale” LES.
For very large filter widths of the size ∆ = 16∆DN S (Fig. 2.6), a very significant part
of turbulent scales are removed by the filtering operator. Hence, the unresolved scales
72

−3

1

x 10

0.8
τ11

0.6
0.4
0.2
0

−2

−1

10

10

0

10

y/δ 0
(a)
−3

4

x 10

τ11

3
2
1
0

−2

−1

10

10

0

10

y/δ 0
(b)
Figure 2.4: SGS stress decomposition at different streamwise locations as obtained from the
filtered DNS data for θ = 7.75o and ∆ = 4∆DN S . (a) X ∗ = −1.1, (b) X ∗ = −0.5. Black
solid line for τ11 , black dashed line for L11 , black dash-dot line for C11 , black dots for R11 .

73

−3

5

x 10

4
τ11

3
2
1
0

−2

−1

10

10

0

10

y/δ 0
(a)
−3

5

x 10

4
τ11

3
2
1
0

−2

−1

10

10

0

10

y/δ 0
(b)
Figure 2.5: SGS stress decomposition at different streamwise locations as obtained from the
filtered DNS data for θ = 7.75o and ∆ = 8∆DN S . (a) X ∗ = 0.25, (b) X ∗ = 1.5. Legend
same as in Fig. 2.4.

74

represented by Rij become very dominant in all the flow regions. At these filter widths,
the LES can be termed as “very-large-scale” LES or VLES. No model currently exists to
accurately capture the behavior of SGS stresses at such large filter widths. For even larger
widths, we would be approaching the realm of RANS computations.
We have seen the effect of filter width on the behavior of individual parts of the SGS
stress. It is evident that with change in filter width, different parts of the SGS stress become
dominant in different regions of the flow. For example, for ∆ = 8∆DN S , the Leonard stress
part is important mostly in the pre-shock region, while the Reynolds stress part is very
important in the shock and post-shock regions. Clearly, choosing the SGS stress model for
simulations of shock-turbulence interactions is not very straightforward and is dependent on
the filter width and shock strength.

2.3.2

SGS terms in the filtered energy equation

The energy conservation for compressible flow simulations can be achieved by solving different
forms of energy equations [47], e.g., the transport equation for internal energy, enthalpy or
the total energy. However, for flows involving shocks, it is often better to use the conservative
form of governing equations which makes the total energy equation an appropriate choice.
The filtered equation for total energy is obtained by applying the filtering operator to the
exact (DNS) total energy equation and subtracting the equation for SGS kinetic energy from
it [38]. The dimensionless form of the filtered total energy equation is given by:
∂ Qj
∂
∂
∂E
+
E + p uj =
σij ui −
− B1 − B2 − B3 + B4 + B5 + B6 − B7 , (2.6)
∂t
∂xj
∂xj
∂xj

75

0.02

τ11

0.015
0.01
0.005
0

−2

−1

10

10

0

10

y/δ 0
(a)
−3

8

x 10

τ11

6
4
2
0

−2

−1

10

10

0

10

y/δ 0
(b)
Figure 2.6: SGS stress decomposition at different streamwise locations as obtained from the
filtered DNS data for θ = 7.75o and ∆ = 16∆DN S . (a) X ∗ = −0.5, (b) X ∗ = 0.25. Legend
same as in Fig. 2.4.

76

where
E=

p
1
+ ρui ui
γ−1 2

Qj =

−µ
∂T
2
(γ − 1)RePrM∞ ∂xj

B1 =

1
∂
(puj − puj )
γ − 1 ∂xj

B2 = p

∂uk
∂u
−p k
∂xk
∂xk

B3 =

∂
(τ u )
∂xj kj k

B4 = τkj

B5 = σkj
B6 =

∂uk
∂xj

∂u
∂uk
− σkj k
∂xj
∂xj

∂
(σ u − σij ui )
∂xj ij i

B7 =

∂
(Q − Qj )
∂xj j

The term B1 in Eq. (3.10) is proportional to the divergence of SGS heat flux term (when
pressure is replaced by density and temperature using the equation of state), B2 is the
pressure-dilatation term, B3 is the SGS transport term, B4 is the SGS production (or sometimes also referred to as the SGS dissipation) and B5 is the SGS contribution to the molecular
dissipation. Terms B6 and B7 arise due to the differences between Favre and regular filtering
operations and were found to be negligible in the present study.
The modeling of B1 can be based on a similarity-type closure,

B1 =

∂
1
ρ uj T − uj T
2 ∂x
γ(γ − 1)M∞ j
77

(2.7)

or a gradient-diffusion type closure,

B1 = −

∂
∂xj

∂T
ρνt
2 ∂x
γ(γ − 1)P rt M∞ j

(2.8)

where γ is the ratio of specific heats, P rt is the SGS turbulent Prandtl number and νt is the
SGS eddy-viscosity. The model for the pressure-dilatation term, B2 , is often neglected in
many practical LES computations. Vreman et al. [48] proposed a similarity model for this
term as:
∂u
∂u
B2 = Cπ pij i − pij i
∂xj
∂xj

(2.9)

The SGS terms B3 and B4 are computed directly using the model for τij . For modeling the
dissipation term, B5 , Vreman et al. [48] proposed several different types of models:

(1)

B5

(2)

B5

= C1 σij

= C2 ρ
(3)

B5

∂ui
∂u
− σ ij i
∂xj
∂xj

,

k 3/2
, k = ui ui − ui ui /2 ,
∆

= C3 ρ

k 3/2
,k =
∆

3/4νt |S|

(2.10)

(2.11)
(2.12)

where k is the modeled SGS kinetic energy. The first model is based on the similarity idea
and the last two are based on the concept of balancing the production and dissipation of
SGS kinetic energy.
The object of discussion here is not to make an appraisal of different SGS models but to
give a quantitative estimate of each SGS term in Eq. 3.10 using the DNS data. In Fig. 2.7
we show various SGS terms in the filtered total energy equation at different streamwise
locations for θ = 7.75o and ∆ = 8∆DN S . The trends are similar for other shock intensities
78

Table 2.1: Color and legend for filtered total energy plots.
Term

Color

Line style

B1
B2
B3
B4
B5

black
black
black
black
blue

solid
dashed
dash-dot
dotted
solid

and filter widths. The divergence of filtered molecular heat flux (∂ Qj /∂xj ) is also plotted
in Fig. 2.7 for showing the relative importance of each SGS term. The choice of ∂ Qj /∂xj
as the reference is arbitrary. The net contribution by all SGS terms in the filtered energy
equation is shown by the dotted red line. This would signify the importance of modeling the
terms in the filtered total energy equation.
At X ∗ = −1.1, the two dominant terms are B1 and B3 which have opposite signs. The
net contribution by B1 , B2 , B3 , B4 , and B5 is very small compared to ∂ Qj /∂xj . This
suggested that for the incoming boundary layer which is a supersonic canonical turbulent
boundary layer, the modeling of B1 , B2 , B3 , B4 , and B5 terms in the total energy equation
can be ignored. This is confirmed in our a posteriori tests which will be discussed later.
At X ∗ = −0.5, which is located in the middle of the interaction region, the peak values
of most terms are comparable to that of ∂ Qj /∂xj . The values of SGS molecular dissipation
(B5 ) and the SGS production terms (B4 ) increase in the interaction region, which is an effect
of shock interacting with the subgrid-level turbulence. The overall contribution of all the
SGS terms has some significance in the inner and outer part of the boundary layer. The
net contribution follows the trend of B5 near the wall and in the outer layer, its peak is

79

0.01

0

−0.01

−0.02

−2

−1

10

10

0

10

y/δ 0
(a)

0.02
0.01
0
−0.01
−0.02

−2

−1

10

10

0

10

y/δ 0
(b)
Figure 2.7: Various SGS terms in the filtered total energy equation at different streamwise
locations as obtained from the filtered DNS data for θ = 7.75o and ∆ = 8∆DN S . (a)
X ∗ = −1.1, (b) X ∗ = −0.5. See Table. 2.1 for legend.

80

comparable to B1 . At X ∗ = 0.25 and X ∗ = 1.5 (not shown), the net contribution of SGS
terms, B1 , B2 , B3 , B4 , and B5 , is negligible in the outer part of the boundary layer and is
only significant near the wall. Hence, the shock affects each SGS term differently in different
flow regions, which is shown by the relative importance of the net contribution of B1 , B2 ,
B3 , B4 , and B5 in different flow regions.
To assess the importance of modeling SGS terms in the total energy equation in practical
LES calculations, we carried out two different a posteriori tests: one in which no models for
SGS terms in the filtered total energy equation were used (LES1) and in the other, models
for B1 , B3 , B4 and B5 were implemented in the LES calculations (LES2). In LES2, the
similarity model for the pressure-dilatation term (B2 ) given by Eq. (2.9) caused numerical
instability and, therefore, was not used. For the term B1 , the gradient-diffusion model
(Eq. (3.24)) is used and B5 was modeled using Eq. (2.11). To model the SGS stresses in
momentum equations, the mixed-time scale model [34] was used. This model has yielded
acceptable results for LES of SBLI and is also relatively inexpensive when compared to the
dynamic Smagorinsky model or the dynamic mixed model. A posteriori tests were performed
for two different shock intensities, θ = 7.75o and θ = 9.25o . The LES grid was four times
coarser than the DNS grid in the x and z directions. To avoid filtering in the wall-normal
direction, the grid in this direction was kept the same as the DNS grid. The setup for
LES calculations were similar to the one described in Ref. [16] where the inflow boundary
condition for both LES and DNS were kept consistent.
In Fig. 2.8, we show the skin-friction predictions by the two LES calculations for θ =
7.75o . The result from filtered DNS is also shown for comparison. As evident from the figure,
there is no significant effect of modeling the terms B1 , B3 , B4 and B5 in the energy equation
except for small differences in the interaction region. It should be mentioned here that by
81

−3

4

x 10

3

Cf

2
1
0
−1
5

10

x

15

20

Figure 2.8: Skin-friction coefficient predictions as obtained from a posteriori LES calculations
for θ = 7.75o and ∆ = 8∆DN S . Solid black line for filtered DNS, red dashed line for LES1,
blue dash-dot line for LES2.

ignoring the model for B5 , we observed some differences in Cf predictions between LES1
and LES2 downstream of the interaction region. Similar results were observed for θ = 9.25o .
Figure 2.9 shows the u u

Reynolds stress profiles for θ = 7.75o at X ∗ = −1.1 and

−0.5. The modeling of (or the lack of) SGS terms in the total energy equation apparently
has no significant effect on the Reynolds stress profiles at X ∗ = −1.1, as also observed in
our a priori analysis. However, there are some differences between the LES calculations
conducted with and without SGS models (LES1 and LES2) in the interaction region and
downstream of the interaction region. In fact, modeling the SGS terms in the total filtered
energy equation slightly deteriorates the results in these regions. This could be because the
82

model for the pressure-dilatation term (B2 ) was ignored, and our a priori results indicated
that this term is not negligible. The same behavior is also observed for the turbulent heat
flux term, u T , shown in Fig. 2.10. Evidently, LES calculations conducted with no models
for B1 , B3 , B4 and B5 in the filtered total energy equation give better results than those
conducted with some of these terms modeled.
The modeling of isotropic part of the SGS stress (τkk ) is often neglected in LES of compressible flows. To assess the importance of this term, we also carried out LES simulations
with and without models for τkk . In one simulation, τkk was assumed to be proportional
to Lkk (similarity-type model) and in another simulation, it was neglected. The results
obtained with the similarity model were similar to those obtained with no model for τkk .

2.3.3

Resolved and residual kinetic energies

The filtered kinetic energy can be decomposed into resolved kinetic energy (Er ) and residual
or SGS kinetic energy (ksgs ).
1
1
1
u i u i = u i u i + (u i u i − u i u i )
2
2
2
Er

(2.13)

ksgs

The ratio ksgs /(ksgs + Er ) provides us with the information on how much SGS energy is
present in the flow in comparison to the total filtered kinetic energy. We have computed
the ratio of ksgs /(ksgs + Er ) from SBLI DNS data for θ = 7.75o and different filter widths.
The results for ∆ = 4∆DN S at different streamwise locations are shown in Fig. 2.11. The
cross symbol on the plots indicates the location of peak values of ksgs . At X ∗ = −1.1, the
contribution of ksgs is only 2.5% near the wall and just 1% at the location of peak ksgs . This

83

0.02

uu

0.015
0.01
0.005
0
0

0.5
y

1

0.5
y

1

(a)

0.06

uu

0.04
0.02
0
−0.02
0

(b)
Figure 2.9: Turbulent Reynolds stress at different streamwise locations as obtained from a
posteriori LES calculations for θ = 7.75o and ∆ = 8∆DN S . (a) X ∗ = −1.1, (b) X ∗ = −0.5.
Solid black line for filtered DNS, red dashed line for LES1, blue dash-dot line for LES2.

84

0.01

uT

0
−0.01
−0.02
−0.03
0

0.5
y

1

0.5
y

1

(a)

0.02

uT

0
−0.02
−0.04
−0.06
0

(b)
Figure 2.10: Turbulent heat flux at different streamwise locations as obtained from a posteriori LES calculations for θ = 7.75o and ∆ = 8∆DN S . (a) X ∗ = −1.1, (b) X ∗ = −0.5.
Legend same as in Fig. 2.9.

85

indicates that at pre-shock locations, the flow is over-resolved by LES with ∆ = 4∆DN S .
At X ∗ = −0.5, the ratio increases to 15% near the wall and 5% near the maximum ksgs
location which are acceptable values for LES. In further downstream locations, the ratio
drops again but is still higher than that of the incoming boundary layer values. With larger
filter widths, the contribution of ksgs to the total filtered kinetic energy increases as expected.
For example, with ∆ = 8∆DN S (Fig. 2.12), the ratio ksgs /(ksgs + Er ) at the location of
peak ksgs increases from 2.5% to 12% as we go from X ∗ = −1.1 to −0.5. For ∆ = 16∆DN S
(Fig. 2.13), the contribution of ksgs to the total filtered kinetic energy in the interaction
region is as high as 60% near the wall. The results in Figs. 2.11, 2.12 and 2.13 elucidate
the significant effect of shock on the distribution of SGS energy in different regions of the
flow. Also note that these results were obtained using ensemble-averaged filtered DNS data.
Instantaneous filtered DNS data can give even higher ratios. With such significant changes
occurring in the distribution of SGS kinetic energy in the flow, it is important for an LES
model to accurately predict the behavior of ksgs .
The effect of shock strength on SGS kinetic energy is evident in the contours of ksgs /(ksgs +
Er ) in Fig. 2.14. The results are shown for ∆ = 8∆DN S . Note that the y-axis is in logarithmic scale for a better illustration of the near-wall results. As observed in Fig. 2.14,
the maximum value for the ratio is approximately 0.4 in the interaction region for all three
shock strengths (The slightly higher values predicted for θ = 6o compared to θ = 7.75o and
θ = 9.25o is possibly due to different grids used for cases with different shock angles). Hence,
the filter width is a more important parameter than the shock strength for the SGS kinetic
energy.
To better understand the behavior of SGS kinetic energy, we examine the transport
equation ksgs below. This would also provide useful information for the one-equation SGS
86

ksgs /(ksgs + Er )

0.025
0.02
0.015
0.01
0.005
0

−2

−1

10

10

0

10

y/δ 0
(a)

ksgs /(ksgs + Er )

0.2
0.15
0.1
0.05
0

−2

−1

10

10

0

10

y/δ 0
(b)
Figure 2.11: The ratio of SGS kinetic energy to the total kinetic energy (ksgs /(ksgs + Er ))
at different streamwise locations as obtained from the filtered DNS data for θ = 7.75o and
∆ = 4∆DN S . (a) X ∗ = −1.1, (b) X ∗ = −0.5.

87

ksgs /(ksgs + Er )

0.08
0.06
0.04
0.02
0

−2

−1

10

10

0

10

y/δ 0
(a)

ksgs /(ksgs + Er )

0.4
0.3
0.2
0.1
0

−2

−1

10

10

0

10

y/δ 0
(b)
Figure 2.12: ksgs /(ksgs + Er ) at different streamwise locations as obtained from the filtered
DNS data for θ = 7.75o and ∆ = 8∆DN S . (a) X ∗ = −1.1, (b) X ∗ = −0.5.

88

ksgs /(ksgs + Er )

0.2
0.15
0.1
0.05
0

−2

−1

10

10

0

10

y/δ 0
(a)

ksgs /(ksgs + Er )

0.8
0.6
0.4
0.2
0

−2

−1

10

10

0

10

y/δ 0
(b)
Figure 2.13: ksgs /(ksgs + Er ) at different streamwise locations as obtained from the filtered
DNS data for θ = 7.75o and ∆ = 16∆DN S . (a) X ∗ = −1.1, (b) X ∗ = −0.5.

89

0

0

10

10

−1

−1

y

10

y

10

−2

−2

10

5

10
10

x

15

20

5

10

(a)

x

15

20

(b)
0

10

−1

y

10

−2

10

5

10

x

15

20

(c)
Figure 2.14: ksgs /(ksgs + Er ) as obtained from the filtered DNS data for different shock
intensities with ∆ = 8∆DN S . (a) θ = 6o , (b) θ = 7.75o , (c) θ = 9.25o . The contours levels
are uniformly spaced between 0 and 0.5. The figures are shown not to scale for clarity.

models which make use of ksgs transport equation. The transport equation for the SGS
kinetic energy is given by:

90

Table 2.2: Color and legend for ksgs and normal SGS stress budget plots.
Term

Color

Line style

convection (I)
turbulent transport (II)
pressure-dilatation (III)
pressure-transport (IV)
viscous diffusion (V)
viscous dissipation (VI)
production (VII)
SGS diffusion (VIII)

black
red
green
blue
black
black
blue
green

solid
solid
solid
transport
dashed
dash-dot
dashed
dashed

∂ρksgs
∂
1 ∂
+
ρuk ksgs +
ρui ui uk − ρui ui uk
∂t
∂xk
2 ∂xk
I

=

II

σik

∂ui
∂u
− σik i
∂xk
∂xk
VI

∂uk
∂u
−p k
∂xk
∂xk

−

III

∂
∂
puk − puk +
σ u − σik ui −
∂xk
∂xk ik i
IV

p

(2.14)

V

− τik Sik +
V II

∂
uτ
∂xk i ik
V III

In this equation, term I is the convection term, term II is the SGS turbulent transport/diffusion, term III is the pressure-dilatation, term IV is the pressure-velocity transport,
term V is the viscous diffusion, term VI is the viscous dissipation, term VII is the SGS production (sometimes also referred to as SGS dissipation) and term VIII is the SGS diffusion.
Notice that terms III, IV, VI, VII and VIII also appear in the filtered total energy equation
(Eq. 3.10) simply because the transport equation for ksgs is used to derive the filtered total
energy equation. We once again use the filtered DNS results for θ = 7.75o and ∆ = 8∆DN S
91

to study the behavior of each term in Eq. (2.14) at selected streamwise locations. . The
budget for ksgs transport equation at various streamwise locations are shown in Fig. 2.15
and Fig. 2.16.
At X ∗ = −1.1, the two most dominant terms are the turbulent transport (II) and the SGS
diffusion (VIII) terms and and as shown in Fig. 2.15 they cancel out each others contributions
to the SGS kinetic energy. Knight et al. [49] proposed the following model for computing
the turbulent transport term (II):

1
ρ ui ui uk − ui ui uk ≈ ui τik
2

(2.15)

which is exactly the same as the SGS diffusion term (VIII). Hence, in the fully developed
supersonic boundary layers, terms II and VIII can be neglected together in the modeling of
transport equation for ksgs . The pressure-dilatation (III) is not negligible and has similar
contribution as the pressure-velocity term (IV) but with an opposite sign. Hence these two
terms also cancel out each others contributions for the incoming boundary layer. Near the
wall, the viscous diffusion term (V) balances the contribution from the viscous dissipation
(VI) term (similar to turbulent kinetic energy). The production (VII) term is responsible
for the generation of ksgs . A positive contribution from the production term implies energy
transfer from the resolved kinetic energy (Er ) or resolved field to the SGS kinetic energy
(ksgs ) or residual field (analogous to the energy transfer from mean and turbulent kinetic
energy fields in RANS) which is known as forward-scatter. It is, however, possible for the
energy to be transferred from the SGS/residual field to the resolved field, referred to as
backscatter. Under such circumstances, there would be a net negative contribution from the
production term (VII) to the ksgs transport equation.

92

−3

x 10
5
0
−5

−2

−1

10

10

0

10

y/δ 0
(a)
−3

x 10
5
0
−5

−2

−1

10

10

0

10

y/δ 0
(b)
Figure 2.15: SGS kinetic energy (ksgs ) budget at different streamwise locations as obtained
from the filtered DNS data for θ = 7.75o and ∆ = 8∆DN S . (a) X ∗ = −1.1, (b) X ∗ = −0.5.
See Table. 2.2 for legend.

93

−3

x 10
5
0
−5

−2

−1

10

10

0

10

y/δ 0
(a)

0.01
0.005
0
−0.005
−0.01
−2

−1

10

10

0

10

y/δ 0
(b)
Figure 2.16: SGS kinetic energy (ksgs ) budget at different streamwise locations as obtained
from the filtered DNS data for θ = 7.75o and ∆ = 8∆DN S . (a) X ∗ = 0.25, (b) X ∗ = 1.5.
See Table. 2.2 for legend.

94

At X ∗ = −0.5 in the interaction region, all the contributing terms to the ksgs are
amplified in their magnitudes. The production term (VII) is the most dominant term in the
interaction region which explains the increase in ksgs in this region (Fig. ??). Once again,
the turbulent transport (II) and the SGS diffusion (VIII) terms balance each other and hence
can be ignored together in the modeling of ksgs equation. The pressure-dilatation term (III)
is significant near the wall, and unlike the incoming boundary layer at X ∗ = −1.1, where it
is balanced by the pressure transport term (IV), the modeling of this term cannot be ignored
in the interaction region. The viscous dissipation term is also significant near the wall, where
it balances the contributions from pressure-dilatation (III) and the viscous-diffusion (VI),
and also away from the wall where it balances the production (VII) term. Hence accurate
modeling of this term is also important.
At X ∗ = 0.25, similar to X ∗ = −0.5, the terms which need modeling are the pressuredilatation (III), pressure-velocity (IV), viscous diffusion (V), viscous dissipation (VI) and the
production (VII) term. The contribution from the pressure-dilatation term is comparable
to the production term near the wall. At X ∗ = 1.5, the profiles of these two terms look
similar to the profiles at X ∗ = −1.1. At all the selected streamwise locations, the turbulent
transport (II) term and the SGS diffusion (VIII) term cancel each others contributions and
hence have a net zero effect on ksgs equation and can be ignored together.
Earlier we had shown the behavior of various SGS stress components in different regions
of the flow. It was observed that τ22 and τ33 SGS increase in magnitude from X ∗ = −0.5
to X ∗ = 0.25 while the τ11 component showed a decrease between these two locations. The
reason for this behavior can be examined by looking at the transport equation for the normal
SGS stresses (τii ). On examination of the transport equation for the normal SGS stresses,
it was observed that at X ∗ = −1.1 the production terms for τ22 and τ33 are negligible and
95

the entire production of ksgs comes from τ11 . At X ∗ = 0.25, the production term for τ22
and τ33 becomes significant and this explains their behavior at X ∗ = 0.25.

2.3.4

Backscatter

As discussed earlier, a positive contribution from the SGS production term to the SGS kinetic
energy equation indicates an energy transfer from the resolved field to the residual field. It
is, however, also possible for the energy to be transferred back from the residual field to
the resolved field (backscatter). The ensemble-averaged values of the SGS production term
indicated no significant backscatter, i.e., the SGS production term always had a positive
contribution to the SGS kinetic energy equation. However, instantaneous flow realizations
show that backscatter does indeed occur in various regions of flow. This can be observed
in Fig. 3.14 where the instantaneous contours of the SGS production term in a y − z plane
at four different streamwise locations are shown. The results in this figure are from the
instantaneous filtered DNS data for θ = 7.75o and ∆ = 8∆DN S . Only regions indicating
backscatter are shown. Figure 3.14 shows that significant backscatter does indeed occur
instantaneously. For example, at X ∗ = −1.1, backscatter occurs in 11% of the selected y − z
plane and at X ∗ = −0.5, it occurs in 16% of that flow region. However, no backscatter is
observed in the ensemble-averaged data. This might explain why constant-coefficient SGS
eddy-viscosity models, which cannot account for the energy backscatter, yield acceptable
results for many steady flow calculations. However, we would like to emphasize here that
shocks do have a very significant effect on the energy transfer between resolved and SGS
scales and generally increase the backscatter, mainly because they increase “small-scale”
turbulence.
To better understand the effect of shock on energy transfer and backscatter, we computed
96

y

2

y

2

1

0
0

1

1

z

2

0
0

3

1

(a)

z

2

3

2

3

(b)

y

2

y

2

1

0
0

1

1

z

2

0
0

3

(c)

1

z

(d)

Figure 2.17: Instantaneous contours of SGS production as obtained from the filtered DNS
data for θ = 7.75o and ∆ = 8∆DN S . Only the backscatter regions are shown. 50 uniformly
spaced contour levels between 0.0 and 0.005. (a) X ∗ = −1.1, (b) X ∗ = −0.5, (c) X ∗ = 0.25,
(d) X ∗ = 1.5.

the probability of backscatter in the simulated flow. Figure 2.18 shows the contours for the
probability of backscatter for θ = 7.75o and ∆ = 8∆DN S . Analogous to the definitions used
by Simpson [30] for identifying flow reversal, a probability greater than 50% would mean
backscatter occurs in the mean flow. As shown in Fig. 2.18 the maximum probability in
97

2.5
0.5

y

2

0.4

1.5

0.3

1

0.2

0.5

0.1

0
−1.5

−1

−0.5

0

0.5

1

0

X∗
Figure 2.18: Probability of backscatter in the x−y symmetry plane as obtained from filtered
DNS data for θ = 7.75o and ∆ = 8∆DN S .

majority of flow regions does not exceed the 40% value and only in the regions close to the
shock, particularly in the expansion fan region, the probability reaches close to 50%. This is
similar to the production term in turbulent kinetic energy equation where the the production
term had negative contribution in the expansion fan region [40]. Figure 2.18 also shows that
the shock increases the probability of backscatter in the interaction region (−1 ≤ X ∗ ≤ 0)
and in the regions downstream of the interaction. Hence the shock affects the dynamics of
energy transfer between the resolved and residual fields.

2.3.5

Directional filtering

In all the results discussed until now, filtering was performed only in x and z directions,
which is a common practice for wall-bounded flows. To examine the effect of shock on SGS
fluctuations in each spatial direction, it would be useful to examine the amount of SGS
98

energy present in each spatial direction. To study this, filtering operations were performed
in only one spatial direction, i.e., x-only and z-only directions. Additionally, filtering was
also performed in y− only direction using the 1-D filter function given by Eq. (2.2). The
filter coefficients were evaluated using the trapezoidal rule for non-uniform grids. It should be
noted that the 1-D filter (Eq. (2.2)) is a symmetric filter and is ill-defined at the wall. Hence,
near the wall, the grid points dictated by the filter width were ignored while computing the
filtered value. Nevertheless, the results discussed here are very useful and help us to better
understand the effects of directional filtering on LES quantities. The results are presented
for θ = 7.75o and ∆ = 8∆DN S and at X ∗ = −1.1, −0.5 and 0.25.
Figures 2.19 and 2.20 show the spatial variations of τ11 , τ22 , τ33 and τ12 SGS stresses
along the wall-normal direction at X ∗ = −1.1. For τ11 , the maximum SGS energy occurs
near the wall in y and z directions, as indicated by the y-only and z-only filtered values.
Very little energy is present in the streamwise direction which is evident in the x-only filtered
results. This indicates a slow growth of SGS turbulence in the streamwise direction. It also
explains why the grids are relatively coarse in the x direction compared to z and y directions
in LES of boundary layers (eg. Ref. [50]). While the overall τ22 stress values are small relative
to those for τ11 , important contributions are observed for the x-only and z-only filters. The
contribution from the y-only filter is negligible which indicates that the SGS turbulence does
not change significantly in the wall-normal direction for τ22 . For τ33 , while the stress values
are small relative to τ11 , they are still significant in all three spatial directions. This suggests
that τ33 is more homogenous than τ11 and τ22 . The peak for τ33 with z-only filtering occurs
near the wall, while the peak is away from the wall for the x-only filtered τ33 . The net effect
of a combined x − z filter results in the peak being shifted away from the wall, as seen in
Fig. 2.1(a). We should also mention here that the τ12 SGS stress is almost entirely due to
99

the SGS fluctuations in the spanwise direction.
In Figs. 2.21 and 2.22, we show the results for τ11 , τ22 , τ33 and τ12 SGS stresses in
the interaction region at X ∗ = −0.5. While the z-only filtered τ11 stress is still the most
dominant one, the contribution from the x-only filtered τ11 stress increases in the interaction
region, indicating that the SGS fluctuations for the streamwise velocity component increase
in the streamwise direction due to the shock. The near wall peak for y-only filtered τ11 stress
disappears at this location because of the separation bubble. The trends for τ22 and τ12 are
similar to those at X ∗ = −1.1. For τ33 , the x-only and z-only filtered stresses are of nearly
same magnitude, indicating that the z component of SGS stress is more homogeneous. At
X ∗ = 0.25 (Figs. 2.23 and 2.24), the trends for τ11 are same as at X ∗ = −0.5 except that
the near-wall peak for y-only filtered τ11 reappears and there is also a second peak away from
the wall. For τ33 , the x-only filtered stress values exceed the z-only filtered values, suggesting
that the streamwise SGS fluctuations are more important. In general, the shock seems to
homogenize the normal components of SGS stresses by increasing the SGS fluctuations in
the streamwise direction.

2.4

Conclusions

A detailed analysis of various subgrid-scale terms in the filtered equations for large-eddy
simulation (LES) of compressible flows is made using accurate DNS database for shockboundary layer interaction (SBLI) that resolves the flow both upstream and downstream
of the interaction region. The SGS stresses were examined in different flow regions and
their behavior were found to be generally similar to those of turbulent Reynolds stresses.
The Leonard, Cross and Reynolds parts of the SGS stresses were also examined in different

100

−3

2.5

x 10

τ 11

2
1.5
1
0.5
0
0

0.2
y /δ 0

0.4

(a)
−4

6

x 10

τ 22

4

2

0
0

0.2
y /δ 0

0.4

(b)
Figure 2.19: Directionally filtered SGS stresses as obtained from the filtered DNS data for
θ = 7.75o and ∆ = 8∆DN S at X ∗ = −1.1. (a) τ11 , (b) τ22 . solid black line for x-only filter,
dashed black line for y-only filter, red dash-dot line for z-only filter

101

−4

4

x 10

τ 33

3
2
1
0
0

0.2
y /δ 0

0.4

(a)
−4

2

x 10

τ 12

0

−2

−4
0

0.2
y /δ 0

0.4

(b)
Figure 2.20: Directionally filtered SGS stresses as obtained from the filtered DNS data for
θ = 7.75o and ∆ = 8∆DN S at X ∗ = −1.1. (a) τ33 , (b) τ12 . solid black line for x-only filter,
dashed black line for y-only filter, red dash-dot line for z-only filter

102

−3

8

x 10

τ 11

6
4
2
0
0

0.2
y /δ 0

0.4

(a)
−3

2.5

x 10

τ 22

2
1.5
1
0.5
0
0

0.2
y /δ 0

0.4

(b)
Figure 2.21: Directionally filtered SGS stresses as obtained from the filtered DNS data for
θ = 7.75o and ∆ = 8∆DN S at X ∗ = −0.5. (a) τ11 , (b) τ22 . solid black line for x-only filter,
dashed black line for y-only filter, red dash-dot line for z-only filter

103

−3

2.5

x 10

τ 33

2
1.5
1
0.5
0
0

0.2
y /δ 0

0.4

(a)
−4

5

x 10

τ 12

0
−5
−10
−15
0

0.2
y /δ 0

0.4

(b)
Figure 2.22: Directionally filtered SGS stresses as obtained from the filtered DNS data for
θ = 7.75o and ∆ = 8∆DN S at X ∗ = −0.5. (a) τ33 , (b) τ12 . solid black line for x-only filter,
dashed black line for y-only filter, red dash-dot line for z-only filter

104

−3

4

x 10

τ 11

3
2
1
0
0

0.2
y /δ 0

0.4

(a)
−3

3

x 10

τ 22

2

1

0
0

0.2
y /δ 0

0.4

(b)
Figure 2.23: Directionally filtered SGS stresses as obtained from the filtered DNS data for
θ = 7.75o and ∆ = 8∆DN S at X ∗ = 0.25. (a) τ11 , (b) τ22 . solid black line for x-only filter,
dashed black line for y-only filter, red dash-dot line for z-only filter

105

−3

2.5

x 10

τ 33

2
1.5
1
0.5
0
0

0.2
y /δ 0

0.4

(a)
−4

5

x 10

τ 12

0
−5
−10
−15
0

0.2
y /δ 0

0.4

(b)
Figure 2.24: Directionally filtered SGS stresses as obtained from the filtered DNS data for
θ = 7.75o and ∆ = 8∆DN S at X ∗ = 0.25. (a) τ33 , (b) τ12 . solid black line for x-only filter,
dashed black line for y-only filter, red dash-dot line for z-only filter

106

regions of the flow for three different filter widths, ∆ = 4∆DN S , 8∆DN S and 16∆DN S
and shock intensities. For the smallest filter width, the contribution of the Leonard stress
component was found to be the most dominant in all the flow regions, suggesting that the
mixed-type models could be appropriate for LES with such filter widths. For intermediate
filter width (∆ = 8∆DN S ), the Leonard stress component was dominant for the incoming
boundary layer while the Reynolds stress component’s contribution became very significant in
the regions downstream of the shock. These suggest that accurate modeling of the Reynolds
stress part of the SGS stress is important close to the shock. For very large filter widths
(∆ = 16∆DN S ), the Reynolds stress part of the SGS stress dwarfs the contributions from
the other two parts.
To further study the shock effect on compressible boundary layer flow, the SGS terms
in the filtered total energy equation were examined in various regions of the flow using
the ensemble-averaged DNS data. It was observed that while the individual unclosed SGS
terms in the filtered energy equation had sizable contributions by themselves, the combined
contribution from all the terms was relatively small when compared to the divergence of
filtered molecular heat flux term. To confirm this, two different a posteriori LES calculations
were performed in which SGS terms in the total energy equation were either neglected or
modeled via suggested models in the literature. The skin-friction coefficient predictions
between the two LES calculations showed no significant differences. Interestingly, comparison
of the turbulent stress and heat flux profiles show that the modeling of SGS terms in the total
energy equation slightly deteriorate the LES predictions in comparison to DNS. This could
be partly because the model for the pressure-dilatation term was ignored as it was found
to cause numerical instability to the LES calculations. Based on these observations, one
could argue that it is probably better to ignore the SGS terms in the total energy equation
107

altogether instead of modeling just a few terms. These results are, however, for the total
energy equation only and might not be valid for the internal energy or enthalpy equations.
The transport equation for SGS kinetic energy was examined using ensemble-averaged
values with the aim providing useful information to aid in the development of one-equation
LES models. The turbulent transport and the SGS diffusion term were found to balance
each others contributions in all the regions of the flow suggesting that modeling of these
two terms can be safely ignored for boundary layers and SBLIs. The contribution from the
pressure-dilatation term was not negligible and became more significant in the interaction
region and downstream of the interaction region. Hence, the modeling of this term cannot be
ignored, which is commonly is done. The ensemble-averaged production term was found to
have a positive contribution at all the selected streamwise locations and also in the majority
of flow region. This suggested there were no backscatter regions in the flow, where energy
is transferred back from the residual field to the resolved field. However, examination of
instantaneous flow samples indicated that backscatter does indeed occur in the flow. To
quantify the occurrence of backscatter, its probability was computed and was found to be
less than 50% in all of the flow regions with higher values in the expansion fan region.
Making an analogy to Simpson’s flow reversal definitions [30], less than 50% probability
indeed confirms that backscatter does not occur in the mean flow in the majority of flow
regions. The high probability of negative energy transfer from residual field to the resolved
filed in the expansion fan region is reminiscent of the behavior of production term in the
turbulent kinetic energy equation, which also has a negative contribution in the expansion
fan region.

108

Chapter 3
SGS Models and their Applications to
LES of Shock-Boundary Layer
Interactions

3.1

Introduction

Shock boundary layer interactions (SBLI) are important to a wide range of applications [51].
The shock wave can significantly alter the structure of supersonic boundary layers and depending on the strength of shock, the boundary layer may even separate in the vicinity of
shock impingement location. There is an amplification of turbulence across the shock and
in the case of boundary layer separation, the shock oscillates at a frequency lower than the
characteristic frequency of the incoming flow [4, 14, 12, 34, 13]. To understand and predict
the complex SBLI phenomena, high-fidelity simulations, i.e., direct numerical simulation
(DNS) and large-eddy simulation (LES), are needed. DNS is limited to simple geometries
and low-to-moderate Reynolds number flows. LES can be conducted at higher Reynolds
numbers and for flows in complex geometries but only computes the large scales and has
some limitations in wall-bounded flows [52]. Although the filtered equations for LES are
derived from Navier-Stokes equations with no major assumptions, there are several unclosed

109

terms appearing in the filtered momentum and energy equations. These terms are often
represented in terms of resolved quantities through subgrid-scale (SGS) closures. The use of
models to include the effects of small scales introduces an uncertainty into the computations.
Also, many of the SGS models were originally developed for incompressible flows. Their application to high speed flows, particularly the shock-turbulence interaction locations where
turbulence and shock effects are both important, has not been fully tested or justified.
The most common type of SGS models used for LES are zero-equation models developed based on the eddy-viscosity and isotropic SGS interaction assumptions which are not
necessarily valid in SBLI. Scale-similarity models [46, 53] and dynamic models [46, 54] are developed to replace or improve the original eddy-viscosity models. However, it is still not clear
how these eddy-viscosity, similarity or mixed models behave in shock-turbulence regions.
The problem of shock wave interacting with a boundary layer has primarily been studied
in two configurations, incident shock interacting with an already fully developed flat-plate
supersonic turbulent boundary layer [8, 55, 5, 4, 33] and shock-boundary layer interaction
in a compression corner [23, 13, 15]. The focus of the present study is on the incident SBLI.
An important characteristic of SBLI is the low-frequency oscillation of the shock system.
Although a few DNS studies [15, 14, 13, 12, 23] of low-frequency shock oscillations have been
reported, the computational cost associated with such simulations is still very significant at
reasonably high Reynolds numbers, Mach numbers and shock intensities.
LES method has also been used, particularly in recent years, for studying the main
characteristics of SBLI. Garnier et al. [33] were the first to conduct LES of SBLI for a
flow Mach number of 2.3 and a flow deflection angle of θ = 8 deg. They compared their
LES results with the available experimental data for similar flow conditions and found a
good agreement between the numerical and experimental results. More recently, Touber and
110

Sandham [34] studied the low-frequency oscillations in SBLI problem using LES with a new
inflow generation technique based on digital filters [56] and the mixed-time-scale closure for
SGS stresses. They also examined the effect of spanwise domain size on the separation bubble
and recommended that the spanwise width be at least the same as the separation bubble size.
Their results were in good agreement with the experiment. Morgan et al. [57] studied the
low-frequency behavior of shock motion using an improved rescaling/recycling method [58].
However, they used no SGS model in their LES studies, arguing that the artificial viscosity
they used for shock capturing mimicked the behavior of an SGS model. Morgan et al. [7]
also carried out extensive parametric studies of SBLI for various wedge angles, Reynolds
numbers, streamwise and spanwise domain sizes. They found that the Reynolds number did
not have any significant effect on the separation bubble size. Increasing the wedge angle
led to an increased size of the interaction region but did not necessarily guarantee higher
probability of reversed flow. Pirozzoli et al. [35] used LES and RANS models to study SBLI
under the conditions of incipient separation. By comparing with experimental results, they
found that both LES and RANS were able to qualitatively predict some basic features of the
interaction.
Given the importance and complexity of shock-turbulence interactions at small scales, the
SGS models are expected to play an important role in the flow development and quantitative
predictions of SBLI under various flow conditions. Our main objectives in this study are to
systematically assess the performance of several widely used SGS models along with a new
dynamic model for SBLI simulations and to evaluate the effects of shock and compressibility
on the models. Comparisons are made, both a priori and a posteriori, with the DNS data
for an incident 34o oblique shock interacting with a flat-plate turbulent boundary layer at
Mach 2. For these selected flow conditions, the boundary layer exhibits a marginal mean
111

flow separation.
The organization of the paper is as follows. In Sections 3.2 and 3.3, the governing
equations and SGS models are presented. In Section 3.4 we describe the numerical method.
The computational setup and parameters and simulation results are given in Section 3.5
followed by concluding remarks in Section 3.7.

3.2

Governing Equations

Written in conservative form, the non-dimesionalized form of governing equations solved in
DNS read as
∂ρ
∂
+
ρuj = 0
∂t ∂xj

(3.1)

∂σij
∂
∂
∂p
(ρui ) +
ρui uj = −
+
∂t
∂xj
∂xi
∂xj

(3.2)

∂Qj
∂E
∂
∂
+
(E + p) uj =
σij ui −
,
∂t
∂xj
∂xj
∂xj

(3.3)

where the viscous stress and heat flux terms are obtained from the following models:

σij =

∂ui ∂uj
+
∂xj
∂xi

2µ 1
Re 2

Qj =

−

1 ∂uk
δ
3 ∂xk ij

−µ
∂T
2 ∂x
(γ − 1)ReP rM∞ j

(3.4)

(3.5)

Here, the temperature dependence of µ is based on Sutherland’s law, µ(T ) = T 1.5 (1 +
C)/(T + C) with C = 0.65. In Eq. (3.3), E is the total energy defined as

E=

p
1
+ ρui ui
γ−1 2
112

(3.6)

Finally, the thermodynamic variables, ρ, p and T are linked by the ideal gas law,

p=

ρT
.
2
γM∞

(3.7)

The filtered equations for LES are obtained by applying the standard filtering and Favre
filtering operators to Eqs. (3.1), (3.2), (3.3) and (3.7). The final dimensionless form of the
filtered continuity, momentum and energy equations are

∂
∂ρ
+
ρuj = 0
∂t ∂xj

(3.8)

∂σij
∂τij
∂
∂
∂
∂p
(ρui ) +
+
−
+
ρui uj = −
σ − σij ,
∂t
∂xj
∂xi
∂xj
∂xj
∂xj ij

(3.9)

∂ Qj
∂
∂
∂E
+
E + p uj =
σij ui −
− B1 − B2 − B3 + B4 + B5 + B6 − B7 , (3.10)
∂t
∂xj
∂xj
∂xj
where the SGS stress term τij and the filtered viscous stresses σij σ ij and are given by

τij = ρ(ui uj − ui uj )

(3.11)

σij =

2µ(T ) 1
Re
2

∂uj
∂ui
+
∂xj
∂xi

−

1 ∂uk
δ
3 ∂xk ij

(3.12)

σ ij =

2µ(T ) 1
Re
2

∂uj
∂ui
+
∂xj
∂xi

−

1 ∂uk
δ
3 ∂xk ij

(3.13)

Qj =

∂T
−µ(T )
2 ∂x
(γ − 1)ReP rM∞ j

(3.14)

The last term on the RHS of Eq. (3.9) appears as a result of non-linearities in the viscous flux
terms due to filtering. Vreman et al. [38] found this term to be negligible for compressible
mixing layer. The same is assumed to be valid in the present study. The choice of total energy
113

equation for energy conservation, although is not unique, is appropriate for formulating the
Navier-Stokes equations in conservative form which is necessary for flows with discontinuities.
In Eq. (3.10), E is the Favre-filtered value of the modified total energy and the subgrid
scale terms B1 , B2 , B3 , B4 , B5 , B6 , and B7 are defined as

B1 =

1
∂
(puj − puj )
γ − 1 ∂xj

B2 = p

∂uk
∂u
−p k
∂xk
∂xk

B3 =

∂
(τ u )
∂xj kj k

∂u
B4 = τkj k
∂xj
∂u
∂u
B5 = σkj k − σkj k
∂xj
∂xj
B6 =

∂
(σ u − σij ui )
∂xj ij i

B7 =

∂
(Q − Qj )
∂xj j

Equations (3.8), (3.9), (3.10) are completed by the filtered state equation.

p=

ρT
2
γM∞

(3.15)

The terms which need modeling are τij , B1 , B2 , B3 , B4 , B5 , B6 , and B7 . The modeling of
these terms is described in the next section.

114

3.3

Subgrid-Scale Models

The filtered equations presented in previous section require closures for the unclosed SGS
terms. Several models have been proposed including the eddy-viscosity type ones in which
the deviatoric part of the SGS stress is modeled as

1
1
τij − τkk δij = −2ρνt
3
2

∂ui ∂uj
+
∂xj
∂xi

−

1 ∂uk
δ
3 ∂xk ij

(3.16)

To compute the isotropic part of the SGS stress, τkk , Yoshizawa [59] proposed the following
model,
2

τkk = 2CI ρ∆ |S|2 ,
where |S| =

(3.17)

2Sij Sij , ∆ is the characteristic filter width and CI is the model constant

which is typically chosen as 0.006 for isotropic turbulence or computed dynamically [54,
60]. However, Eq. (3.17) correlates poorly with the DNS data for compressible isotropic
turbulence [61].
In this study, four different SGS stress models, namely, the mixed-time-scale model
(MTSM), the dynamic Smagorinsky model (DSM), the dynamic mixed model (DMM) and
a new dynamic model, termed the compressible serial decomposition model (CSDM) are
examined. A brief description of each model is given below.
1. Mixed-Time-Scale Model (MTSM): The mixed time scale model has been used by Touber and Sandham [34] for LES study of low-frequency unsteadiness in SBLI. Although
this model gives an incorrect asymptotic behavior of eddy-viscosity near the wall, it
has been found to yield acceptable results for SBLI. The model details are given in the

115

following set of equations.
νt = CM T S kes TS

kes = ui − ui
−1
TS

=

∆
√
kes
ui =

ui − ui

−1

CT

+

(3.18)

−1

|S|

ρui
ρ

CM T S = 0.03, CT = 10

2. Dynamic Smagorinsky Model (DSM): In the dynamic Smagorinky model, the Smagorinsky coefficient, Cs , is computed dynamically [62]. The compressible form of the model
is given by,
2

νt = Cs ∆ |S|
Cs = −

(3.19)

1
(Gij − 3 Gkk δij )Mij
2
2∆ Mij Mij

(3.20)

In Eq. (3.20), ‘ ’ denotes the averaging in homogeneous direction and the tensors Gij
and Mij are computed by the following relations.

Gij = ρui uj −

Mij =

∆
∆

2

1
ρ|S| S ij − S kk δij
3

ρui ρuj
ρ
1
− ρ|S| Sij − Skk δij
3

The charactersitic width of the test filter, ∆, is chosen to be twice of the grid filter
width, ∆.

116

3. Dynamic Mixed Model (DMM): One of the known drawbacks of eddy-viscosity models is the inherent assumption that SGS stress and resolved strain rate tensors are
‘aligned’. This assumption is generally not correct, particularly in wall-bounded flows
where there is significant anisotropy. The model predictions also correlates poorly with
the DNS data. The dynamic mixed model (DMM) addresses these issues by including
the Leonard part of stress tensor in the model. This results in ‘lesser’ modeling since
Leonard stress is respresented only by the resolved quantities and a better representation of energy back-scatter to the resolved scales. Salvetti and Banerjee [46] and
Vreman [53] used the DMM model in their a priori analysis and a posteriori simulations of compressible (but subsonic) homogeneous turbulence and planar mixing layer.
The model is described in the following equations.

1
τij − τkk δij =
3

1
Lij − Lkk δij
3

∂uj
∂ui
+
∂xj
∂xi

1
2

− 2ρνt

−

1 ∂uk
δ
3 ∂xk ij

(3.21)

Lij = ρ ui uj − ui uj
2

νt = Cm ∆ |S|
Cm = −

1
Gij − 1 Gkk δij − Hij − 3 Hkk δij
3

Mij

2

2∆ Mij Mij
Hij = ρui uj −

ρui ρuj
ρ

4. Compressible Serial Decomposition Model (CSDM): The SGS stress tensor τij can be
decomposed into modified Leonard stress, modified Cross stress and modified Reynolds
stress components (Eq. 3.25). In the dynamic mixed model, the Cross term and the
117

Reynolds term are modeled together using an eddy-viscosity type model [46]. In the
compressible serial decomposition model, they are modeled using a similarity type
model following the same dynamic procedure used for DSM or DMM. For the sake
of brevity, we refrain from describing the details of the model in this paper. Further
details about this model can be found in Ref. [63, 64]. The basic formulation for CSDM
is given below.
τij = Lij + Cα Fij

(3.22)

Fij = ρ vi vj − vi v j
Cα =

Gij − Hij Nij
Nij Nij

Nij = ρv i vj −

ρv i ρv j
ρ

vi = ui − ui

In the filtered total energy equation, the term B1 is modeled as

B1 = −

∂
∂xj

ρνt
∂T
2 ∂x
γ(γ − 1)P rt M∞ j

(3.23)

for the eddy-viscosity type models. For the mixed type models which includes the dynamic
mixed model and the new compressible serial decomposition model, B1 is modeled as
∂
B1 = −
∂xj

ρνt
∂T
2 ∂x
γ(γ − 1)P rt M∞ j
L
qj = ρ uj T − uj T

118

L
∂qj
1
+
2
γ(γ − 1)M∞ ∂xj

(3.24)

For CSDM, we use MTSM to compute the eddy-viscosity used in modeling B1 . The turbulent
Prandtl number, P rt , is assumed to have a constant value of 0.89 for the simulated boundary
layer. Once τij is known, B3 and B4 are computed from their definitions. The other unclosed
terms, B2 , B5 , B6 , and B7 were neglected by Garnier et al. [33] in their LES of SBLI and
the same is done here.

3.4

Numerical Method

Numerical simulation of shock-turbulence interactions require high-order numerical schemes
for accurate computations of both the shock wave and turbulence. Since it is very difficult
to resolve the shock structure, the use of a shock-fitting scheme that treats the shock as a
discontinuity becomes necessary. In this study, we use a recently developed seventh-order
Monotonicity-Preserving (MP) [11, 10] scheme to discretize the Euler fluxes. The MP scheme
has proved to be less dissipative and more accurate than the standard WENO scheme [11]
for shock and turbulence calculations. To discretize the viscous fluxes, we use a sixth-order
compact scheme [20] with a fourth-order boundary formulation. The viscous terms are
discretized in conservative form. However, for the sake of comparison, we also performed
LES computations with the viscous terms discretized in Laplacian form [22] and found no
significant differences between the two forms of discretization. The governing equations
are integrated in time by a third-order low-storage Runge-Kutta scheme. The details of
numerical method are presented in Ref. [11] and also in Appendix A.

119

3.5
3.5.1

Results
Computational Setup and Parameters

The computational parameters for DNS and LES are listed in tables 3.1 and 3.2, respectively.

Table 3.1: Simulation Parameters for DNS of Mach 2 SBLI
Parameter

Value

M∞
Pr
Reδ

2.0
0.72
13500

Reδ2
Tw
Lx
Ly
Lz
Nx
Ny
Nz
∆x+
+
∆yw
∆z +
β
θ
XI
Cf

1475
1.717T∞
20δref
13δref
3δref
1333
252
201
3.89
0.65
3.89
34o
4.67o
10δref
0.0031

ref

In these tables, Reδ

ref

Comment

at x = 6

at x = 6
at x = 6
at x = 6

at x = 6

is the Reynolds number based on the reference boundary layer

thickness, δref (computed from DNS data as described below), and the freestream velocity,
density and viscosity. Reδ2 is the Reynolds number based on momentum thickness (δ2 ) and
the freestream quantities. Also in tables 3.1 and 3.2, x, y and z denote the streamwise,
wall-normal and spanwise directions, respectively. The grid is uniform in x and z directions.
120

Table 3.2: Simulation Parameters for LESof Mach 2 SBLI
Parameter

Value

M∞
Pr
Reδ

2.0
0.72
13500

Tw
Lx
Ly
Lz
Nx
Ny
Nz
β
θ
XI

1.717T∞
20δref
13δref
3δref
333
252
51
34o
4.67o
10δref

Comment

ref

1/4th of DNS
1/4th of DNS

In the wall-normal direction, a hyperbolic sine function was used upto y ≈ 5δref and then
the grid was smoothly varied to a uniform spacing up to the top boundary. In this work, the
wall-normal grid resolution is identical for both DNS and LES since filtering is performed
only in x and z directions. A sample x − y slice of the DNS grid is shown in Fig. 3.1. For
clarity, we only show every other 20th grid point in the x direction and every other 4th grid
point in the y direction. The grid in the blue region after x = 20 is coarse and indicates
sponge region.
When comparing LES results with DNS results, particularly when making model assessments as being done in the present study, it is important to minimize the effects of inlet
boundary condition on simulations. In order to keep the inflow boundary condition consistent for both DNS and LES, we first carried out an “auxiliary” DNS of a canonical turbulent
boundary layer with a rescaling-recycling inflow condition similar to that used by Pirrozzoli
et al. [25]. Once the flow reached a statistically steady condition, a time-history of primitive
121

Figure 3.1: An x − y slice of the sample DNS grid.
variables was collected at a y − z plane located at a distance of 23 boundary layer thicknesses from the inlet (this was also the location of recycling plane). A total of 5000 samples
were collected for a period of T = 75δref /U∞ which corresponded to approximately 3 flowthrough times and were periodically used as the inflow boundary condition for both DNS
and LES of SBLI. The inflow turbulence generation method employed in this work is similar
to that previously used in other SBLI simulations [23, 65, 36]. For LES, the inflow DNS
quantities were filtered in x and z directions using a box filter in physical space and then
down-sampled onto the LES grid. It should be noted here that the streamwise flow varibles
at the inflow boundary are obtained from a time-dependent periodic box data. Thus, for
filtering in x-direction, Taylor’s hypothesis was used to convert the temporal coordinate into
a spatial coordinate using a convection velocity of U∞ . A schematic of the computational
procedure is shown in Fig. 3.2.
At the outlet, a sponge layer [24] is used and supersonic exit condition is applied to
all variables. At the wall, no-slip isothermal wall with wall temperature (Tw ) close to the
adiabatic wall temperature is used. The wall pressure is computed by assuming that the
122

Figure 3.2: A schematic of the computational procedure for DNS and LES of SBLI.
pressure remains constant across the boundary layer, i.e., ∂p/∂y = 0. Symmetry boundary
condition is applied to all variables at the freestream and at the spanwise boundaries periodic
boundary condition is used. The shock is introduced at the inlet boundary by imposing
discontinuous conditions according to the Rankine-Hugoniot relations.
In order to evaluate the performance of SGS models for SBLI simulations, the flow
statistics at three different streamwise locations, (i) x = 6, upstream of the mean separation
point, (ii) x = 10, downstream but very close to the mean reattachment point and (iii)
x = 15, further downstream of the mean reattachment point, are compared. Here, all the
spatial coordinates are non-dimensionalized by the reference boundary layer thickness δref
unless otherwise stated.
Figure 3.3 illustrates a typical impinging shock-turbulent boundary layer interaction.
The shock is identified using the modified Ducros shock sensor [8] and the vortex structures
are identified using λ2 criterion [2]. The vortex structures are shown through isosurfaces
of λ2 = −3 and colored based on the local velocity magnitude. The effect of shock on
turbulence is clearly very significant.
123

Figure 3.3: Flow structure in a typical shock/boundary layer interaction configuration. Vortex structures identified using λ2 criterion [2] and colored based on the local values of instantaneous velocity magnitude.
To demonstrate the validity of our DNS results, Van Driest transformed mean streamwise
velocity profile of the incoming boundary layer is compared in Fig. 3.4 with zero-pressure
gradient Mach 2 DNS results of Pirozzoli and Bernardini [1] for two different Reynolds
+
numbers. The mean DNS velocity profile shows a linear relationship, UV D = y + , until

y + ≈ 6 and a logarithmic behavior (Eq. (1.14)) between y + ≈ 30 and y + ≈ 100. The
density-scaled profiles of Reynolds stresses are shown in Fig. 3.5 in wall units. Evidently,
there is a good comparison between the density-scaled profiles of the current DNS with those
of Pirozzoli and Bernardini. [1]. To further establish the numerical accuracy of current DNS

124

25
20

+
UV D

15
10
5
0 0
10

1

2

10

10

3

10

y+
Figure 3.4: Van Driest transformed mean streamwise velocity profile of the incoming boundary layer in inner scaling calculated from DNS data at x = 6. Black dashed line represents
+
UV D = y + and Eq. (1.14) is represented by black dash-dot line.
data, we computed various terms in the turbulent kinetic energy equation at different regions
of the flow and found them to be in balance.

3.5.2

A Priori Assessment

The four models described in Sec. 3.3 are tested a priori for the SBLI problem using a single
realization of instantaneous DNS data. Although a priori testing of SGS models provides
only a partial assessment of LES, it is still a good direct method of evaluating the models.

125

8
R11

6

Rij

4
R33
2
R22

0

R12

−2
0

10

1

2

10

10

3

10

y+
Figure 3.5: Density-scaled profiles of Reynolds stresses obtained by DNS at x = 6.
To make useful conclusions, filtered values from the DNS data are down-sampled onto an
LES grid four times coarser than the DNS grid in the streamwise and spanwise directions.
The modeled SGS fluxes are then computed on the LES grid and compared with the true
SGS fluxes. The filter function used for filtering the DNS data is a box filter in physical
space and filtering is performed only in x and z directions.
Following Germano’s [45] redefinition of stresses, the SGS stress tensor can be decom-

126

posed into modified Leonard’s stress, Cross stress and Reynolds stress components.

τij = Lij + Cij + Rij

(3.25)

where

Lij = ρ ui uj − ui uj
Cij = ρ ui uj + ui uj − ui uj − ui uj
Rij = ρ ui uj − ui uj
Figures 3.6 and 3.7 show the three parts of SGS stress tensor τ11 for two different filter widths,
∆ = 4∆DN S (Fig. 3.6) and ∆ = 8∆DN S (Fig. 3.7), where ∆DN S is the characteristic size of
DNS grid. Note that since the LES grid is four times coarser than the DNS grid in x and z
directions, ∆LES = 4∆DN S . Hence, ∆ = 4∆DN S would imply that the characteristic filter
width and the characteristic LES grid size are the same, i.e., ∆ = ∆LES . As expected, the
magnitude of subgrid-scale tensor τ11 increases with filter width. The modified Leonard’s
stress component (Lij ) and the modified Reynolds stress component (Rij ) have the largest
and the smallest contributions, respectively, when ∆ = 4∆DN S . However, the contribution
of Rij increases significantly and becomes comparable to Lij near the wall and to Cij away
from the wall when the filter width is doubled, i.e., ∆ = 8∆DN S . This is expected since
the Reynolds part of the total stress represents the unresolved scales and as the filter width
is increased, the contribution of unresolved scales also increases. Henceforth, in the results
and discussions of a priori assessment of SGS models we use ∆ = 8∆DN S since it presents
a greater challenge for SGS models to correctly capture the behavior of SGS stresses and

127

−3

2

x 10

SGS Stresses

1.5

1

0.5

0
0

0.2

0.4

y

0.6

0.8

1

Figure 3.6: Leonard, Cross and Reynolds parts of the spanwise-averaged τ11 SGS stress
component at x = 6 for ∆ = 4∆DN S . total stress - black solid line; L11 - dashed black line;
C11 - dash-dot red line; R11 - blue dotted line.
fluxes.
Figure 3.8 shows the comparison between the true and the modeled values of spanwiseaveraged trace-free τ11 SGS stress tensor. The reported τ11 values are normalized by the
local mean wall shear stress and are plotted against the wall-normal distance normalized by
viscous length scale. Since the eddy-viscosity models assume an alignment between the stress
and the strain-rate tensors and do not take into account the Leonard stress component, they
compare very poorly with the true stress. On the other hand, the DMM and CSDM account

128

−3

4

x 10

3.5

SGS Stresses

3
2.5
2
1.5
1
0.5
0
0

0.2

0.4

y

0.6

0.8

1

Figure 3.7: Leonard, Cross and Reynolds parts of the spanwise-averaged τ11 SGS stress
component at x = 6 for ∆ = 8∆DN S . total stress - black solid line; L11 - dashed black line;
C11 - dash-dot red line; R11 - blue dotted line.
for the Leonard’s term and, understandably, show a better agreement with the true stress
values. In fact, CSDM predicts the trend of true stress better than the DMM. These trends
are reiterated in Fig. 3.9 where we compare the spanwise-averaged contours of τ11 . For the
sake of brevity, we only compare DSM, DMM and CSDM predictions with the exact values
in this figure. It is evident that the shape and levels of isocontours predicted by DMM and
CSDM are similar to those considered to be exact while the DSM’s predictions seem to be
very poor. Between DMM and CSDM, the contours of CSDM show a greater similarity to

129

the contours of true τ11 stress.
Figure 3.10 shows the spanwise-averaged plots of τ12 SGS stress tensor. The DSM gives
a marginally improved prediction for τ12 as compared to τ11 predictions. The eddy-viscosity
models assume the stress components to be proportional to the strain-rate components. Near
the wall, the gradients are generally very high in the wall-normal direction as compared to
the gradients parallel to the wall. This makes the x−y component of strain-rate tensor more
significant and hence, a better prediction of τ12 . The predictions of MTSM are the poorest
since the turbulent viscosity computed by the model is small compared to the turbulent
viscosity computed by DSM. The DMM and CSDM again show better agreement with the
true stress in comparison to other models. Again, the predictions by CSDM are qualitatively
better than those of DMM.
A priori analysis of SGS models indicate that DSM and DMM correctly predict the
Smagorinsky coefficient, Cs , to be zero at the wall. However, at other locations Cs values
predicted by the DMM are found to be smaller and fluctuate less than those of DSM. This is
consistent with the results previously reported in the literature [46]. The dynamic coefficient
of CSDM is not of the same order as the dynamic coefficients of DSM or DMM. Its value is
generally an order of magnitude greater than DMM or DSM. However, the fluctuations in
CSDM’s coefficient are smaller than those of DMM.
Figure 3.11 shows the spatial variations of the streamwise component of SGS heat flux
qk = ρ uk T − uk T

along the wall-normal direction. For eddy-viscosity type models, qk is

modeled as
qk = −

ρνt ∂ T
P rt ∂xk

130

(3.26)

3

τ11 /ρ w u2
τ

2
1
0
−1
0

50

100

100

y

y

150

150

+

+

(a)

4

τ11 /ρ w u2
τ

3
2
1
0
−1
0

50
(b)

Figure 3.8: Spanwise-averaged values of the trace-free τ11 component of SGS stress tensor
at (a) x = 6, (b) x = 10. Legend : open triangles = True stresses, dotted black line =
Modeled-MTSM, dash-dot green line = Modeled-DSM, dashed blue line = Modeled-DMM,
solid red line = Modeled-CSDM.

131

−4

−3

x 10

1

x 10

1

5
0.8

4

0.4

2

y

3

6

0.6

y

0.6

8

0.8

0.4
0.2

1

0.2

4
2

0
0
5

10
x

0
5

15

10
x

(a)

15

0

(b)
−3

−3

x 10

1

x 10

1

2.5

0.8

6

0.6

1.5

0.6

4

0.4

1

y

2

y

0.8

0.5

0.2
0
5

0.4
0.2

0
10
x

2
0

0
5

15
(c)

10
x

15
(d)

Figure 3.9: Spanwise-averaged isocontours of trace-free τ11 SGS stress component. (a) True
stress, (b) Modeled-DSM, (c) Modeled-DMM, (d) Modeled-CSDM.

132

1

τ12 /ρ w u2
τ

0
−1
−2
−3
−4
0

50

100

100

y

150

150

+

(a)

τ12 /ρ w u2
τ

0.5

0

−0.5

−1
0

50
y

+

(b)
Figure 3.10: Spanwise-averaged values of the τ12 SGS stress component at (a) x = 10, (b)
x = 15. See Fig. 3.8 for legend.

133

For the dynamic mixed model, it is modeled as,

qk = −

ρνt ∂ T
L
+ qk
P rt ∂xk

(3.27)

L
where qk is the Leonard part of the heat flux vector computed as,

L
qk = ρ uk T − uk T

For the CSDM, the SGS heat flux is modeled as,

L
F
qk = qk + Cβ qk

(3.28)

F
qk = ρ vk θ − v k θ

Cβ =

[Gk − Hk ] Nk
Nk Nk

Gk = ρuk T −

ρuk ρT
ρ

Hk = ρuk T −

ρuk ρT
ρ

Nk = ρv k θ −

ρv k ρθ
ρ

vk = uk − uk
θ =T −T
Similar to τ11 and τ12 results shown in Figs. 3.8 and 3.10, the best predictions are given by
CSDM. Similar trends are observed for other components of SGS stresses and heat fluxes.
134

−3

1

x 10

q1

0
−1
−2
−3
0

50

100

100

y

150

150

+

(a)
−3

1

x 10

q1

0
−1
−2
−3
0

50
y

+

(b)
Figure 3.11: Spanwise-averaged values of q1 SGS heat flux component at (a) x = 10, (b)
x = 15. See Fig. 3.8 for legend.

135

To better illustrate the performance of models, we have plotted in Fig. 3.12 and Fig. 3.13
the correlation coefficients between the true and modeled SGS stress and heat flux components. The correlations are plotted against a non-dimensional parameter termed ψ which is
a measure of the filtered density gradient. The parameter ψ is similar to the parameter NS
used by Wu and Martin [13] in their numerical schlieren visualizations and is defined as

ψ = 1.0 − exp −10.0 φ − φmin / φmax − φmin

,

(3.29)

where φ = | ρ| is the magnitude of the filtered density gradient. This transformation amplifies the small density gradients in the flow. Here, we chose the data in the domain between
[5.0, 0.0, 0.0] and [20.0, 1.0, 3.0] for computing the correlation coefficients. The maximum
and the minimum values of ψ are then found within this domain and the interval between
the minimum and maximum values is divided into 10 equal sub-intervals. The correlation
coefficients for τ11 and q1 corresponding to each sub-interval are shown for DSM, DMM and
CSDM in Fig 3.12 and Fig 3.13. Evidently, the correlation coefficients are higher for DMM
and CSDM, which is consistent with other a priori results we have shown in this paper. Between DMM and CSDM, the CSDM gives higher correlations, which is again consistent with
other a priori results shown earlier. Interestingly, for all models the general trend is that
the correlation coefficient for SGS heat fluxes, and to a lesser extent for the SGS stresses,
increases with an increase in the density gradient magnitude. A possible but not definite
explanation for this trend is provided in Fig. 3.14 where the amount of backscatter or reverse
energy transfer from subgrid field to the resolved resolved field for various levels of ψ are
shown. The backscatter occurs in the the regions where the SGS dissipation, εsgs ≡ τij Sij
is positive, i.e., τij Sij > 0. Figure 3.14 shows that the amount of backscatter generally

136

1
0.8

Φτ11

0.6
0.4
0.2
0

0

0.2

0.4

0.6

0.8

1

ψ
Figure 3.12: Correlation coefficients between the true and the modeled values of τ11 SGS
stress for different values of density gradient parameter ψ. Legend : open circles = DSM,
open triangles = DMM , open squares = CSDM.
decreases with an increase in density gradient. The SGS models tested here are expected to
perform better when the transfer of energy from the unresolved to resolved scales is minimal.
This is corroborated in Fig. 3.15 where it is shown that the correlation coefficient between
the true and modeled SGS dissipation decreases appreciably for CSDM in the backscatter
regions. The same trend is also observed for DSM and DMM, although the correlation
coefficient for CSDM in the backscatter regions is still higher than other models.
Analogous to the SGS kinetic energy dissipation, εsgs , we have the SGS scalar ‘dissipation’, εφsgs , which appears in the SGS or residual scalar variance equation. For temperature,

137

1
0.8

Φq1

0.6
0.4
0.2
0

0

0.2

0.4

0.6

0.8

1

ψ
Figure 3.13: Correlation coefficients between the true and the modeled values of q1 heat flux
component for different values of density gradient parameter ψ. Legend : open circles =
DSM, open triangles = DMM , open squares = CSDM.
the SGS scalar dissipation is defined as,

∂T
εφsgs = qk
∂xk

(3.30)

The ‘backscatter’ for temperature occurs when the temperature variance values transfer from
the unresolved to resolved fields (i.e.,εφsgs > 0). Similar to that for εsgs , the backscatter for
εφsgs shown in Fig. 3.14 decreases with increasing density gradient parameter ψ, although
the decrease in εφsgs with ψ is much more pronounced than that for εsgs . This again might
explain why the correlation coefficients for the SGS heat fluxes increase with increasing ψ.
138

100
Backscatter (percentage)

sgs
φsgs

80
60
40
20
0
0

0.2

0.4

0.6

0.8

1

ψ
Figure 3.14: Variations of backscatter with the density gradient parameter ψ.

3.5.3

A Posteriori Assessment

In the a posteriori assessment of SGS models, the Reynolds-averaged statistics (represented
by

) are obtained by averaging in time over a period of 60δref /U∞ and in the homogeneous

spanwise direction. The time period used for statistical averaging was found to be sufficient
for the convergence of second-order statistics. The LES statistics are compared with those
obtained from the filtered DNS (FDNS) data. The FDNS data was obtained by filtering and
averaging DNS data on the fly. The characteristic width of the filter function employed for
LES computations and for filtering the DNS data was chosen to be twice the characteristic

139

size of the LES grid, with filtering being performed only in x and z directions, i.e., ∆ =
√
2 ∆x∆z. In LES of compressible flows, the isotropic part of SGS stress (τkk ) cannot be
neglected. Here, to model τkk in MTSM, DSM and DMM, we use the trace of modified
Leonard stress, i.e., τkk ≈ Lkk .
Figure 3.16 shows the turbulent kinetic energy profiles at x = 6. In addition to FDNS and
LES results with different SGS models, we also show the profiles of DNS and LES without
SGS model (No-Model). The effect of SGS model is quite clear in the turbulent kinetic
energy profiles. Since in all LES calculations we use the same filtered inflow velocity, the

Correlation Coefficient

1
0.8
0.6
0.4
0.2
0
0

whole domain
backscatter regions
forwardscatter regions
0.2

0.4

0.6

0.8

1

ψ
Figure 3.15: Correlation coefficients for the SGS dissipation, εsgs , as a function of the density
gradient parameter ψ as predicted by CSDM.
140

0.012
0.01

TKE

0.008
0.006
0.004
0.002
0
0

0.2

0.4

y

0.6

0.8

1

Figure 3.16: Turbulent kinetic energy profiles at x = 6. Legend : open circles = DNS,
open triangles = Filtered DNS, crosses = LES-No-Model, dotted black line = LES-MTSM,
dash-dot green line = LES-DSM, dashed blue line = LES-DMM, solid red line = LES-CSDM.
profiles for turbulent kinetic energy are the same at inlet. However at x = 6, as shown in
Fig. 3.16, the LES predictions with different SGS models significantly differ from each other.
The prediction by LES with No-Model is comparable to DNS but not FDNS. In fact, LES
with No-Model even slightly over-predicts the peak DNS value. On coarser grids like the ones
used in our LES, without an SGS model there is an energy build-up at smaller scales which
cannot be dissipated properly and leads to an over-prediction of turbulent kinetic energy
and perhaps even a failure of the simulations at higher Reynolds numbers. This clearly

141

demonstrates the necessity of SGS models. In Fig. 3.17, the 1D-energy spectra obtained by
LES with No-Model, MTSM and CSDM is compared with that of FDNS. The spectra for LES
with No-Model and MTSM show an over-prediction at all wavenumbers, which is consistent
with the over-prediction of turbulent kinetic energy in Fig. 3.16. At higher wavenumbers, the
MTSM provides a slightly better comparison than No-Model. The CSDM, on the other hand,
shows a very good agreement with the FDNS spectra, particularly at higher wavenumbers
where it correctly predicts the decay of energy. Comparing the spectra for FDNS and LES
with No-model suggests that our numerical scheme is not excessively dissipative.

−1

10

−2

10

−3

E(k)

10

−4

10

FDNS
LES−NoModel
LES−MTSM
LES−CSDM

−5

10

−6

10

−1

0

10

10
k

1

10

Figure 3.17: One-dimensional energy spectra in spanwise direction at x = 6 and y ≈ 0.07.

142

4

x 10

−3

3

2
Cf

−3

4

x 10

1
2
0
0
0

−1
5

10
x

5
15

Figure 3.18: Evolution of mean skin-friction coefficient, Cf . See Fig. 3.16 for legend.
Figure 3.18 shows the skin-friction coefficient, Cf obtained from FDNS and LES data.
The skin-friction coefficient is defined as

Cf =

1 ρw
2 ρ∞

uτ 2
U∞

(3.31)

Figure 3.18 shows that the shock is strong enough for creating a small recirculation region
in the flow. The predictions by various models can be compared in three different flow
regions: upstream of the separation point, in the interaction region and downstream of the
reattachment point. The Cf predictions by the new model (CSDM) are very good upstream
of the separation point. For example, at x = 6 the difference between FDNS and CSDM
143

results is only +1.5% (positive sign indicates over-prediction by the model). All other SGS
models under-predict Cf with DSM giving the poorest prediction. At x = 6, the errors
of No-Model, MTSM, DSM and DMM are −8%, −10%, −17% and −7%, respectively. It
should be noted that the Cf values of FDNS and LES are identical at the inflow boundary.
The inset in Fig. 3.18 shows the evolution of Cf downstream of the inflow boundary. The
under-prediction of Cf by MTSM, DSM and DMM downstream of inflow boundary indicates
the inability of these models to accurately model the SGS stresses. On the other hand, the
new model properly accounts for the SGS effects as LES results with CSDM are very close
to those of FDNS.
Table 3.3: Separation region sizes
FDNS
Lsep

0.54δref

No-Model MTSM
0.58δref

0.6δref

DSM

DMM

CSDM

0.94δref

0.62δref

0.36δref

In the interaction region, none of the SGS models are accurate and can completely capture
the FDNS results. Table 3.3 shows the lengths of separation bubble obtained from FDNS
and LES data. While LES with No-Model, MTSM, DSM and DMM over-predict the size
of separation region with varying degrees, the FDNS results are under-predicted by CSDM.
LES with No-Model, MTSM and DMM predict the onset of separation well but compared
to FDNS, the mean separation point occurs earlier when DSM is employed or delayed when
CSDM is used. Downstream of the reattachment point, the predictions by MTSM and CSDM
are in good agreement with FDNS, although it should be noted that MTSM under-predicted
Cf upstream of the mean separation point. Similar to results shown upstream of separation,
LES with No-Model, DSM or DMM under-predict Cf .
Figure 3.19 shows the downstream variations of mean wall pressure normalized by the

144

1.6
1.5

Pw

1.4
1.3
1.2
1.1
1
5

10
x

15

Figure 3.19: Evolution of mean wall-pressure, Pw . See Fig. 3.16 for legend.
upstream wall pressure. It can be observed in this figure that the wall pressure begins to
rise at x = 7.2 even though the nominal impingement location is at x = 10 which is due to
the upstream propagation of the disturbances through the subsonic region of the boundary
layer. In general, all models predict the FDNS data satisfactorily. The CSDM predicts the
pressure rise at a slightly later location in comparison to other models and also in comparison
to FDNS.
Figure 3.20 compares the Van Driest transformed mean streamwise velocity profiles in
inner scaling at different streamwise locations. At x = 6 (not shown), the FDNS profile
exhibits a log-law profile (Eq. 1.14) between 30 < y + < 100 . At this location, the boundary

145

layer is unaffected by the shock and exhibits the profile of a canonical turbulent boundary
layer. Also, at this location LES results obtained by CSDM are in excellent agreement
with those of FDNS. However, the results obtained by LES with No-Model, MTSM and
DMM are different than those of FDNS. DSM seems to be the least accurate model. Since
the velocity profiles are normalized by friction velocity, uτ , the variations in the slope of
logarithmic function for various models is directly related to the Cf predictions. At x = 10
(Fig. 3.20(a)), just downstream of the mean reattachment point, the velocity profiles do not
exhibit any logarithmic behavior although the universal near-wall U + = y + behavior is still
obeyed. Once again, CSDM shows the best agreement with FDNS profile in comparison
to other models. At x = 15 (Fig. 3.20(b)), the boundary layer starts to relax back to
a canonical boundary layer profile although it has not fully recovered from the effect of
adverse pressure gradient induced by the shock. At this location, the profile exhibits a
logarithmic-like behavior albeit with different constants and a stronger wake region. This is
very similar to a separated incompressible boundary layer [31]. All models show reasonably
good agreement with FDNS at this location. Overall, the predictions by CSDM are more
accurate than other tested models for mean streamwise velocity at all three streamwise
locations.
Figure 3.21 compares the R11 Reynolds stress component obtained by DNS and LES
at different locations. At x = 6 (Fig. 3.21(a)), LES with any model or with No-Model
over-predicts R11 with No-Model giving the largest values. However, the best agreement
with FDNS profile is obtained by the CSDM. At x = 10 (Fig. 3.21(b)), the peak value of
R11 increases due to the shock and this trend is well captured by LES. Again, the best
agreement with FDNS profile is obtained by CSDM. Interestingly, the LES results obtained
by No-Model and those with eddy-viscosity models are very similar at this location. Similar
146

50

+
UV D

40
30
20
10
0 0
10

2

10
y+
(a)

25

+
UV D

20
15
10
5
0 0
10

2

10
y+
(b)

Figure 3.20: Van-Driest transformed mean streamwise velocity profiles in inner scaling. (a)
x = 10, (b) x = 15. Black dashed line represents law of the wall. See Fig. 3.16 for legend.

147

trends are observed at x = 15 (not shown), although the predicted LES trends do not quite
match with those of FDNS.
For the R22 component of Reynolds stress, Fig. 3.22 again shows a general trend of overpredicting the FDNS stresses by LES employing No-model, MTSM, DMM or even CSDM at
all selected locations, even though the best prediction is by CSDM. At x = 15, the predictions
by MTSM and No-Model are similar. Similar behaviors are observed for the R33 and R12
component of Reynolds stress (not shown) with the best predictions given by CSDM at all
three selected locations.

3.6

3D effects of lateral walls on SBLI

Up until now, majority of the SBLI simulations have been carried out with an assumption
of homogeneity in the spanwise directions. Several experimental results indicate that the
lateral walls of wind tunnel influence the physics of SBLI. A recent numerical study by
Garnier [66] has shown that the corner separation generated due to the presence of lateral
walls reduces the effective span of the wind tunnel, strengthens the interaction and the
strongest fluctuations are found at the corners. This problem has been gaining a steady
interest because of its practical relevance, for example, at the inlet of a scramjet.
To examine the effect of lateral walls on SBLI, some preliminary calculations have been
performed. A unique approach was adopted to overcome the specific challenges to this
problem. The first challenge was to generate a realistic turbulent state at the interaction region. The presence of lateral walls meant standard inflow turbulence generation methods for
boundary layers could not be used. To generate a realistic turbulent state in the interaction
region, an auxiliary LES simulation was first carried out. In the auxiliary LES, a laminar

148

0.02

u1 u1

0.015
0.01
0.005
0
0

0.5
y

1

0.5
y

1

(a)

0.03

u1 u1

0.02

0.01

0
0
(b)

Figure 3.21: Comparisons of R11 Reynolds stress component obtained by DNS and LES. (a)
x = 6, (b) x = 10. See Fig. 3.16 for legend.

149

−3

8

x 10

u2 u2

6
4
2
0
0

0.5
y

1

0.5
y

1

(a)
−3

5

x 10

u2 u2

4
3
2
1
0
0
(b)
Figure 3.22: Comparisons of R22 Reynolds stress component obtained by DNS and LES. (a)
x = 10, (b) x = 15. See Fig. 3.16 for legend.

150

boundary layer profile was specified at the inflow and the boundary layer was perturbed by
a blowing-suction forcing at the wall [67]. This is achieved by computing the wall-normal
velocity component at the wall in a region between a ≤ x ≤ b using the following relation:
v(x, z, t) = U∞ Af (x)g(z)h(t)

(3.32)

where
√
f (x) = 4sinθ(1 − cosθ)/ 27
10

g(z) =

Zl sin(2πl(z/Lz + φl ))
l=0

10

Zl = 1, Zl = 1.25Zl+1
l=0
5

h(t) =

Tm sin(2πm(βt + φm ))
m=0

10

Tm = 1, Tm = 1.25Tm+1
l=0

The amplitude for the perturbations was chosen equal to A = 0.1 and the forcing frequency,
β = 0.1. The domain size for the auxiliary LES was, x = 4.5in, y = 1.3in and z = 0.5in.
The blowing and suction strip was placed between 0.5in ≤ x ≤ 1in. The freestream Mach
number was 2.75. The Smagorinsky SGS model was used for the computations. Once the
flow reached a realistic turbulent state downstream, a time history of primitive variables was
stored from a selected downstream plane. This ‘inflow box’ was then used periodically for the
main SBLI simulation. The second challenge was the generation of an impinging shock. A
compression ramp with an inviscid/slip wall boundary condition at the freestream was used

151

to generate the impinging shock. The wedge angle was equal to 6o . In Fig. 3.23 we show the
y−z slices for instantaneous temperature contours. The nominal shock impingement location
was at x = 3 and as can be seen in the figure, the mixing is enhanced in the interaction
region and also downstream of the interaction region and there is a very significant influence
of the boundary layers developing on the lateral walls. It should be noted that the spanwise
width for these computations was only 0.5in which is very narrow and hence, 3D effects are
expected to be very strong for narrow spanwise domain widths. A further detailed analysis
is required to make a direct comparison between the cases with and without lateral walls.

Figure 3.23: Instantaneous slices of temperature contours at different streamwise locations
for 3DSBLI.

152

3.7

Conclusions

A priori and a posteriori assessments of SGS models are made for the problem of 34o
incident shock wave interacting with a Mach 2 turbulent boundary layer via direct numerical
simulation (DNS) and large-eddy simulation (LES) data. In both DNS and LES, the Euler
terms of the governing equations are solved by a new seventh-order monotonicity preserving
scheme while the diffusion terms are solved by a sixth-order compact finite difference scheme.
The SGS stress and scalar flux models tested are the mixed time scale model (MTSM),
the dynamic Smagorinsky model (DSM), the dynamic mixed model (DSM) and the new
compressible serial decomposition model (CSDM).
A priori testing of the models indicate the superiority of similarity-type models like
DMM and CSDM over eddy-viscosity type closures like DSM. Between DMM and CSDM,
the predictions by CSDM are more accurate. The decomposition of SGS fluxes into Leonard,
Cross and Reynolds components shows the importance of taking into account the Leonard
component since it is a major contributor toward the total SGS stress or heat flux. For
larger filter widths, the contribution of Reynolds component also becomes important since it
represents the direct effect of unresolved scales. Hence, a proper modeling of this term is also
necessary for accurate LES computations. The compressible serial decomposition model uses
a similarity concept to model this term and it is more accurate than the eddy-viscosity type
model employed by the dynamic mixed model. This is clearly reflected in the correlation
coefficients computed between the true and modeled SGS stresses and heat fluxes. The
correlation coefficients computed for the three dynamic models are found to increase with
increasing the density gradient magnitude. A possible but by no means definite explanation
for this trend can be the amount of backscatter which is shown to decrease with increasing

153

density gradient magnitudes.
A posteriori comparisons made for mean velocity profiles and filtered Reynolds stresses
indicate the significance of SGS model effect and are generally consistent with our a priori
testing of SGS models. For a posteriori assessment of the models, special attention was made
to keep the inflow boundary condition consistent in DNS and LES. Also, the characteristic
filter width for a posteriori assessment was chosen to be ∆ = 2∆LES since, (a) it presents
a greater challenge for the SGS models to correctly capture SGS terms behavior, and (b)
it gives a better indication of model performance for higher Reynolds number flows where
the grid requirements are more demanding. The skin-friction coefficient is shown to be best
predicted by CSDM in comparison to other SGS models for the incoming boundary layer. In
the interaction region, however, the predictions by CSDM is not very accurate, which might
be due to CSDM not being sufficiently dissipative. The performance of other models is also
not adequate. Overall, the LES results obtained by CSDM are found to be comparable with
the FDNS values, despite the fact that the SGS terms B2 , B5 , B6 , and B7 are neglected
in our LES calculations as is done in most LES studies of compressible flows. These terms
represent the interactions between the resolved and unresolved fields and non-linear effects
in the filtered energy equation. It is possible that these terms could be significant in shockturbulence interactions at higher shock intensities and/or flow Mach numbers, something
that has to be tested in future studies. In general, the LES results obtained by CSDM
appear to be very promising. However, certain features of the model need to be improved.
Also, its predictions for other compressible turbulent flows needs to be further examined.

154

BIBLIOGRAPHY

155

BIBLIOGRAPHY
[1] Pirozzoli, S. and Bernardini, M., “Turbulence in Supersonic Boundary Layers at Moderate Reynolds Number,” Journal of Fluid Mechanics, Vol. 688, 2011, pp. 120–168.
[2] Jeong, J. and Hussain, F., “On the Identification of a Vortex,” Journal of Fluid Mechanics, Vol. 285, 1995, pp. 69–94.
[3] Green, J. E., “Reflexion of an Oblique Shock Wave by a Turbulent Boundary Layer,”
Journal of Fluid Mechanics, Vol. 40, No. 1, 1970, pp. 81–95.
[4] Dupont, P., Haddad, C., and Debieve, J. F., “Space and Time Organization in a
Shock-Induced Separated Boundary Layer,” Journal of Fluid Mechanics, Vol. 559,
2006, pp. 255–277.
[5] Souverein, L. J., Dupont, P., Debieve, J. F., Dussauge, J. P., van Oudheusden, B. W.,
and Scarano, F., “Effect of Interaction Strength on Unsteadiness in Turbulent ShockWave-Induced Separations,” AIAA Journal , Vol. 48, No. 7, 2010, pp. 1480–1493.
[6] Touber, E. and Sandham, N. D., “Comparison of Three Large-Eddy Simulations of
Shock-Induced Turbulent Separation Bubbles,” Shock Waves, Vol. 19, 2009, pp. 469–
478.
[7] Morgan, B., Kawai, S., and Lele, S. K., “A Parametric Investigation of Oblique Shockwave/Turbulent Boundary Layer Interaction Using LES,” AIAA Paper No. 2011-3430 ,
June 2011.
[8] Pirozzoli, S. and Bernardini, M., “Direct Numerical Simulation Database for Impinging
Shock Wave/Turbulent Boundary-Layer Interaction,” AIAA Journal , Vol. 49, No. 6,
2011, pp. 1307–1312.
[9] Pirozzoli, S., “Generalized Conservative Approximations of Split Convective Derivative
Operators,” Journal of Computational Physics, Vol. 2010, 2010, pp. 7180–7190.
[10] Suresh, A. and Huynh, H. T., “Accurate Monotonicity-Preserving Schemes with
Runge-Kutta Time Stepping,” Journal of Computational Physics, Vol. 136, No. 1,
1997, pp. 83–99.

156

[11] Li, Z. and Jaberi, F. A., “High-Order Finite Difference Schemes for Numerical Simulations of Supersonic Turbulent Flows,” International Journal for Numerical Methods
in Fluids, Vol. 68, No. 6, 2012, pp. 740–766.
[12] Pirozzoli, S. and Grasso, F., “Direct Numerical Simulation of Impinging Shock
Wave/Turbulent Boundary Layer Interaction at M = 2.25,” Physics of Fluids, Vol. 18,
No. 6, 2006, pp. 065113.
[13] Wu, M. and Martin, M. P., “Direct Numerical Simulation of Supersonic Turbulent
Boundary Layer over a Compression Ramp,” AIAA Journal , Vol. 45, No. 4, 2007,
pp. 879–889.
[14] Priebe, S., Wu, M., and Martin, M. P., “Direct Numerical Simulation of a ReflectedShock-Wave/Turbulent-Boundary-Layer Interaction,” AIAA Journal , Vol. 47, No. 5,
2009, pp. 1173–1185.
[15] Priebe, S. and Martin, M. P., “Low-Frequency Unsteadiness in Shock Wave-Turbulent
Boundary Layer Interaction,” Journal of Fluid Mechanics, Vol. 699, 2012, pp. 1–49.
[16] Jammalamadaka, A., Li, Z., and Jaberi, F. A., “Subgrid Scale Models for Large-Eddy
Simulations of Shock/Turbulent Boundary Layer Interactions,” AIAA Journal , Vol. 51,
No. 5, 2013, pp. 1174–1188.
[17] Larsson, J. and Lele, S. K., “Direct Numerical Simulation of Canonical
Shock/Turbulence Interaction,” Physics of Fluids, Vol. 21, 2009, pp. 126101.
[18] Bernardini, M. and Pirozzoli, S., “The Wall Pressure Signature of Transonic
Shock/Boundary Layer Interaction,” Journal of Fluid Mechanics, Vol. 671, March
2011, pp. 288–312.
[19] Kennedy, C. A. and Gruber, A., “Reduced Aliasing Formulations of the Convective
Terms Within the Navier-Stokes Equations for a Compressible Fluid,” Journal of Computational Physics, Vol. 227, 2008, pp. 1676–1700.
[20] Lele, S. K., “Compact Finite Difference Schemes with Spectral-Like Resolution,” Journal of Computational Physics, Vol. 103, No. 1, 1992, pp. 16–42.
[21] Pirozzoli, S., “Numerical Methods for High-Speed Flows,” Annual Review of Fluid
Mechanics, Vol. 43, 2011, pp. 163–194.

157

[22] Sandham, N. D., Li, Q., and Yee, H. C., “Entropy Splitting for High-Order Simulation
of Compressible Turbulence,” Journal of Computational Physics, Vol. 178, No. 2, 2002,
pp. 307–322.
[23] Adams, N. A., “Direct Simulation of the Turbulent Boundary Layer Along a Compression Ramp at M = 3 and Reθ = 1685,” Journal of Fluid Mechanics, Vol. 420, 2000,
pp. 47–83.
[24] Adams, N. A., “Direct Numerical Simulation of Turbulent Compression Ramp Flow,”
Theoretical Computational Fluid Dynamics, Vol. 12, No. 2, 1998, pp. 109–129.
[25] Pirozzoli, S., Bernardini, M., and Grasso, F., “Direct Numerical Simulation of Transonic Shock/Boundary Layer Interaction Under Conditions of Incipient Separation,”
Journal of Fluid Mechanics, Vol. 657, 2010, pp. 361–393.
[26] Thompson, K. W., “Time Dependent Boundary Conditions for Hyperbolic Systems,
II,” Journal of Computational Physics, Vol. 89, No. 2, 1990, pp. 439–461.
[27] Poinsot, T. J. and Lele, S. K., “Boundary Conditions for Direct Simulations of Compressible Viscous Flows,” Journal of Computational Physics, Vol. 101, No. 1, 1992,
pp. 104–129.
[28] Ganapathisubramani, B., Clemens, N. T., and Dolling, D. S., “Effects of Upstream
Boundary Layer on the Unsteadiness of Shock-Induced Separation,” Journal of Fluid
Mechanics, Vol. 585, 2007, pp. 369–394.
[29] Huang, P. G., Coleman, G. N., and Bradshaw, P., “Compressible Turbulent Channel
Flows: DNS Results and Modelling,” Journal of Fluid Mechanics, Vol. 305, 1995,
pp. 185–218.
[30] Simpson, R. L., “Turbulent Boundary Layer Separation,” Annual Review of Fluid
Mechanics, Vol. 21, 1989, pp. 205–234.
[31] Na, Y. and Moin, P., “Direct Numerical Simulation of a Separated Turbulent Boundary
Layer,” Journal of Fluid Mechanics, Vol. 374, No. 1, 1998, pp. 379–405.
[32] Lele, S. K., “Compressibility Effects on Turbulence,” Annual Review of Fluid Mechanics, Vol. 26, 1994, pp. 211–254.
[33] Garnier, E., Sagaut, P., and Deville, M., “Large Eddy Simulation of Shock/Boundary
Layer Interaction,” AIAA Journal , Vol. 40, No. 10, 2002, pp. 1935–1944.
158

[34] Touber, E. and Sandham, N. D., “Large-Eddy Simulation of Low-Frequency Unsteadiness in a Turbulent Shock-Induced Separation Bubble,” Theoretical Computational
Fluid Dynamics, Vol. 23, No. 2, 2009, pp. 79–107.
[35] Pirozzoli, S., Beer, A., Bernardini, M., and Grasso, F., “Computational Analysis of
Impinging Shock-Wave Boundary Layer Interaction Under Conditions of Incipient Separation,” Shock Waves, Vol. 19, No. 6, 2009, pp. 487–497.
[36] Loginov, M. S., Adams, N. A., and Zheltovdov, A. A., “Large Eddy Simulation
of Shock-Wave/Turbulent-Boundary-Layer Interaction,” Journal of Fluid Mechanics,
Vol. 565, 2006, pp. 135–169.
[37] Stolz, S., Adams, N. A., and Kleiser, L., “The Approximate Deconvolution Model
for Large-Eddy Simulations of Compressible Flows and its Application to ShockTurbulent-Boundary-Layer Interaction,” Physics of Fluids, Vol. 13, No. 10, 2001,
pp. 2985–3001.
[38] Vreman, B., Geurts, B., and Kuerten, H., “A Priori Tests of Large Eddy Simulation
of Compressible Mixing Layer,” Journal of Engineering Mathematics, Vol. 29, No. 5,
1995, pp. 299–327.
[39] Martin, M. P., Piomelli, U., and Candler, G. V., “Subgrid-Scale Models for Compressible Large-Eddy Simulations,” Theoretical Computational Fluid Dynamics, Vol. 13,
2000, pp. 361–376.
[40] Jammalamadaka, A., Li, Z., and Jaberi, F. A., “Numerical Simulations of Shock Wave
Interactions with a Supersonic Turbulent Boundary Layer,” submitted to AIAA Journal.
[41] Ducros, F., Laporte, F., Souleres, T., Guinot, V., Moinat, P., and Caruelle, B., “HighOrder Fluxes for Conservative Skew-Symmetric-like Schemes in Structured Meshes:
Application to Compressible Flows,” Journal of Computational Physics, Vol. 161,
No. 1, 2000, pp. 114–139.
[42] Roe, P. L., “Approximate Riemann Solvers, Parameter Vectors, and Difference
Schemes,” Journal of Computational Physics, Vol. 135, 1997, pp. 250–258.
[43] Huynh, H. T., “Accurate Monotone Cubic Interpolation,” SIAM Journal on Numerical
Analysis, Vol. 30, No. 1, 1993, pp. 57–100.
[44] Liu, X. D., Osher, S., and Chan, T., “Weighted Essentially Non-Oscillatory Schemes,”
Journal of Computational Physics, Vol. 115, 1994, pp. 200–212.
159

[45] Germano, M., “A Proposal for a Redefinition of the Turbulent Stresses in the Filtered
Navier-Stokes Equations,” Physics of Fluids, Vol. 29, No. 7, 1986, pp. 2323–2324.
[46] Salvetti, M. V. and Banerjee, S., “A Priori Tests of a New Dynamic Subgrid-Scale
Model for Finite-Difference Large-Eddy Simulations,” Physics of Fluids, Vol. 7, No. 11,
1995, pp. 2831–2847.
[47] Piomelli, U., “Large-Eddy Simulation: Achievements and Challenges,” Progress in
Aerospace Sciences, Vol. 35, No. 4, 1999, pp. 335–363.
[48] Vreman, B., Geurts, B., and Kuerten, H., “Subgrid-Modelling in LES of Compressible
Flow,” Applied Scientific Research, Vol. 54, 1995, pp. 191–203.
[49] Knight, D., Zhou, G., Okong’o, N., and Shukla, V., “Compressible Large Eddy Simulation using Unstructured Grids,” AIAA Paper No. 98-0535 , 1998.
[50] Stolz, S. and Adams, N. A., “Large-Eddy Simulation of High-Reynolds Number Supersonic Boundary Layers using the Approximate Deconvolution Model and a Rescaling
and Recycling Technique,” Physics of Fluids, Vol. 15, No. 8, 2003, pp. 2398–2412.
[51] Dolling, D. S., “Fifty Years of Shock-Wave/Boundary-Layer Interaction Research:
What Next?” AIAA Journal , Vol. 39, No. 8, 2001, pp. 1517–1531.
[52] Pope, S. B., Turbulent Flows, Cambridge University Press, 3rd ed., 2003.
[53] Vreman, B., Direct and Large-Eddy Simulation of the Temporal Turbulent Compressible Mixing Layer , Ph.D. thesis, Department of Applied Mathematics, University of
Twente, Enschede, Netherlands, 1995.
[54] Moin, P., Squires, K., Cabot, W., and Lele, S., “A Dynamic Subgrid-Scale Model for
Compressible Turbulence and Scalar Transport,” Physics of Fluids, Vol. 3, No. 11,
1991, pp. 2746–2757.
[55] Dupont, P., Piponniau, S., Sidorenko, A., and Debieve, J. F., “Investigation by Particle
Image Velocimetry Measurements of Oblique Shock Reflection with Separation,” AIAA
Journal , Vol. 46, No. 6, 2008, pp. 1365–1370.
[56] Klein, M., Sadiki, A., and Janicka, J., “A Digital Filter Based Generation of Inflow
Data for Spatially Developing Direct Numerical or Large Eddy Simulations,” Journal
of Computational Physics, Vol. 186, 2003, pp. 652–665.

160

[57] Morgan, B., Kawai, S., and Lele, S. K., “Large-Eddy Simulation of an Oblique Shock
Impingement on a Turbulent Boundary Layer,” AIAA Paper No. 2010-4467 , July
2010.
[58] Morgan, B., Larsson, J., Kawai, S., and Lele, S. K., “Improving Low-Frequency
Characteristics of Recycling/Rescaling Inflow Turbulence Generation,” AIAA Journal , Vol. 49, No. 3, 2011, pp. 582–597.
[59] Yoshizawa, A., “Statistical Theory for Compressible Turbulent Shear Flows, with the
Application to Subgrid Modeling,” Physics of Fluids, Vol. 29, No. 7, 1986, pp. 2152–
2164.
[60] Spyropoulous, E. and Blaisdell, G. A., “Large-Eddy Simulation of a Spatially Evolving
Turbulent Boundary-Layer Flow,” AIAA Journal , Vol. 36, No. 11, 1998, pp. 1983–
1990.
[61] Speziale, C. G., Erlebacher, G., Zang, T. A., and Hussaini, M. Y., “The SubgridScale Modeling of Compressible Turbulence,” Physics of Fluids, Vol. 31, No. 4, 1988,
pp. 940–942.
[62] Lilly, D. K., “A Proposed Modification of the Germano Subgrid-Scale Closure
Method,” Physics of Fluids, Vol. 4, No. 3, 1992, pp. 633–635.
[63] Jaberi, F. A. and Colucci, P. J., “Large Eddy Simulation of Heat and Mass Transport
in Turbulent Flows. Part 1: Velocity Field,” International Journal for Heat and Mass
Transfer , Vol. 46, No. 10, 2003, pp. 1811–1825.
[64] Li, Z. and Jaberi, F. A., “Large-Eddy Simulations of Shock-Turbulence Interactions
and Compressible Scalar Mixing,” (to be published).
[65] Rizzetta, D. P., Visbal, M. R., and Gaitonde, D. V., “Large-Eddy Simulation of Supersonic Compression-Ramp Flow by High-Order Method,” AIAA Journal , Vol. 39,
No. 12, 2001, pp. 2283–2292.
[66] Garnier, E., “Stimulated Detached Eddy Simulation of Three-Dimensional
Shock/Boundary Layer Interaction,” Shock Waves, Vol. 19, 2009, pp. 479–486.
[67] Pirozzoli, S., Grasso, F., and Gatski, T. B., “Direct Numerical Simulation and Analysis
of a Spatially Evolving Turbulent Boundary Layer at M=2.25,” Physics of Fluids,
Vol. 16, No. 7, 2004, pp. 530–545.

161

[68] Pirozzoli, S., Bernardini, M., and Grasso, F., “On the Dynamical Relevance of Coherent Vortical Structures in the Turbulent Boundary Layers,” Journal of Fluid Mechanics, Vol. 648, 2010, pp. 325–349.
[69] Jammalamadaka, A. and Jaberi, F. A., “A Priori Estimates of Subgrid-Scale Terms
for Shock-Boundary Layer Interactions,” submitted to Physics of Fluids.
[70] Bernardini, M. and Pirozzoli, S., “Wall Pressure Fluctuations Beneath Supersonic
Turbulent Boundary Layers,” Physics of Fluids, Vol. 23, 2011, pp. 085102.
[71] Bernardini, M. and Pirozzoli, S., “Inner/Outer Layer Interactions in Turbulent Boundary Layers: A Refined Measure for the Large-Scale Amplitude Modulation Mechanism,” Physics of Fluids, Vol. 23, 2011, pp. 061701.
[72] Pirozzoli, S., “Conservative Hybrid Compact-WENO Schemes for Shock-Turbulence
Interaction,” Journal of Computational Physics, Vol. 178, No. 1, 2002, pp. 81–117.
[73] Touber, E. and Sandham, N. D., “Low-Order Stochastic Modelling of Low-Frequency
Motions in Reflected Shock-Wave/Boundary Layer Interactions,” Journal of Fluid Mechanics, Vol. 671, 2011, pp. 417–465.
[74] Touber, E. and Sandham, N. D., “Large-Eddy Simulations of an Oblique Shock Impinging on a Turbulent Boundary Layer: Effect of Spanwise Confinement on the LowFrequency Oscillations,” Turbulence and Interactions, edited by M. Deville, T.-H. Le,
and P. Sagaut, Vol. 110 of Notes on Numerical Fluid Mechanics and Multidisciplinary
Design, Springer Berlin/Heidelberg, 2010, pp. 361–367.
[75] Johnsen, E., Larsson, J., Bhagatwala, A. V., Cabot, W. H., Moin, P., Olson, B. J.,
Rawat, P. S. Shankar, S. K., Sjogreen, B., Yee, H. C., Zhong, X., and Lele, S. K.,
“Assessment of High-Resolution Methods for Numerical Simulations of Compressible
Turbulence with Shock Waves,” Journal of Computational Physics, Vol. 229, 2010,
pp. 1213–1237.
[76] Li, N., Balaras, E., and Piomelli, U., “Inflow Conditions for Large-Eddy Simulations
of Mixing Layers,” Physics of Fluids, Vol. 12, No. 4, 2000, pp. 935–938.
[77] Jarrin, N., Benhamadouche, S., Laurence, D., and Prosser, R., “A Synthetic-EddyMethod for Generating Inflow Conditions for Large-Eddy Simulations,” International
Journal of Heat and Fluid Flow , Vol. 27, 2006, pp. 585–593.
[78] di Mare, L., Klein, M., Jones, W. P., and Janicka, J., “Synthetic Turbulence Inflow
Conditions for Large-Eddy Simulation,” Physics of Fluids, Vol. 18, 2006, pp. 025107.
162

[79] Xu, S. and Martin, M. P., “Assessment of Inflow Boundary Conditions for Compressible
Turbulent Boundary Layers,” Physics of Fluids, Vol. 16, No. 7, 2004, pp. 2623–2639.
[80] Katzer, E., “On the Lengthscales of Laminar Shock/Boundary-Layer Interaction,”
Journal of Fluid Mechanics, Vol. 206, 1989, pp. 477–496.
[81] Piponniau, S., Dussauge, J. P., Debieve, J. F., and Dupont, P., “A Simple Model for
Low-Frequency Unsteadiness in Shock-Induced Separation,” Journal of Fluid Mechanics, Vol. 629, 2009, pp. 87–108.
[82] Park, N., Yoo, J. Y., and Choi, H., “Discretization Errors in Large Eddy Simulations:
On the Suitablity of Centered and Upwind-Biased Compact Difference Schemes,” Journal of Computational Physics, Vol. 198, 2004, pp. 580–616.
[83] Maeder, T., Adams, N. A., and Kleiser, L., “Direct Simulation of Turbulent Supersonic
Boundary Layers by an Extended Temporal Approach,” Journal of Fluid Mechanics,
Vol. 429, 2001, pp. 187–216.
[84] Gatski, T. B. and Erlebacher, G., “Numerical Simulation of a Spatially Evolving Supersonic Turbulent Boundary Layer,” Tech. rep., NASA, 2002, NASA Tech. Memo
2002-211934.
[85] Guarini, S. E., Moser, R. D., Shariff, K., and Wray, A., “Direct Numerical Simulation
of a Supersonic Turbulent Boundary Layer at Mach 2.5,” Journal of Fluid Mechanics,
Vol. 414, 2000, pp. 1–33.
[86] Spina, E. F., Smits, A. J., and Robinson, S. K., “The Physics of Supersonic Turbulent
Boundary Layers,” Annual Review of Fluid Mechanics, Vol. 26, 1994, pp. 287–319.
[87] Grasso, F. and Pirozzoli, S., “Shock-Wave-Vortex Interactions: Shock and Vortex
Deformations, and Sound Production,” Theoretical Computational Fluid Dynamics,
Vol. 13, 2000, pp. 421–456.
[88] Thomas, F. O., Putnam, C. M., and Chu, H. C., “On the Mechanism of Unsteady
Shock Oscillation in Shock Wave/Turbulent Boundary layer Interactions,” Experiments
in Fluids, Vol. 18, 1994, pp. 69–81.
[89] Andreopoulous, J. and Muck, K. C., “Some New Aspects of the Shock-Wave/Boundary
Layer Interaction in Compression Ramp Flows,” Journal of Fluid Mechanics, Vol. 180,
1987, pp. 405–428.

163

[90] Lee, S., Lele, S. K., and Moin, P., “Direct Numerical Simulation of Isotropic Turbulence
Interacting With a Weak Shock Wave,” Journal of Fluid Mechanics, Vol. 251, 1993,
pp. 533–562.
[91] Dussauge, J. P. and Piponniau, “Shock/Boundary-Layer Interactions: Possible Sources
of Unsteadiness,” Journal of Fluids and Structures, Vol. 24, 2008, pp. 1166–1175.
[92] Dussauge, J. P., Dupont, P., and Debieve, J. F., “Unsteadiness in Shock Wave Boundary Layer Interactions,” Aerospace Science and Technology, Vol. 10, 2006, pp. 85–91.
[93] Wu, M. and Martin, M. P., “Analysis of Shock Motion in Shockwave and Turbulent
Boundary Layer Interaction Using Direct Numerical Simulation Data,” Journal of
Fluid Mechanics, Vol. 594, 2008, pp. 71–83.
[94] Urbin, G. and Knight, D., “Large-Eddy Simulation of a Supersonic Boundary Layer
Using an Unstructured Grid,” AIAA Journal , Vol. 39, No. 7, 2001, pp. 1288–1295.
[95] Suman, M. and Mahesh, K., “DNS of Unsteady Shock Boundary Layer Interaction,”
49th AIAA Aerospace Sciences Meeting and Exhibit, Orlando, Florida, January 2011.
[96] Cook, A. W. and Riley, J. J., “Direct Numerical Simulation of a Turbulent Reactive
Plume on a Parallel Computer,” Journal of Computational Physics, Vol. 129, 1996,
pp. 263–283.
[97] Visbal, M. R. and Gaitonde, D. V., “On the Use of Higher-Order Finite-Difference
Schemes on Curvilinear and Deforming Meshes,” Journal of Computational Physics,
Vol. 181, 2002, pp. 155–185.
[98] Lee, S., Lele, S. K., and Moin, P., “Interaction of Isotropic Turbulence with Shock
Waves: Effect of Shock Strength,” Journal of Fluid Mechanics, Vol. 340, 1997, pp. 225–
247.
[99] Jammalamadaka, A., Li, Z., and Jaberi, F. A., “Numerical Study of Shock-Boundary
Layer Interactions for Various Shock and Flow Conditions,” 49th AIAA Aerospace
Sciences Meeting and Exhibit, Orlando, Florida, January 2011.
[100] Jammalamadaka, A., Li, Z., and Jaberi, F. A., “Large-Eddy Simulation of Turbulent
Boundary Layer Interaction with an Oblique Shock Wave,” 48th AIAA Aerospace
Sciences Meeting and Exhibit, Orlando, Florida, January 2010.
[101] Smits, A. J. and Dussauge, J. P., Turbulent Shear Layers in Supersonic Flows,
Springer, 2nd ed., 2006.
164

[102] Thompson, K. W., “Time Dependent Boundary Conditions for Hyperbolic Systems,”
Journal of Computational Physics, Vol. 68, No. 1, 1987, pp. 1–24.
[103] Sagaut, P., Garnier, E., Tromeur, E., Larcheveque, L., and Labourasse, E., “Turbulent
Inflow Conditions for Large-Eddy Simulation of Compressible Wall-Bounded Flows,”
AIAA Journal , Vol. 42, No. 3, 2004, pp. 469–477.
[104] Spalart, P. R., Strelets, M., and Travin, A., “Direct Numerical Simulation of LargeEddy-Break-Up Devices in a Boundary Layer,” International Journal of Heat and Fluid
Flow , Vol. 27, 2006, pp. 902–910.
[105] Martin, M. P., “Direct Numerical Simulation of Hypersonic Turbulent Boundary Layers. Part 1. Initialization and Comparison with Experiments,” Journal of Fluid Mechanics, Vol. 570, 2007, pp. 347–364.
[106] Rudy, D. H. and Strikwerda, J. C., “A Nonreflecting Outflow Boundary Condition
for Subsonic Navier-Stokes Calculations,” Journal of Computational Physics, Vol. 36,
1980, pp. 55–70.
[107] Ghosal, S. and Moin, P., “The Basic Equations for the Large Eddy Simulation of
Turbulent Flows in Complex Geometries,” Journal of Computational Physics, Vol. 118,
1995, pp. 24–37.
[108] Horiuti, K., “A New Dynamic Two-Parameter Mixed Model for Large-Eddy Simulation,” Physics of Fluids, Vol. 9, No. 11, 1997, pp. 3443–3464.
[109] Vreman, B., Geurts, B., and Kuerten, H., “On the Formulation of the Dynamic Mixed
Subgrid-Scale Model,” Physics of Fluids, Vol. 6, No. 12, 1994, pp. 4057–4059.
[110] Sagaut, P. and Grohens, R., “Discrete Filters for Large Eddy Simulation,” International Journal for Numerical Methods in Fluids, Vol. 31, 1999, pp. 1195–1220.
[111] Germano, M., Piomelli, U., Moin, P., and Cabot, W. H., “A Dynamic Subgrid-Scale
Eddy Viscosity Model,” Physics of Fluids, Vol. 3, No. 7, 1991, pp. 1760–1765.
[112] Erlebacher, G., Hussaini, M. Y., Speziale, C. G., and Zang, T. A., “Toward the
Large-Eddy Simulation of Compressible Turbulent Flows,” Journal of Fluid Mechanics, Vol. 238, 1992, pp. 155–185.
[113] Wu, X. and Moin, P., “Direct Numerical Simulation of Turbulence in a Nominally ZeroPressure Gradient Flat Plate Boundary Layer,” Journal of Fluid Mechanics, Vol. 630,
2009, pp. 5–41.
165

[114] Li, Z. and Jaberi, F. A., “Numerical Investigations of Shock-Turbulence Interactions in
a Planar Mixing Layer,” 48th AIAA Aerospace Sciences Meeting and Exhibit, Orlando,
Florida, January 2010.
[115] Li, Z. and Jaberi, F. A., “Large Scale Simulations of High-Speed Turbulent Flows,”
47th AIAA Aerospace Sciences Meeting and Exhibit, Orlando, Florida, January 2009.
[116] Banaeizadeh, A., Li, Z., and Jaberi, F. A., “Large-Eddy Simulations of Compressible Turbulent Reacting Flows,” 48th AIAA Aerospace Sciences Meeting and Exhibit,
Orlando, Florida, January 2010.
[117] Li, Q. and Coleman, G., “DNS of an Oblique Shock Wave Impinging Upon a Turbulent
Boundary Layer,” Direct and Large-Eddy Simulation V (ERCOFTAC Series 9), edited
by R. Friedrich, B. J. Geurts, and O. Metais, 2003, pp. 387–396.
[118] Harten, A., Engquist, B., Osher, S., and Chakravarthy, S. R., “Uniformly High-Order
Accurate Essentially Nonoscillatory Schemes, III,” Journal of Computational Physics,
Vol. 71, 1987, pp. 231–303.
[119] Shu, C. W. and Osher, S., “Efficient Implementation of Essentially Non-Oscillatory
Shock-Capturing Schemes,” Journal of Computational Physics, Vol. 77, 1988, pp. 439–
471.
[120] Zang, Y., Street, R. L., and Koseff, J. R., “A Dynamic Mixed Subgrid-Scale Model and
its Application to Turbulent Recirculating Flows,” Physics of Fluids, Vol. 5, No. 12,
1993, pp. 3186–3196.
[121] Jiang, G. and Shu, C. W., “Efficient Implementation of Weighted ENO Schemes,”
Journal of Computational Physics, Vol. 126, 1996, pp. 202–228.
[122] Jaberi, F. A., Colucci, P. J., James, S., Givi, P., and Pope, S. B., “Filtered Mass
Density Function for Large-Eddy Simulation of Turbulent Reacting Flows,” Journal
of Fluid Mechanics, Vol. 401, 1999, pp. 85–122.
[123] Lenormand, E., Sagaut, P., and Ta Phuoc, L., “Large-Eddy Simulation of Subsonic
and Supersonic Channel Flow at Moderate Reynolds Number,” International Journal
for Numerical Methods in Fluids, Vol. 32, 2000, pp. 369–406.
[124] Pirozzoli, S., Bernardini, M., and Grasso, F., “Characterization of Coherent Vortical
Structures in a Supersonic Turbulent Boundary Layer,” Journal of Fluid Mechanics,
Vol. 613, 2008, pp. 205–231.

166