Resolvent-based study of compressibility effects on supersonic turbulent boundary layers

TL;DR

This paper uses resolvent analysis to examine compressibility effects on supersonic turbulent boundary layers, validating Morkovin's hypothesis and revealing universal scaling laws and flow structures across Mach numbers 2 to 4.

Contribution

It applies resolvent analysis to supersonic boundary layers, demonstrating the validity of Morkovin's hypothesis and identifying universal scaling behaviors of turbulence modes.

Findings

Resolvent modes exhibit Mach wave radiation in the supersonic region.

Semi-local scaling reveals universal and self-similar mode structures.

Principal modes effectively capture energy distribution between momentum and thermodynamics.

Abstract

The resolvent formulation of McKeon & Sharma (2010) is applied to supersonic turbulent boundary layers to study the validity of Morkovin's hypothesis, which postulates that high-speed turbulence structures in zero pressure-gradient turbulent boundary layers remain largely the same as its incompressible counterpart. Supersonic zero-pressure-gradient turbulent boundary layers with adiabatic wall boundary conditions at Mach numbers ranging from 2 to 4 are considered. Resolvent analysis highlights two distinct regions of the supersonic turbulent boundary layer in the wave parameter space: the relatively supersonic region and the relatively subsonic region. In the relatively supersonic region, where the flow is supersonic relative to the freestream, resolvent modes display structures consistent with Mach wave radiation that are absent in the incompressible regime. In the relatively subsonic region, we show that the low-rank approximation of the resolvent operator is an effective approximation of the full system and that the response modes predicted by the model exhibit universal and geometrically self-similar behaviour via a transformation given by the semi-local scaling. Moreover, with the semi-local scaling, we show that the resolvent modes follow the same scaling law as their incompressible counterparts in this region, which has implications for modelling and the prediction of turbulent high-speed wall-bounded flows. We also show that the thermodynamic variables exhibit similar mode shapes to the streamwise velocity modes, supporting the strong Reynolds analogy. Finally, we demonstrate that the principal resolvent modes can be used to capture the energy distribution between momentum and thermodynamic fluctuations.

Tables

Table 1. Table 1: The freestream Mach number M ∞ subscript 𝑀 M_{\infty} , friction Reynolds number \Rey τ subscript \Rey 𝜏 \Rey_{\tau} , momentum thickness Reynolds number based on the wall dynamic viscosity \Rey δ 2 subscript \Rey subscript 𝛿 2 \Rey_{\delta_{2}} , and grid resolutions for the different cases. N i subscript 𝑁 𝑖 N_{i} is the number of grid points in the x i subscript 𝑥 𝑖 x_{i} direction and N c subscript 𝑁 𝑐 N_{c} is the number of grid points for the wave speed c = ω / κ 1 𝑐 𝜔 subscript 𝜅 1 c=\omega/\kappa_{1} . Δ ​ λ i Δ subscript 𝜆 𝑖 \Delta\lambda_{i} and Δ ​ x 2 Δ subscript 𝑥 2 \Delta x_{2} are the grid resolutions in the wall-parallel and wall-normal directions, respectively.

$M_{\infty }$	${\text{Rey}}_{\tau }$	${\text{Rey}}_{{\delta }_2}$	$N_1$	$N_2$	$N_3$	$N_c$	$\log \left(\Delta {\lambda }_1/\delta \right)$	$\log \left(\Delta {\lambda }_3/\delta \right)$	$\Delta x_{2,\min }/\delta$
0	445.5	1100.0
2	447.7	1327.4
2	898.5	3027.2	81	401	81	25	0.05	0.05	$3.08\times {10}^{-5}$
3	502.0	1815.7
4	504.6	2129.6

Table 2. Table 2: Expected length scales for the universal modes of the resolvent operator for the turbulent boundary layer.

Class	${\kappa }_1$-scale	${\kappa }_3$-scale	$x_2$-scale	$c$-scale	${\left(u_1\right)}_1,{\left(\rho \right)}_1,{\left(T\right)}_1$-scale
Inner	$\frac{{\kappa }_1\delta }{{\text{Rey}}_{\tau }^{\ast }}$	$\frac{{\kappa }_3\delta }{{\text{Rey}}_{\tau }^{\ast }}$	$x_2^{\ast }$	$c^{\ast }$	${\widetilde{\left(q_i\right)}}_1\frac{1}{\sqrt{{\text{Rey}}_{\tau }^{\ast }}}$
Outer	${\kappa }_1\delta {\text{Rey}}_{\tau }^{\ast }$	${\kappa }_3\delta$	$\frac{x_2}{\delta }$	$\frac{c}{{\overline{u}}_{1,\infty }^{\star }}$	${\widetilde{\left(q_i\right)}}_1$
Logarithmic	${\kappa }_1{x}_2^c{x}_2^{c\ast }$	${\kappa }_3{x}_2^c$	$\frac{x_2}{x_2^{\text{c}}}$	–	${\widetilde{\left(q_i\right)}}_1\sqrt{\frac{x_2^c}{\delta }}$

Equations

$$\rho \left(\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}\right)=$$

$$\frac{\partial \rho }{\partial t}+u_j\frac{\partial \rho }{\partial x_j}=$$

$$\rho \left(\frac{\partial T}{\partial t}+u_j\frac{\partial T}{\partial x_j}\right)=$$

$$+\gamma (\gamma -1)\frac{M^2}{\text{Rey}}\mu \left[\frac{\partial u_i}{\partial x_j}\frac{\partial u_i}{\partial x_j}+\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_i}+\lambda {\left(\frac{\partial u_k}{\partial x_k}\right)}^2\right],$$

$$M=\frac{\breve{u}}{\sqrt{\gamma \mathcal{R}\breve{T}}},\text{Rey}=\frac{\breve{\rho }\breve{u}\breve{l}}{\breve{\mu }},\text{Pran}=\frac{\breve{\mu }{c}_p}{\breve{k}},$$

$$\mu (T)=\frac{T^{3/2}(1+C)}{T+C},$$

$$\mathbf{q}(x_1,x_2,x_3,t)={\iiint }_{-\infty }^{\infty }\hat{\mathbf{q}}(x_2;{\kappa }_1,{\kappa }_3,\omega )e^{\mathrm{i}({\kappa }_1{x}_1+{\kappa }_3{x}_3-\omega t)}\mathrm{d}{\kappa }_1\mathrm{d}{\kappa }_3\mathrm{d}\omega ,$$

$$-\mathrm{i}\omega {\hat{u}}_i+{\bar{u}}_1{\partial }_1{\hat{u}}_i+{\hat{u}}_2{\partial }_2{\bar{u}}_i=-\frac{1}{\gamma M^2}\left({\partial }_i\hat{T}+{\bar{T}}^2{\partial }_i\hat{\rho }+\hat{\rho }\bar{T}{\partial }_i\bar{T}+\bar{T}\hat{T}{\partial }_i\bar{\rho }\right)$$

$$+\frac{\bar{T}}{\text{Rey}}\left[\bar{\mu }{\partial }_j\left({\partial }_j{\hat{u}}_i+{\partial }_i{\hat{u}}_j\right)+\bar{\lambda }{\partial }_i\left({\partial }_j{\hat{u}}_j\right)+\frac{\partial \bar{\mu }}{\partial T}{\partial }_j\hat{T}\left({\partial }_j{\bar{u}}_i+{\partial }_i{\bar{u}}_j\right)\right]+{\hat{f}}_i$$

$$-\mathrm{i}\omega \hat{\rho }+{\bar{u}}_1{\partial }_1\hat{\rho }+{\hat{u}}_2{\partial }_2\bar{\rho }=-\bar{\rho }{\partial }_i{\hat{u}}_i+{\hat{f}}_4$$

$$-\mathrm{i}\omega \hat{T}+{\bar{u}}_1{\partial }_1\hat{T}+{\hat{u}}_2{\partial }_2\bar{T}=-(\gamma -1)\bar{T}{\partial }_i{\hat{u}}_i$$

$$+\frac{\gamma \bar{T}}{\text{Pran}\text{Rey}}\left[\bar{\mu }{\partial }_j{\partial }_j\hat{T}+\frac{{\partial }^2\bar{\mu }}{\partial T^2}({\partial }_2\bar{T}{)}^2\hat{T}+2\frac{\partial \bar{\mu }}{\partial T}{\partial }_2\bar{T}{\partial }_2\hat{T}+\frac{\partial \bar{\mu }}{\partial T}{\partial }_2^2\bar{T}\hat{T}\right]$$

