This draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics
1
Spectral analysis of jet turbulence
Oliver T. Schmidt
1
†
, Aaron Towne
2
, Georgios Rigas
1
, Tim Colonius
1
,
and Guillaume A. Br`es
3
1
California Institute of Technology, Pasadena, CA 91125, USA
2
Center for Turbulence Research, Stanford University, Stanford, CA 94305, USA
3
Cascade Technologies Inc., Palo Alto, CA 94303, USA
(Received xx; revised xx; accepted xx)
Informed by LES data and resolvent analysis of the mean flow, we examine the struc-
ture of turbulence in jets in the subsonic, transonic, and supersonic regimes. Spectral
(frequency-space) proper orthogonal decomposition is used to extract energy spectra
and decompose the flow into energy-ranked coherent structures. The educed structures
are generally well predicted by the resolvent analysis. Over a range of low frequencies
and the first few azimuthal mode numbers, these jets exhibit a low-rank response
characterized by Kelvin-Helmholtz (KH) type wavepackets associated with the annular
shear layer up to the end of the potential core and that are excited by forcing in the
very-near-nozzle shear layer. These modes too the have been experimentally observed
before and predicted by quasi-parallel stability theory and other approximations–they
comprise a considerable portion of the total turbulent energy. At still lower frequencies,
particularly for the axisymmetric mode, and again at high frequencies for all azimuthal
wavenumbers, the response is not low rank, but consists of a family of similarly amplified
modes. These modes, which are primarily active downstream of the potential core, are
associated with the Orr mechanism. They occur also as sub-dominant modes in the range
of frequencies dominated by the KH response. Our global analysis helps tie together
previous observations based on local spatial stability theory, and explains why quasi-
parallel predictions were successful at some frequencies and azimuthal wavenumbers, but
failed at others.
Key words:
1. Introduction
Large-scale coherent structures in the form of wavepackets play an important role in
the dynamics and acoustics of turbulent jets. In particular, their spatial coherence makes
wavepackets efficient sources of sound (Crighton & Huerre 1990; Jordan & Colonius 2013).
They are most easily observed in forced experiments, where periodic exertion establishes
a phase-reference. The measurements of a periodically forced turbulent jet by Crow &
Champagne (1971) served as a reference case for wavepacket models. Early examples of
such models include the studies by Michalke (1971) and Crighton & Gaster (1976), who
established the idea that the coherent structures can be interpreted as linear instability
waves evolving around the turbulent mean flow.
Wavepackets in unforced jets exhibit intermittent behavior (Cavalieri
et al.
2011) and
†
Email address for correspondence: oschmidt@caltech.edu
arXiv:1711.06296v1 [physics.flu-dyn] 16 Nov 2017
2
O. Schmidt, A. Towne, G. Rigas, T. Colonius, G. Br`es
as a result are best understood as statistical objects that emerge from the stochastic
turbulent flow. We use spectral proper orthogonal decomposition (Lumley 1970; Towne
et al.
2017
c
) to extract these structures from the turbulent flow. SPOD has been applied
to a range of jets using data from both experimental (Glauser
et al.
1987; Arndt
et al.
1997; Citriniti & George 2000; Suzuki & Colonius 2006; Gudmundsson & Colonius 2011),
and numerical studies (Sinha
et al.
2014; Towne
et al.
2015; Schmidt
et al.
2017).
Wavepackets have been extensively studied, mainly because of their prominent role
in the production of jet noise, and models based on the parabolized stability equations
(PSE) have proven to be successful at modeling them (Gudmundsson & Colonius 2011;
Cavalieri
et al.
2013; Sinha
et al.
2014). This agreement between SPOD modes and PSE
solution breaks down at low frequencies and for some azimuthal wavenumbers, and in
general downstream of the potential core. In the present study, we show that the SPOD
eigenspectra unveil low-rank dynamics, and we inspect the corresponding modes and
compare them with predictions based on resolvent analysis. The results show that two
different mechanisms are active in turbulent jets, and they explain the success and failure
of linear and PSE models.
Resolvent analysis of the turbulent mean flow field seeks sets of forcing and response
modes that are optimal with respect to the energetic gain between them. When applied
to the mean of a fully turbulent flow, the resolvent forcing modes can be associated with
nonlinear modal interactions (McKeon & Sharma 2010) as well as stochastic inputs to the
flow, for example the turbulent boundary layer in the nozzle that feeds the jet. Garnaud
et al.
(2013) interpreted the results of their resolvent analysis of a turbulent jet in the light
of experimental studies of forced jets by Moore (1977) and Crow & Champagne (1971).
They found that the frequency of the largest gain approximately corresponds to what the
experimentalists referred to as the
preferred frequency
, i.e. the frequency where external
harmonic forcing in the experiments triggered the largest response. In the context of jet
aeroacoustics, Jeun
et al.
(2016) restricted the optimal forcing to vortical perturbations
close to the jet axis, and the optimal responses to the far-field pressure. Their results
show that suboptimal modes have to be considered in resolvent-based jet noise models.
In this paper, we use resolvent analysis to model and explain the low-rank behavior
of turbulent jets revealed by SPOD. Recent theoretical connections between SPOD and
resolvent analysis (Towne
et al.
2015; Semeraro
et al.
2016; Towne
et al.
2017
c
) make
the latter a natural tool for this endeavor and provide a framework for interpreting our
results.
The remainder of the paper is organized as follows. The three large-eddy simulation
databases used to study different Mach number regimes are introduced in
§
2. At first,
the focus is on the lowest Mach number case in
§§
3-5. We start by analyzing the data of
this case using (mainly) SPOD in
§
3, followed by the resolvent analysis in
§
4. The results
of the SPOD and the resolvent model are compared in
§
5.
§
6 addresses Mach number
effects observed in the remaining two cases representing the transonic and the supersonic
regime. Finally, the results are discussed in
§
7. For completeness, we report in appendix
A all results for the higher Mach number cases that were omitted in
§
6 for brevity.
2. Large eddy simulation
The unstructured flow solver “Charles” (Br`es
et al.
2017
b
) is used to perform large-
eddy simulations of three turbulent jets at Mach numbers
M
j
=
U
j
/a
j
of 0.4, 0.9,
and 1.5. All jets are isothermal and the supersonic jet is ideally expanded. The nozzle
geometry is included in the computational domain, and synthetic turbulence combined
with a wall model is applied inside the nozzle to obtain a fully turbulent boundary
Guidelines for authors
3
case
M
j
Re
p
0
p
∞
T
0
T
∞
n
cells
d
tc
∞
D
t
sim
c
∞
D
subsonic
0.4 0
.
45
·
10
6
1.117 1.03 15
.
9
·
10
6
1
·
10
−
3
2000
transsonic 0.9 1
.
01
·
10
6
1.7
1.15 15
.
9
·
10
6
1
·
10
−
3
2000
supersonic 1.5 1
.
76
·
10
6
3.67 1.45 31
·
10
6
2
.
5
·
10
−
4
500
Table 1: Parameters of the large-eddy simulations.
layer inside the nozzle. The jets are further characterized by their Reynolds numbers
Re
=
ρ
j
U
j
D/μ
j
, which correspond approximately to the laboratory values in a set of
companion experiments. The reader is referred to Br`es
et al.
(2017
b
) for further details on
the numerical method, meshing strategy and subgrid-scale model. A detailed validation
of the
M
j
= 0
.
9 jet including the nozzle-interior turbulence modeling (i.e., synthetic
turbulence, wall model) can be found in Br`es
et al.
(2017
a
). The subscripts
j
and
∞
refer
to jet and free-stream conditions,
a
is the speed of sound,
ρ
density,
D
nozzle diameter,
μ
dynamic viscosity,
T
temperature, and
U
j
to the axial jet velocity on the centerline of
the nozzle exit, respectively. Throughout this paper, the flow is non-dimensionalized by
its nozzle exit values, pressure by
ρ
j
U
2
j
, lengths by
D
, and time by
D/U
j
. Frequencies
are reported in terms of the Strouhal number
St
=
ω/
(2
πM
j
), where
ω
is the non-
dimensional angular frequency. The parameters of the three simulations are listed in
table 1:
p
0
/p
∞
is the nozzle pressure ratio,
T
0
/T
∞
the nozzle temperature ratio,
n
cells
the number of control volumes, d
tc
∞
/D
the computational time step, and
t
sim
/c
∞
D
the
total simulation time after the flow became stationary, i.e. without initial transients. The
unstructured LES data is first interpolated onto a
n
x
×
n
r
×
n
θ
structured cylindrical
grid spanning
x,r,θ
∈
[0
,
30]
×
[0
,
6]
×
[0
,
2
π
]. Snapshots are saved with a temporal
separation
∆tc
∞
/D
(in acoustical units). For both the spectral analysis in
§
3, and the
resolvent model in
§
4, we Reynolds decompose a flow quantity
q
into the long-time mean
denoted by
(
·
) and the fluctuating part (
·
)
′
as
q
(
x,r,θ,t
) =
q
(
x,r,θ
) +
q
′
(
x,r,θ,t
)
.
(2.1)
A visualization of the subsonic jet is shown in figure 1. Only the domain of interest for
this study is shown (the full LES domain is much larger and includes flow within the
nozzle as well as far-field sponge regions).
The mean centerline velocity of the three jets is compared in figure 2(a). The plateau
close to the nozzle characterizes the potential core, whose length increases with Mach
number. A weak residual shock pattern is observed for
x
.
10 in the supersonic case.
The streamwise development of turbulent jets is usually described in terms of an
initial
development region
(0
.
x/D
.
25), and a
self-similar region
(
x/D
&
25), see e.g. Pope
(2000). In the latter, the jet is fully described by a self-similar velocity profile, and a
constant spreading and velocity-decay rate. For a dynamical description of the flow,
we further divide the initial development region into two parts. The
initial shear-layer
region
extends up to the end of the potential core (0
.
x/D
.
5 for
M
= 0
.
4) and is
characterized by a constant velocity jump over the shear-layer, and a linearly increasing
shear-layer thickness. The
developing jet region
lies between the end of the potential core
and the start of the self-similar region (5
.
x/D
.
25). In this region, the centerline
velocity transitions rapidly to its asymptotic decaying behavior. The mean radial velocity
profile, on the other hand, has not yet converged to its self-similar downstream solution.
The 1
/x
-scaling of the centerline velocity with streamwise distance (Pope 2000) is
indicated for the subsonic case. The momentum thickness shown in figure 2(b) increases
4
O. Schmidt, A. Towne, G. Rigas, T. Colonius, G. Br`es
Figure 1: Instantaneous flow field of the subsonic jet: (a) streamwise cross-section along
the jet axis; (b-e) transverse planes at different streamwise locations
x
. The streamwise
velocity fluctuation
u
′
x
is shown. The potential core and the jet width are indicated as
lines of constant ̄
u
x
at 99% and 5% of the jet velocity
U
j
, respectively.
Figure 2: Mean flow characteristics for the three jets: (a) centerline velocity
u
c
=
u
x
(
x,r
= 0); (b) momentum thickness
δ
θ
=
∫
r
1
0
̄
ρ
̄
u
x
ρ
c
u
c
(
1
−
̄
u
x
u
c
)
d
r
, where the integral
is taken from the centerline to
r
1
defined as ̄
u
x
(
r
1
) = ̄
u
∞
+ 0
.
01 ̄
u
j
(the ̄
u
∞
term is
included to account for a small coflow included in the simulations). The centerline velocity
becomes inversely proportional to the axial distance shortly after the potential core, and
the momentum thickness increases approximately linearly.
approximately linear in all three regions in all cases. These characteristic velocities and
length scales in each region imply different frequency scalings that will become important
later.
Guidelines for authors
5
Interpolated database
SPOD
case
∆tc
∞
D
n
x
n
r
n
θ
x
1
r
1
n
t
n
freq
n
ovlp
n
blk
subsonic
0.2
656 138 128 30 6 10,000
256
128
78
transsonic
0.2
656 138 128 30 6 10,000
256
128
78
supersonic
0.1
698 136 128 30 6 5,000
256
128
39
Table 2: Parameters of the structured cylindrical grid of the interpolated database (left),
and the SPOD parameters (right).
In
§§
3-5, we will focus on the
M
= 0
.
4 jet as the main conclusions are similar for all
three Mach number regimes, and we address Mach number effects in detail in
§
6.
3. Spectral analysis of the LES data
In this section, we use spectral proper orthogonal decomposition (SPOD) to identify
coherent structures within the three turbulent jets. This form of proper orthogonal
decomposition (POD) identifies energy-ranked modes that each oscillate at a single
frequency, are orthogonal to other modes at the same frequency, and, as a set, optimally
represent the space-time flow statistics. SPOD was introduced by Lumley (1967, 1970)
but has been used sparingly compared to the common spatial form of POD (Sirovich
1987; Aubry 1991) and dynamic mode decomposition (Schmid 2010). However, recent
work by Towne
et al.
(2017
c
) showed that SPOD combines the advantages of these other
two methods–SPOD modes represent coherent structures that are dynamically significant
and optimally account for the random nature of turbulent flows. This makes SPOD an
ideal tool for identifying coherent structures within the turbulent jets considered in this
paper.
We seek modes that are orthogonal in the inner product
〈
q
1
,
q
2
〉
E
=
∫∫∫
q
∗
1
diag
(
T
γ
ρM
2
,
ρ,
ρ,
ρ,
ρ
γ
(
γ
−
1)
TM
2
)
q
2
r
d
x
d
r
d
θ,
(3.1)
which are optimal in an induced compressible energy norm
〈·
,
·〉
E
(Chu 1965). The
energy weights are defined for the state vector
q
= [
ρ,u
x
,u
r
,u
θ
,T
]
T
(
x,r,θ,t
) of primitive
variables density
ρ
, cylindrical velocity components
u
x
,
u
r
and
u
θ
, and temperature
T
.
The notation (
·
)
∗
indicates the Hermitian transpose. We discretize the inner product
defined by equation (3.1) as
〈
q
1
,
q
2
〉
E
=
q
∗
1
W
q
2
,
(3.2)
where the weight matrix
W
accounts for both numerical quadrature weights and the
energy weights. Since the jet is stationary and symmetric with respect to rotation about
the jet axis, it can be decomposed into azimuthal Fourier modes
ˆ
(
·
)
m
of azimuthal
wavenumber
m
, temporal Fourier modes
ˆ
(
·
)
ω
of angular frequency
ω
, or combined spatio-
temporal Fourier modes Fourier modes
ˆ
(
·
)
mω
as
q
(
x,r,θ,t
) =
∑
m
ˆ
q
m
(
x,r,t
)e
i
mθ
=
∑
ω
ˆ
q
ω
(
x,r,θ
)e
i
ωt
=
∑
m
∑
ω
̃
q
mω
(
x,r
)e
i
mθ
e
i
ωt
.
(3.3)
To calculate the SPOD, the data is first segmented into sequences
Q
=
[
q
(1)
q
(2)
···
q
(
n
freq
)
]
each containing
n
freq
instantaneous snapshots of
q
which are considered to be
statistically independent realizations of the flow under the ergodic hypothesis.
6
O. Schmidt, A. Towne, G. Rigas, T. Colonius, G. Br`es
Details on the interpolated databases and spectral estimation parameters are
listed in table 2. The spatio-temporal Fourier transform of the
l
-th block yields
ˆ
Q
(
l
)
mω
k
=
[
ˆ
q
(
l
)
mω
1
ˆ
q
(
l
)
mω
2
···
ˆ
q
(
l
)
mω
n
freq
]
, where
ˆ
q
(
l
)
mω
k
is the
l
-th realization of the Fourier
transform at the
k
-th discrete frequency. A periodic Hann window is used to minimize
spectral leakage. The ensemble of
n
blk
Fourier realizations of the flow at a given
frequency
ω
k
and azimuthal wavenumber
m
are now collected into a data matrix
ˆ
Q
mω
k
=
[
ˆ
q
(1)
mω
k
ˆ
q
(2)
mω
k
···
ˆ
q
(
n
blk
)
mω
k
]
. For a particular
m
and
ω
k
, the SPOD modes are
found as the eigenvectors
Ψ
mω
k
=
[
ψ
(1)
mω
k
ψ
(2)
mω
k
···
ψ
(
n
blk
)
mω
l
]
, and the modal energy as
the corresponding eigenvalues
Λ
mω
k
= diag
(
λ
(1)
mω
k
,λ
(2)
mω
k
,
···
,λ
(
n
blk
)
mω
k
)
of the weighted
cross-spectral density matrix
ˆ
S
mω
k
=
ˆ
Q
mω
k
ˆ
Q
∗
mω
k
as
ˆ
S
mω
k
WΨ
mω
k
=
Ψ
mω
k
Λ
mω
k
.
(3.4)
The modes are sorted by decreasing energy, i.e.
λ
(1)
mω
k
>
λ
(2)
mω
k
>
···
>
λ
(
n
blk
)
mω
k
. This
formulation guarantees that the SPOD modes have the desired orthonormality property
〈
ψ
(
i
)
mω
k
,
ψ
(
j
)
mω
k
〉
E
=
δ
ij
, where
δ
ij
is the Kronecker delta function, and are optimal in
terms of modal energy in the norm induced by equation (3.2). For brevity, we denote the
l
-th eigenvalue and the pressure component of the corresponding eigenmode as
λ
l
and
ψ
(
l
)
p
, respectively. The dependence on a specific azimuthal wavenumber and frequency
is implied and given in the description. Since the SPOD modes are optimal in terms of
energy, we sometimes refer to the first SPOD mode as the
leading
or
optimal
mode and
to the subsequent lower-energy modes as
suboptimal
modes.
Since we wish to express the data in terms of modes that oscillate at real and positive
frequencies, we take the temporal Fourier transform in equation (3.3) first. Statistical
homogeneity in
θ
implies that averaged quantities are the same for any
±
m
. After
verifying statistical convergence, we add the contributions of positive and negative
m
.
The distribution of power into the frequency components a signal is comprised of is
referred to as its power spectrum. Power spectra are commonly expressed in terms of
the power spectral density (PSD), which we will introduce later in equation (3.5). In
the context of SPOD, we are interested in finding a graphical representation that can
be interpreted in a similar way. In each frequency bin, the discrete SPOD spectrum is
represented by the decreasing energy levels of the corresponding set of eigenfunctions.
There is no obvious continuity in the modal structure of the most energetic mode (or
any other) from one frequency bin to the next. Nevertheless, it is instructive to examine
how the modal energy changes as a function of frequency, so that in what follows we plot
the SPOD eigenvalues for each mode
l
as functions
λ
l
(
St
) of frequency and refer to the
resulting set of curves as the
SPOD eigenvalue spectrum
.
The SPOD eigenvalue spectrum of the axisymmetric component of the subsonic jet is
shown in figure 3(a). The red shaded area highlights the separation between the first and
second modes. A large separation indicates that the leading mode is significantly more
energetic than the others. When this occurs, the physical mechanism associated with the
first mode is prevalent, and we say that the flow exhibits
low-rank behavior
.
The low-rank behavior is apparent over the frequency band 0
.
2
.
St
.
2, and peaks
at
St
≈
0
.
6. It is most pronounced for
m
= 0 and
m
= 1, shown in panel 3(a) and 3(b),
respectively. With increasing azimuthal wavenumber, the low-rank behaviour becomes
less and less pronounced. For
m
= 1, the dominance of the optimal gain persists to low
frequencies, whereas it cuts off below
St
≈
0
.
2 for
m
= 0. For
m
= 2 in figure 3(c), even