Linear response analysis of supersonic turbulent channel flows with a large parameter space

Abstract

In this work, the linear responses of turbulent mean flow to both harmonic and stochastic forcing are investigated for supersonic channel flow. Well-established universal relations are utilized to obtain efficiently the mean profiles with a large parameter space, with the bulk Mach number up to 5 and the friction Reynolds number up to $10^4$, so a systematic parameter study is feasible. The most amplified structure takes the form of streamwise velocity and temperature streaks forced optimally by the streamwise vortices. The outer peak of the pre-multiplied energy amplification corresponds to the large-scale motion, whose spanwise wavelength ($\lambda _z^+$) is very insensitive to compressibility effects. In contrast, the classic inner peak representing small-scale near-wall motions disappears for the stochastic response with increasing Mach number. Meanwhile, the small-scale motions become much less coherent. A decomposition of the forcing identifies different effects of the incompressible counterpart and the thermodynamic components. Wall-cooling effects, arising with high Mach number, increase the spacing of the most amplified near-wall streaks; the spacing becomes nearly invariant with Mach number if expressed in semi-local units. Meanwhile, the coherence of stochastic response with $\lambda _z^+>90$ is enhanced, but on the other hand, with $\lambda _z^+<90$ it is decreased. The geometrical self-similarity of the response in the mid-$\lambda _z$ range is still roughly satisfied, insensitive to Mach number. Finally, theoretical analyses of the perturbation equations are presented to help with understanding the scaling of energy amplification.

1. Introduction

Understanding and modelling high-speed turbulent flows are of fundamental importance in aerodynamic applications (Spina, Smits & Robinson 1994). The turbulent boundary layer influences surface drag and heat transfer dramatically. Thereby, developing accurate prediction models has long been the pursuit for reliable vehicle design. Compared with the incompressible case, our knowledge of compressible turbulent flows is rather limited, for two main reasons. First, the compressible flow suffers from coupling effects of various factors, which, in addition to the Reynolds number, include the Mach number, surface heat transfer, shock waves and high-enthalpy effects (Coleman, Kim & Moser 1995; Duan & Martín 2011; Fu et al. 2021; Fu, Bose & Moin 2022). Second, the available experimental and numerical databases are rather limited because of the restricted facilities and computational resources.

With the advancement of modern supercomputers, direct numerical simulations (DNS) (see e.g. Moin & Mahesh 1998) are increasingly powerful for the research of compressible turbulent flows. For example, recent DNS data from Zhang, Duan & Choudhari (2018) and Huang, Duan & Choudhari (2022) cover a freestream Mach number ${{Ma}}_\infty$ range of up to 14 for flat-plate boundary layers with adiabatic and cold walls. In addition, the maximum friction Reynolds number ${{Re}}_\tau$ reaches as high as 4000 for the boundary layer (Pirozzoli & Bernardini 2013, with ${{Ma}}_\infty =2$) and 2000 for the channel flow (Yao & Hussain 2020, with a bulk Mach number ${{Ma}}_b=1.5$). Utilizing these DNS data, the effects of Mach number, wall cooling, pressure gradient, high-enthalpy, and so on, on turbulent scalings and structures, can be isolated and studied. For example, Duan, Beekman & Martín (2011) found that many of the scaling relations, though derived based on the Morkovin hypothesis, remained valid for boundary-layer flows with ${{Ma}}_\infty$ up to 12. In addition, recent work of Di Renzo, Fu & Urzay (2020), Di Renzo & Urzay (2021) and Passiatore et al. (2022) demonstrated the applicability of these scaling relations on high-enthalpy boundary layers with thermochemical non-equilibrium effects.

Despite much progress, the parameter ranges of compressible turbulent DNS databases are still quite limited, especially for channel flows. The highest ${{Ma}}_b$ reported for a channel flow is approximately 4 (Trettel & Larsson 2016). The ${{Re}}_\tau$ of that case is near 1000, but the centreline friction Reynolds number in semi-local units, ${{Re}}_{\tau,c}^*$, is only approximately 200. As mentioned above, the highest ${{Re}}_{\tau,c}^*$ of a channel flow case is reported by Yao & Hussain (2020), while the ${{Ma}}_b$ of that case is limited to 1.5. Any effort to extend current parameter ranges is extremely resource-demanding. For example, the grid number in the case of Yao & Hussain (2020) is over 8 billion. If ${{Ma}}_b$ is increased to 3 at the same ${{Re}}_\tau$, then the required grid number is evaluated to be dauntingly over 20 billion. The substantial computational cost prohibits a systematic parameter study of various influencing factors.

Compared to instantaneous DNS data, the calculation of only turbulent mean flow can be substantially cheaper, such as through the Reynolds-averaged Navier–Stokes (RANS) or large-eddy simulations (LES) approaches. With only mean profiles in hand, one can still investigate some important turbulence characteristics using linear models, which will be introduced at length in § 1.2. Therefore, the combination of cheap reliable mean-flow calculation and the mean-flow-based analysis can be an effective way to explore the turbulence at high Mach and Reynolds numbers beyond the current capability of DNS. As a demonstration, Cossu, Pujals & Depardon (2009) and Hwang & Cossu (2010b) used a curve-fitted incompressible mean flow ${{Re}}_\tau >10^{4}$, and analysed the transient growth and the responses to forcing for boundary-layer and channel flows. Pickering et al. (2021) proposed a resolvent-based model with the mean flow from LES and captured successfully the characteristics of turbulent coherent structures in the shear flow. This type of method combination will be deployed in this work for canonical compressible channel flows. The focus is on the linear response of mean profiles and coherent structures, as an effort to extend the present knowledge frontier on the turbulence characteristics in high Mach and Reynolds number regimes beyond the scope of existing DNS databases. More specifically, the first aim is to utilize well-established universal relations to obtain turbulent mean profiles with a large parameter space. The second is to develop the linear response solvers to harmonic and stochastic forcing for compressible turbulent flows. The third is to investigate the response characteristics with varying Mach and Reynolds numbers, and provide supporting theoretical analysis.

1.1. Universal relations for compressible wall-bounded turbulent flow

It is well known that incompressible wall-bounded turbulent flows exhibit nearly universal mean streamwise velocity profiles versus the wall-normal coordinate when normalized with the wall viscous units. Therefore, a universal curve-fitted expression is available, and a prevailing one for the channel flow is from Cess (1958). This universal curve is also pursued in compressible turbulent flows for the mean velocity and temperature. For instance, the velocity transformation is designed to transform the compressible streamwise velocity profile to the incompressible counterpart. The pioneering work was accomplished by van Driest (1951), and subsequent improved versions were proposed by Zhang et al. (2012), Trettel & Larsson (2016), Volpiani et al. (2020), among others. Very recently, the total-stress-based transformation proposed by Griffin, Fu & Moin (2021b) collapses successfully nearly all available DNS data, with and without heat transfer, for the boundary-layer, channel and pipe flows. The relative errors are typically less than 3 % when ${{Re}}_\tau > 200$. Note that the performance of this transformation for non-canonical turbulent boundary-layer flows has been reported further in Bai, Griffin & Fu (2022). In terms of the temperature profile, Crocco (1932) and Busemann (1931) derived a temperature–velocity relation showing that the mean temperature was nearly a quadratic function of the mean streamwise velocity. Later, a less-restricted relation was proposed by Walz (1969) through introducing the recovery temperature. Furthermore, Duan & Martín (2011) made an improvement by introducing the mean enthalpy and provided more universal fitting coefficients for the quadratic expression based on the data of wide parameter ranges. This relation is shown to be sufficiently accurate for both boundary-layer and channel flows, even with high-enthalpy effects (Passiatore et al. 2022). Also, Zhang et al. (2014) proposed a general Reynolds analogy for compressible wall-bounded turbulent flows by introducing a general recovery factor.

The success of these universal relations makes it possible to obtain the mean profiles by solving an inverse problem under proper flow conditions in the outer region. This idea was adopted recently by Griffin, Fu & Moin (2022) to realize a highly accurate wall-modelled LES framework and to estimate the grid points required for DNS and LES (Griffin, Fu & Moin 2021a). For the channel flow with symmetric isothermal boundaries, Song et al. (2022) proposed an empirical relation between the centreline mean temperature and velocity. This can serve as an outer boundary condition to close the inverse problem, such that the mean-flow profiles can be obtained by solving a simple ordinary differential equation (ODE).

1.2. Linear analysis of turbulent wall-bounded flows

The linear models for analysing turbulent flows have received increasing interest, especially over the last decade (see McKeon 2017). This method of analysis shares some fundamental procedures with linear stability theory, and borrows some critical thoughts from control theory (Taira et al. 2020; Jovanović 2021). The starting point is the linearized Navier–Stokes equation around a mean state for the Fourier-decomposed perturbation. The model equation can be written as

(1.1)\begin{equation} \mathcal{L}_{mn}\,{\hat{\boldsymbol{q}}}_{mn} = {\hat{\boldsymbol{f}}}_{mn}, \end{equation}

where $\mathcal {L}$ is a linear operator related to the mean flow as well as the temporal and spatial wavenumbers, the subscripts $m$ and $n$ are the Fourier indexes, and ${\hat {\boldsymbol {q}}}$ denotes the perturbation shape function. The nonlinear term ${\hat {\boldsymbol {f}}}$ was ignored in early works, and all the eigenmodes of $\mathcal {L}$ were found to be stable for the incompressible turbulent profile (Malkus 1956). Instead of discarding it totally, Reynolds & Hussain (1972) modelled the nonlinear term into the linear operator by using the Boussinesq approximation and the eddy viscosity $\mu _t$. These eddy-viscosity-enhanced linear models have been adopted widely ever since. Better behaviour of artificially induced waves was reproduced in their results by considering $\mu _t$, but the mean profile was still stable with ${{Re}}_\tau \sim 1000$.

As no unstable modes were found, linear model research was quiet for a while until two sets of mathematical tools were applied successfully. The first is the non-modal instability theory with pioneering work from Trefethen et al. (1993). Specifically, as the eigenvectors of $\mathcal {L}$ are not orthogonal to each other, ${\hat {\boldsymbol {q}}}$ can experience strong transient growth even though each mode is asymptotically stable. This was identified by Butler & Farrell (1993) for incompressible channel flows. Furthermore, del Álamo & Jiménez (2006) (also Pujals et al. 2009) found that by using the eddy-viscosity-enhanced model, two peaks of the transient growth with different spanwise wavenumbers were well related to the near-wall motions and the outer-layer large-scale motions, respectively. The second methodology originates from control theory (Schmid & Henningson 2001; Kim & Bewley 2006; Zare, Georgiou & Jovanović 2020), and is called the input–output or resolvent analysis. By rewriting (1.1) as ${\hat {\boldsymbol {q}}}_{mn} = \mathcal {L}_{mn}^{-1}\,{\hat {\boldsymbol {f}}}_{mn}$ (this is an oversimplified interpretation and the theory details are presented in § 3), ${\hat {\boldsymbol {f}}}$ is regarded as the input or forcing, and ${\hat {\boldsymbol {q}}}_{mn}$ is the output or response. Their relation is established through the transfer matrix $\mathcal {L}_{mn}^{-1}$, known as the resolvent. Based on artificially selected forcing, the response characteristics of turbulent mean flow can be obtained through singular-value and Karhunen–Loève (KL) decompositions. Hwang & Cossu (2010a,b) analysed the linear responses to both harmonic and stochastic forcing for incompressible Couette and channel flows. Similar to the transient growth, the two most amplified structures are the near-wall and large-scale motions of infinite streamwise wavelength. As the most amplified modes can occupy a large portion of the response energy, they represent the characteristic flow structures, which can be deployed for reduced-order modelling. For example, through analysing the most amplified response, McKeon & Sharma (2010) emphasized the vital role of the critical layer for incompressible flows, where the phase velocity of the response equalled the mean velocity. The viscosity effects are significant in the critical layer, so their influence domain can extend well beyond the immediate vicinity of the wall. In the follow-up work, Moarref et al. (2013) designed a low-rank model based on the principal response mode. Different classes of scalings for the streamwise turbulence intensity were found.

The results of the above non-modal growth and resolvent analyses are also closely related. The strong resemblance was demonstrated among the amplified structures of transient growth, optimal harmonic and stochastic responses (Hwang & Cossu 2010b). Meanwhile, the linearly amplified modes all exhibit self-similarity in the mid-wavelength range related to the logarithmic region (del Álamo & Jiménez 2006; Pujals et al. 2009; Moarref et al. 2013; McKeon 2019; Vadarevu et al. 2019), which provides mathematical evidence to the well-known attached-eddy model (AEM) initially proposed by Townsend (1976) (also see the recent review of Marusic & Monty 2019).

In most of the research mentioned above, the nonlinear forcing term ${\hat {\boldsymbol {f}}}$ is modelled ideally, instead of an extraction from real flow data. Specifically, ${\hat {\boldsymbol {f}}}$ is either assumed to be a stochastic signal or taken in an optimal form allowing for the strongest perturbation growth. Consequently, the results using the artificially designed forcing are more qualitative than quantitative, compared with the analysis based on DNS data (Zare, Jovanović & Georgiou 2017). The ‘colour’ of the forcing can lead to different weights of the amplified modes, potentially changing the relative importance of different modes. The weights of modes can be obtained by partly matching the DNS data. The resulting energy spectra and other statistics from the resolvent analysis are well consistent with those from DNS (Moarref et al. 2014; McMullen, Rosenberg & McKeon 2020). As another way of improvement, different optimization processes were designed by Zare et al. (2017) and Hwang & Eckhardt (2020) to obtain the cospectrum of the forcing signal, which is closer to realistic flows. To gain a complete understanding of the forcing statistics for incompressible turbulence, considerable efforts have been made recently to analyse the forcing based on the DNS data (Amaral et al. 2021; Morra et al. 2021; Nogueira et al. 2021). It is revealed that the forcing exhibits a relatively high level of spatial–temporal coherence. The low-rank nature of the forcing will benefit the construction of more reliable reduced-order models of turbulent flow. In addition, Morra et al. (2021) clarified why the introduction of $\mu _t$ in the linear operator leads to an improved response prediction compared to that without $\mu _t$. That is because the weights or the combination coefficients of the principal resolvent modes are closer to those from DNS. In this sense, the usage of $\mu _t$ is also a partial modelling of the forcing colour, consistent with its low-rank nature.