$$+\gamma (\gamma -1)\frac{M^2\bar{T}}{\text{Rey}}\left[2\bar{\mu }{\partial }_2{\bar{u}}_1{\partial }_2{\hat{u}}_1+2\bar{\mu }{\partial }_2\bar{u}{\partial }_1{\hat{u}}_2+\frac{\partial \bar{\mu }}{\partial T}({\partial }_2{\bar{u}}_1{)}^2\hat{T}\right]+{\hat{f}}_5,$$

$$\frac{\mathrm{d}}{\mathrm{d}{x}_2}\left(\bar{\mu }\frac{\mathrm{d}{\bar{u}}_1}{\mathrm{d}{x}_2}\right)=0,\frac{\mathrm{d}}{\mathrm{d}{x}_2}\left(\frac{\bar{\mu }}{\text{Pran}}\frac{\mathrm{d}\bar{T}}{\mathrm{d}{x}_2}\right)+(\gamma -1)M^2\bar{\mu }\frac{\mathrm{d}{\bar{u}}_1}{\mathrm{d}{x}_2}=0,\bar{\rho }\bar{T}=1.$$

$$\hat{\mathbf{q}}(x_2;{\kappa }_1,{\kappa }_3,\omega )={\left[-\mathrm{i}\omega \text{mathsfbi}I+\text{mathsfbi}L({\kappa }_1,{\kappa }_3,\omega )\right]}^{-1}\hat{\mathbf{f}}(x_2;{\kappa }_1,{\kappa }_3,\omega ),$$

$$2K=(\mathbf{q},\mathbf{q}{)}_K=\Vert \mathbf{q}{\Vert }_K^2={\int }_0^{\infty }\bar{\rho }{u}_i^{†}{u}_i\mathrm{d}{x}_2.$$

$$2E=(\mathbf{q},\mathbf{q}{)}_E=\Vert \mathbf{q}{\Vert }_E^2$$

$$\text{mathsfbi}H=\sum _{j=1}^{\infty }{\mathbf{\psi }}_j(x_2;{\kappa }_1,{\kappa }_3,\omega ){\sigma }_j({\kappa }_1,{\kappa }_3,\omega ){\mathbf{\varphi }}_j^{†}(x_2;{\kappa }_1,{\kappa }_3,\omega ),$$

$$({\mathbf{\psi }}_i(x_2;{\kappa }_1,{\kappa }_3,\omega ),{\mathbf{\psi }}_j(x_2;{\kappa }_1,{\kappa }_3,\omega ){)}_E={\delta }_{ij},$$

$$({\mathbf{\varphi }}_i(x_2;{\kappa }_1,{\kappa }_3,\omega ),{\mathbf{\varphi }}_j(x_2;{\kappa }_1,{\kappa }_3,\omega ){)}_E={\delta }_{ij},$$

$$\hat{\mathbf{f}}(x_2;{\kappa }_1,{\kappa }_3,\omega )$$

$$\hat{\mathbf{q}}(x_2;{\kappa }_1,{\kappa }_3,\omega )$$

$$u_i(x_2=0)=0,(T-\bar{T})(x_2=0)=0.$$

$$(u_i-{\bar{u}}_i)(x_2\to \infty ),(\rho -\bar{\rho })(x_2\to \infty ),(T-\bar{T})(x_2\to \infty )<\infty ,$$

$$\mathbf{\chi }(x_2)=\sum _{n=0}^{N_2-1}{b}_n{T}_n\left(\frac{x_2-a}{x_2+a}\right)=\sum _{n=0}^{N_2-1}{b}_n{T}_n(x_2^{\prime }),$$

$$\overline{M}(x_2)=\frac{({\kappa }_1{\bar{u}}_1(x_2)-\omega )M_{\infty }}{{\left({\kappa }_1^2+{\kappa }_3^2\right)}^{1/2}\bar{T}(x_2{)}^{1/2}},$$

$$E_{qq}(x_2,\omega )=\sum _{i=1}^5\frac{\iint {\kappa }_1^2{\kappa }_3{\left[{\sigma }_1({\kappa }_1,{\kappa }_3,\omega )|(q_i{)}_1|({\kappa }_1,x_2,{\kappa }_3,\omega )\right]}^2\mathrm{d}\log {\kappa }_1\mathrm{d}\log {\kappa }_3}{\max \iint {\kappa }_1^2{\kappa }_3{\sigma }_1({\kappa }_1,{\kappa }_3,\omega {)}^2\mathrm{d}\log {\kappa }_1\mathrm{d}\log {\kappa }_3}.$$

$$\overline{M}=\frac{({\kappa }_1{\bar{u}}_1-{\kappa }_1\bar{c})M_{\infty }}{{\left({\kappa }_1^2+{\kappa }_3^2\right)}^{1/2}{\bar{T}}^{1/2}}=1$$

$$x_2^{\ast }$$

$${\bar{u}}_1^{\ast }$$

$${\text{Rey}}_{\tau }^{\ast }(x_2)=\frac{\bar{\rho }({\tau }_w/\bar{\rho }{)}^{1/2}\delta }{\bar{\mu }}$$

$${\bar{u}}_1^{\star }={\bar{u}}_1^{+}{\left(\frac{\bar{\rho }}{{\bar{\rho }}_w}\right)}^{1/2},$$

$$\frac{\bar{T}}{{\bar{T}}_{\infty }}=1+r\frac{\gamma -1}{2}{M}_{\infty }^2\left[1-{\left(\frac{{\bar{u}}_1}{{\bar{u}}_{1,\infty }}\right)}^2\right],$$

$${\text{Rey}}_{\tau }^{\ast }={\text{Rey}}_{\tau }\frac{\bar{\mu }/{\bar{\mu }}_w}{\sqrt{\bar{\rho }/{\bar{\rho }}_w}}={\text{Rey}}_{\tau }\frac{{\left(\bar{T}/{\bar{T}}_w\right)}^{3/2}+C_w}{{\left(\bar{T}/{\bar{T}}_w\right)}^{3/2}\left(1+C_w\right)},$$

$${\text{Rey}}_{\tau }^{\ast }(x_2\to \infty )={\text{Rey}}_{\tau }\frac{1+C_w{\left(1+r(\gamma -1)M_{\infty }^2/2\right)}^{3/2}}{1+C_w}$$

Full text

1 Introduction

The prediction and modelling of turbulent high-speed wall-bounded flows remain an active field of study for their tremendous technological importance in the aerospace industry with respect to high-speed vehicles. The turbulent boundary layers determine the aerodynamic drag and heat transfer, which are quantities of interest for accurate performance assessment. For supersonic and hypersonic turbulent boundary layers, a major role has been historically played by experiments, with direct numerical simulations (DNS) becoming more and more common in the last decade.

Experimental investigations of supersonic and hypersonic turbulent boundary layers have been conducted historically with hot-wire anemometry (Kistler, 1959; Laderman & Demetriades, 1974; Owen et al., 1975; Spina & Smits, 1987; Konrad & Smits, 1998) (see also Roy & Blottner, 2006, for a review), which suffers from uncertainties associated with mixed-mode sensitivity (Kovasznay, 1953) and was later shown to suffer from poor frequency response and spatial resolution (Williams et al., 2018). In addition to hot-wire anemometry, direct measurements of spatially varying velocity fields of high-speed turbulent boundary layers have been attempted using particle image velocimetry (Ekoto et al., 2008; Tichenor et al., 2013; Peltier et al., 2016; Williams et al., 2018), which range up to Mach number of 7.5 for a flat-plate turbulent boundary layer.

Complementary to experiments, DNS of high-speed turbulent boundary layers have been conducted to overcome the experimental difficulties and gain access to the full three-dimensional structure of the turbulent flow field. Several DNS have been conducted with an emphasis on studying Morkovin’s scaling in turbulent boundary layers at moderate freestream Mach numbers (Guarini et al., 2000; Maeder, 2000; Pirozzoli et al., 2004; Martín, 2007; Shahab et al., 2011; Bernardini & Pirozzoli, 2011; Pirozzoli & Bernardini, 2011; Hadjadj et al., 2015; Poggie et al., 2015; Trettel & Larsson, 2016; Modesti & Pirozzoli, 2016) for both adiabatic and isothermal wall boundary conditions. Hypersonic studies for turbulent boundary layers at higher Mach numbers have also recently been available (Duan et al., 2010, 2011; Duan & Martin, 2011; Lagha et al., 2011; Zhang et al., 2014, 2018) with data sets of up to Mach number of 20.

In spite of the recent developments in numerical experiments, simulations of supersonic turbulent boundary layers still remain a daunting challenge. DNS in the incompressible regime (Simens et al., 2009), as well as earlier experiments (Erm & Joubert, 1991), show that a fully developed state of the boundary layer requires the use of an extremely long computational domain, which makes accurate numerical simulations computationally demanding. Simulations in the supersonic regime are further slowed down by the inherently larger computational effort and the possible occurrence of ‘eddy shocklets’. This creates a critical demand for model-based approaches that describe and predict the behaviour of turbulent flows at technologically relevant high Reynolds numbers in the supersonic regime. In particular, assessing Morkovin’s hypothesis (Morkovin, 1962) has been the focus of many high-speed turbulent boundary layer studies. Morkovin (1962) concluded from the analysis of supersonic boundary layer data available at the time that for moderate Mach numbers “the essential dynamics of these shear flows will follow the incompressible pattern.” The hypothesis was used and reformulated by Bradshaw (1974) to indicate that high-speed boundary layers can be computed using the same model as that of low speeds as long as the density fluctuations are weak. Another consequence of Morkovin’s hypothesis is the analogy between the temperature and velocity fields that leads to velocity-temperature relations such as the classical Walz formula (Walz, 1969) and the strong Reynolds analogy and its variants (Morkovin, 1962; Gaviglio, 1987; Huang et al., 1995; Zhang et al., 2014). It also motivates the so-called ‘compressibility transformations’ that transform the mean velocity and Reynolds stress profiles from a compressible boundary layer to equivalent incompressible profiles by accounting for mean property variations across the thickness of the boundary layer (Van Driest, 1951; Brun et al., 2008; Zhang et al., 2012; Trettel & Larsson, 2016; Yang & Lv, 2018). Moreover, previous studies of coherent structures in compressible boundary layers show that the structure of the flow associated with the near-wall region shares several features with the incompressible boundary layer as well as the so-called internal shear layers (Pirozzoli et al., 2008; Ringuette et al., 2008; Lagha et al., 2011), further supporting the hypothesis.

Resolvent analysis provides a tool to understand and predict flows using a low-rank approximation of the linear sub-system with the nonlinear interactions treated as a forcing term. Previous studies on incompressible flows show that resolvent analysis can capture a range of phenomena already observed in wall turbulence from near-wall streaks and quasi-streamwise vortices (McKeon & Sharma, 2010), hairpin vortices (Sharma & McKeon, 2013) and the corresponding pressure signature (Luhar et al., 2014), very large scale motions (McKeon & Sharma, 2010), to scaling of the statistics (Moarref et al., 2013). A similar approach, developed by Jovanović & Bamieh (2005), was used to study the roles of various disturbances in laminar boundary layer transition and optimal transient growth in laminar and turbulent boundary layers (Monokrousos et al., 2011; Sipp & Marquet, 2013; Alizard et al., 2015). As such, resolvent analysis yields an efficient basis for the flow, which can be used to provide a low order representation of the key dynamical processes of turbulence. Several studies have investigated more sophisticated means of shaping the nonlinear forcing from data and analytical considerations (Moarref et al., 2014; Zare et al., 2017; Illingworth et al., 2018; Towne et al., 2019; Morra et al., 2019). For turbulent channel flows (and consequently for turbulent boundary layers), the resolvent velocity modes have a universal scaling with Reynolds number, and in the overlap region of the mean velocity, admit geometric self-similarity (Moarref et al., 2013, 2014). This observation has strong implications for modelling: resolvent modes throughout the entire region can be described in terms of modes assessed at one wall-normal plane. For example, under the correct scaling of the wall-parallel wavenumbers, the self-similar hierarchies of response modes give rise to self-similar families of vortical structures. However, there has been less work applying this operator based decomposition to fully-developed compressible flows, with the exception of the recent studies on compressible jet flows (Jeun et al., 2016; Towne et al., 2018; Schmidt et al., 2018), subsonic aerofoils (Yeh & Taira, 2019), turbulent Couette flows (Dawson & McKeon, 2019), and hypersonic boundary layer transitions (Cook et al., 2018; Dwivedi et al., 2018).

The main goal of this work is to apply the resolvent analysis framework for compressible flows to supersonic turbulent boundary layers. The resolvent framework allows the decomposition of the governing equations in the wavenumber-frequency space, permitting an in-depth comparison of the underlying mechanisms in the flow. We utilise this tool to compare the mechanisms driving the incompressible and compressible boundary layer which will allow us to assess Morkovin’s hypothesis on a mode-by-mode basis.

The paper is organised as follows. We first introduce the compressible Navier-Stokes equations and the resolvent operator in §2, where we discuss the relevant resolvent norm, boundary conditions, and computational methods. In §3, we discuss the characteristics of the resolvent modes for the supersonic turbulent boundary layer and define the relatively supersonic and subsonic region. We show that the response modes in the relatively supersonic region display Mach waves. We also highlight the low-rank behaviour of the resolvent operator and discuss the necessary conditions for the universality of the resolvent modes in the relatively subsonic region. In §4, we then provide the Reynolds and Mach number scaling for the principal resolvent modes and amplification factor for the inner, logarithmic and outer region of the boundary layer and demonstrate that the leading response mode is enough to predict the energy distribution between momentum and thermodynamic fluctuations. Finally, we conclude the paper in §5.

2 Resolvent formulation of compressible zero-pressure-gradient turbulent boundary layer

2.1 Compressible Navier-Stokes equations

The nondimensional compressible Navier-Stokes equations for a perfect gas are given by

[TABLE]

[TABLE]

where $\rho$, $p$, $u_{i}$, $T$ are, respectively, density, pressure, velocity components, and temperature. Variables $\mu$ and $\lambda$ are the coefficients of first and second viscosity, respectively, $k$ is the coefficient of thermal conductivity, $\gamma=c_{p}/c_{v}$ is the ratio of specific heats, and $\delta_{ij}$ is the Kronecker delta. We formulate the equations in nondimensionalised form using the Mach, Reynolds, and Prandtl numbers, given respectively by

[TABLE]

where $\breve{(\cdot)}$ denotes reference (dimensional) quantities, $l$ is a length scale, and $\mathcal{R}$ is the universal gas constant. The system is closed with the equation of state, $p=\rho T$.

Here, we assume constant specific heat coefficients and constant Prandtl number, $\Pran=0.72$, and we set $\gamma=1.4$ (diatomic gas). Furthermore, we assume that viscosity varies with temperature according to the Sutherland formula

[TABLE]

with $C=S/\breve{T}$, where $S=110.4$K, and that the second coefficient of viscosity $\lambda$ follows the Stokes’ assumption $\lambda=-2/3\mu$.

2.2 Resolvent operator

Assuming a fully developed, locally parallel flow with the directions $x_{1}$, $x_{2}$ and $x_{3}$ signifying the streamwise, wall-normal, and spanwise directions, respectively, the state variable $\boldsymbol{q}=[q_{1},q_{2},q_{3},q_{4},q_{5}]^{\intercal}=[u_{1},u_{2},u_{3},\rho,T]^{\intercal}$ is decomposed using the Fourier transform in homogeneous directions and time,

[TABLE]

where $\hat{(\cdot)}$ denotes variables in the transformed domain, $\mathrm{i}=\sqrt{-1}$, and the triplet $(\kappa_{1},\kappa_{3},\omega)$ identifies the streamwise and spanwise wavenumbers and the temporal frequency, respectively. Here, the superscript $\intercal$ denotes transpose and $\dagger$ will denote conjugate transpose.

The mean turbulent state, $\bar{\boldsymbol{q}}(x_{2})=[\bar{u}_{1}(x_{2}),0,0,\bar{\rho}(x_{2}),\bar{T}(x_{2})]^{\intercal}$, corresponds to $(\kappa_{1},\kappa_{3},\omega)=(0,0,0)$. Furthermore, with the parallel-flow assumption, which is reasonable as the base flow variations of the zero-pressure-gradient turbulent boundary layer depend on the viscous time scale compared to the much faster convective time scale for fluctuations, the mean momentum equation (1) gives a constant $\bar{p}(x_{2})$. In the remainder of the paper, we scale the pressure such that $\bar{p}=1$ for simplicity. Note that the spatially-developing turbulent boundary layer can still be studied using resolvent analysis, and the associated increase in computational effort is not prohibitive. However, the interpretation of the underlying physical mechanisms is significantly more straightforward for the quasi-parallel, one-dimensional mean.

Following a similar approach to McKeon & Sharma (2010), the governing equations (1)–(3) can be rewritten in the Fourier domain for each $(\kappa_{1},\kappa_{3},\omega)\neq(0,0,0)$ as