The above results are mainly for incompressible turbulence. The extension and application of linear models to compressible turbulent flows are still limited. Alizard et al. (2015) first presented an analysis of linear transient growth on the turbulent boundary layers with ${{Ma}}_\infty$ up to 4. The inner- and outer-layer modes, analogous to the incompressible ones, are identified, and the weak effects of compressibility are revealed on these characteristic structures. In addition to the similarities, an emerging problem is the dynamics of the suggested streaky temperature perturbation and its coupling with the velocity components. Recently, Bae, Dawson & McKeon (2020a,b) applied the resolvent formulation to supersonic turbulent boundary layers with both adiabatic and cold walls. Their results highlight distinct features below and above the relative sonic line of the response. Moreover, there are implications for the geometrically self-similar structure of response modes and the analogous rules for the incompressible counterpart. Dawson & McKeon (2020) investigated the same cases with emphasis on predicting the shapes of the resolvent mode at different Mach numbers. Up to now, to the authors’ knowledge, the analysis framework of stochastic forcing and response has not been applied to supersonic turbulent flows, although there were successful attempts on the stability and receptivity of compressible laminar flows (Alizard et al. 2022; Madhusudanan & McKeon 2022), which demonstrates the feasibility and value of this linear framework in interpreting the physics of compressible flows. Therefore, a systematic study on the linear responses of supersonic turbulent flows to harmonic and stochastic forcing is required, which is the main focus of the present work.

1.3. Streaky motions and compressibility effects

Streamwise streaks are the most common coherent structures across different layers in wall-bounded turbulent flows. Their behaviours in incompressible flows have been studied in great detail, and are introduced briefly below (see e.g. Marusic et al. 2010b; Jiménez 2013, for more details). The streaks are narrow regions elongated in the streamwise direction, containing mainly the fluctuations of streamwise velocity. In the near-wall region, their mean streamwise length and spanwise spacing are approximately 1000 and 100 wall units in the buffer layer, first measured in experiments (Kline et al. 1967; Smith & Metzler 1983). Furthermore, they are dominant in carrying turbulence energy. The near-wall streaks are also known to be able to self-sustain along with the streamwise vortices, where the linear mechanisms of transient growth and instability play an essential role (Hamilton, Kim & Waleffe 1995; Schoppa & Hussain 2002). The streaky motions also exist at larger scales in the outer region, and are more often termed the large-scale or very-large-scale motions (Kim & Adrian 1999; del Álamo & Jiménez 2003). For the channel flow, the peak of the pre-multiplied spectra of streamwise velocity corresponds to characteristic streamwise and spanwise scales $5h\unicode{x2013}6h$ and $h\unicode{x2013}2h$, respectively (del Álamo & Jiménez 2003; Lee & Moser 2015), where $h$ is the channel half-height. These large-scale motions are extremely energetic. Meanwhile, they are least affected by the presence of the wall; on the other hand, they have a modulation effect on the near-wall turbulence intensity (Hutchins & Marusic 2007; Marusic, Mathis & Hutchins 2010a).

In compressible flows, additional temperature and density streaks are present due to the wall-normal gradients of the mean flow (see e.g. Coleman et al. 1995). Superficial similarities of near-wall streaks were demonstrated in boundary-layer flows for the low, moderate and high Mach number (up to 12) cases (Duan et al. 2011). Meanwhile, two crucial features of compressibility effects on the near-wall motion are noted. The first is the variation of the spacing of streaks, on which different trends were reported. For example, wall cooling is shown to increase the streak spacing (Coleman et al. 1995; Duan, Beekman & Martin 2010), while the increase of ${{Ma}}$ has a decreasing effect with an adiabatic wall condition (Duan et al. 2011). Further analysis shows that the spacing variation is consistent with the change of ${{Re}}_\tau ^*$. Therefore, the spacing experiences much less variation in semi-local units (Morinishi, Tamano & Nakabayashi 2004; Patel, Boersma & Pecnik 2016; Yao & Hussain 2020). The second feature is that wall cooling enhances the coherence of the near-wall motions. Meanwhile, the vortical structures become less chaotic (Coleman et al. 1995; Morinishi et al. 2004; Duan et al. 2010; Zhang et al. 2022). In comparison, the large-scale motions in the outer layer are less influenced by the compressibility effects from both experimental and numerical observations (Spina et al. 1994; Pirozzoli, Bernardini & Grasso 2008; Williams et al. 2018; Cheng & Fu 2022a).

The linear models of turbulence have long been used to investigate streaky motions. They are especially suitable for studying the structure, coherence and instability of streaks (see e.g. Schoppa & Hussain 2002; del Álamo & Jiménez 2006; Hwang & Cossu 2010b; McKeon & Sharma 2010). Nevertheless, as introduced in § 1.2, these linear models have rarely been deployed for compressible turbulent flows. The combination of the ODE mean-flow solver and the linear response analysis developed in this work enables a systematic study of the compressibility effects on the streaky motions in the channel flow. The parameter ranges considered are ${{Ma}}_b$ from 0 to 5, and ${{Re}}_\tau$ up to $10^4$, both of which are far beyond the current DNS database.

The remainder of the article is organized as follows. Section 2 describes the governing equations and the approaches to calculating the mean turbulent flow. Section 3 provides the formulation of the linear response analysis to both harmonic and stochastic forcing. The main results are discussed in § 4, and some theoretical proofs are provided in § 5. Finally, the work is summarized in § 6.

2. Governing equations and mean-flow calculations

2.1. Compressible Navier–Stokes equations

We consider a canonical compressible turbulent channel flow with symmetric isothermal boundaries. The flow illustration and coordinate set-up are shown in figure 1. With the assumption of a calorically perfect gas, the non-dimensional Navier–Stokes equations are written as

(2.1a)$$\begin{gather} \frac{\partial \rho}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\rho {\boldsymbol{u}}\right) = 0, \end{gather}$$

(2.1b)$$\begin{gather}\rho \left(\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}\right) ={-}\boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot}\left[\frac{\mu}{{{Re}}_b} \left( {\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} {{\boldsymbol{u}}^{\rm T}}} \right) \right] - \frac{2}{3}\,\boldsymbol{\nabla} \left( \frac{\mu}{{{Re}}_b}\,{\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}} \right), \end{gather}$$

(2.1c)\begin{gather} \rho \left(\frac{\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} T \right) - {{Ec}} \left(\frac{\partial p}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} p \right) = \frac{{{Ec}}\,\mu}{{{Re}}_b} \left[ \boldsymbol{\nabla} \boldsymbol{u} \boldsymbol{:} \left( {\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} {{\boldsymbol{u}}^{\rm T}}} \right) - \frac{2}{3} \left( {\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}} \right)^2 \right]\nonumber\\ + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\kappa}{{{Re}}_b}\, \boldsymbol{\nabla} T\right), \end{gather}

(2.1d)\begin{gather} p = \frac{\rho T}{\gamma {{Ma}}_b^2}, \end{gather}

where $\rho$, $\boldsymbol {u}=[u,v,w]^{\rm T}$, $p$ and $T$ are the fluid density, velocity, pressure and temperature, respectively, and $\mu$ and $\kappa$ are the viscosity and conductivity. For deriving (2.1), the following non-dimensionalizations are employed:

(2.2)\begin{equation} \left.\begin{gathered} {\boldsymbol{x}} = \dfrac{{\boldsymbol{x}}^d}{h^d}, \quad \boldsymbol{u} = \dfrac{\boldsymbol{u}^d}{U_b^d}, \quad t = \dfrac{t^d U_b^d}{h^d}, \quad \rho = \dfrac{\rho^d}{\rho_b^d}, \quad p = \dfrac{p^d}{\rho_b^d U_b^{d2}}, \quad T = \dfrac{T^d}{T_w^d}, \\ \mu = \dfrac{\mu^d}{\mu_w^d}, \quad \kappa = \dfrac{\kappa^d}{\mu_w^d c_p^d}. \end{gathered}\right\} \end{equation}

The resulting Mach, Reynolds, Prandtl and Eckert numbers are defined as

(2.3a–d)\begin{align} {{Ma}}_b = \dfrac{U_b^d}{a_w^d}, \quad {{Re}}_b = \dfrac{\rho_b^d U_b^d h^d}{\mu_w^d}, \quad {{Pr}} = \dfrac{\mu^d c_p^d}{\kappa^d}=0.72, \quad {{Ec}} = \dfrac{U_b^{d2}}{c_p^d T_w^d} = (\gamma-1)\,{{Ma}}_b^2, \end{align}

where the superscript $d$ stands for dimensional quantities, and the subscript $w$ denotes quantities at the wall ($y=\pm 1$). Also, $h$ is the channel half-height, and $c_p$ is the isobaric specific heat. The speed of sound is $a^d=(\gamma R^d T^d)^{1/2}$, with $R$ the gas constant, and the specific heat ratio $\gamma =1.4$. The bulk density $\rho _b^d$ and bulk velocity $U_b^d$ are selected to meet the requirement

(2.4)\begin{equation} \int_{{-}1}^1 {\bar{\rho} \, {\mathrm{d}} y} = \int_{{-}1}^1 {\bar{\rho} \bar{u} \, {\mathrm{d}} y} = 2, \end{equation}

where the overbar denotes a mean variable. The viscosity is obtained through Sutherland's law, where the fitting constant is 110.4 K and the reference temperature is 293.15 K.

Figure 1. Schematic of the coordinate set-up for the channel flow.

In wall-bounded turbulent flows, it is common to express the variables in wall viscous units with a superscript $+$ as

(2.5a–e)\begin{equation} {\boldsymbol{x}}^+ = \dfrac{{\boldsymbol{x}}^d}{\delta_\nu^d}, \quad \boldsymbol{u}^+ = \dfrac{\boldsymbol{u}^d}{u_\tau^d}, \quad \rho^+ = \dfrac{\rho^d}{\rho_w^d}, \quad T^+ = \dfrac{T^d}{T_w^d} = T, \quad \mu^+ = \dfrac{\mu^d}{\mu_w^d} = \mu, \end{equation}

where $u_\tau ^d = \sqrt {\tau _w^d/\rho _w^d}$ is the friction velocity, $\tau _w$ is the wall shear, and $\delta _\nu ^d=\mu _w^d/(\rho _w^d u_\tau ^d)$ is the length unit. The corresponding friction Reynolds number is ${{Re}}_\tau =h^d/\delta _v^d$. Furthermore, Huang, Coleman & Bradshaw (1995) suggested a set of semi-local units for a better collapse of the turbulence statistics with the incompressible counterparts. The velocity and length units are based on the local mean density and viscosity as $u_\tau ^{*d}=\sqrt {\tau _w^d/\bar {\rho }^d}$ and $\delta _\nu ^{*d}=\bar {\mu }^d/(\bar {\rho }^d u_\tau ^{*d})$, hence $y^*=y_d/\delta _v^{*d}$ and ${{Re}}_\tau ^*=h^d/\delta _v^{*d}$. Hereinafter, the quantities using semi-local units are expressed with a superscript *.

The basic variable of (2.1) is $\boldsymbol {q}=[\rho,u,v,w,T]^{\textrm {T}}$, which is decomposed into a mean part and a perturbed part as

(2.6)\begin{equation} \boldsymbol{q} \left(x,y,z,t\right) = \bar{\boldsymbol{q}} \left(y\right) + \tilde{\boldsymbol{q}} \left(x,y,z,t\right), \end{equation}

where the mean component in the channel flow is taken as ${\bar {\boldsymbol {q}}}=[\bar {\rho }(y), \bar {U}(y),0,0,\bar {T}(y)]^{\textrm {T}}$. The calculation of ${\bar {\boldsymbol {q}}}$ is described in § 2.2, and the operations with $\tilde {\boldsymbol {q}}$ will be elucidated in § 3.

2.2. Calculation of turbulent mean flow

For the incompressible channel flow, a semi-empirical expression proposed by Cess (1958) is used widely to describe the mean velocity profile $\bar {U}_{inc}^+$ at different ${{Re}}_\tau$ as

(2.7) \begin{equation} \left. \begin{array}{l@{}} \displaystyle \bar{U}_{inc}^+ (y^+;{{Re}}_\tau) = \int_{0}^{y^+} {\dfrac{1 - y^+{/}{{Re}}_\tau}{1 + \nu_t}\,{\mathrm{d}} y^+} ,\\ \displaystyle \nu_t = \dfrac{1}{2} \sqrt{1 + f_{\nu_t}^2} - \dfrac{1}{2}, \quad f_{\nu_t} = \dfrac{\kappa_c\,{{Re}}_\tau}{3}\,(1-y^2) (1+2y^2) \left[1 - \exp \dfrac{{{Re}}_\tau\,(|y|-1)}{A_{\nu_t}} \right] , \end{array} \right\} \end{equation}

where $\nu _t=\mu _t/\bar {\rho }$ is the kinematic eddy viscosity, and $\kappa _c$ and $A_{\nu _t}$ are two fitting parameters. Based on the mean profile of ${{Re}}_\tau =2003$, the values recommended by del Álamo & Jiménez (2006) are $\kappa _c=0.426$ and $A_{\nu _t}=25.4$, respectively. By including the DNS dataset of ${{Re}}_\tau$ up to 5200, the values adopted in the present work are $\kappa _c=0.418$ and $A_{\nu _t}=25.1$ for all ${{Re}}_\tau$ (see more details in Appendix A). The analytical equation (2.7) has been utilized widely in the linear and resolvent analyses of incompressible flows over wide parameter ranges (see e.g. Butler & Farrell 1993; del Álamo & Jiménez 2006; Moarref et al. 2013). However, no such expressions are currently available for the compressible channel flow.