[TABLE]

where $\hat{\boldsymbol{f}}$ contains the nonlinear terms and $(\partial_{1},\partial_{2},\partial_{3})=(\mathrm{i}\kappa_{1},\mathrm{d}/\mathrm{d}x_{2},\mathrm{i}\kappa_{3})$ (see Mack, 1984, for details on the formation of the linear operator), and the mean state equations as

[TABLE]

The full set of equations characterised by the mean and the fluctuating equations are unclosed; however, by assuming that the mean velocity profile is known, the shape of the mean profile acts as a constraint on the full closure of the nonlinear terms for the fluctuating equations. The Fourier domain equations can then be equivalently expressed as

[TABLE]

where $\mathsfbi{I}$ is the identity matrix and $\mathsfbi{L}$ is the linearised operator of the governing equations around the supersonic turbulent mean profile. The operator $\mathsfbi{H}=\left[-\mathrm{i}\omega\mathsfbi{I}+\mathsfbi{L}(\kappa_{1},\kappa_{3},\omega)\right]^{-1}$ is called the resolvent operator and exists if there are no eigenvalues of $\mathsfbi{L}$ equal to $\mathrm{i}\omega$.

2.3 Choice of resolvent norm

From (11), we wish to find a decomposition of the resolvent operator that enables us to identify high gain input and output modes with respect to the linear operator. For the resolvent analysis, the decomposition is given by the Schmidt decomposition (called the singular value decomposition for the discrete case). However, this decomposition must be accompanied by a choice of inner product and the corresponding norm. In the case of the incompressible resolvent operator, the natural and physically meaningful norm is the kinetic energy norm, which is defined as

[TABLE]

Unfortunately, there is no obvious choice for the compressible case and the standard incompressible kinetic energy norm becomes a seminorm on this space; however, Chu (1965) introduced a norm that eliminates pressure-related energy transfer terms (compression work),

[TABLE]

which has been used in numerous other studies of compressible flows where the definition of an inner product is required (e.g., Hanifi et al., 1996; Malik et al., 2006; Özgen & Kırcalı, 2008; Malik et al., 2008; de Pando et al., 2014; Bitter & Shepherd, 2014; Dawson & McKeon, 2019), and this norm will be used for the remainder of the paper. Discussion of other possible inner products and assumptions for compressible flows are given in Rowley et al. (2004). A study on the sensitivity of the resolvent modes with respect to the compressible inner product and the standard incompressible kinetic energy inner product was performed by Dawson & McKeon (2019) for planar Couette flow, and the differences in mode shapes were attributed to changes in mean profiles and the compressible fluctuation equations rather than the choice of the inner product.

We take the Schmidt decomposition of the resolvent, namely,

[TABLE]

with an orthogonality condition

[TABLE]

and $\sigma_{j}\geq\sigma_{j+1}\geq 0$. The $\boldsymbol{\phi}_{j}$ and $\boldsymbol{\psi}_{j}$ form the right and left Schmidt bases (singular vectors) for the forcing and response fields, and the real $\sigma_{j}$ are the singular values. This decomposition is unique up to a pre-multiplying unitary complex factor on both bases corresponding to a phase shift (Young, 1988).

This basis pair can then be used to decompose any arbitrary forcing and the resulting state vector at a particular Fourier component such that

[TABLE]

Clearly the forcing shape that gives the largest energy is given by $a_{j}=\delta_{1j}$, i.e. when the forcing is aligned with the principal singular vector. Moreover, we later show that the resolvent operator is low-rank, i.e. $\sigma_{1}\gg\sigma_{2}$, where the flow is most energetic. Thus, in this paper, we focus on the principal Schmidt vectors (singular vectors), i.e. the principal forcing mode $\boldsymbol{\phi}_{1}$ and the principal response mode $\boldsymbol{\psi}_{1}=[(q_{1})_{1},(q_{2})_{1},(q_{3})_{1},(q_{4})_{1},(q_{5})_{1}]^{\intercal}=[(u_{1})_{1},(u_{2})_{1},(u_{3})_{1},(\rho)_{1},(T)_{1}]^{\intercal}$, and the principal singular value $\sigma_{1}$.

2.4 Boundary conditions for the resolvent operator

For the compressible boundary layer, the boundary conditions at the wall are given by

[TABLE]

The boundary conditions on the velocity fluctuations are the usual no-slip conditions, and the boundary condition on the temperature fluctuation is consistent for a gas flowing over a solid wall.

The boundary conditions at the freestream are given by

[TABLE]

which are less restrictive than requiring all fluctuations to be zero at infinity. However, in supersonic flow, waves may propagate to infinity and this boundary condition allows the waves with constant amplitude to be included. In practice, this is achieved by imposing a free boundary condition on the top boundary.

2.5 Computational methods

Following the Schmidt decomposition in (14), there are formally an infinite number of singular values. However, we solve the discrete equations using a spectral collocation method with the number of points in the wall-normal direction given by $N_{2}$, limiting the number of singular values to $5N_{2}$, the size of the state vector $\boldsymbol{q}$.

The discrete points in $x_{2}$ are given by a rational transformation of the Chebyshev collocation points. The Chebyshev collocation points are defined as $x_{2}^{\prime}=\cos({\pi j}/({N_{2}-1}))$ for $j=0,1,\dots,N_{2}-1$ in the domain $-1\leq x_{2}^{\prime}\leq 1$. The rational transformation $x_{2}=a{(1+x_{2}^{\prime})}/{(1-x_{2}^{\prime})}$ maps $x_{2}^{\prime}$ to the semi-infinite domain, where $a/\delta=2$ is the wall-normal location containing $(N_{2}-1)/2$ points (Grosch & Orszag, 1977; Christov, 1982) and $\delta$ is the boundary layer thickness corresponding to the location where mean velocity reaches $99\%$ of the freestream velocity. With this transformation, the value of a function $\boldsymbol{\chi}(x_{2})$ can be expressed as

[TABLE]

where $T_{n}$ is the $n$th order Chebyshev polynomial and $b_{n}$ is the coefficient for the corresponding Chebyshev polynomial.

The rational transformation of the Chebyshev collocation points allows the use of spectral methods for the semi-infinite domain. The stretching of the collocation points such that the majority of the points lie within $x_{2}/\delta\leq 2$ is appropriate for most of the energy-containing modes that are located within the boundary layer. However, in supersonic flows, waves may propagate to infinity and thus become under-resolved due to the fact the grid resolution $\Delta x_{2}\rightarrow\infty$ as $x_{2}\rightarrow\infty$. Nonetheless, the limitation stems from resolving a semi-infinite domain with a finite number of discrete points, and thus any choice of decomposition in the semi-infinite domain suffers from this limitation. The resolvent modes considered for this paper are primarily located within the boundary layer, and thus we are not affected by this limitation.

The turbulent mean profiles for the supersonic case are obtained from Bernardini & Pirozzoli (2011) and Pirozzoli & Bernardini (2011), where DNS of a spatially evolving zero-pressure-gradient supersonic turbulent boundary layer with the wall temperature set to its nominally adiabatic value are computed. The results of the resolvent analysis for the compressible boundary layer are compared against the results from the resolvent analysis of the incompressible boundary layer with the mean velocity profile obtained from Jiménez et al. (2010). The cases are chosen such that $\Rey_{\tau}$ is similar to the Reynolds number for the incompressible case ($\Rey_{\tau}\approx 450$). While $\Rey_{\tau}$ is not the ideal Reynolds number scaling for the outer region and a better comparison would be the momentum thickness Reynolds number $\Rey_{\delta_{2}}$, we make use of the friction Reynolds number regardless. An additional case with $M_{\infty}=2$ and $\Rey_{\tau}=900$ was chosen to observe Reynolds number effects. The specific cases used are listed in table 1 along with their respective resolution in numerical computations, where $N_{2}$ is the number of collocation grid points in the $x_{2}$ direction, $N_{1}$ and $N_{3}$ are the number of spatial frequencies for $\kappa_{1}$ and $\kappa_{3}$, and $N_{c}$ is the number of grid points for the wave speed $c=\omega/\kappa_{1}$.

3 Characteristics of resolvent modes for compressible turbulent boundary layer

We first examine the energy contained in the principal response mode. The energy contribution of $\boldsymbol{\psi}_{k}$ to the total response subject to broadband forcing in the wall-normal direction can be quantified by $\sigma_{k}^{2}/(\sum_{j}\sigma_{j}^{2})$. Figure 1 shows the principal energy contribution from the principal response mode $\boldsymbol{\psi}_{1}$ for the incompressible case and the $M_{\infty}=4$ case at two wall-normal locations $x_{2}^{+}=15$ and $x_{2}/\delta=0.2$, where the superscript $+$ denotes wall units defined in terms of $\bar{\rho}$ and $\mu$ at the wall and the friction velocity $u_{\tau}$.

The results from the incompressible and compressible turbulent boundary layers show similarities in the region where the principal energy contribution of the incompressible boundary layer is concentrated, i.e., where the low-rank approximation is valid for the incompressible regime. This region coincides with the most energetic wavenumbers from DNS of incompressible channel flows identified by the premultiplied energy spectra given by the contour lines in the figure. While the premultiplied energy spectra shown here are from a DNS of channel flow, many studies in the past show similarities in higher-order streamwise statistics for internal and external flows (Mochizuki & Nieuwstadt, 1996; Metzger & Klewicki, 2001; Del Alamo et al., 2004; Jiménez & Hoyas, 2008; Monty et al., 2009), especially in the near-wall region. Moreover, while comparisons between energy spectra of incompressible and compressible boundary layers are not available, first-order statistics of the two cases for adiabatic wall boundary conditions are known to collapse given a correct velocity transformation (Pirozzoli & Bernardini, 2011), making the comparison still meaningful. While data such as two-dimensional energy spectra at various wall-normal heights for incompressible and compressible boundary layers are not readily obtainable, access to such data could facilitate the analysis in the future.

The most notable difference between the incompressible and compressible results is the triangular region marked by freestream relative Mach number (Mack, 1984), $\overline{M}_{\infty}=\overline{M}(x_{2}\rightarrow\infty)$, greater than unity (figure 1(b)). The relative Mach number, defined as

[TABLE]

can be understood as the local Mach number of the mean flow in the direction of the wavenumber vector, $[\kappa_{1},\kappa_{3}]^{\intercal}$, relative to the wave speed at a given wall-normal location $x_{2}$. A three-dimensional depiction of the principal energy contribution as a function of streamwise and spanwise wavenumbers and wave speeds is given in figure 2 for the incompressible case and the compressible case at three different Mach numbers. It is clear that the region with $\overline{M}_{\infty}>1$, i.e. the relatively supersonic region, increases with Mach number and grows from the wall towards the freestream. In linear stability theory, $\overline{M}_{\infty}$ has been used to classify disturbances as subsonic, sonic, or supersonic (Mack, 1984; Schmid & Henningson, 2000). Moreover, it has been shown that if $\overline{M}_{\infty}>1$, a compressible boundary layer is unstable to inviscid waves regardless of any other feature of the velocity and temperature profiles (Mack, 1984). Considering that the family of modes with $\overline{M}_{\infty}>1$ does not have any counterpart in incompressible boundary layers, it is expected that the largest deviation between the compressible and incompressible boundary layers occurs in this regime. In particular, the irregular low-rank behaviour present in the relatively supersonic region in figure 1(b) and figure 2(b)–(d) is due to the discrete acoustic eigenmodes of the system approaching the wave speed $c$ (Appendix A) and thus giving resonant amplification of the resolvent operator (Dawson & McKeon, 2019). This is a consequence of discretisation, and in the continuous case, this resonant amplification effect will be present everywhere in the relatively supersonic region. The extent of the relatively supersonic region, where the incompressible and compressible boundary layers show a significant difference, increase with increasing Mach number, and may be an indicator of why Morkovin’s hypothesis fails for high Mach numbers.

The second, although less significant, difference can be seen for the two spectra at $x_{2}/\delta=0.2$ shown in figure 1(c) and (d), where the region of high energy contribution covers a much wider range of wall-parallel wavelengths for the $M_{\infty}=4$ case, with the ‘nose’ of the spectrum located at a smaller $(\lambda_{1}/\delta,\lambda_{3}/\delta)$.

The first observation regarding the relative Mach number of unity can be explained by the formation of Mach waves in the relatively supersonic region. The second can be comprehended in terms of the correct scaling required for the compressible boundary layer. These main two observations will be discussed in the remainder of this section.

3.1 Relatively supersonic region and Mach waves

As displayed by the energy contribution in figure 1, the resolvent operator is shown to exhibit low-rank behaviour for the supersonic turbulent boundary layer as well as for the incompressible case. Moarref et al. (2013) showed that significant understanding of the scaling of wall turbulence can be obtained by using the simple rank-1 model. Here, we employ the same rank-1 model by only keeping the most energetic forcing and response directions corresponding to $\sigma_{1}$ and compute the premultiplied one-dimensional energy density of the principal response of $\mathsfbi{H}$ defined as

[TABLE]

In figure 3(a) and (b), we plot the energy density as a function of wave speed $c$ and $x_{2}$ for both the incompressible and supersonic ($M_{\infty}=2$, $\Rey_{\tau}=450$) turbulent boundary layers. Unlike the incompressible case, where the energy density is localised around the mean velocity profile (solid line), the compressible case shows a second region which is displaced from the mean velocity. The overlaid relative sonic line (dashed line), where the velocity profile corresponds to relative streamwise Mach number of unity at each wall-normal location, shows that the displacement of the areas of high energy density is indicative of Mach waves. The relative sonic line $\bar{c}$ is given by solving

[TABLE]

at each $x_{2}$ for $\kappa_{3}=0$, i.e. $\bar{c}=\bar{u}_{1}-\bar{T}^{1/2}/M_{\infty}$, and indicates the minimum streamwise velocity at each $x_{2}$ where a relatively supersonic region exists. From conditionally sampling the energy intensity for $\overline{M}_{\infty}<1$ or $\overline{M}_{\infty}>1$ by evaluating (24) with the integrand of the numerator multiplied by an indicator function for each case (as shown in figure 3(c) and (d)), the two phenomena can clearly be separated, and the second region of high energy density can be attributed entirely to the relatively supersonic region of $\overline{M}_{\infty}>1$ and the existence of modes resembling Mach waves.

In figure 4, we plot a few of the principal streamwise velocity response modes $(u_{1})_{1}$ for both the incompressible and compressible ($M_{\infty}=2,\Rey_{\tau}=450)$ cases at $(\lambda_{1}/\delta,\lambda_{3}/\delta)=(0.01,10)$, which lie within the region $\overline{M}_{\infty}>1$ for the compressible case for a variety of wave speeds. Note that the modes under consideration are essentially two-dimensional, as we need $\lambda_{1}\ll\lambda_{3}$ in order for the relatively supersonic region to exist for a wide range of wall-normal locations. This aspect ratio has not been studied in the past in the context of unforced turbulent boundary layers and the mode shapes that occur here are different from the nominal three-dimensional case.

Examining the principal response modes, we see the modes for the incompressible case (figure 4(a)) are centred at the critical layer $x_{2}^{c}$, where $\bar{u}_{1}(x_{2}^{c})=c$, but the modes corresponding to $\overline{M}_{\infty}>1$ for the compressible case (figure 4(b)) are centred at the sonic layer, $x_{2}^{s}$, where $\overline{M}(x_{2}^{s})=1$. Although not shown, the other velocity, density, and temperature modes, and consequently the pressure modes, are also centred in the sonic layer. These response modes are consistent with Mach waves that propagate towards the freestream and the concept of eddy shocklets (Phillips, 1960; Ffowcs Williams & Maidanik, 1965), where the instantaneous supersonic events cause local shock-like structures in the boundary layer. However, the analysis of the formation of eddy shocklets from superimposed Mach waves requires additional knowledge of phase for each wave parameter and is not considered here.

3.2 Relatively subsonic region and universality of resolvent

modes

In order for the resolvent modes to exhibit universal behaviour for different $\Rey_{\tau}$, the modes must have a narrow footprint such that the support of the modes in the wall-normal direction is localised. This is because, otherwise, resolvent modes are affected by the mean velocity at various wall-normal locations that scale differently with $\Rey_{\tau}$ (Moarref et al., 2013). Moreover, the necessary condition for the existence of geometrically self-similar resolvent modes is the presence of a logarithmic region in the turbulent mean velocity profile.