As introduced in § 1.1, there are already adequate universal relations for solving the inverse problem. Consequently, the mean profile can be obtained quickly by solving an ODE. The algorithm involves four universal relations: (1) the incompressible velocity profile as shown in (2.7); (2) the compressible velocity transformation from Trettel & Larsson (2016) or Griffin et al. (2021b); (3) the algebraic relation from Duan & Martín (2011) between the temperature and velocity distributions; (4) the outer boundary condition from Song et al. (2022) connecting the centreline temperature and velocity. Only two parameters are required as the inputs, i.e. ${{Ma}}_b$ and ${{Re}}_\tau$ (or equivalently, ${{Re}}_b$, ${{Re}}_\tau ^*$ or others), and the solution takes less than one second on a laptop. The algorithm and the solver verification are elaborated in Appendix A.

With this efficient solver in hand, the mean turbulent flow of wide ranges of ${{Re}}_\tau$ and ${{Ma}}_b$ can be obtained. As a demonstration, figure 2 depicts the contours of ${{Re}}_{\tau,c}^*$ and $\bar {T}_c$ as functions of ${{Re}}_\tau$ and ${{Ma}}_b$, with subscript $c$ denoting the channel centre. The black contour with ${{Re}}_{\tau,c}^*=200$ denotes the lower available boundary of the solver. As can be seen, at a fixed ${{Re}}_\tau$, ${{Re}}_{\tau,c}^*$ drops rapidly with the increase of ${{Ma}}_b$, which heavily restricts the range of ${{Re}}_{\tau,c}^*$ of the DNS database. The present solver enables us to explore the flow characteristics far beyond the parameter range of the DNS database. Although its accuracy cannot be directly proved rigorously, the following three facts further strengthen our confidence in high ${{Ma}}_b$ and ${{Re}}_\tau$ regimes.

1. (i) The velocity transformation (see (A2)) has been verified for flows with ${{Ma}}_b$ up to 14, and tends to be more accurate at higher ${{Re}}_{\tau,c}^*$ (Trettel & Larsson 2016; Griffin et al. 2021b).

2. (ii) The accuracy of the velocity–temperature relation (see (A3)) has been proved for flows with ${{Ma}}_b$ up to 14 (Duan & Martín 2011).

3. (iii) The $\bar {T}_c$ shown in figure 2(b) is insensitive to ${{Re}}_\tau$ when ${{Re}}_{\tau,c}^*>200$. Meanwhile, $\bar {T}_c$ increases at a comparable speed with the adiabatic wall temperature (Song et al. 2022).

Figure 2. Contours of (a) the friction Reynolds number in semi-local units ${{Re}}_{\tau,c}^*$, and (b) the centreline temperature. The black line in (b) is the contour line of ${{Re}}_{\tau,c}^*=200$.

It is worth mentioning that the ratio of $\bar {T}_c/T_w$ increases quickly with ${{Ma}}_b$ rise. For example, $\bar {T}_c/T_w$ is over 5.2 at ${{Ma}}_b=5$. This results in a relatively large temperature gradient, and hence strong heat transfer at the wall. Therefore, the supersonic channel flow case can be considered as a highly cooled wall case in analogy to the boundary-layer flow. In this sense, the increase of ${{Ma}}_b$ in the channel flow leads to a wall-cooling effect simultaneously, which will be discussed further in § 4.3. It is also important to note that since $T_w<\bar {T}_c$, ${{Re}}_{\tau }^*$ decreases with ${{Ma}}_b$ rise, as shown in figure 2(a), contrary to the trend of an adiabatic or moderately cooled wall boundary-layer case with $T_w>T_\infty$.

3. Formulation of the linear response analysis

3.1. Linearized Navier–Stokes equation

To study the characteristics of $\tilde {\boldsymbol {q}}$, (2.1) is linearized around the mean state $\bar {\boldsymbol {q}}$. The resulting equation of $\tilde {\boldsymbol {q}}$ is written into a matrix form as

(3.1)\begin{align} &\boldsymbol{\mathsf{F}}\,\frac{\partial \tilde{\boldsymbol{q}}}{\partial t} + \boldsymbol{\mathsf{A}}\,\frac{\partial \tilde{\boldsymbol{q}}}{\partial x} + \boldsymbol{\mathsf{B}}\,\frac{\partial \tilde{\boldsymbol{q}}}{\partial y} + \boldsymbol{\mathsf{C}}\,\frac{\partial \tilde{\boldsymbol{q}}}{\partial z} + \boldsymbol{\mathsf{D}} \tilde{\boldsymbol{q}} \nonumber\\ &\quad= \boldsymbol{\mathsf{H}}_{xx}\,\frac{{\partial ^2}\tilde{\boldsymbol{q}}}{\partial x^2} + \boldsymbol{\mathsf{H}}_{yy}\,\frac{{\partial ^2} \tilde{\boldsymbol{q}}}{\partial y^2} +\boldsymbol{\mathsf{H}}_{zz}\,\frac{{\partial ^2}\tilde{\boldsymbol{q}}}{\partial z^2} + \boldsymbol{\mathsf{H}}_{xy}\,\frac{{\partial ^2}\tilde{\boldsymbol{q}}}{\partial x\,\partial y} + \boldsymbol{\mathsf{H}}_{yz}\,\frac{{\partial ^2}\tilde{\boldsymbol{q}}}{\partial y\,\partial z} + \boldsymbol{\mathsf{H}}_{xz}\,\frac{{\partial ^2}\tilde{\boldsymbol{q}}}{\partial x\,\partial z} + {\tilde{\boldsymbol{N}}}, \end{align}

where the first eleven terms are linear terms of $\tilde {\boldsymbol {q}}$, and ${\tilde {\boldsymbol {N}}}$ collects all the nonlinear terms. The coefficients $\boldsymbol{\mathsf{F}}$, $\boldsymbol{\mathsf{A}}$, $\boldsymbol{\mathsf{B}}$, $\boldsymbol{\mathsf{C}}$, $\boldsymbol{\mathsf{D}}$ and $\boldsymbol{\mathsf{H}}$ are all $5\times 5$ matrices related only to ${\bar {\boldsymbol {q}}}$. Their specific expressions are listed in Appendix C.

In a turbulent flow, $\tilde {\boldsymbol {q}}$ and ${\tilde {\boldsymbol {N}}}$ can have large amplitudes. To facilitate the prediction capability of the linear model, turbulence models are used for the linearization of ${\tilde {\boldsymbol {N}}}$, such as the widely used eddy-viscosity-enhanced linear model. The introduction of $\mu _t$ leads to better results in predicting the response growth, because it partially models the ‘colour’ of the nonlinear forcing term (Moarref & Jovanović 2012; Morra et al. 2019, 2021). In compressible flows, however, the expression of ${\tilde {\boldsymbol {N}}}$ is far more complicated than its incompressible counterpart. Following Alizard et al. (2015) (also Pickering et al. 2021), two substitutions in (3.1) are made for the turbulence modelling, i.e.

(3.2a,b)\begin{equation} \bar{\mu} \rightarrow \bar{\mu} + \mu_t, \quad \bar{\kappa} \rightarrow \bar{\kappa} + {\kappa}_t = \frac{\bar{\mu}}{{{Pr}}} + \frac{\mu_t}{{{Pr}}_t}. \end{equation}

Here, $\kappa _t$ and ${{Pr}}_t$ are the turbulent heat conductivity and Prandtl number, and $\mu _t$ and $\kappa _t$ are non-dimensionalized in the same way as $\mu$ and $\kappa$ in (2.2). The assumptions behind (3.2a,b) are the Boussinesq assumption and the classic Reynolds analogy. The usage of the former is analogous to the incompressible counterpart. The deployment of the latter assumes implicitly that subject to the modelled part of the forcing, the diffusion of the internal energy of the response is proportional to the Reynolds stress part. We give more remarks on the usage of (3.2a,b). The two assumptions above may not be very accurate throughout the flow field. Nevertheless, the introduction of $\mu _t$ and $\kappa _t$ can improve the prediction performance of the linear model as long as the modelling points to the ‘right direction’. This is satisfied here for the canonical channel flow, which is away from the destruction of separation bubbles, shock/boundary-layer interactions, shedding vortex, etc. In this relatively ‘clean’ turbulence, the RANS models have proven their capability extensively (Wilcox 2006). Furthermore, as noted by Cossu et al. (2009) and will be shown in § 3.4, the results of the most amplified response are not very sensitive to the shapes of $\mu _t$ and ${{Pr}}_t$.

3.2. Responses to harmonic and stochastic forcing

As the mean flow is homogeneous in the streamwise and spanwise directions, the following Fourier decomposition is applied on $\tilde {\boldsymbol {q}}$ as

(3.3)\begin{equation} \tilde{\boldsymbol{q}} (x,y,z,t) = \iint_{-\infty}^{\infty} {{\hat{\boldsymbol{q}}}(y, t) \exp \left[{\mathrm{i}} (k_x x + k_z z)\right] {\mathrm{d}} k_x \,{\mathrm{d}} k_z}, \end{equation}

where $k_x$ and $k_z$ are the streamwise and spanwise wavenumbers, respectively, and ${\hat {\boldsymbol {q}}}$ is the shape function. The same decomposition is performed on ${\tilde {\boldsymbol {N}}}$. The resulting ${\hat {\boldsymbol {N}}}$ involves the convolution of ${\hat {\boldsymbol {q}}}$ components due to the nonlinearity.

After substituting (3.3) into (3.1), the following equation is arrived at for a single mode with $k_x\ne 0$ or $k_z\ne 0$:

(3.4)\begin{equation} \boldsymbol{\mathsf{F}}\,\frac{\partial {\hat{\boldsymbol{q}}}}{\partial t} + \hat{\boldsymbol{D}} {\hat{\boldsymbol{q}}} + {\hat{\boldsymbol{B}}}\,\frac{\partial {\hat{\boldsymbol{q}}}}{\partial y} + {\hat{\boldsymbol{C}}}\,\frac{\partial^2 {\hat{\boldsymbol{q}}}}{\partial y^2} = {\hat{\boldsymbol{N}}}, \end{equation}

where the rearranged matrices are

(3.5) \begin{equation} \left. \begin{array}{l@{}} {\hat{\boldsymbol{D}}} = \boldsymbol{\mathsf{D}} + {\mathrm{i}} k_x \boldsymbol{\mathsf{A}} + {\mathrm{i}} k_z \boldsymbol{\mathsf{C}} + k_x^2 \boldsymbol{\mathsf{H}}_{xx} + k_x k_z {\boldsymbol{\mathsf{H}}_{xz}} + k_z^2\boldsymbol{\mathsf{H}}_{zz},\\ {\hat{\boldsymbol{B}}} = \boldsymbol{\mathsf{B}} - {\mathrm{i}} k_x \boldsymbol{\mathsf{H}}_{xy} - {\mathrm{i}} k_z \boldsymbol{\mathsf{H}}_{yz},\\ {\hat{\boldsymbol{C}}} ={-}\boldsymbol{\mathsf{H}}_{yy}. \end{array} \right\} \end{equation}

Equation (3.4) can be rewritten into a standard form as

(3.6a–c)\begin{equation} \frac{\partial {\hat{\boldsymbol{q}}}}{\partial t} = \boldsymbol{L} {\hat{\boldsymbol{q}}} + {\hat{\boldsymbol{f}}}, \quad \boldsymbol{L} ={-}\boldsymbol{\mathsf{F}}^{{-}1}\hat{\boldsymbol{D}} - \boldsymbol{\mathsf{F}}^{{-}1}\boldsymbol{\hat B}\, \frac{\partial}{\partial y} - \boldsymbol{\mathsf{F}}^{{-}1}\hat{\boldsymbol{C}}\,\frac{\partial^2}{\partial y^2}, \quad {\hat{\boldsymbol{f}}} = \boldsymbol{\mathsf{F}}^{{-}1} \hat{\boldsymbol{N}}, \end{equation}

where $\boldsymbol {L}$ is a linear operator, and ${\hat {\boldsymbol {f}}}=[\hat {f}_\rho,\hat {f}_u,\hat {f}_v,\hat {f}_w,\hat {f}_T]^{\textrm {T}}$ is the forcing term.

Before discussing the energy amplification of responses, an energy norm of ${\hat {\boldsymbol {q}}}$ (E) needs to be defined. For compressible flows, a widely used form was proposed by Chu (1965) as

(3.7)\begin{align} \|{\hat{\boldsymbol{q}}}\|^2 = ( {\hat{\boldsymbol{q}}}, {\hat{\boldsymbol{q}}} )_E = \int_{{-}1}^{1} { \left( \bar{\rho} {\hat{\boldsymbol{u}}}^{H} {\hat{\boldsymbol{u}}} + \frac{R \bar{T}}{\bar{\rho}}\,\hat{\rho}^{\dagger} \hat{\rho} + \frac{\bar{\rho} c_v}{{{Ec}}\, \bar{T}}\,\hat{T}^{\dagger} \hat{T} \right) {\mathrm{d}} y} = \int_{{-}1}^{1} {{\hat{\boldsymbol{q}}}^{H} \boldsymbol{M} {\hat{\boldsymbol{q}}} \, {\mathrm{d}} y} \equiv 2 \hat{E}, \end{align}

under the current definition of non-dimensionalization. Here, the superscript ${\dagger}$ stands for the complex conjugate, $c_v=c_v^d/c_p^d=1/\gamma$, $\boldsymbol {M}$ is a diagonal matrix, and $\hat {E}$ is the energy norm. Compared to the kinetic energy used in incompressible flows, (3.7) additionally considers the contribution from the acoustic and entropy components of ${\hat {\boldsymbol {q}}}$ (see e.g. Chen, Wang & Fu 2022). More importantly, from the numerical point of view, $\boldsymbol {M}$ is strictly positive definite, which guarantees the invertibility of the transfer matrix (see (3.9a,b)). For later use, we re-express (3.7) after the wall-normal discretization (see § 3.3 for details) in a vector form as