In order for the resolvent modes to be universal for the supersonic boundary layer, not only do the modes have to be localised in $x_{2}$, but the mean velocity profile must have a scaling law similar to that of the incompressible case such that different regions of the mean profiles collapse for various $M_{\infty}$ and $\Rey_{\tau}$. In compressible flows, viscous heating causes non-uniform mean density and viscosity, which results in a mean velocity profile that no longer satisfies the scaling of its incompressible counterpart. Many attempts have been made to recover the scalings in this regime (Wilson, 1950; Van Driest, 1951; Coles, 1964; Zhang et al., 2012; Trettel & Larsson, 2016, among others), with particular emphasis on the logarithmic region. Most of these attempts have been made by seeking a transformation of $\bar{u}_{1}$ and $x_{2}$ such that the compressible velocity profile maps onto an equivalent incompressible profile. The most recent of these approaches given by Trettel & Larsson (2016) utilises a semi-local scaling in $x_{2}$ and an integrated stress-balance condition, which assumes that the sum of viscous and Reynolds stresses in both the raw and transformed states must be equal, for the scaling of $\bar{u}_{1}$ such that

[TABLE]

Here, the subscript $w$ indicates quantities evaluated at the wall, and $\tau_{w}$ is the wall shear stress. The results of this transformation are illustrated in figure 5(a) and (b), and an improved collapse of the mean velocity profile in the inner and logarithmic region for the various Mach numbers is achieved. The semi-local scaling $x_{2}^{*}$ used here was introduced by Lobb et al. (1955), revisited by Huang et al. (1995) and Coleman et al. (1995) and generalised by Trettel & Larsson (2016). This scaling gives rise to a semi-local Reynolds number (Cebeci & Bradshaw, 2012; Patel et al., 2015) at each wall-normal location

[TABLE]

such that $x_{2}^{*}=(x_{2}/\delta)\Rey_{\tau}^{*}$. While the transformation proposed by Trettel & Larsson (2016) works well for the inner and logarithmic region, the collapse is not as good for the outer region. We find that the best collapse for the outer region is achieved with the transformation

[TABLE]

which is equivalent to scaling the velocity with the semi-local $u_{\tau}^{*}=\sqrt{\tau_{w}/\bar{\rho}}$ instead of $u_{\tau}=\sqrt{\tau_{w}/\bar{\rho}_{w}}$, and the results are given in figure 5(c) and (d). A different transformation for the outer layer is expected since the transformation given in (27) is based on the idea that momentum conservation can be achieved by satisfying the stress balance condition, which only holds in the inner layer of nearly parallel shear flow at reasonable turbulence Mach numbers. Note that despite the better scaling in the outer region, the collapse is not perfect, which is a known issue for low Reynolds number boundary layer flows. Still, a universal mean velocity profile for the inner, logarithmic and outer layer can be achieved by utilising the semi-local scaling.

Given that the transformation of the turbulent supersonic mean velocity profile in (26)–(29) produces a reasonable match to the incompressible profile in the inner, outer and logarithmic regions, the relatively subsonic ($\overline{M}_{\infty}<1$) resolvent modes are expected to have a universal behaviour in both Reynolds and Mach number, as the one-dimensional energy density conditioned to $\overline{M}_{\infty}<1$ in figure 3(c) shows a localisation with respect to $x_{2}$. Moreover, the existence of a logarithmic region in the transformed mean streamwise velocity profile satisfies the necessary condition for the geometrically self-similar modes to be present. Due to the transformation in both $\bar{u}_{1}$ and $x_{2}$, the scaling of the resolvent modes should be with respect to the semi-local variables, $x_{2}^{*}$ and $\Rey_{\tau}^{*}$. The detailed scaling laws for response modes of the inner, outer, and logarithmic region are given in §4.1.

The semi-local scaling also explains the discrepancy between the principal energy contribution in the incompressible boundary layer (figure 1(c)) and the supersonic one (figure 1(d)). For the supersonic case of $M_{\infty}=4$, the semi-local Reynolds number at $x_{2}/\delta=0.2$ is $\Rey_{\tau}^{*}=1020$, which is significantly larger than the Reynolds number of the incompressible case ($\Rey_{\tau}=450$). In figure 6, we instead compare against the results from the incompressible turbulent boundary at $\Rey_{\tau}=1040$ with the mean velocity profile from Schlatter & Örlü (2010) and the premultiplied energy spectra from turbulent channel flow at $\Rey_{\tau}=950$ (Del Alamo et al., 2004), which show a better qualitative comparison between the two principal energy contribution spectra than the comparison given in figure 1.

4 Scaling of the principal singular value and resolvent modes

4.1 Scaling of the principal response mode

As discussed in §3.2, the mean streamwise velocity profile with the semi-local scaling collapses to the incompressible boundary layer profile. This implies that the same scaling used in Moarref et al. (2013) for the incompressible channel flow can be extended to the compressible boundary layer by using the length-scale $x_{2}^{*}$ and Reynolds number $\Rey_{\tau}^{*}$.

Note that $\Rey_{\tau}^{*}$ depends on Mach number through the variation in density and temperature (see figure 7). For an adiabatic wall, the mean temperature profile is given by

[TABLE]

where $r=\Pran^{1/3}$ is the recovery factor (Walz, 1969). Given the definition of the semi-local Reynolds number (28) and the mean equation of state $\bar{p}=\bar{\rho}\bar{T}=1$, we have

[TABLE]

where $C_{w}=S/T_{w}$. When $x_{2}=0$, we have that $\Rey_{\tau}^{*}=\Rey_{\tau}$ as expected. In the limit $x_{2}\rightarrow\infty$, $\bar{T}/\bar{T}_{w}=1/\left(1+r(\gamma-1)M_{\infty}^{2}/2\right)$ and we have

[TABLE]

giving approximately a $M_{\infty}^{3}$ dependence for high Mach numbers.

In order for a fair comparison between the incompressible and compressible response modes, the compressible velocity response modes must be normalised by the kinetic energy content in the response modes due to the orthonormality constraint of the singular vectors and the different norms used for the two cases. We define turbulent kinetic energy and turbulent thermodynamic energy as

[TABLE]

respectively. We normalise the velocity, density and temperature modes such that

[TABLE]

The relationship between $E_{K}$ and $E_{T}$ is discussed later in §4.3. The proposed scaling of the resolvent for the different classes of the compressible boundary layer is summarised in table 2. Note that the combined scaling of the normalised response modes, the singular values and the ratio of $E_{K}$ to $E_{T}$ is of importance when examining the response mode of the resolvent operator.

4.1.1 Outer region

Because the difference between $\Rey_{\tau}$ and $\Rey_{\tau}^{*}$ grows as $x_{2}$ increases (see figure 7), the difference between the two scalings is expected to be most pronounced in the outer region, defined as $(\bar{u}^{\star}_{1,\infty}-\bar{u}_{1}^{\star})\lesssim 6$ with $\kappa_{3}/\kappa_{1}\gtrsim\epsilon\Rey_{\tau}^{*}/\Rey_{\tau,\min}^{*}$, $\epsilon^{2}\approx 10$ (see Appendix B and Moarref et al., 2013, for details). In this region, the $x_{2}$ dependent coefficients of $\mathsfbi{H}$ such as $\bar{u}_{1}-c$ (and $\bar{T}$ and $\bar{\rho}$) are independent of $\Rey_{\tau}$. And in the incompressible case, the balance between the viscous dissipation term and the mean advection term in the resolvent requires scaling of the spanwise coordinate in $\delta$ and the streamwise coordinate with $\delta^{+}$ (Moarref et al., 2013). Thus, the universal modes in the outer region for the incompressible case are given by the wave parameters

[TABLE]

in wall units. Using semi-local variables, we expect the wave parameters for the compressible case to be

[TABLE]

Due to the scaling of the wall-normal distance in outer units as well as the orthonormality condition for the response modes, i.e. $((\boldsymbol{q})_{1},(\boldsymbol{q})_{1})_{E}=1$, the mode height is expected to scale in outer units. In figure 8, we plot the normalised and scaled principal streamwise velocity, density and temperature response modes, with the reference parameters $\Rey_{\tau}^{\text{ref}}=445.5$, $\kappa_{1}^{\text{ref}}\delta=1$, $\kappa_{3}^{\text{ref}}\delta=10$ and $c^{\text{ref}}/\bar{u}_{1,\infty}=0.98$. These reference parameters correspond to the very large scale structures for the incompressible case. The collapse among different Mach numbers is excellent for the semi-local scaling, and the streamwise velocity modes for the supersonic cases are indistinguishable from the incompressible response modes. By comparison, the modes do not collapse for various Mach numbers when scaled in wall units, with the modes becoming narrower in $x_{2}$ and the peak of the modes shifting closer to the wall with increasing Mach number. Although not shown, the wall-normal velocity response modes exhibit a similar collapse when scaled with $\Rey_{\tau}^{*}$ as opposed to $\Rey_{\tau}$, the proposed scaling in Sharma et al. (2017) for the incompressible case. The collapse is not as good for the spanwise velocity response modes (not shown), although the mode shape is still very similar among different Mach numbers. Furthermore, the mode shapes are comparable for the streamwise velocity, density, and temperature. The resemblance of the density and temperature response modes can be explained by the fact that from the equation of state, the two are linked with the pressure fluctuations. Moreover, the similarity between the streamwise velocity and thermal modes reinforces the strong Reynolds analogy, which links the transport of momentum and heat transfer, and it can be concluded that the velocity and temperature profiles are correlated.

4.1.2 Inner region

In the case of the inner region, the deviation of $\Rey_{\tau}$ and $\Rey_{\tau}^{*}$ is not as significant, as can be seen in figure 7, leading to a similar result for both the wall-unit scaling and the semi-local scaling. For the universal modes in the inner region, where the streamwise and spanwise coordinates are given in inner units, we expect the universal wave parameters to be

[TABLE]

for $x_{2}^{\text{ref}*}<30$ and $\kappa_{3}/\kappa_{1}\gtrsim\epsilon$, and analogously defined for $(q_{i})_{1}^{\mathrm{in},+}$. In figure 9, we plot the normalised and scaled principal streamwise velocity, density and temperature response modes for reference Reynolds number, $\Rey_{\tau}^{\text{ref}}=445.5$ and wave parameters $\kappa_{1}^{\text{ref}}\delta=1$, $\kappa_{3}^{\text{ref}}\delta=10$ and $x_{2}^{\text{ref}*}=10$ (and $x_{2}^{\text{ref}+}=10$ for comparison). Here, the scaling of the wall-normal distance is in semi-local (or wall) units. Thus, the orthonormality condition for the velocity and thermodynamic modes requires that the response mode height is expected to scale as $1/\sqrt{\Rey_{\tau}^{*}}$ (or $1/\sqrt{\Rey_{\tau}}$). Again, the collapse of the response modes among different Mach numbers is good for both cases, and the streamwise velocity response modes also collapse with the incompressible case. The wall-normal and spanwise velocity response modes also collapse for the various Mach numbers to the incompressible case but are not shown for brevity. The mode shapes for density and temperature are almost identical for the same reasons discussed above. The temperature and density modes scaled in wall units are not shown, but they are almost identical to the ones given in semi-local units.

4.1.3 Logarithmic region

Finally, while the $\Rey_{\tau}$ for the cases under consideration is too small for a clearly defined logarithmic region, we consider the self-similar response modes. The wave parameters in this region are given by

[TABLE]

for wave speeds $c$ in the logarithmic region ($30/\Rey_{\tau}^{*}<x_{2}/\delta<0.15$) (Tennekes et al., 1972), where $x_{2}^{\text{ref}}$ denotes the critical layer for $c^{\text{ref}}$. The wave parameters in wall units are given analogously for $(q_{i})_{1}^{\mathrm{log},+}$. In figure 10, we plot the normalised and scaled principal streamwise velocity, density and temperature response modes for reference Reynolds number, $\Rey_{\tau}^{\text{ref}}=445.5$, and reference wave parameters $\kappa_{1}^{\text{ref}}\delta=1$, $\kappa_{3}^{\text{ref}}\delta=10$, and $c^{\text{ref}}/\bar{u}_{1,\infty}=0.5$. The orthonormality condition for the velocity and thermodynamic quantities give the response mode height scaling of $\sqrt{x_{2}^{c}/\delta}$. Similar to the outer region, the semi-local scaling gives a better collapse among the supersonic response modes for various Mach numbers. While the agreement with the incompressible case is not perfect, some improvement is made from the use of the semi-local scale. As mentioned earlier, the $\Rey_{\tau}$ for cases under consideration is too low for an actual logarithmic layer and a better collapse is expected for higher Reynolds numbers. Additionally, as in the case of the inner and outer region, the temperature and density modes are identical and are similar to the streamwise velocity mode shape as well.

4.2 Scaling of the principal singular values

In addition to the universality and self-similarity of the response modes, the scaling of the principal singular values is also investigated. The expected scaling of the singular values for the incompressible case is given by $1/\Rey_{\tau}$ in the inner region, $x_{2}^{+2}x_{2}/\delta$ in the logarithmic region and $\Rey_{\tau}^{2}$ in the outer region (Moarref et al., 2013; Sharma et al., 2017). The scaling of the principal singular values can be found by performing a scaling analysis on the resolvent operator $\mathsfbi{H}$ by assessing the Reynolds number dependency of the terms in the linearised operator $\mathsfbi{L}$ (see Appendix B for details). The elements of the resolvent operator matrix, and thus the leading singular value, can be shown to follow $1/\sqrt{\Rey_{\tau}^{*}\Rey_{\tau}^{\star}}$ in the inner region, $x_{2}^{\star}x_{2}^{*}x_{2}/\delta$ in the logarithmic region, and $\Rey_{\tau}^{*}\Rey_{\tau}^{\star}$ in the outer region, where $\Rey_{\tau}^{\star}=\bar{\rho}u_{\tau}\delta/\bar{\mu}$ and $x_{2}^{\star}=\Rey_{\tau}^{\star}x_{2}/\delta$. This is due to the presence of the factor $\bar{T}\bar{\mu}/\Rey$ (or equivalently $\bar{\mu}/(\Rey\bar{\rho})$) in the governing equations combined with the semi-local scaling of $\kappa_{1}$, $\kappa_{3}$ and $x_{2}$. The proposed scaling for the principal singular values is given in figure 11. Note that this scaling will converge to the scaling given in Moarref et al. (2013) in the incompressible limit.

4.3 Scaling of the kinetic and thermodynamic energy ratio

Due to the orthonormality constraint of the resolvent modes, the comparison between the compressible and incompressible resolvent modes in the previous section was for normalised response modes $\widetilde{(q_{i})}_{1}$. However, it is equally important to assess the distribution of energy among the kinetic and thermodynamic variables for the supersonic cases.

In figure 12(a), (b) and (c), we plot the spectra of the ratio of thermodynamic energy and kinetic energy $E_{T}/E_{K}$ for the resolvent modes for $M_{\infty}=2$ and $\Rey_{\tau}=450$ for two different wall-normal heights. For the incompressible case, this ratio would be uniformly zero. Close to the wall, where a relatively supersonic region is present, the thermodynamic energy clearly dominates in the relatively supersonic region as expected. Further away from the wall at $x_{2}/\delta=0.2$ and $x_{2}/\delta=0.5$, the thermodynamic energy still dominates in a smaller region with $\lambda_{3}>\lambda_{1}$ with the strongest amplification centred around $\kappa_{1}\delta\approx 0.5$. This may be explained by the observation from Pirozzoli & Bernardini (2011) that in this region, the thermal streaks are more elongated in the spanwise direction compared to the velocity streaks. It is also consistent with the observation that large-scale pressure-carrying eddies or wavepackets, which are correlated with thermodynamic fluctuations, are more coherent in the spanwise direction for both incompressible (Sillero et al., 2014) and hypersonic (Duan et al., 2016) boundary layers.

We also plot the ratio of turbulent kinetic energy to the sum of the mean-square density and temperature fluctuations obtained from DNS (Bernardini & Pirozzoli, 2011; Pirozzoli & Bernardini, 2011) as a function of $x_{2}$ in figure 12(d). In particular, we plot

[TABLE]

where $rms$ denotes the root-mean-squared fluctuations from DNS. For all wall-normal locations, the ratio increases with Mach number. Moreover, the ratio increases as a function of $x_{2}$ for a fixed Mach number. In order to compare the results from the resolvent modes to DNS, we define the energy ratio of the most energetic mode as

[TABLE]

where $\Phi_{u_{1}u_{1}}$ is the premultiplied streamwise energy spectra for the channel flow at $\Rey_{\tau}=550$ obtained a priori from Del Alamo et al. (2004). The $(E_{T}/E_{K})^{\text{res}}$ given by these wave parameters is plotted in figure 12(d). The agreement between the ratio of kinetic and thermodynamic energy given by the most energetic principal response mode of the resolvent analysis and the DNS is excellent in the logarithmic region. While DNS data are used to identify the wave parameters of the most energetic mode, note that the spectra are for incompressible channel flow rather than the compressible turbulent boundary layer. The discrepancy in the outer region, especially for the higher Mach numbers, may be due to the fact that the difference in the spectra of the boundary layer and channel flow, which was used to choose the wave parameters for the most energetic modes, becomes more pronounced in this region. This may possibly be mitigated by utilising energy spectra from boundary layer simulations. In the inner region, the estimated energy ratio plateaus, deviating from the DNS profile. This could be due to the increased contributions from relatively supersonic region, which is more prevalent in the near-wall region (see figures 1 and 2). Additionally, it is shown in LeHew et al. (2011) that the energetic contribution of structures with convection velocities less than $10u_{\tau}$ is negligible in real turbulent flows, which corresponds to the region where the mismatch is pronounced.