(3.8a,b)\begin{equation} ( {\hat{\boldsymbol{q}}}, {\hat{\boldsymbol{q}}} )_E = \hat{\boldsymbol{Q}}^{H} \boldsymbol{G}_{MI}\,\hat{\boldsymbol{Q}} \equiv \hat{\boldsymbol{Q}}_{MI}^{H}\,\hat{\boldsymbol{Q}}_{MI}, \quad \hat{\boldsymbol{Q}}_{MI} = \boldsymbol{G}_{MI}^{1/2}\,\hat{\boldsymbol{Q}}, \end{equation}

which is convenient for applying standard matrix functions. Here, the global vector $\hat {\boldsymbol {Q}}=[{\hat {\boldsymbol {q}}}_1,{\hat {\boldsymbol {q}}}_2,\ldots,{\hat {\boldsymbol {q}}}_{N_y}]^{\textrm {T}}$ contains all the wall-normal components, the global matrix $\boldsymbol {G}_{MI} = \textrm {blkdiag}\{ c_1 \boldsymbol {M}_1, \ldots, c_{N_y} \boldsymbol {M}_{N_y} \}$, with $c_j$ the coefficient for numerical integration, and $\boldsymbol {G}_{MI}^{1/2}$ is from the Cholesky decomposition of $\boldsymbol {G}_{MI}$.

In the following, the responses to harmonic and stochastic forcing are accounted for, in turn. A harmonic forcing takes the form ${\hat {\boldsymbol {f}}}(y,t)={\check {\boldsymbol {f}}}(y) \exp (-{\mathrm {i}} \omega t)$, where $\omega$ is the circular frequency. After a sufficiently long time, the response is also harmonic as ${\hat {\boldsymbol {q}}}(y,t)={\check {\boldsymbol {q}}}(y) \exp (-{\mathrm {i}} \omega t)$ if the system is linearly stable. From (3.6a–c), the transfer matrix $\boldsymbol {H}$ between ${\check {\boldsymbol {f}}}$ and ${\check {\boldsymbol {q}}}$ is

(3.9a,b)\begin{equation} {\check{\boldsymbol{q}}} = \boldsymbol{H} {\check{\boldsymbol{f}}}, \quad \boldsymbol{H} =(-{\mathrm{i}}\omega \boldsymbol{I} - \boldsymbol{\mathsf{L}})^{{-}1}, \end{equation}

where $\boldsymbol {H}$ is also known as the resolvent. The energy amplification factor between the forcing and the response is defined as $R(\omega,k_x,k_z)=\|\check {\boldsymbol {q}}\|^2 / \|\check {\boldsymbol {f}}\|^2$; then the optimal response is

(3.10)\begin{equation} R_{max}(k_x,k_z) = R(\omega_{opt};k_x,k_z) = \max_{\omega} \frac{\|\check{\boldsymbol{q}}\|^2}{\|\check{\boldsymbol{f}}\|^2} = \|\boldsymbol{H}\|_\infty^2, \end{equation}

where $\omega _{opt}$ is the optimal frequency. Note that $R_{max}$ is also called the $H_\infty$ norm of the transfer matrix (Zhou, Doyle & Glover 1996). Mathematically, $\|\boldsymbol {H}\|_\infty$ is obtained through a singular value decomposition of $\boldsymbol {H}$, thus $R_{max}=\sigma _1^2$, where $\sigma _1$ is the maximum singular value. Consequently, the relation between $\check {\boldsymbol {q}}$ and $\check {\boldsymbol {f}}$ is

(3.11)\begin{equation} \check{\boldsymbol{q}} = \sum_j { \sigma_j \boldsymbol{\cdot} ({\check{\boldsymbol{\phi}}}_j, \check{\boldsymbol{f}} )_E \boldsymbol{\cdot} {\check{\boldsymbol{\psi}}}_j}, \end{equation}

where $\sigma _1\ge \sigma _j\ge 0$ are the singular values in descending order, and ${\check {\boldsymbol {\psi }}}_j$ and ${\check {\boldsymbol {\phi }}}_j$ are the corresponding response and forcing modes (McKeon & Sharma 2010), respectively, which satisfy the orthogonal relations

(3.12a,b)\begin{equation} (\check{\boldsymbol{\psi}}_i,\check{\boldsymbol{\psi}}_j)_E = \delta_{ij}, \quad (\check{\boldsymbol{\phi}}_i,\check{\boldsymbol{\phi}}_j)_E = \delta_{ij}, \end{equation}

where $\delta _{ij}$ denotes the Kronecker delta. For the optimal response, the unit forcing is $\check {\boldsymbol {\phi }}_1$ and the response is $\check {\boldsymbol {\psi }}_1$. As the channel mean flow is symmetric around the centreline, $\check {\boldsymbol {\psi }}_j$ and $\check {\boldsymbol {\phi }}_j$ are either symmetric or antisymmetric. For the mode of small wavelength, the singular values even appear in pairs (i.e. $\sigma _{2j}=\sigma _{2j-1}$), so the shape function is not unique. For uniqueness, the shape function is normalized such that $\check {u}_1$ is real at its peak amplitude (Moarref et al. 2013).

When the system (3.6a–c) is driven by a stochastic forcing, the response is also stochastic. Note that the colour of the nonlinear forcing has been modelled partially by introducing $\mu _t$ and ${{Pr}}_t$, so the remaining part of the forcing, i.e. ${\hat {\boldsymbol {f}}}$ here, is assumed to be white-in-time. Following previous works (e.g. Farrell & Ioannou 1993; Hwang & Cossu 2010b), ${\hat {\boldsymbol {f}}}$ is a $\delta$-correlated Gaussian white noise process with zero mean:

(3.13a–c)\begin{equation} {\hat{\boldsymbol{f}}} = \boldsymbol{B} {\hat{\boldsymbol{f}}}_0, \quad \langle {\hat{\boldsymbol{f}}}_0 \rangle = 0, \quad \langle {\hat{\boldsymbol{f}}}_0(t_1)\, {\hat{\boldsymbol{f}}}_0^{H}(t_2) \rangle = \boldsymbol{I}\,\delta_d (t_1 - t_2), \end{equation}

where $\delta _d$ is the Dirac function, ${\hat {\boldsymbol {f}}}_0$ is the unmasked forcing, and $\boldsymbol {B}$ is a mask matrix allowing for the calculation of separate responses to different forcing components (Jovanović & Bamieh 2005; Madhusudanan & McKeon 2022). Specifically, $\boldsymbol {B}$ is an identity matrix in normal cases. If, for example, only $\hat {f}_T$ – i.e. the temperature component of the forcing – is considered, then $\boldsymbol {B}=\boldsymbol {B}_T=\textrm {diag}([0\ 0\ 0\ 0\ 1])$.

The variance of the response $V=\langle \|{\hat {\boldsymbol {q}}}\|^2\rangle$ is focused on as a measure of energy amplification, where $\langle {\cdot }\rangle$ denotes the ensemble average. Mathematically, $V$ is also referred to as the $H_2$ norm of the transfer matrix, representing a temporal integral of $\textrm {tr}(\boldsymbol {H}\boldsymbol {H}^{H})$ in terms of $\omega$ (Zhou et al. 1996), where $\textrm {tr}({\cdot })$ is the matrix trace. Equivalently, we have

(3.14a,b)\begin{equation} V = \langle\hat{\boldsymbol{Q}}_{MI}^{H}\,\hat{\boldsymbol{Q}}_{MI}\rangle = \text{tr}(\boldsymbol{X}_{MI}), \quad \boldsymbol{X}_{MI} = \langle\hat{\boldsymbol{Q}}_{MI}\,\hat{\boldsymbol{Q}}_{MI}^{H}\rangle. \end{equation}

Here, the global covariance tensor $\boldsymbol {X}_{MI}$ is obtained by solving the algebraic Lyapunov equation

(3.15a,b)\begin{align} \boldsymbol{G}_{L,MI}\,\boldsymbol{X}_{MI} + \boldsymbol{X}_{MI}\, \boldsymbol{G}_{L,MI}^{H} + \boldsymbol{G}_{B,MI}\,\boldsymbol{G}_{B,MI}^{H} = {\boldsymbol{0}}, \quad \langle \hat{\boldsymbol{Q}}\,\hat{\boldsymbol{Q}}^{H} \rangle = \boldsymbol{G}_{MI}^{{-}1/2}\,\boldsymbol{X}_{MI}\,\boldsymbol{G}_{MI}^{{-}1/2{H}}, \end{align}

where $\boldsymbol {G}_{L,MI}=\boldsymbol {G}_{MI}^{1/2}\,\boldsymbol {G}_{L}\, \boldsymbol {G}_{MI}^{-1/2}$ and $\boldsymbol {G}_{L}$ is the global matrix of $\boldsymbol {L}$, and a similar definition for $\boldsymbol {G}_{B,MI}$. As a result, $V$ is obtained using (3.14a,b) without calculating the temporal integral. We can proceed to solve the eigenvalue problem of $\langle \hat {\boldsymbol {Q}} \hat {\boldsymbol {Q}}^{H} \rangle$ as

(3.16)\begin{equation} \langle \hat{\boldsymbol{Q}} \hat{\boldsymbol{Q}}^{H} \rangle\,\boldsymbol{G}_{MI} {\hat{\boldsymbol{Y}}}_j = \theta_j {\hat{\boldsymbol{Y}}}_j, \end{equation}

where $\theta _j$ and ${\hat {\boldsymbol {Y}}}_j$ are the eigenvalue and eigenfunction, and ${\hat {\boldsymbol {Y}}}_j$ satisfies the same orthogonal relation as (3.12a,b). Note that ${\hat {\boldsymbol {Y}}}_j$ is also called the KL decomposition or the ‘proper orthogonal decomposition’ mode, representing the leading energy-containing structure of the response. The forcing mode is solved using the back Lyapunov equation. More details can be found in Farrell & Ioannou (1993), Zhou et al. (1996) and Bamieh & Dahleh (2001).

The shape function of the harmonic forcing will be denoted as ${\hat {\boldsymbol {q}}}$ not $\check {\boldsymbol {q}}$ below, where there is no ambiguity. This is intended for consistency with the notation of the stochastic forcing and thus convenience for writing.

3.3. Numerical consideration

The procedures in § 3.2 require a discretization in the wall-normal direction. The differential matrices are from the Chebyshev collocation point method (Trefethen 2001). By default, $N_y=501$ points are used, which is abundant to ensure grid independence (Bae et al. 2020b). It is worth mentioning that the $c_j$ for numerical integration come from those of Clenshaw & Curtis (1960), rather than the direct inverse of the differential matrix, to ensure the positive definiteness of $\boldsymbol {G}_{MI}$.

The problems of singular value decomposition, Lyapunov equation and eigenvalue, as described in (3.11), (3.15a,b) and (3.16), respectively, are solved using the svd, eig and lyap functions of the Matlab software. In terms of the boundary condition, a no-slip and isothermal wall is assumed on both sides. Therefore, the perturbations at the wall satisfy

(3.17)\begin{equation} \hat{u}_w = \hat{v}_w = \hat{w}_w = \hat{T}_w = 0 \quad \text{at}\ y={\pm} 1, \end{equation}

and $\hat {\rho }_w$ is solved through the perturbed continuity equation. For channel flows, $\hat {\rho }_w$ has the same finite amplitudes at both wall sides, and oscillating unphysical modes may appear. These unphysical modes are recognized and removed in the program as they cannot converge with increasing $N_y$ (Bewley & Liu 1998). Verification of the linear response solvers with previous results is provided in Appendix B.

3.4. Eddy viscosity and thermal conductivity

If the turbulent statistics are not available, then $\mu _t$ can be evaluated from the mean streamwise momentum (2.1b) as

(3.18)\begin{equation} 0 ={-}\frac{{\mathrm{d}} \bar{P}}{\mathrm{d}\kern0.06em x} + \frac{{\mathrm{d}}}{{\mathrm{d}} y} \left[ (\bar{\mu} + \mu_t)\, \frac{{\mathrm{d}} \bar{U}}{{\mathrm{d}} y} \right] \quad \rightarrow \quad \mu_t = \frac{1 - y^+{/}{{Re}}_\tau}{{{\mathrm{d}} \bar{U}^+}/{{\mathrm{d}} y^+}} - \bar{\mu}, \end{equation}

the same as that in (2.7). For ${{Pr}}_t$, its expression by definition is

(3.19)\begin{equation} {{Pr}}_t = \frac{\overline{\tilde{u}\tilde{v}}\,\partial \bar{T} / \partial y}{\overline{\tilde{v}\tilde{T}} \, \partial \bar{U} / \partial y}, \end{equation}

which is shown to be approximately unity in most regions in the channel flow (Huang et al. 1995; Modesti & Pirozzoli 2016). Therefore, ${{Pr}}_t$ can be assumed to be simply a constant. The linear response results of different ${{Pr}}_t$ are compared below.

To highlight the influence of ${{Pr}}_t$, the case of a relatively high ${{Ma}}_b$ of 4 is selected, and ${{Re}}_\tau ={6000}$. Figure 3(a) gives the distribution of $\kappa _t$ with ${{Pr}}_t=0.8$, 1.0 and 1.2, respectively. Despite the relatively large difference in ${{Pr}}_t$, the amplification factor of the optimal harmonic response is little affected, as shown in figure 3(b). A visible difference exists only in the region $k_zh<1$. The same conclusion applies to the stochastic response (not shown here). For the convenience of calculation, ${{Pr}}_t=1$ is selected for later use.

Figure 3. (a) Eddy thermal conductivity and (b) energy amplification factor of the optimal harmonic response ($k_x=0$) based on the mean flow of different turbulent Prandtl number. The case parameters are ${{Ma}}_b=4$ and ${{Re}}_\tau ={6000}$.

3.5. Linear response calculation based on the ODE solver

Although the ODE solver can provide an accurate mean flow, the corresponding linear response results may have large deviations from those using the DNS mean flow since there are first- and second-order derivatives of the mean flow in (3.4). Here, the reliability of the ODE-based mean flow in predicting the linear response is examined further.

Two representative cases are selected. The first is the incompressible ${{Re}}_\tau =5200$ case from Lee & Moser (2015), whose ${{Re}}_\tau$ is very high for a channel flow. This helps to examine the present solver when extended to the high-Reynolds-number regime. Note that ${{Ma}}_b$ is set to 0.01 in the present solver for this case, to avoid singularities. The second case is the ${{Ma}}_b=1.5$, ${{Re}}_\tau =1910$ case from Yao & Hussain (2020) (no. (15) in table 2 in Appendix A), which has both a high ${{Re}}_\tau$ and obvious compressibility effects. The energy amplification factors of both optimal harmonic and stochastic responses for the two cases are plotted in figure 4. As can be seen, the results based on the mean flow from the ODE solver match well with those using the DNS mean flow. A small visible difference exists near $k_zh\sim 10^{2}$ in figure 4(b). This is attributed to the difference of $\bar {U}$ at the junction of the buffer and logarithmic layers, which originates from the deviation of the incompressible curve fit in (2.7) (see Appendix A for more details). In short, figure 4 demonstrates the reliability of the present ODE solver in calculating the linear response of the turbulent mean flow.

Figure 4. Energy amplification factors of the (a,b) optimal harmonic and (c,d) stochastic forcing based on the mean flows from DNS and the ODE solver. Two cases are (a,c) the ${{Re}}_\tau =5200$ case from Lee & Moser (2015), and (b,d) the ${{Ma}}_b=1.5$, ${{Re}}_\tau =1910$ case from Yao & Hussain (2020).

4. Results and discussions

4.1. Basic characteristics of responses

The basic characteristics of the responses to the harmonic and stochastic forcing are analysed first. The case parameters are selected to be ${{Ma}}_b=2.0$ and ${{Re}}_\tau =6000$, which is a benchmark case in this work. A parameter study will be conducted in § 4.3.

First, the curves of $R_{max}$ and $V$ with different $k_x$ and $k_z$ are plotted in figures 5(a) and 5(b). An overall feature is that the $R_\textrm {max}$ and $V$ values in the large-scale region (in terms of the spanwise wavelength $\lambda _z=2{\rm \pi} /k_x$) are orders of magnitude higher than those in the small-scale region, reflecting the energetic nature of the former (Kim & Adrian 1999). Meanwhile, $R_{max}$ and $V$ of the large-scale motion are strongly dependent on $k_x$ when $\lambda _z>0.3h$. The modes with the largest $R_{max}$ and $V$ both have an infinite streamwise wavelength (i.e. $k_x=0$); for the harmonic forcing, the most amplified mode is also steady ($\omega _{opt}=0$). In contrast, $R_{max}$ and $V$ in the small-scale region are very insensitive to $k_x$, where the curves of different $k_x$ nearly collapse with each other.

Figure 5. (a,b) Energy amplification factors, and (c,d) corresponding pre-multiplied factors of the (a,c) optimal harmonic and (b,d) stochastic forcing in the ${{Ma}}_b=2$, ${{Re}}_\tau =6000$ case.

Hwang & Cossu (2010b) found, in their incompressible flow case, that there exist a $k_z^{-2}$ scaling and a $k_z^{-1}$ scaling in the mid-$k_z$ range for $R_{max}$ and $V$, respectively. We also apply these scalings for the compressible case here for examination. A theoretical analysis on the linear operator related to these scalings will be provided in § 5. The pre-multiplied energy amplification factors $k_z^2 R_{max}$ and $k_z V$ are plotted in figures 5(c) and 5(d). The variation of $k_z^2 R_{max}$ bears a strong resemblance to that in the incompressible case, where the classic bimodal structure is observed (del Álamo & Jiménez 2006; Cossu et al. 2009). These two peaks correspond to the large-scale motions in the outer layer (termed ‘outer peak’) and the small-scale motions in the near-wall region (termed ‘inner peak’), respectively, as also noted by Alizard et al. (2015) using the transient growth analysis. The $\lambda _z^+$ of the inner peak, denoted $\lambda _{z,i}^+$, is approximately 100, while the $\lambda _z$ of the outer peak, denoted $\lambda _{z,o}$, decreases with $k_x$ rise and has maximum value approximately $3.7h$ at $k_x=0$. For $k_z V$, it also peaks near $\lambda _{z,o}=3.8h$ in the large-scale region, but exhibits distinct features from the incompressible case in the small-scale region (see figure 19(b) in Appendix B for reference). The slope of $V$ does not decrease but increases at $k_zh>400$, so the inner peak is missing in $k_z V$. For reference, the variations of $\theta _1$ and $k_z \theta _1$ are also plotted in figures 5(b) and 5(d). As stated in § 3.2, $\theta _1$ is the largest eigenvalue of the covariance matrix, and represents the mode that receives the largest energy gain. As can be seen, $k_z \theta _1$ has the same trend of variation as $k_z V$ in the large-scale region, but exhibits the inner peak in the small-scale region. In addition, the inner peak also locates at approximately $\lambda _z^+=100$, very close to that of $k_z^2 R_{max}$. The different behaviour of $k_z V$ here from the incompressible case indicates that the compressibility effects can be the cause. This is discussed further below.

The mask matrix $\boldsymbol {B}$ (see (3.13a–c)) is used to study the responses to different forcing components. Therefore, $\boldsymbol {B}$ is separated into two parts, i.e. the kinematic part and the thermodynamic part (Madhusudanan & McKeon 2022), as

(4.1a–c)\begin{equation} \boldsymbol{B} = \boldsymbol{B}_{\boldsymbol{u}} + \boldsymbol{B}_{\rho T}, \quad \boldsymbol{B}_{\boldsymbol{u}}=\text{diag}([0\ 1\ 1\ 1\ 0]), \quad \boldsymbol{B}_{\rho T}=\text{diag}([1\ 0\ 0\ 0\ 1]). \end{equation}

The variances of the responses to the two forcing parts are denoted $V_{\boldsymbol {u}}$ and $V_{\rho T}$, respectively. By definition, $V_{\boldsymbol {u}}$ results from the forcing applied only to three velocity equations, i.e. ${\hat {f}}_\rho ={\hat {f}}_T=0$, the same as that in an incompressible flow. Likewise, $V_{\rho T}$ results from the forcing components present only in the density and temperature equations, i.e. ${\hat {f}}_u={\hat {f}}_v={\hat {f}}_w=0$. It is worth mentioning that ${\hat {f}}_T=0$ does not mean $\hat {T}=0$ because of the variable coupling. Figure 6(a) gives the distribution of $V$, $V_{\boldsymbol {u}}$ and $V_{\rho T}$. Approximate linear superposition $V \approx V_{\boldsymbol {u}} + V_{\rho T}$ is satisfied with relative error less than 2 %. The shape of $k_z V_{\boldsymbol {u}}$ is quite similar to that in the incompressible flow (see figure 19 for reference). It tends to decrease at $k_zh>400$, so the inner peak related to the near-wall motion can be recognized, also located at $\lambda _z^+\approx 100$. In comparison, $k_z V_{\rho T}$ increases monotonically with $k_z$. It is one order of magnitude smaller than $k_z V_{\boldsymbol {u}}$ around the outer peak, but rises to be the larger one with $k_z$ higher than the inner peak ($\lambda _z^+<100$), which results simultaneously in the quick increase of $k_z V$. Therefore, it is the thermodynamic part of the forcing – i.e. the components in the equations of $\hat {\rho }$ and $\hat {T}$ – that leads to the different trend of $k_z V$ in the compressible flow case.

Figure 6. (a) Pre-multiplied energy growth factors, and (b) contribution of the two most energetic KL modes of the stochastic response to different forcing components ($k_x=0$).

Furthermore, the KL modes of the three covariance tensors corresponding to $V$, $V_{\boldsymbol {u}}$ and $V_{\rho T}$ are computed to see the contribution from the most energetic modes. The contribution is measured by the ratio $r_V=(\theta _1+\theta _2)/V$. The second mode ($\theta _2$) is also included because in the small-$\lambda _z$ region, the $\theta _j$ appear in pairs, i.e. $\theta _1=\theta _2$. As shown in figure 6(b), the two maxima of $r_{V_{\boldsymbol {u}}}$ are 89 % and 42 %, near the outer and inner peaks, respectively, indicating that a large portion of $V_{\boldsymbol {u}}$ is contributed by one mode pair (or two modes). This energetic mode is representative in describing the local turbulent characteristics. In comparison, $r_{V_{\rho T}}$ decreases continuously with $k_z$, and is less than 3 % at $k_zh>380$ ($\lambda _z^+<100$). This suggests that no single mode dominates in terms of the energy contribution, i.e. the response is of little coherence in the small-$\lambda _z$ region, though $V_{\rho T}$ is large. As $r_{V_{\boldsymbol {u}}}$ is much larger than $r_{V_{\rho T}}$, $\theta _{1,\boldsymbol {u}}$ dominates $\theta _{1,\rho T}$ in most regions. Therefore, $\theta _1$ nearly equals $\theta _{1,\boldsymbol {u}}$, or in other words, the most energetic mode to the total variance is contributed mainly by $\theta _{1,\boldsymbol {u}}$. The $r_V$ at the outer and inner peaks are 81 % and 21 %, respectively, pulled down by $r_{V_{\rho T}}$. The above comparison demonstrates clearly the distinct feature between the responses to the kinematic and thermodynamic components of the forcing. The former is analogous to the incompressible counterpart, exhibiting characteristic large- and small-scale motions. The latter significantly affects the small-scale motion, and is of little coherence.

4.2. Response structure

The shape functions of the modes at the outer and inner peaks ($k_x=0$) in figure 5 are plotted in figure 7, where both the response and forcing are provided. The wall-normal and spanwise velocities are the two largest components in the input signal, and so are the streamwise velocity and temperature components in the output signal. One significant feature is that the shape functions of the optimal harmonic and stochastic forcing/response are nearly identical, in both shape and amplitude, as also noted by Hwang & Cossu (2010b) in the incompressible case. This is because ${\hat {\boldsymbol {q}}}$ of the harmonic response is very insensitive to $\omega$ at $k_x=0$, so ${\hat {\boldsymbol {q}}}$ is nearly unchanged from a harmonic forcing to a broadband stochastic one. From a physical point of view, the convection term ${\bar {\boldsymbol {U}}}\boldsymbol {\cdot }\boldsymbol {\nabla }{{\hat {\boldsymbol {q}}}}$ vanishes with $k_x=0$, hence the response shape is determined primarily by the local mean-flow gradient. In terms of different components, the peak locations of $|\hat {v}|$ and $|\hat {w}|$ are different, forming a streamwise vortex to be shown later. The maxima of $|\hat {u}|$ and $|\hat {T}|$ of the inner-peak mode are at $y^+=11$ and $y^+=10$, both of which lie within the buffer layer. For the outer-peak mode, both $\hat {u}$ and $\hat {T}$ penetrate deep into and below the logarithmic region, as labelled by the dashed line at $y^+=100$ in figure 5(d) (recall that ${{Re}}_\tau =6000$). This is indicative of an amplitude modulation effect of the outer-layer motions on the near-wall velocity and temperature (Marusic et al. 2010a; Cheng, Shyy & Fu 2022; Yu & Xu 2022). Specifically, the forcing imposed in the outer region (especially $\hat {v}$) induces large amplitudes of streamwise velocity and temperature near the wall as the response. Meanwhile, $\hat {T}$ penetrates deeper towards the wall than $\hat {u}$, owing to the relatively large temperature gradient near the wall.

Figure 7. Shape functions of different components and the energy norm of the (a,b) forcing and (c,d) response at $k_x=0$ for the mode at the (a,c) inner peak ($\lambda _z^+=100$) and (b,d) outer peak ($\lambda _z=3.7h$) in figure 5.

As discussed in figure 6, the most energetic modes at the outer and inner peaks result mainly from the kinematic parts of the forcing. In this sense, $\hat {T}$ behaves more like a passive scalar (noted by Coleman et al. 1995), resulting from the interaction between velocity components of the forcing and the mean temperature gradient. This is demonstrated in figure 8, where the contours of $\tilde {u}$ and $\tilde {T}$ in the $y$–$z$ plane are depicted for optimal harmonic response. The shape of the stochastic response is nearly the same (see figure 7) and is thus not shown here. The velocity vector $[\tilde {w},\tilde {v}]$ of the forcing is also plotted to identify the in-plane motion. As can be seen, the classic combination of quasi-streamwise vortices and streamwise streaks is recognized for the modes of both inner and outer peaks. Moreover, the streamwise vortex as the input strongly promotes the streaks of streamwise velocity and temperature as the output through the lift-up effect (Ellingsen & Palm 1975; Schoppa & Hussain 2002). For both modes here, the temperature streaks look more stuck to the wall, as already shown in figure 7.

Figure 8. Contours of the (a,c) streamwise velocity and (b,d) temperature of the optimal harmonic response for the modes of the (a,b) inner peak ($\lambda _z^+=100$) and (c,d) outer peak ($\lambda _z=3.7h$) in the $y$–$z$ plane. The velocity vector is based on the spanwise and wall-normal velocities of the forcing.