5 Conclusions

We have applied resolvent analysis for the compressible Navier-Stokes equations to a supersonic zero-pressure-gradient turbulent boundary layer. From the low-rank approximation formulated for individual wall-parallel wavenumbers and frequencies, we have identified two distinct regions in the wave parameter space: the relatively supersonic region and the relatively subsonic region.

In the relatively supersonic region, marked by relative Mach number greater than unity, we show that the resolvent modes are centred around the relative sonic line rather than the critical layer and that the majority of the energy is carried by the thermodynamic fluctuations. These response modes are consistent with acoustic Mach waves propagating towards the freestream and the idea of eddy shocklets, where the instantaneous supersonic events cause local shock-like structures in the boundary layer. Additionally, the modes in this region are also shown to follow a modal amplification mechanism. The range of wave parameters corresponding to the relatively supersonic region, where the compressibility effects are concentrated, grows with Mach number, which might be an indicator of why Morkovin’s hypothesis fails for high Mach numbers.

In the relatively subsonic region, we show that the principal response modes are localised around the critical layer corresponding to the mean velocity profile. Furthermore, with the semi-local scaling proposed by Trettel & Larsson (2016), the mean velocity profiles collapse for various Reynolds and Mach numbers. This provides the necessary condition for the resolvent modes to exhibit universality and geometrical self-similarity when scaled with the semi-local Reynolds number and wall-normal distance. We show that the principal response modes are indeed universal and self-similar and that they follow the same scaling laws as the incompressible boundary layer when normalised by the fluctuating kinetic energy. This validates the notion of Morkovin’s hypothesis for the relatively subsonic region on a mode-by-mode basis. Moreover, the velocity modes and the temperature and density modes are qualitatively similar, consistent with the strong Reynolds analogy. We also provide scaling laws for the amplification factor of the principal response mode.

Finally, we show that the energy distribution between the velocity fluctuations and the thermodynamic fluctuations can be predicted from the energy distribution of the most energetic response mode. Coupled with the universality and self-similarity of the resolvent modes in the relatively subsonic region, this has implications in modelling and prediction of high-speed turbulence. As in the incompressible case, the self-similar resolvent modes facilitate analytical developments in the logarithmic region of the boundary layer. Additionally, this allows prediction tools developed for resolvent analysis of incompressible fluids to be applied to supersonic boundary layers.

The results show that the main difference between the compressible Navier-Stokes equations and the incompressible equations are due to a different scaling law originating from density variations in the wall-normal direction and the acoustic contribution in the relatively supersonic region. The full nonlinear closure then propagates the deviation through triadic interactions, which results in variations in the mean velocity profile. Future efforts will be focused on studying the effect of the feedback loop by incorporating limited self-interactions to estimate the forcing term as in Rosenberg et al. (2019) for the incompressible Navier-Stokes equations. Also, further efforts are necessary to study the effects of higher Mach numbers and different wall boundary conditions such as cooled walls.

The authors acknowledge support from the Air Force Office of Scientific Research grant FA9550-16-1-0232. The authors also thank Dr. Adrián Lozano-Durán and Dr. Minjeong Cho for their insightful comments on the manuscript.

Appendix A Modal amplification mechanism in the relatively supersonic

region

The spectrum of the linear operator $\mathsfbi{L}$ is obtained from solving the eigenvalue problem $\mathsfbi{L}\boldsymbol{q}=\mathrm{i}\omega\boldsymbol{q}$. Here, we plot the results for the $M_{\infty}=4$ case in terms of the complex wave speed $c=\Re(c)+\mathrm{i}\Im(c)$ in figure 13 for two different grid resolutions in $x_{2}$, $N_{2}=401$ and $N_{2}=801$. We also tested finite domain heights, compared to the semi-infinite domain, and the resulting spectra are robust to this change (not shown). We see that besides the discrete modes, the spectra consist of continuous branches that are represented as discrete points which either converge towards a vertical line, as expected for the viscous modes in both incompressible and compressible boundary layers, or form a horizontal line close to the real axis for the acoustic modes. For the spectrum of $(\kappa_{1}\delta,\kappa_{3}\delta)=(0.2,2)$, relative Mach number, $\overline{M}_{\infty}=0.39$, is subsonic, and thus relatively supersonic modes are absent in this case. However, the spectrum for $(\kappa_{1}\delta,\kappa_{3}\delta)=(1,2)$ (figure 13(b)), where the relatively supersonic region is present, shows an additional feature indicative of a continuous branch of ‘acoustic’ eigenmodes close to $\Im(c)=0$. In the relatively supersonic region, the resolvent norm shows peaks in both the leading singular value and the energy contained in the principal resolvent mode due to spectral amplification, i.e., the resonant effects from the operator becoming close to normal as the wave speed approaches the discrete acoustic eigenvalues. This shows that the irregular patterns in figure 1(b) are the consequence of discretisation, and that the resonant effects will be universally present in the continuous case.

For a given wave speed $c$, we can plot the distance $d_{\Lambda}$, under the norm given in (13), from $\omega=c\kappa_{1}$ to the closest discrete acoustic eigenvalue as a function of wavenumbers $\kappa_{1}$ and $\kappa_{3}$. In figure 14, the inverse of the minimum eigenvalue distance is shown to correlate well with the leading singular value, which leads to correlation with the principal energy contribution $\sigma_{1}^{2}/(\sum_{j}\sigma_{j}^{2})$. As mentioned before, this is due to the resonant amplification of the resolvent operator through modal amplification mechanisms. The change in the discrete location of the acoustic eigenvalues as a function of $\kappa_{1}$ and $\kappa_{3}$ explains the irregular patterns in figure 1(b) in the relatively supersonic region. However, the patterns are a purely numerical phenomenon and are indicative of the fact that the resolvent operator does not exist.

Appendix B Reynolds number scaling of the principal singular value

In the inner region where the collapse of the mean velocity profile is achieved through the semi-local scaling, the coordinates scale as

[TABLE]

We then analyse the scaling of the linear operator $\mathsfbi{H}$ by assessing the Reynolds number dependency of the terms in the linearised operator $\mathsfbi{L}$. In this process, we assume that the Mach number plays a secondary role, since we know $\Rey_{\tau}^{*}$ (or $\Rey_{\tau}^{\star}$) reflects the Mach number dependency as seen in §4.1. Doing so, the linear operator $\mathsfbi{L}$ approximately scales as

[TABLE]

As $\Rey_{\tau}^{*}\approx\Rey_{\tau}^{\star}$ in this region, the scaling can be seen as $\mathsfbi{L}\sim\sqrt{\Rey_{\tau}^{*}\Rey_{\tau}^{\star}}$, giving $\mathsfbi{H}=(\mathrm{i}\omega\mathsfbi{I}+\mathsfbi{L})^{-1}\sim 1/\sqrt{\Rey_{\tau}^{*}\Rey_{\tau}^{\star}}$. Thus, the leading singular value is expected to be proportional to $1/\sqrt{\Rey_{\tau}^{*}\Rey_{\tau}^{\star}}$.

In the outer region, the coordinates scale as

[TABLE]

Additionally, if we assume that $\kappa_{3}^{2}$ dominates $\kappa_{1}^{2}$ for all values of $\Rey_{\tau}^{*}$, such that for the streamwise wavenumber in the outer coordinate given by $\kappa_{1}^{-}=\Rey_{\tau}^{*}\kappa_{1}$ we have $\kappa_{3}/\kappa_{1}^{-}>\epsilon/\Rey_{\tau}^{*}$, the linear operator $\mathsfbi{L}$ scales as

[TABLE]

The resolvent $\mathsfbi{H}$ then scales as $\Rey_{\tau}^{*}\Rey_{\tau}^{\star}$, which gives the scaling for the leading singular value $\sigma_{1}$.

Finally, for the logarithmic region, the differential operators are scaled as

[TABLE]

We also assume that the spanwise coordinate dominates the spanwise coordinate, i.e. $(\kappa_{3}/\kappa_{1})>\epsilon$ with a conservative estimation of $\epsilon\approx\sqrt{10}$, and arrive at

[TABLE]

And a similar analysis as the inner and outer regions reveals that the singular values scale with $x_{2}^{*}x_{2}^{\star}x_{2}$.