More discussion is provided on the structure of the outer-peak mode. As shown in figures 8(c) and 8(d), the contour of the response almost fills in the whole spatial domain, and extends to as long as $\lambda _{z,o}=3.7h$ in the spanwise direction. This value is very close to that in incompressible channel flows (e.g. Cossu et al. 2009; Pujals et al. 2009). Actually, as will be shown in § 4.3, $\lambda _{z,o}$ is nearly independent of ${{Ma}}_b$ and ${{Re}}_\tau$, standing as a pretty robust feature of the optimal response in the outer layer. However, as noted by previous researchers, such a large-scale structure has not been observed in the DNS data or the experiment. The outer peak of the pre-multiplied energy spectra from DNS indicates a spanwise scale $\lambda _z\approx h- 2h$ for both incompressible (Lee & Moser 2015) and compressible (Yao & Hussain 2020) channel flows. Part of the reason is that the structures in figure 8 correspond to infinite $\lambda _x$ ($k_x=0$), but in actual flows it is hard to extend to that long. According to figure 5(c), $\lambda _{z,o}$ decreases at a lower $\lambda _x$, to be closer to the DNS data. In addition, nonlinear effects of the response growth can also lead to the difference (Cossu et al. 2009).

In the incompressible flow case, the $\hat {u}$ response with different $k_z$ values between those of the inner and outer peaks exhibits a geometrically self-similar structure (Pujals et al. 2009; Hwang & Cossu 2010b; McKeon 2019). Specifically, the normalized $|\hat {u}|$ of different $k_z$ can be collapsed approximately using a wall-normal coordinate $y^+/\lambda _z^+$. The amplitude peaks lie within the logarithmic layer, so this self-similar structure stands as a supporting proof for the well-known AEM (Townsend 1976). In addition, the outer-layer large-scale motions and near-wall small-scale motions can be regarded as the largest and smallest attached eddies, respectively (Hwang 2015; Cheng et al. 2019, 2022). Although abundant evidence has appeared to support the AEM in the incompressible turbulence, its behaviour and validity in compressible turbulence are far from clear. Here, we provide some understanding of the compressible AEM by examining the geometrical self-similarity of the response.

The streamwise velocity, temperature and energy norm of the optimal harmonic response are plotted in figure 9, with $\lambda _z^+$ ranging from 800 to 4000 ($\lambda _z/h=0.13- 0.67$). They all exhibit a streamwise-elongated streaky structure, as observed in the high-speed experiment (Ganapathisubramani, Clemens & Dolling 2006). The peaks of $|\hat {u}|$ of these modes lie within the logarithmic layer, but those of $|\hat {T}|$ are much closer to the wall due to large mean temperature gradient there. To collapse the highly scattering distribution of the shape function, the wall-normal coordinate $y^+$ is normalized by $\lambda _z^+$. As the amplitude of ${\hat {\boldsymbol {q}}}$ is determined by (3.12a,b) based on the integral form in (3.7), the amplitude is normalized correspondingly using $\lambda _z^+$. As shown in figures 9(d)–9(f), the normalized curves of different $\lambda _z^+$ are much more concentrated, and are collapsed to a first approximation. Also, the maximum of the normalized $|\hat {T}|$ experiences larger variations with $\lambda _z^+$ than those of $|\hat {u}|$ and $\hat {E}$. For comparison, the shapes of the inner-peak mode ($\lambda _z^+=100$), as shown in figure 7, are also plotted in black dashed lines. As can be seen, the curve of $\lambda _z^+=100$ apparently deviates from others of higher $\lambda _z^+$. This is reasonable because its amplitude peak is already within the buffer layer (see figure 7), out of the logarithmic region, which violates the assumptions of the AEM. Another wall-normal coordinate widely used in compressible turbulence is the semi-local one, $y^*$. But here, the ratio $y^*/\lambda _z^*$ equals $y^+/\lambda _z^+$, thus leading to no difference in collapsing curves. As the unit $\delta _\nu ^*$ and thus $\lambda _z^*$ are not invariant at different $y^+$, other attempts at normalization can be deployed to reflect the difference between $y^+$ and $y^*$. For example, $\lambda _z^*(y^+)$ can be replaced with $\lambda _z^*(y_m^+)$, where $y_m^+$ is a fixed point for each $\lambda _z^+$. As a choice of $y_m^+$, the peak location of $|\hat {u}|$ (or $|\hat {T}|$, $\hat {E}$), denoted $y_{U_m}^+$, can be regarded as a characteristic wall-normal scale of the response. Nevertheless, using $y^*/\lambda _z^*(y_{U_m}^+)$ as the wall-normal coordinate results in no visible improvement in collapsing curves (hence not shown here).

Figure 9. Wall-normal distribution of the (a,d) streamwise velocity, (b,e) temperature, and (c,f) energy norm for the optimal harmonic response ($k_x=0$) with different spanwise wavelengths. The black dashed lines are the corresponding results of the $\lambda _z^+=100$ mode. The coordinates in (d–f) are normalized by the spanwise wavelength.

Finally, the response modes of $\lambda _z^+<100$ are focused on because they have different behaviours from the incompressible case (see figure 6). Figures 10(a) and 10(b) display two representative cospectra, $\langle \hat {u}\hat {T}^{\dagger} \rangle$ and $\langle \hat {\rho }\hat {T}^{\dagger} \rangle$, for the modes of $\lambda _z^+=100$ and $\lambda _z^+=20$, respectively. The cospectrum is obtained by extracting the sub-elements of $\langle {\hat {\boldsymbol {q}}}{\hat {\boldsymbol {q}}}^{H} \rangle$ in (3.15a,b). Moreover, to study the contributors to these cospectra, the responses to the kinematic and thermodynamic parts of the stochastic forcing are calculated separately, denoted as $\langle {\hat {\boldsymbol {q}}}{\hat {\boldsymbol {q}}}^{H} \rangle _{\boldsymbol {u}}$ and $\langle {\hat {\boldsymbol {q}}}{\hat {\boldsymbol {q}}}^{H} \rangle _{\rho T}$, respectively. Their relative contributions to the cospectrum are plotted in figures 10(c) and 10(d) in percentages. As can be seen, over 95 % of $\langle \hat {u}\hat {T}^{\dagger} \rangle$ is contributed by $\langle \hat {u}\hat {T}^{\dagger} \rangle _{\boldsymbol {u}}$ for both modes, indicating the strong coupling between the streaks of streamwise velocity and temperature induced by the streamwise vortices as shown in figure 8. On the other hand, the contribution from the thermodynamic part of the forcing is negligible. With $\lambda _z^+$ decreasing, the amplitude of $\langle \hat {u}\hat {T}^{\dagger} \rangle$ quickly diminishes (see figure 6), to less than half of $\langle \hat {\rho }\hat {T}^{\dagger} \rangle$ at $\lambda _z^+=20$. This is reasonable because as the streak motion is closer to the wall, the $\hat {v}$ component is restricted by the wall boundary, thus the streamwise vortex is suppressed. Compared with $\langle \hat {u}\hat {T}^{\dagger} \rangle$, $\langle \hat {\rho }\hat {T}^{\dagger} \rangle$ has two hot zones. One is around the buffer layer as well, resulting from the coupling of the temperature and density streaks. The other is attached to the wall within the viscous layer, contributed mainly by $\langle \hat {\rho }\hat {T}^{\dagger} \rangle _{\rho T}$. As shown in figure 6, the response of the second hot zone is of little coherence, and is connected to the large mean density gradient there through the perturbed continuity equation. With $\lambda _z^+$ decreasing from 100 to 20, the peak amplitude of $\langle \hat {\rho }\hat {T}^{\dagger} \rangle$ shifts to the wall, reflecting the strengthening of the relative importance of $\langle {\hat {\boldsymbol {q}}}{\hat {\boldsymbol {q}}}^{H} \rangle _{\rho T}$.

Figure 10. (a,b) Component cospectrum in the $y^+$–$y^+$ plane, and (c,d) the relative contribution from the response to different components of the forcing. The two $y^+$ in the horizontal and vertical coordinates in (a,b) are for the two variables in $\langle {\cdot }\rangle$, respectively. The response mode is $k_x=0$, with (a,c) $\lambda _z^+=100$ and (b,d) $\lambda _z^+=20$. Each pair of plots (contour and bar) in (a,c) and (b,d) represents the same variable as labelled.

4.3. Effects of Reynolds and Mach numbers

The cases of different ${{Ma}}_b$ and ${{Re}}_\tau$ are investigated further, with an emphasis on compressibility effects. Figure 11 depicts the energy amplification factors of the optimal harmonic responses for different cases. Both the outer and inner units are deployed to study the characteristics of the outer and inner peaks. For the outer-region motion, one significant feature from figures 11(a) and 11(c) is that the outer peaks all correspond to $\lambda _{z,o}\approx 3.6h$ for different ${{Re}}_\tau$ and ${{Ma}}_b$. Here, ${{Re}}_\tau$ is from $10^{3}$ to $10^{4}$, and ${{Ma}}_b$ is from 0.1 to 5, both covering a relatively large range. The weak dependence of $\lambda _{z,o}$ on ${{Re}}_\tau$ has already been reported in previous research (see e.g. Pujals et al. 2009), attributed to the fact that the large-scale motion in the outer layer cannot directly ‘feel’ the presence of the wall, and is hence least affected (Pirozzoli 2014). Furthermore, the weak dependence of $\lambda _{z,o}$ on ${{Ma}}_b$ is demonstrated in the present work. This trend is consistent with the observation in experiment (Spina et al. 1994; Williams et al. 2018) and from DNS of the energy spectrum at different ${{Ma}}_b$ (Modesti & Pirozzoli 2016; Yao & Hussain 2020; Cogo et al. 2022; Cheng & Fu 2022b). As discussed in figure 6, the large-scale structure corresponding to the outer peak results mainly from the kinematic part of the forcing, not the thermodynamic part. Thereby, its motion nature can partly explain the weak dependence on ${{Ma}}_b$.

Figure 11. Pre-multiplied energy amplification factors of the optimal harmonic responses ($k_x=0$) for the cases with (a,b) different ${{Re}}_\tau$ (${{Ma}}_b=2$), and (c,d) different ${{Ma}}_b$ (${{Re}}_\tau =6000$). Plots (a,c) are in the outer units, and (b,d) are in the inner units. Note that the inset in (b) gives the same results as (b) except that ${{Ma}}_b=4$.

For the mode of the inner peak, it is observed that the pre-multiplied energy growths of different ${{Re}}_\tau$ are nearly collapsed using the inner units for the Mach 2 case (see figure 11b), where $\lambda _{z,i}^+\approx 100$. This collapse is also achieved when ${{Ma}}_b$ is increased to 4, as displayed in the inset, except that $\lambda _{z,i}^+$ rises to 180. Actually, $\lambda _{z,i}^+$ indeed increases with ${{Ma}}_b$, as shown in figure 11(d). Meanwhile, with the increase of ${{Ma}}_b$, the $k_z^{-2}$ scaling of $R_{max}$ in the mid-$k_z$ range has more significant deviations, and $k_z^2 R_{max}$ of the inner peak rises in both outer and inner units. This trend will be discussed further in § 5 by analysing the perturbation equations. Through a comparison with the $k_z^2 R_{max}$ of the outer peak, the relative importance of the near-wall motion is enhanced, in terms of amplifying the response energy.

The variation of $\lambda _{z,i}^+$ with ${{Ma}}_b$ is further quantified, as depicted in figure 12(a). The $\lambda _{z,i}^+$ evaluated from $k_z^+ \theta _1^+$ of the stochastic response is also plotted. As discussed in figure 5, the variation of $k_z \theta _1$ bears a strong resemblance to $k_z^2 R_{max}$. It is observed that $\lambda _{z,i}^+$ increases by over three times with ${{Ma}}_b$ lifted up from 0.1 to 5, whether evaluated from $k_z^+ \theta _1^+$ or $k_z^{2+} R_{max}^+$. As noted in § 2.2, the increase of ${{Ma}}_b$ results simultaneously in a wall-cooling effect for the channel flow. Wall cooling is shown in the DNS results to rapidly increase the spanwise spacing of the near-wall streaks for both boundary-layer and channel flows (Morinishi et al. 2004; Duan et al. 2010; Zhang et al. 2022), consistent with the trend in figure 12(a). For more quantitative comparison, available DNS results of channel flows are added. The data are from Yao & Hussain (2020) with a maximum ${{Ma}}_b$ of 1.5. The $\lambda _{z,i}^+$ plotted are calculated from the local maximum of their pre-multiplied energy spectra of the streamwise velocity at a fixed $y^*=15$. As the spectra data are not publicly available, the points are evaluated by the present authors from their published figures, so error bars are added to display the potential error. There is a systematic difference between the $\lambda _{z,i}^+$ from the linear response analysis and the DNS results, which has also been noted and interpreted in previous research, and will be discussed at length later.

Figure 12. Variation of the spanwise wavelength of the inner-peak mode with ${{Ma}}_b$, expressed using the (a) inner units and (b) semi-local units for both the optimal harmonic and stochastic responses. The results evaluated from the DNS energy spectra at $y^*=15$ are also plotted (Yao & Hussain 2020).

As demonstrated by Yao & Hussain (2020), if expressed using semi-local units, $\lambda _{z,i}^*$ values are all approximately 120 for the three cases with varying ${{Ma}}_b$. To apply this scaling, the wall-normal height to determine $\delta _\nu ^*$ is selected to be $y_{E_m}^+$, i.e. the peak location of $\hat {E}$, as introduced in § 4.2. The resulting variation of $\lambda _{z,i}^*$ with ${{Ma}}_b$ is plotted in figure 12(b). It is observed that both $\lambda _{z,i}^*$ from the optimal and stochastic responses are nearly invariant with ${{Ma}}_b$ up to 5, consistent with the trend in the DNS data in terms of the spanwise spacing of the near-wall streaks (Morinishi et al. 2004; Patel et al. 2016; Zhang et al. 2022). Nevertheless, the current DNS database of channel flows cannot reach a ${{Ma}}_b$ as high as 5. In this sense, the result in figure 12(b) is predictive, waiting for more reference points at higher ${{Ma}}_b$ for further examination. From the variation of $\delta _v^*$, it is conjectured that the increase of $\lambda _{z,i}^+$ with ${{Ma}}_b$ rise is due mainly to the $\bar {T}/T_w$ increase, or equivalently the $\bar {\rho }/\rho _w$ decrease.

More discussion is provided on the systematic deviation of $\lambda _{z,i}^+$ and $\lambda _{z,i}^*$ in figure 12 from the DNS data. The first reason is that only the principal mode (pair) is considered here. Though it explains a large portion of the total energy amplification (see figure 6), secondary modes may also play crucial roles. An approach for improvement is to consider several leading modes, such as the work by McMullen et al. (2020), which impressively reproduced the energy spectra in an incompressible case. However, the determination of weights for different modes involves the usage of instantaneous DNS data, which limits the prediction capability, especially for compressible flows. Second, the energy amplification factor in figure 11 at a given $(k_x, k_z, \omega )$ is a global measure of the response at all wall-normal heights, so it is not a direct analogy to the energy spectra in the DNS, which is fixed at a specific height. As noted by Cossu et al. (2009), the $\lambda _{z,i}^+$ of the inner peak denotes the ‘most probable’ value of the spanwise spacing, while that from DNS is more like a ‘mean’ value.

The energy norm of the inner-peak mode is displayed in figure 13 for the cases of varying ${{Re}}_\tau$ and ${{Ma}}_b$. The usage of $\hat {E}$ is a universal measure of the response structure rather than $\hat {u}$ or $\hat {T}$ alone across different ${{Ma}}_b$. As discussed in figure 12, $y^*$ is used for $\hat {E}$ of various ${{Ma}}_b$, and hence also for the cases of varying ${{Re}}_\tau$ for consistency. As can be seen, the $\hat {E}$ profiles of different ${{Re}}_\tau$ are nearly collapsed in figures 13(b) and 13(c). In comparison, the $\hat {E}$ profiles of different ${{Ma}}_b$ exhibit visible differences; $y_{E_m}^+$ decreases continuously with ${{Ma}}_b$ rises, which is due to the increasingly large temperature gradient of the mean flow.

Figure 13. Wall-normal distributions of the normalized energy norm of the inner-peak mode for the cases (a) ${{Re}}_\tau =6000$ and varying ${{Ma}}_b$, (b) ${{Ma}}_b=2$ and varying ${{Re}}_\tau$, and (c) ${{Ma}}_b=4$ and varying ${{Re}}_\tau$.

The geometrical self-similarity of the response in the mid-$k_z$ range has been discussed in § 4.2 for the case ${{Ma}}_b=2$. Here, the self-similar characteristics are examined further for the cases of different ${{Ma}}_b$. Figure 14 gives the normalized energy norms of different $\lambda _z^+$ with ${{Ma}}_b=0.1$, 3 and 5, in the same style as in figure 9. The minimum $\lambda _z^+$ displayed increases from 500 to 800 with ${{Ma}}_b$ rise because, as shown in figure 9, the mode with $\lambda _z^+$ close to $\lambda _{z,i}^+$ has obvious deviations, and $\lambda _{z,i}^+$ increases at a higher ${{Ma}}_b$ (see figure 12a). As can be seen, the energy norms of different $\lambda _z$ in the mid-$\lambda _z$ range all approximately collapsed with ${{Ma}}_b$ ranging from 0.1 to 5. The above observation demonstrates that the geometrical self-similarity of the response is little influenced by the compressibility effect, which is helpful in understanding the AEM for supersonic flows.

Figure 14. Wall-normal distributions of the normalized energy norm for the cases with ${{Re}}_\tau ={6000}$ and (a) ${{Ma}}_b=0.1$, (b) ${{Ma}}_b=3$, and (c) ${{Ma}}_b=5$, for the optimal harmonic response ($k_x=0$) with different spanwise wavelengths.

As shown in figure 6, $V_{\boldsymbol {u}}$ and $V_{\rho T}$ have totally different behaviours at different $k_z$. Their variations with ${{Ma}}_b$ are investigated further, as plotted in figure 15. The variation of $k_z V_{\boldsymbol {u}}$ with ${{Ma}}_b$ is strongly analogous to $k_z^2 R_{max}$ (see figure 11d). Both the inner and outer peaks are present in the curves of different ${{Ma}}_b$. Furthermore, the $\lambda _{z,o}^+$ of $V_{\boldsymbol {u}}$ is nearly unchanged corresponding to $\lambda _{z,o}=3.6h$, while $\lambda _{z,i}^+$ increases continuously with ${{Ma}}_b$, as already described in figure 12(b). The $k^{-1}$ scaling in the mid-$k_z$ range gradually deviates at high ${{Ma}}_b$, and the relative importance of the inner peak is enhanced in amplifying the response. In comparison, $V_{\rho T}$ experiences a much stronger variation with ${{Ma}}_b$ than $V_{\boldsymbol {u}}$. As ${{Ma}}_b$ increases, $k_z^+ V_{\rho T}^+$ in the large-scale region tends to diminish, while that in the small-scale region increases dramatically. Consequently, for the case ${{Ma}}_b>0.5$, $k_z^+ V_{\rho T}^+$ starts to monotonically increase with the rise of $k_z^+$. It is also observed that both $V_{\boldsymbol {u}}$ and $V_{\rho T}$ are subject to the $k^{-1}$ scaling when ${{Ma}}_b=0.1$ in the mid-$k_z$ range (approximately $k_z^+$ from 0.002 to 0.1), as labelled by the horizontal dashed lines in figure 15. With the increase of ${{Ma}}_b$, this scaling gradually deviates, especially for $V_{\rho T}$. In short, the motion due to the kinematic part of the forcing is much less sensitive to the compressibility effect than the thermodynamic part. With ${{Ma}}_b$ increasing, the thermodynamic components begin to dominate the total variance of the stochastic response in the small-scale region (roughly $\lambda _z^+<60$).

Figure 15. Pre-multiplied variances of the stochastic responses to the forcing of only the (a) velocity components and (b) thermodynamic components in the cases of varying ${{Ma}}_b$ (${{Re}}_\tau =6000$).

As discussed in § 4.1, the response corresponding to $V_{\rho T}$ is much less coherent. Therefore, for the response to full forcing ($V$), the increase of $V_{\rho T}$ with ${{Ma}}_b$ can decrease the dominance of the principal mode pair and thus the coherence of the flow. This is examined by plotting the ratio $r_V=(\theta _1+\theta _2)/V$ for the cases of different ${{Ma}}_b$, as displayed in figure 16. The ratio corresponding to the optimal harmonic response, $r_\sigma =(\sigma _1^2+\sigma _2^2)/\sum \sigma _j^2$, is also plotted. As can be seen, $r_V$ and $r_\sigma$ are not affected by ${{Ma}}_b$ in the large-scale region near the outer peak. The motion there, with $r_V$ and $r_\sigma$ both larger than 80 %, is dominated by that induced by the kinematic part of the forcing, also insensitive to ${{Ma}}_b$ and ${{Re}}_\tau$ (see figure 11). In the small-scale region, $r_V$ experiences different trends of variation with ${{Ma}}_b$ in the ranges larger and smaller than $\lambda _z^+=90$, respectively. It is the same for $r_\sigma$ except that the dividing line is at $\lambda _z^+=130$. Specifically, $r_V$ (also $r_\sigma$) increases with ${{Ma}}_b$ for $\lambda _z^+>90$, i.e. the principal mode pair has an increasingly large contribution to the total energy growth, which is attributed to the increase of $V_{\boldsymbol {u}}^+$. The DNS results also indicate that wall cooling tends to increase the coherence of the near-wall motion in this wavelength range. The interpretation is that the anisotropic motion is suppressed by wall cooling due to the reduction in energy exchange (Coleman et al. 1995; Morinishi et al. 2004; Duan et al. 2010; Zhang et al. 2022). In contrast, it is observed here that $r_V$ (also $r_\sigma$) quickly decreases with ${{Ma}}_b$ at $\lambda _z^+<90$, indicating a reduction in coherence and also the dominance of the principal mode pair for the total energy growth. This trend is obviously due to the enhancement of the response induced by the thermodynamic part of the forcing, as shown in figure 15(b).

Figure 16. Contribution from the principal mode pair to the total energy amplification (in %) for the (a) stochastic and (b) optimal harmonic responses in the cases of varying ${{Ma}}_b$.

5. Theoretical analysis of the linear operator

It is observed in figures 11(c) and 15 that the pre-multiplied energy amplification of the optimal harmonic and stochastic response gradually deviates from the $k_z^{-2}$ and $k_z^{-1}$ scalings in the mid-$k_z$ range at higher ${{Ma}}_b$. We appeal to a theoretical analysis on (3.1) to support this observation.

To simplify the form of the equations, the power law $\mu =T^n$ is adopted temporarily for derivation use instead of the original Sutherland's law. This results in only a minor difference because $\mu$ is much smaller than $\mu _t$. The mode of $k_x=0$ is considered, and the terms related directly to $\bar {\rho }$ and $\bar {T}$ are combined to cluster the effect of compressibility in the equation. For example, the continuity and streamwise momentum equations are written as

(5.1) \begin{equation} \left. \begin{array}{l@{}} \displaystyle \dfrac{\partial}{\partial t} \dfrac{\hat{\rho}}{\bar{\rho}} - \dfrac{\bar{T}_y}{\bar{T}}\,\hat{v} + \left( \hat{v}_y + {\mathrm{i}} k_z \hat{w} \right) = \dfrac{\hat{f}_\rho}{\bar{\rho}} ,\\ \displaystyle \dfrac{\partial \hat{u}}{\partial t} + \bar{U}_y \hat{v} - \dfrac{1}{{{Re}}_b} \left[ \dfrac{\tilde{\mu}_{t,y}}{\bar{\rho}}\,\hat{u}_y + \dfrac{\tilde{\mu}_t}{\bar{\rho}} \left(\hat{u}_{yy} - k_z^2 \hat{u}\right) \right]\\ \displaystyle \quad{}- \dfrac{n}{{{Re}}_b}\,\dfrac{\bar{\mu}}{\bar{\rho}} \left[ \left(\bar{U}_{yy} + (n-1) \bar{U}_y\,\dfrac{\bar{T}_y}{\bar{T}}\right) \dfrac{\hat{T}}{\bar{T}} + \bar{U}_y\,\dfrac{\hat{T}_y}{\bar{T}} \right] = \hat{f}_u , \end{array} \right\} \end{equation}

where the subscript $y$ denotes the wall-normal derivative, and $\tilde {\mu }_t=\mu _t+\bar {\mu }$ is the total viscosity. As can be seen, the mean-flow coefficients are $\bar {U}$, $\bar {T}$, $\tilde {\mu }_t$, $\bar {\mu }$ and their wall-normal derivatives. The other three equations can be rearranged in the same way; their expressions are specified in Appendix C. Therefore, if the basic variable $\hat{\boldsymbol{q}}$ is replaced by $\hat{\boldsymbol{q}}_0$, then the operator form of the equation is

(5.2a,b)\begin{equation} \frac{\partial \hat{\boldsymbol{q}}_0}{\partial t} = \mathcal{L} \hat{\boldsymbol{q}}_0 + \hat{\boldsymbol{f}}_0, \quad \hat{\boldsymbol{q}}_0 = \left[\frac{\hat{\rho}}{\bar{\rho}}, \hat{u},\hat{v},\hat{w},\frac{\hat{T}}{\bar{T}}\right], \end{equation}

and the linear operator is

(5.3)\begin{equation} \mathcal{L} = \mathcal{L} \left( \bar{U}_y,\bar{U}_{yy}, \frac{\bar{T}_y}{\bar{T}},\frac{\bar{T}_{yy}}{\bar{T}}, \frac{\bar{T}}{{{Ma}}_b^2}, \frac{\tilde{\mu}_t}{\bar{\rho}}, \frac{\tilde{\mu}_{t,y}}{\bar{\rho}}, \frac{\bar{\mu}}{\bar{\rho}} \right). \end{equation}

The variation of these mean-flow coefficients reflects the influencing path of compressibility on the response behaviours. Note that $\kappa _t$ and $\bar {\kappa }$ are related to viscosity through the constants ${{Pr}}$ and ${{Pr}}_t$, so they are not included explicitly in the bracket. In the following, the behaviours of these mean-flow coefficients are analysed for the cases of different $k_z$.

The self-similarity of the shape function is in terms of the normalized coordinate $y^+/\lambda _z^+$, so a rescaled coordinate is introduced as $\check {y}=\bar {k}_z (y+1)$, where $\bar {k}_z=k_z/\check {k}_z$, and $\check {k}_z$ is an arbitrarily selected reference value. Note that the overhead check in this section represents the rescaled variable, different from that defined in § 3.2 for a shape function The rescaled form of (5.2a,b) is written as

(5.4)\begin{equation} \frac{\partial \hat{\boldsymbol{q}}_0}{\partial \check{t}} = \tilde{\mathcal{L}} \hat{\boldsymbol{q}}_0 + \frac{1}{\bar{k}_z}\, \hat{\boldsymbol{f}}_0, \end{equation}

where $\check {t}=\bar {k}_z t$. If the rescaled operator $\tilde {\mathcal {L}}$ is independent of $\bar {k}_z$, i.e. $\tilde {\mathcal {L}}=\check {\mathcal {L}}$, then Hwang & Cossu (2010b) proved that there are a $k_z^{-2}$ scaling for $R_{max}$ and a $k_z^{-1}$ scaling for $V$ in the mid-$k_z$ range (see figure 5). Conversely, the variable that cannot be made independent of $\bar {k}_z$ is responsible for the failure of the scaling.

As $\hat {\boldsymbol {q}}_0$ is a solution defined on the $\check {y}$-coordinate, its wall-normal derivative follows

(5.5)\begin{equation} \frac{\partial}{\partial y} = \bar{k}_z\,\frac{\partial}{\partial \check{y}}. \end{equation}

However, this chain rule is not applicable for the derivatives of mean-flow quantities because they are solutions on the $y$-coordinate. Taking the streamwise momentum equation as an example, the rescaled form is written as

(5.6) \begin{equation} \left. \begin{array}{l@{}} \displaystyle \dfrac{\partial \hat{u}}{\partial \check{t}} + \dfrac{\bar{U}_y}{\bar{k}_z}\,\hat{v} - \dfrac{1}{{{Re}}_b} \left[ \dfrac{\tilde{\mu}_{t,y}}{\bar{\rho}}\,\hat{u}_y + \dfrac{\bar{k}_z \tilde{\mu}_t}{\bar{\rho}} \left(\hat{u}_{\check{y}\check{y}} - \check{k}_z^2 \hat{u}\right) \right] \\ \displaystyle \quad{}- \dfrac{n}{{{Re}}_b}\,\dfrac{\bar{k}_z \bar{\mu}}{\bar{\rho}} \left[ \left(\dfrac{\bar{U}_{yy}}{\bar{k}_z^2} + (n-1)\,\dfrac{\bar{U}_y}{\bar{k}_z}\,\dfrac{\bar{T}_y}{\bar{k}_z \bar{T}}\right) \dfrac{\hat{T}}{\bar{T}} + \dfrac{\bar{U}_y}{\bar{k}_z}\, \dfrac{\hat{T}_{\check y}}{\bar{T}} \right] = \dfrac{\hat{f}_u}{\bar{k}_z}. \end{array} \right\} \end{equation}

The mean-flow coefficients are analysed as follows. The peaks of the perturbation in figure 9 are within the logarithmic region, where the streamwise velocity from the current velocity transformation (see (A2b)) satisfies

(5.7)\begin{equation} \bar{\mu}\,\frac{\partial \bar{U}^+}{\partial y^+} = \frac{{\mathrm{d}} \bar{U}_{TL}^+}{{\mathrm{d}} y^*} = \frac{1}{\kappa_c y^*} \quad \rightarrow \quad \bar{U}_{y} = \frac{u_\tau}{\kappa_c}\,\frac{\sqrt{\bar{T}}}{y+1}, \end{equation}

where $\kappa _c$ is the Kármán constant. Therefore, the rescaled velocity gradient is

(5.8)\begin{equation} \frac{\bar{U}_y(y)}{\bar{k}_z} = \frac{u_\tau}{\kappa_c}\, \frac{\sqrt{\bar{T}(y)}}{\check{y}} = \check{\bar U}_{\check{y}}\, \frac{\sqrt{\bar{T}(y)}}{\sqrt{\bar{T}(\check{y})}} \equiv K_1(y,\check{y})\,\check{\bar U}_{\check{y}}. \end{equation}

As can be seen, only if $K_1=1$, i.e. $\bar {T}(y)$ is constant, is ${\bar {U}_y}/{\bar {k}_z}=\check {\bar U}_{\check {y}}$ independent of $\bar {k}_z$ in the logarithmic region. Similarly, the temperature gradient from (A3) is

(5.9)\begin{equation} \frac{\bar{T}_y}{\bar{T}} = \frac{ C_1 \left[2 (C_2 - C_3) \bar{U} + C_3 \bar{U}_c \right]}{1 + C_1 \bar{U} \left[(C_2 - C_3) \bar{U} + C_3 \bar{U}_c \right]}\,\bar{U}_y \equiv C_T(y)\,\bar{U}_y, \end{equation}

so the rescaled temperature gradient is

(5.10)\begin{equation} \frac{\bar{T}_y(y)}{\bar{k}_z\,\bar{T}(y)} = C_T(y)\, \frac{\bar{U}_y(y)}{\bar{k}_z} = K_2\,\frac{\check{\bar T}_{\check y}(\check{y})}{\check{\bar T}(\check{y})}, \quad \text{where } K_2(y,\check{y}) = \frac{\sqrt{\bar{T}(y)}\, C_T(y)}{\sqrt{\bar{T}(\check{y})}\,C_T(\check{y})}. \end{equation}

For the eddy viscosity, (3.18) in the logarithmic region gives

(5.11)\begin{equation} \frac{\tilde{\mu}_t}{\bar{\rho}} = \frac{\kappa_c\,{{Re}}_\tau}{\rho_w}\, \sqrt{\bar{T}}\,(y+1) \left[ 1 - (y+1) \right], \end{equation}

so the rescaled version is

(5.12)\begin{equation} \frac{\bar{k}_z\,\tilde{\mu}_t(y)}{\bar{\rho}(y)} \approx \frac{\kappa_c\, {{Re}}_\tau}{\rho_w}\,\sqrt{\bar{T}}\,\check{y} \left( 1 - \check{y} \right) = K_1\,\frac{\check{\tilde{\mu}}_t(\check{y})}{\check{\bar{\rho}}(\check{y})}, \end{equation}

where the approximation is due to the assumption that $(y+1)\ll 1$ (recall that the wall is at $y=-1$). The other two mean-flow coefficients in (5.7) are rescaled as

(5.13) \begin{align} \left. \begin{array}{c@{}} \displaystyle \dfrac{\bar{U}_{yy}(y)}{\bar{k}_z^2} = \Bigg(\dfrac{\bar{T}_y}{2\bar{k}_z \bar{T}} - \dfrac{1}{\bar{k}_z (y+1)}\Bigg) \dfrac{\bar{U}_y(y)}{\bar{k}_z} = K_1 K_2 \Bigg[ \check{\bar U}_{\check{y}\check{y}}(\check{y}) + \Bigg(1-\dfrac{1}{K_2}\Bigg) \dfrac{\check{\bar{U}}_{\check{y}}(\check{y})}{\check{y}} \Bigg], \\ \displaystyle \dfrac{\tilde{\mu}_t(y)}{\bar{\rho}(y)} \approx K_1 K_2 \Bigg[ \dfrac{\check{\tilde{\mu}}_{t,\check{y}}(\check{y})}{\check{\bar{\rho}}} + \Bigg(1-\dfrac{1}{K_2}\Bigg) \dfrac{\check{\tilde{\mu}}_t(\check{y})}{\check{y}\, \check{\bar{\rho}}(\check{y})} \Bigg]. \end{array} \right\} \end{align}

Therefore, if $K_1=K_2=1$, then ${\bar {U}_{yy}}/{\bar {k}_z^2}$ and ${\tilde {\mu }_t}/{\bar {\rho }}$ are independent of $\bar {k}_z$. The rescaled expressions in (5.8), (5.10), (5.12) and (5.13) are all explicitly irrelevant of $\bar {k}_z$, and dependent on $K_1$ and $K_2$ instead. The remaining mean-flow coefficient in (5.7), $\bar {k}_z \bar {\mu }/\bar {\rho }$, is related only to $\bar {T}$ and cannot be effectively made independent of $\bar {k}_z$. The other four equations can be analysed similarly. As a result, all the mean-flow coefficients can be categorized into two groups: (i) the ones that can be rescaled using $K_1$ and $K_2$ in the logarithmic region, including $\bar {U}_y$, $\bar {U}_{yy}$, $\bar {T}_y/\bar {T}$, $\tilde {\mu }_t/\bar {\rho }$ and $\tilde {\mu }_{t,y}/\bar {\rho }$ (see (5.3)); (ii) the ones related only to $\bar {T}$ that cannot be made explicitly irrelevant to $\bar {k}_z$, including $\bar {\mu }/\bar {\rho }$ and $\bar {T}/{{Ma}}_b^2$.

In summary, three factors are primarily responsible for the potential failure of the $k_z^{-2}$ and $k_z^{-1}$ scalings. The first is related to the logarithmic region. If a large portion of the shape function is out of the logarithmic region, then (5.7) and equations following all fail. This is why the scaling applies only in the mid-$\lambda _z$ range and is consistent with the assumption of the AEM. The second factor is the substantial variations of $K_1$ and $K_2$ at different wall-normal heights. The third factor is the enormous contribution from the terms related to $\bar {\mu }/\bar {\rho }$ and $\bar {T}/{{Ma}}_b^2$. The latter two are absent in the incompressible case, both related directly to the distribution of mean temperature.

To see the variations of $K_1$ and $K_2$, the distributions of the functions $\sqrt {\bar {T}}$ and $\sqrt {\bar {T}}\,C_T$ (see (5.8) and (5.10)) are plotted in figure 17 for the cases of different ${{Ma}}_b$ and ${{Re}}_\tau =6000$. As can be seen, with the rise of ${{Ma}}_b$, the two functions, and thus $K_1$ and $K_2$ at fixed $y$ and $\check {y}$, experience much larger variations. Meanwhile, the other coefficients related to $\bar {T}$ also change rapidly with $y^+$. This explains the trend observed in figures 11(c) and 15 that the $k^{-2}$ and $k^{-1}$ scaling for $R_{max}$ and $V$ have larger deviations at higher ${{Ma}}_b$. In short, for the high-${{Ma}}_b$ case of the channel flow, the sizeable mean temperature gradient near the wall destroys the $k_z$ power scaling of the response in the mid-$\lambda _z$ range.

Figure 17. Wall-normal distributions of the two mean-flow-related coefficients corresponding to (a) $K_1$ and (b) $K_2$, for the cases of different ${{Ma}}_b$ (${{Re}}_\tau =6000$).

The above analysis should also have implications for the geometrical self-similarity of the response as discussed before. For example, $\sqrt {\bar {T}}$ and $\sqrt {\bar {T}}\,C_T$ in figure 17 change more slowly in the logarithmic region, especially for the former, which benefits the self-similarity of shape functions. However, the operation of the matrix inversion in (3.9a,b) obscures further analytical derivation on the perturbation equations and the transfer matrix. In fact, the results in figure 14 indicate the insensitivity of the self-similarity of the energy norm to ${{Ma}}_b$.

6. Summary

In this work, the linear responses of turbulent mean flow to both harmonic and stochastic forcing are investigated for supersonic channel flow. Well-established universal relations are utilized to obtain efficiently the mean profiles with a large parameter space, with ${{Ma}}_b$ up to 5 and ${{Re}}_\tau$ up to $10^{4}$. As a result, a systematic parameter study is feasible. Comparisons with the DNS data demonstrate the reliability of the present solver in both calculating the mean flow and serving as the base flow for the linear response analysis.

Both the optimal harmonic and stochastic forcing are considered. The response of $k_x=0$ has the largest $R_{max}$ and $V$. The pre-multiplied energy growth $k_z^2 R_{max}$ exhibits the classic bi-modal structure in terms of $k_z$. The two peaks correspond to the large-scale motion in the outer layer, and the near-wall small-scale motion in the buffer layer, respectively. For $k_z V$, however, the inner peak disappears when ${{Ma}}_b>0.5$; $k_z V$ increases monotonically with $k_z$. A further decomposition of the forcing shows that $k_z V_{\boldsymbol {u}}$, resulting from the kinematic part of the forcing, is quite similar to the incompressible counterpart, exhibiting the inner peak. However, the thermodynamic part, $k_z V_{\rho T}$, increases quickly at a larger $k_z$, and the resulting motions are much less coherent than the former. Specifically, the contribution of the most energetic mode of $V_{\rho T}$ to the energy growth is less than 3 % at $\lambda _z^+<100$ (${{Ma}}_b=2$).

The spanwise wavenumber of the outer peak is $\lambda _{z,o}\approx 3.6h$, nearly unchanged at different ${{Ma}}_b$ and ${{Re}}_\tau$. The motion takes the form that the input is the streamwise vortex, and the output is the streamwise streak of velocity and temperature. The temperature streak penetrates deeper towards the wall than the streamwise velocity. The spanwise wavenumber of the inner peak $\lambda _{z,i}^+$ is nearly independent of ${{Re}}_\tau$, but increases quickly with ${{Ma}}_b$, due mainly to the wall-cooling effects in terms of the smaller ratio $T_w/\bar {T}_c$. From ${{Ma}}_b=0.1$ to 5, $\lambda _{z,i}^+$ is lifted by over three times. If using semi-local units, $\lambda _{z,i}^*$ is approximately unchanged with ${{Ma}}_b$. Wall-cooling effects enhance the coherence of the flow for the stochastic response of $\lambda _{z,m}^+>90$ and the optimal harmonic response of $\lambda _{z,m}^+>130$, consistent with the DNS data. On the contrary, it is found that the coherence of the response with $\lambda _z^+<\lambda _{z,m}^+$ decreases with ${{Ma}}_b$. This is due to the strengthening of the less-coherent motion related to $V_{\rho T}$. As a result, the cospectrum of $\langle \hat {\rho }\hat {T}^{\dagger} \rangle$ has a hot zone attached to the wall within the viscous layer.

The response modes with $\lambda _z$ between the inner and outer peaks exhibit a geometrical self-similarity to a first approximation, with ${{Ma}}_b$ from 0.1 to 5. The normalized energy norm is approximately collapsed using the wall-normal coordinate $y^+/\lambda ^+_z$. The self-similarity is related to the logarithmic region, relatively insensitive to ${{Ma}}_b$ and thus the compressibility effects. Finally, a theoretical analysis is conducted on the perturbation equations. The behaviour of the mean-flow coefficients is analysed within the logarithmic region. More importantly, two crucial parameters are summarized related to the dependence of the linear operator on $k_z$. The two parameters are related directly to the mean temperature, experiencing stronger variations in the wall-normal direction with ${{Ma}}_b$ increase, which explains the trend that the $k_z^{-2}$ and $k_z^{-1}$ scalings of $R_{max}$ and $V$ in the mid-$k_z$ range have increasingly large deviations with ${{Ma}}_b$ rise.

Compared with the DNS data, the present results are more qualitative than quantitative. Some interpretations of the differences have been provided earlier in the text. Future attention will be paid to the improvement of the modelling of nonlinear forcing term, and the development of a prediction model on the energy spectra and space–time properties of the fluctuations.