Study of the linear models in estimating coherent velocity and temperature structures for compressible turbulent channel flows

Abstract

Linear models, based on stochastically forced linearized equations, are deployed for spectral linear stochastic estimation (SLSE) of the velocity and temperature fluctuations in compressible turbulent channel flows with a bulk Mach number of 1.5. Through comparison with the direct numerical simulation (DNS) data, an eddy-viscosity-enhanced model (eLNS) outperforms the one not enhanced (LNS) in computing the coherence and amplitude ratio of streamwise velocity at different wall-normal heights, but they both largely deviate from DNS regarding the temperature prediction. For further investigation, the eigenspectra and pseudospectra of the linear operators are scrutinized. The eddy viscosity is shown to stabilize the eigenmodes and decrease the non-normality of the vortical modes. Consequently, the relative importance of acoustic and entropy modes increases, and they can contribute 20 % to 55 % of the response growth, which is not supported by DNS. Hence, it is an intrinsic defect of the eLNS model introduced by turbulence modelling. After a procedure of cospectrum decomposition, the contributions of acoustic and entropy components are filtered out. The resulting SLSE quantities for velocity, temperature and their coupling are basically agreeable with DNS, demonstrating that the coherent temperature fluctuation is dominated by advection and other vortical motions, instead of the compressibility effects. Moreover, a parameter study of Reynolds and Mach numbers (from 0.3 to 4) is conducted. The semi-local units are shown to well collapse the velocity SLSE quantities to the incompressible case for streamwise-elongated structures of high coherence.

1. Introduction

Understanding and modelling compressible turbulent flows are fundamentally significant for aerodynamic applications (Bradshaw 1977). Turbulent boundary layers strongly affect the surface drag and heat transfer, so accurate predictive models are important for reliable vehicle design and flow control (Gatski & Bonnet 2013). Compared with the incompressible one, additional thermodynamic processes are present in high-speed turbulent flows, such as heat transfer, acoustic fluctuations, dilatational work and high-enthalpy effects (Di Renzo, Fu & Urzay 2020; Di Renzo & Urzay 2021; Fu et al. 2021; Fu, Bose & Moin 2022). Due to the theoretical complexity and practical significance, compressible wall-bounded turbulence has become a hot topic.

The mean flow properties of compressible turbulence have been extensively studied. The hypothesis of Morkovin (1962) is widely verified, stating that at moderate free stream Mach numbers (${\textit {Ma}}_\infty \lesssim 5$), the dilatation effect is small, and any differences from incompressible turbulence can be accounted for by mean variations of fluid properties (Coleman, Kim & Moser 1995; Pirozzoli, Grasso & Gatski 2004; Duan, Beekman & Martin 2010). Lagha et al. (2011) even show that the hypothesis is applicable with ${\textit {Ma}}_\infty$ up to 20 for flat-plate boundary layers. Thereby, velocity transformation can be designed using only the mean flow. The transformed streamwise velocity profile can match the incompressible one within and below the logarithmic region with very high accuracy (van Driest 1951; Trettel & Larsson 2016; Griffin, Fu & Moin 2021; Bai, Griffin & Fu 2022). Meanwhile, the mean temperature is shown to be nearly a quadratic function of the mean velocity (Walz 1969; Duan & Martín 2011; Zhang et al. 2014). Strong Reynolds analogy (SRA) and its extensions are proposed, connecting the temperature and velocity fluctuations (Morkovin 1962; Huang, Coleman & Bradshaw 1995; Zhang et al. 2014). In a recent work, Cheng & Fu (2023a) note that the profile of turbulent Prandtl number ${\textit {Pr}}_t$ in supersonic channels is very close to the incompressible case with heat transfer (Antonia, Abe & Kawamura 2009). They suggest that the coherent temperature fluctuations behave like a passive scalar, and the coupling between the velocity and temperature fields results largely from the advection effect.

In addition to the mean flow, much progress has been made in understanding the structures of velocity and temperature fluctuations, owing to the advancement of experimental and numerical techniques, especially the direct numerical simulation (DNS, Moin & Mahesh 1998). In the near-wall region, small-scale temperature streaks are observed elongated in the streamwise direction, analogous to the streamwise velocity streaks (Coleman et al. 1995; Duan et al. 2010). The spanwise spacing of streaks is shown to decrease with rising ${\textit {Ma}}$ on adiabatic walls, and increase by wall cooling; it experiences much less variation if expressed in semi-local units (Morinishi, Tamano & Nakabayashi 2004). Within and above the logarithmic layer, energy-containing large-scale and very-large-scale motions (LSMs and VLSMs), as first reported in incompressible flows (Kim & Adrian 1999), are identified in supersonic boundary layers experimentally (Ganapathisubramani, Clemens & Dolling 2006; Bross, Scharnowski & Kähler 2021) and numerically (Ringuette, Wu & Martín 2008; Huang, Duan & Choudhari 2022). Unlike the near-wall motions, the characteristic scales of LSMs and VLSMs are not very sensitive to compressibility effects (Pirozzoli, Bernardini & Grasso 2008; Williams et al. 2018; Cheng & Fu 2022b). The LSMs populating the logarithmic and outer regions are known to exert footprints, i.e. have effects of velocity modulation and superposition, on the near-wall motions (Abe, Kawamura & Choi 2004; Hutchins & Marusic 2007), which is the basis of the inner–outer interaction model (IOIM, Marusic, Mathis & Hutchins 2010). Using the DNS of supersonic boundary layers, Bernardini & Pirozzoli (2011) and Pirozzoli & Bernardini (2011) demonstrate that the temperature fluctuations also experience inner–outer interaction, though the modulation effect seems weaker than the velocity. Yu & Xu (2022) extend IOIM to realize a predictive model for both velocity and temperature fields in compressible boundary layers, considering the modified SRA and the superposition and modulation from LSMs. Meanwhile, the inner–outer interaction of temperature fluctuations indicates that they can also be regarded as a wall-attached quantity, hence described by the well-known attached-eddy model (AEM, Townsend 1976). This is supported by accumulating evidence concerning the geometrical self-similarity of temperature in the logarithmic region (Cheng & Fu 2022b; Yu et al. 2022; Chen et al. 2023). Moreover, Cheng & Fu (2023a) point out that in the logarithmic region, only the scales corresponding to the attached eddies and VLSMs are firmly coupled; the Reynolds number $Re$ acts as the crucial similarity parameter in constructing the coupling, rather than the Mach number.

To estimate and model coherent large-scale structures, the spectral linear stochastic estimation (SLSE) approach is helpful. SLSE originates from the extensively used stochastic estimation (Adrian 1979; Adrian & Moin 1988), which can provide an estimated velocity signal given a measurement at another point from experiments or simulations. Considering the multi-scale feature of turbulence, SLSE is performed in the Fourier space, which takes advantage of the coherence of large-scale structures and suppresses small-scale random noise (Tinney et al. 2006; Baars, Hutchins & Marusic 2016). A kernel function $H_L$ is designed to relate the Fourier components at different wall-normal heights, based on the linear coherence spectrum (LCS). Over the years, SLSE has been widely deployed to study the multi-scale structures in incompressible flows and to help understand and improve AEM and IOIM (Baars, Hutchins & Marusic 2017; Encinar & Jiménez 2019; Cheng & Fu 2022a, 2023b; Cheng, Shyy & Fu 2022). Nevertheless, its application in compressible flows is still limited. Very recently, Cheng & Fu (2022b, 2023a) extend the SLSE framework for compressible flows. The self-similarity of thermodynamic quantities and the coupling effects between velocity and temperature are thus investigated. This extended SLSE framework is considered in this work for compressible channel flows.

By definition, the LCS and $H_L$ in SLSE involve ensemble-averaged cospectra, so their calculation requires a series of instantaneous measurements or numerical data, which are more challenging to obtain than the mean flow. In an alternative method, Madhusudanan, Illingworth & Marusic (2019) obtain $H_L$ using the linearized incompressible Navier–Stokes (NS) equations, which requires only the mean flow as the input. Their stochastically forced, eddy-viscosity-enhanced linear model gives agreeable LCS with the DNS data, except in the near-wall region. For improvement, Gupta et al. (2021) model the amplitude distribution of the stochastic forcing using mean flow quantities. The resulting LCS and $H_L$ are quite close to the DNS data, though the prediction can still struggle in the near-wall region. Notably, the linearized NS equations have long been utilized in turbulence research, as will be introduced at length later. Thereby, it is a natural thought to develop a linear model for compressible turbulent flows, so that the LCS and $H_L$ can be obtained using only mean flow quantities. In this way, the coherence of velocity and temperature fluctuations, as well as their coupling effects, can be further investigated. Such a linear model is of particular significance for compressible flows, because the instantaneous experimental and DNS data are much more limited than the incompressible counterparts due to larger parameter space and higher facility requirements (Gatski & Bonnet 2013). It is the objective of this work to develop a linear model for SLSE in compressible turbulent channel flows. We find that a direct extension of the incompressible model to the compressible case does not provide satisfactory SLSE results, especially for temperature fluctuations. Thereby, we will scrutinize the mathematical properties and physical relevance of different linear models, and assess carefully their behaviours.

To set more grounds for the present work, the approaches based on the linearized NS equations for turbulent flows are introduced in more detail. Early works find that the incompressible turbulent mean flow (one dimensional) is globally stable, i.e. there are no unstable modes through linear stability analysis (Malkus 1956; Reynolds & Hussain 1972). For this globally stable system, non-modal instability theory and the resolvent-based input–output analysis, first applied for laminar flows (Trefethen et al. 1993; Schmid & Henningson 2001), are extended to study the transient behaviour of turbulent mean flows. These linear approaches are successful in studying the multi-scale characteristic motions (del Álamo & Jiménez 2006; Hwang & Cossu 2010; Abe, Antonia & Toh 2018), constructing low-rank predictive models (McKeon & Sharma 2010; Moarref et al. 2013; Illingworth, Monty & Marusic 2018), designing flow control strategies (Moarref & Jovanović 2012; Ran, Zare & Jovanović 2020) and so on. Many of these works adopt the eddy-viscosity-enhanced models, where the Reynolds stress fluctuation is linearized using the eddy viscosity $\mu _t$, to partly model the colour of the forcing. The eddy-viscosity-enhanced models are shown to perform better than those without using $\mu _t$, especially for estimating fluctuations between different heights (Reynolds & Hussain 1972; Illingworth et al. 2018; Madhusudanan et al. 2019; Morra et al. 2021; Symon, Illingworth & Marusic 2021). The linearized-equation-based approaches have also been deployed for compressible turbulent flows. Alizard et al. (2015) perform the transient growth analysis on turbulent boundary layers with ${\textit {Ma}}_\infty$ up to 4. The inner and outer peak modes are identified from the energy growth curves with different spanwise wavenumbers, analogous to the incompressible cases. The former mode is related to the near-wall streaky motions, and the latter represents the outer-layer LSMs (and VLSMs). Bae, Dawson & McKeon (2020) extend the resolvent analysis to supersonic boundary layers and highlight the distinct features of the relatively supersonic region in the wavenumber space, which is not present in incompressible flows. Also, Chen et al. (2023) develop the linear response analyses subject to both optimal harmonic and stochastic forcing for supersonic channel flows, and analyse the response characteristics over wide ranges of ${\textit {Ma}}$ and $Re$. To model the linearized Reynolds stress and turbulent heat flux, $\mu _t$ and SRA are introduced in some of these works (Alizard et al. 2015; Pickering et al. 2021; Chen et al. 2023). Nevertheless, no quantitative comparisons have been made with DNS or experiments in these works, so the accuracy of these linear models remains unclarified. It is also a task of the present work to make such a quantitative assessment regarding SLSE using the DNS database constructed in the present authors’ group.

In addition, there are some theoretical issues that require to be addressed regarding the linearized equations for compressible turbulent flows. For example, the singular value of the resolvent by Bae et al. (2020) (their figure 1) has relatively strong oscillations when the response mode is relatively supersonic. Madhusudanan & McKeon (2022) (and also our own experience) note that relatively supersonic modes have amplified acoustic components in the free stream, and the results can be dependent on the height of the computational domain (which is artificially selected) for boundary layers. These facts suggest that there can be some fundamental differences in the linear models between incompressible and compressible flows, in addition to the widely addressed similarities. However, the differences have not been fully resolved. Also, in deriving the compressible eddy-viscosity-enhanced models, multiple options are available to linearize different terms, leading to numerous combinations. They deserve a careful discussion through inter-model comparisons.

The primary attention of this work is on the following three aspects, which are also its value. First, we carefully discuss the derivation and mathematical properties of the linearized equations for compressible turbulent flows, which is instructive for interpreting subsequent results (§§ 2–4). Second, we utilize the linear model to perform SLSE for streamwise velocity, temperature, and their coupling, aiming to predict the coherent velocity and temperature fluctuations based on the measurement at outer locations. Direct comparison with DNS is conducted to assess their behaviours (§ 5). Third, we discuss the physics inferred from SLSE results, especially for temperature, and conduct a parameter study to see the effects of Reynolds and Mach numbers (§§ 5 and 6). Some remarks are provided on the scope of this work. For incompressible flows, many linear models have been developed to provide continuously improved results approaching the DNS and experimental data (see the references above). These models have different treatments in linearizing the Reynolds stress, modelling the forcing, utilizing or not instantaneous DNS data and so on. Consequently, it is nearly impossible for the present work to extend and compare among all these models, especially considering the still larger parameter space in compressible flows. Thereby, we hope that the linear model in this work is both simple and instructive. Some possible extensions and future works are discussed in §§ 7 and 8.

2. Governing equations and dataset

2.1. Compressible Navier–Stokes equations

The canonical compressible turbulent channel flow is considered, as illustrated in figure 1. Assuming a calorically perfect gas, the governing NS equations are

(2.1a)\begin{gather} \frac{\partial \rho}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u}) = 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[ \mu ( {\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} {\boldsymbol{u}}^{\text{T}}}) \right] - \frac{2}{3} \boldsymbol{\nabla} ( \mu \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} ), \end{gather}

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

where $\rho$, $\boldsymbol {u}=[u,v,w]^ {\text {T}}$, $T$ and $p=\rho RT$ are the fluid density, velocity, temperature and pressure, respectively; $R$ and $c_p$ are the gas constant and isobaric specific heat; and $\mu$ and $\kappa$ are the molecular viscosity and thermal conductivity. Characteristic non-dimensional parameters are the Mach, Reynolds, Prandtl and Eckert numbers,

(2.2a–d)\begin{gather} {\textit{Ma}}_b = \dfrac{U_b}{a_w}, \quad Re_b = \dfrac{\rho_b U_b h}{\mu_w}, \quad {\textit{Pr}} = \dfrac{\mu c_p}{\kappa}=0.72, \quad {\textit{Ec}} = \dfrac{U_b^2}{c_p T_w} = (\gamma_0-1) {\textit{Ma}}_b^2, \end{gather}

where the subscript $w$ denotes quantities at the wall ($y={0,2h}$) and $h$ is the channel half-height. The speed of sound is $a=(\gamma _0 R T)^{1/2}$ with a constant specific heat ratio $\gamma _0=c_p/c_v=1.4$. The subscript 0 here is to distinguish from the LCS notation in § 2.3. The bulk density $\rho _b$ and bulk velocity $U_b$ are defined as $\rho _b = \int _{0}^{2h} {\bar {\rho } \, {\mathrm {d}} y}/(2h)$ and $\rho _b U_b = \int _{0}^{2h} {\overline {\rho u} \, {\mathrm {d}} y}/(2h)$, where the overbar denotes mean variables. The viscosity is calculated through Sutherland's law, where the fitting constant is 110.4 K and the reference temperature is 293.15 K. No-slip and isothermal walls are set on both sides as the boundary condition.

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

The basic variable set in (2.1) is ${\boldsymbol q}=[\rho,u,v,w,T]^ {\text {T}}$. Both the Reynolds average ($\bar {\varphi }$ with $\varphi$ a time-dependent variable) and Favre average ($\tilde {\varphi }=\overline {\rho \varphi }/\bar {\rho }$, $\varphi \ne \rho$) are deployed. The resulting two fluctuations are denoted as $\varphi ^ {\prime }$ and $\varphi ^ {{\prime \prime }}$, respectively. For wall-bounded turbulence, it is common to use wall viscous units with a superscript $+$, as $\boldsymbol {x}^+={\boldsymbol {x}}/{\delta _\nu }$, $\rho ^+=\rho /\rho _w$, $\boldsymbol {u}^+={\boldsymbol {u}}/{u_\tau }$ and $T^+={T}/{T_\tau }$. Here the viscous length unit is $\delta _\nu =\mu _w/(\rho _w u_\tau )$, the friction velocity and temperature are $u_\tau =(\tau _w/\rho _w)^{1/2}$ and $T_\tau =Q_w/(\rho _w c_p u_\tau )$, where $\tau _w$ and $Q_w$ are the wall mean shear and heat flux. The friction Reynolds number is $Re_\tau =h/\delta _v$. Furthermore, semi-local units are adopted, expressed with a superscript $*$ as $u_\tau ^*=(\tau _w/\bar {\rho })^{1/2}$, $\delta _\nu ^*=\bar {\mu }/(\bar {\rho } u_\tau ^*)$, so $y^*=y/\delta _v^*$ and $Re_\tau ^*=h/\delta _v^*$.

2.2. DNS dataset and mean flow calculation

A series of DNSs for compressible turbulent channel flows has been conducted by the present authors, as reported at length (Cheng & Fu 2022b). The calculations adopt two ${\textit {Ma}}_b$ of 0.8 and 1.5, and cover a $Re_b$ range of 7667–20 020 ($Re_\tau ={{436\unicode{x2013}1150}}$). The domains are all rectangular, in the same sizes $L_x\times L_z\times L_y=4{\rm \pi} h\times 2{\rm \pi} h\times 2h$. Previous studies have verified that this set-up can capture most of the large-scale motions in the outer region (Agostini & Leschziner 2014). Readers can refer to Cheng & Fu (2022b, 2023a) for more analyses on the DNS statistics. The case ${\textit {Ma}}_b=1.5$, $Re_b={20\ 020}$ ($Re_\tau ={1150}$, $Re_\tau ^*={780}$) is selected as the benchmark case for this study, which is of the highest ${\textit {Ma}}_b$ and $Re$ in our dataset. The valuable DNS data allow a quantitative evaluation of the linear model results.

In addition to DNS, well-established universal relations can be used to obtain efficiently the turbulent mean flows, by solving a set of ordinary differential equations (ODEs) (Griffin, Fu & Moin 2022, 2023; Chen et al. 2023; Song, Zhang & Xia 2023). The ODE-based method of Chen et al. (2023) is used here, which requires only the values of ${\textit {Ma}}_b$ and $Re_b$ as the input. We will demonstrate in § 6 that the ODE-based mean flow can be used in the linear models for SLSE, which will serve the parameter study.

2.3. Spectral linear stochastic estimation

As motivated in § 1, SLSE predicts the coherent portion of fluctuations at different heights as a reflection of the superposition effects. For consistency with the equations in § 3.1, the SLSE quantities are defined based on the Favre average. As the mean flow is homogeneous in the wall-parallel directions, Fourier decomposition is applied on ${\boldsymbol q}^ {{\prime \prime }}$ as

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

where $k_x$ and $k_z$ are the streamwise and spanwise wavenumbers, and $\hat {\boldsymbol q}^ {{\prime \prime }}$ is the shape function. In SLSE, the spectral signals $\hat {u}^ {{\prime \prime }}$ between the locations of measurement $y_m$ and prediction $y_p$ are connected through a complex-valued kernel function $H_L$ (Tinney et al. 2006),

(2.4a)\begin{gather} \left.\begin{array}{c@{}} \hat{u}_p^{{\prime\prime}} (y_p, t; k_x, k_z) = H_{L,uu} (y_p,y_m; k_x,k_z) \hat{u}^{{\prime\prime}} (y_m,t; k_x,k_z),\\ H_{L,uu} (y_p, y_m; k_x, k_z) = \dfrac{\langle \hat{u}^{{\prime\prime}} (y_p, t; k_x, k_z) \hat{u}^{{{\prime\prime}}{\dagger}} (y_m, t; k_x, k_z) \rangle}{\langle \hat{u}^{{\prime\prime}} (y_m, t; k_x, k_z) \hat{u}^{{{\prime\prime}}{\dagger}} (y_m, t; k_x, k_z) \rangle}, \end{array}\right\} \end{gather}

where $\langle {\cdot } \rangle$ stands for ensemble average and superscript ${\dagger}$ denotes the complex conjugate. If the time-independent $H_{L,uu}$ is obtained, then $\hat {u}^ {{\prime \prime }}(y_p,t;k_x,k_z)$ and $u^ {{\prime \prime }}(x,y_p,z,t)$ at $y_p$ can be estimated ($u_p^ {{\prime \prime }}$ or $\hat {u}_p^ {{\prime \prime }}$) given an instantaneous field ($u^ {{\prime \prime }}$ or $\hat {u}^ {{\prime \prime }}$) at $y_m$ considering all coherent motions (Baars et al. 2016). In the following, $k_x$ and $k_z$ in the brackets will be omitted if there is no ambiguity. It is worth mentioning that $H_{L,uu}$ for compressible flows is defined differently among the literature. The $u^ {{\prime \prime }}$ in (2.4a) can be replaced with $u^ {\prime }$ or other density-weighting forms. Cheng & Fu (2022b) show that different density-weighting approaches do not alter the energy distribution among scales for the benchmark case here. Also, with or without density weighting negligibly affects the streamwise and spanwise length scales of the pre-multiplied spectra peaks. Therefore, the definition in (2.4a) is used throughout. Similar to that for $u^ {{\prime \prime }}$, SLSE can be deployed to estimate $T^ {{\prime \prime }}(y_p)$ given $T^ {{\prime \prime }}(y_m)$ or $u^ {{\prime \prime }}(y_m)$. The latter combination concerns the coupling effects between $u^ {{\prime \prime }}$ and $T^ {{\prime \prime }}$ (Cheng & Fu 2023a). The two kernel functions are defined as

(2.4b)\begin{equation} H_{L,TT} (y_p, y_m) = \dfrac{\langle \hat{T}^{{\prime\prime}} (y_p, t) \hat{T}^{{{\prime\prime}}{\dagger}} (y_m, t) \rangle}{\langle \hat{T}^{{\prime\prime}} (y_m, t) \hat{T}^{{{\prime\prime}}{\dagger}} (y_m, t) \rangle}, \quad H_{L,uT} (y_p, y_m) = \dfrac{\langle \hat{T}^{{\prime\prime}} (y_p, t) \hat{u}^{{{\prime\prime}}{\dagger}} (y_m, t) \rangle}{\langle \hat{u}^{{\prime\prime}} (y_m, t) \hat{u}^{{{\prime\prime}}{\dagger}} (y_m, t) \rangle}. \end{equation}

Following Tinney et al. (2006), the amplitude of $H_{L,uu}$ (also for $H_{L,TT}$, $H_{L,uT}$) is further decomposed into two parts, as

(2.5)\begin{equation} |H_{L,uu}(y_p, y_m)|^2 = A_{pm,uu}^2 (y_p, y_m) \gamma_{uu}^2 (y_p, y_m), \end{equation}

where $A_{pm}(y_p,y_m;k_x,k_z)$ measures the fluctuation amplitude ratio at a specific length scale, and $\gamma ^2(y_p,y_m;k_x,k_z)$ is the 2-D LCS; their expressions are

(2.6a)\begin{equation} A_{pm,uu}(y_p,y_m) = \sqrt{\frac{\langle|\hat{u}^{{\prime\prime}}(y_p)|^2 \rangle}{\langle|\hat{u}^{{\prime\prime}}(y_m)|^2\rangle}}, \quad \gamma_{uu}^2 (y_p, y_m) = \frac{|\langle \hat{u}^{{\prime\prime}} (y_p, t) \hat{u}^{{{\prime\prime}}{\dagger}} (y_m, t) \rangle|^2}{\langle |\hat{u}^{{\prime\prime}} (y_p, t)|^2 \rangle \langle |\hat{u}^{{\prime\prime}} (y_m, t)|^2 \rangle}. \end{equation}

The denominator of $\gamma _{uu}^2$ comprises two individual 2-D energy spectra of $u^ {{\prime \prime }}$ at $y_p$ and $y_m$, while the numerator is the absolute value of the $u^ {{\prime \prime }}$ cospectrum between the two heights. Therefore, $\gamma _{uu}^2$ reflects the maximum correlation coefficient at each Fourier scale, hence instructive for understanding structural coherence. By definition, $0\le \gamma ^2 \le 1$; $\gamma ^2=1$ represents perfect coherence while $\gamma ^2=0$ means no coherence. The phase of $H_{L,uu}$ is not reflected in (2.5), which also matters as a measure of the structural angle difference in the $x$–$z$ plane (Baars et al. 2016). Similar to (2.6a), $\gamma _{TT}^2$ and $\gamma _{uT}^2$ are defined as

(2.6b)\begin{equation} \gamma_{TT}^2 (y_p, y_m) = \dfrac{|\langle \hat{T}^{{\prime\prime}} (y_p, t) \hat{T}^{{{\prime\prime}}{\dagger}} (y_m, t) \rangle|^2}{\langle |\hat{T}^{{\prime\prime}} (y_p, t)|^2 \rangle \langle |\hat{T}^{{\prime\prime}} (y_m, t)|^2 \rangle}, \quad \gamma_{uT}^2 (y_p, y_m) = \dfrac{|\langle \hat{T}^{{\prime\prime}} (y_p, t) \hat{u}^{{{\prime\prime}}{\dagger}} (y_m, t) \rangle|^2}{\langle |\hat{T}^{{\prime\prime}} (y_p, t)|^2 \rangle \langle |\hat{u}^{{\prime\prime}} (y_m, t)|^2 \rangle}. \end{equation}

Calculating $H_L$, $\gamma ^2$ and $A_{pm}^2$ requires the 2-D cospectrum of $\hat {\boldsymbol q}^ {{\prime \prime }}$,

(2.7)\begin{equation} \boldsymbol{\varPhi}(y,y^{\prime};k_x,k_z) = \langle \hat{\boldsymbol q}^{{\prime\prime}} (y,t;k_x,k_z)\hat{\boldsymbol q}^{{{\prime\prime}}{H}} (y^{\prime},t;k_x,k_z) \rangle, \end{equation}

where $H$ denotes the Hermite transpose. The cospectrum is computed using a time series of instantaneous measurements or numerical data. As an alternative, it can be obtained through a linear model based on the linearized (2.1), as introduced in § 1. The construction of the linear model is discussed below.

3. Linear models for compressible flows

Although the linear models in incompressible flows have been extensively addressed, their compressible counterparts deserve a careful discussion.

3.1. Mean flow and fluctuation equations

To establish a linear model, the governing equation for ${\boldsymbol q}^ {{\prime \prime }}$ is derived first. Following standard procedures, the time-averaged (2.1) is

(3.1a)\begin{gather} \frac{\partial \bar{\rho}}{\partial t} + \tilde{u}_j \frac{\partial \bar{\rho}}{\partial x_j} + \bar{\rho} \frac{\partial \tilde{u}_j}{\partial x_j} = 0, \end{gather}

(3.1b)\begin{equation}\bar{\rho} \left( \frac{\partial \tilde{u}_i}{\partial t} + \tilde{u}_j \frac{\partial \tilde{u}_i}{\partial x_j} \right) ={-}\frac{\partial \bar{p}}{\partial x_i} + \frac{\partial \bar{\tau}_{ij}}{\partial x_j} + \frac{\partial}{\partial x_j} \underbrace{(-\bar{\rho} {\widetilde{u_i^{{\prime\prime}}} u_j^{\prime\prime}} )}_\text{Reynolds stress}, \quad i=1,2,3, \end{equation}

(3.1c) \begin{align} \bar{\rho} c_v \left( \frac{\partial \tilde{T}}{\partial t} + \tilde{u}_j \frac{\partial \tilde{T}}{\partial x_j} \right) &= \frac{\partial}{\partial x_j} \left(\tilde{\kappa} \frac{\partial \tilde{T}}{\partial x_j}\right) - \bar{p} \frac{\partial \tilde{u}_j}{\partial x_j} + \bar{\tau}_{ij} \frac{\partial \tilde{u}_i}{\partial x_j} + \frac{\partial}{\partial x_j} \underbrace{(-\bar{\rho} c_v {\widetilde {u_j^{{\prime\prime}} T^{{\prime\prime}}}})}_{\substack{\text{turbulent} \\ \text{heat flux}}} \nonumber\\ &\quad - \underbrace{\overline{p^{\prime} \frac{\partial u_j^{{\prime\prime}}}{\partial x_j}}}_{\substack{\text{pressure} \\ \text{dilatation}}} + \underbrace{\overline{\tau_{ij}^{\prime} \frac{\partial u_i^{{\prime\prime}}}{\partial x_j}}}_{\substack{\text{viscous} \\ \text{dissipation}}}. \end{align}

Here the mean pressure is $\bar {p}=\bar {\rho }R\tilde {T}$, and the viscous stress (and also for $\tilde {\kappa }$) is approximated, following the suggestion by Gatski & Bonnet (2013), as

(3.2)\begin{equation} \left. \begin{array}{c@{}} \bar{\tau}_{ij} \approx \tilde{\mu} \left( \dfrac{\partial \tilde{u}_i}{\partial x_j} + \dfrac{\partial \tilde{u}_j}{\partial x_i} - \dfrac{2}{3} \dfrac{\partial \tilde{u}_k}{\partial x_k} \delta_{ij} \right), \\ \tau_{ij}^{\prime} \approx \tilde{\mu} \left( \dfrac{\partial u_i^{{\prime\prime}}}{\partial x_j} + \dfrac{\partial u_j^{{\prime\prime}}}{\partial x_i} - \dfrac{2}{3} \dfrac{\partial u_k^{{\prime\prime}}}{\partial x_k} \delta_{ij} \right) + \mu^{{\prime\prime}} \left( \dfrac{\partial \tilde{u}_i}{\partial x_j} + \dfrac{\partial \tilde{u}_j}{\partial x_i} - \dfrac{2}{3} \dfrac{\partial \tilde{u}_k}{\partial x_k} \delta_{ij} \right), \end{array} \right\} \end{equation}

where $\tilde {\mu }\approx \mu (\tilde {T})$, $\mu ^ {{\prime \prime }}\approx ({\partial \mu }/{\partial T}) T^ {{\prime \prime }}$, and $\delta _{ij}$ is the Kronecker delta. Equation (3.1) uses the Favre-averaged quantities, so it can be regarded as a governing equation for the mean flow $\tilde {\boldsymbol q}=[\bar {\rho },\tilde {u}_i,\tilde {T}]^ {\textrm {T}}$, though there are unclosed terms regarding ${\boldsymbol q}^ {{\prime \prime }}=[{\rho }^ {\prime },u_i^ {{\prime \prime }},T^ {{\prime \prime }}]^ {\textrm {T}}$. The physical meanings of these unclosed terms, as underbraced, have been elaborated before (e.g. Huang et al. 1995). In particular, the last two terms in (3.1c) are crucial for the energy transfer between the internal energy and turbulent kinetic energy $k\equiv {\widetilde {u_j^ {{\prime \prime }} u_j^ {{\prime \prime }}}}/2$ (Lele 1994). For reference, the governing equation for $k$ is provided as

(3.3)\begin{equation} \bar{\rho} \left({ \frac{\partial k}{\partial t} + \tilde{u}_j \frac{\partial k}{\partial x_j} }\right) ={-}\bar{\rho} {\widetilde{u_i^{{\prime\prime}} u_j^{{\prime\prime}}}} \frac{\partial \tilde{u}_i}{\partial x_j} - \frac{\partial}{\partial x_j} (\bar{\rho} {\widetilde{u_j^{{\prime\prime}}} k^{{\prime\prime}}} + \overline{p^{\prime} u_j^{{\prime\prime}}} - \overline{\tau_{ij}^{\prime} u_i^{{\prime\prime}}}) + \overline{p^{\prime} \frac{\partial u_j^{{\prime\prime}}}{\partial x_j}} - \overline{\tau_{ij}^{\prime} \frac{\partial u_i^{{\prime\prime}}}{\partial x_j}}. \end{equation}

The governing equations for ${\boldsymbol q}^ {{\prime \prime }}$ are obtained after a subtraction between (2.1) and (3.1). The fluctuating continuity equation is

(3.4a) \begin{equation} \left(\frac{\partial \rho^{\prime}}{\partial t} + \tilde{u}_j \frac{\partial \rho^{\prime}}{\partial x_j}\right) + \left( \bar{\rho} \frac{\partial u_j^{{\prime\prime}}}{\partial x_j} + \frac{\partial \tilde{u}_j}{\partial x_j} \rho^{\prime} + \frac{\partial \bar{\rho}}{\partial x_j} u_j^{{\prime\prime}} \right) = \underbrace{ - \frac{\partial \rho^{\prime} u_j^{{\prime\prime}}}{\partial x_j}}_{N_\rho^{{\prime\prime}}}. \end{equation}

The equation is rearranged so that the terms on the left-hand side are linear terms of $\boldsymbol {q}^ {{\prime \prime }}$, while the one on the right-hand side is nonlinear, denoted as $N_\rho ^ {{\prime \prime }}$. Similarly, the fluctuating momentum equation ($i=1,2,3$) takes the form of

(3.4b) \begin{align} & \bar{\rho} \left( \frac{\partial u_i^{{\prime\prime}}}{\partial t} + \tilde{u}_j \frac{\partial u_i^{{\prime\prime}}}{\partial x_j} \right) + ( \rho^{\prime} \tilde{u}_j + \bar{\rho} u_j^{{\prime\prime}}) \frac{\partial \tilde{u}_i}{\partial x_j} + \frac{\partial p^{\prime}}{\partial x_i} - \frac{\partial \tau_{ij}^{\prime}}{\partial x_j}\nonumber\\ & \quad ={-} \frac{\partial}{\partial x_j} \underbrace{[\bar{\rho} (u_i^{{\prime\prime}} u_j^{{\prime\prime}} - {\widetilde{u_i^{{\prime\prime}} u_j^{{\prime\prime}}}}) ]}_{\text{Reynolds stress fluctuation}} \underbrace{- \left( \frac{\partial \rho^{\prime} u_i^{{\prime\prime}}}{\partial t} + \frac{\partial \rho^{\prime} u_i^{{\prime\prime}} \tilde{u}_j}{\partial x_j} + \rho^{\prime} u_j^{{\prime\prime}} \frac{\partial \tilde{u}_i}{\partial x_j} \right)}_{N_{u_i}^{{\prime\prime}}}. \end{align}

Same as (3.4a), the linear and nonlinear terms are placed on the two sides. The nonlinear terms are classified into two groups. The first group is the fluctuation of the Reynolds stress tensor, and the second group collects all the second- and higher-order terms related to $\rho ^ {\prime }$. The classification in (3.4b) is not unique. For example, one can expand $N_{u_i}^ {{\prime \prime }}$ using (3.4a). Different possibilities will be discussed in § 7. Notably, when DNS is conducted, there is a temporarily-varying body force term in the streamwise momentum equation to fix the mass flux (e.g. Yao & Hussain 2020). This term is spatially uniform in the current DNS data, so it appears only at the scale $k_x=k_z=0$, hence not included in (3.4b). The fluctuating internal energy equation is

(3.4c) \begin{align} & \bar{\rho} c_v \left( \frac{\partial T^{{\prime\prime}}}{\partial t} + \tilde{u}_j \frac{\partial T^{{\prime\prime}}}{\partial x_j} \right) + c_v ( \rho^{\prime} \tilde{u}_j + \bar{\rho} u_j^{{\prime\prime}}) \frac{\partial \tilde{T}}{\partial x_j} - \frac{\partial}{\partial x_j} \left(\tilde{\kappa} \frac{\partial T^{{\prime\prime}}}{\partial x_j} + \kappa^{{\prime\prime}} \frac{\partial \tilde{T}}{\partial x_j}\right) \nonumber\\ & \quad + \left( \bar{p} \frac{\partial u_j^{{\prime\prime}}}{\partial x_j} + p^{\prime} \frac{\partial \tilde{u}_j}{\partial x_j} \right) - \left( \bar{\tau}_{ij} \frac{\partial u_i^{{\prime\prime}}}{\partial x_j} + \tau_{ij}^{\prime} \frac{\partial \tilde{u}_i}{\partial x_j} \right) ={-} \frac{\partial}{\partial x_j}\!\underbrace{[\bar{\rho} c_v ( u_j^{{\prime\prime}} T^{{\prime\prime}} - {\widetilde{u_j^{{\prime\prime}} T^{{\prime\prime}}}} ) ]}_{\substack{\text{turbulent heat flux} \\ \text{fluctuation}}} \nonumber\\ & \quad - \underbrace{\left(p^{\prime} \frac{\partial u_j^{{\prime\prime}}}{\partial x_j} - \overline{p^{\prime} \frac{\partial u_j^{{\prime\prime}}}{\partial x_j}} \right)}_{\substack{\text{pressure dilatation} \\ \text{fluctuation}}} + \underbrace{\left(\tau_{ij}^{\prime} \frac{\partial u_i^{{\prime\prime}}}{\partial x_j} - \overline{\tau_{ij}^{\prime} \frac{\partial u_i^{{\prime\prime}}}{\partial x_j}} \right)}_{\substack{\text{dissipation rate} \\ \text{fluctuation}}}\underbrace{- c_v \left( \frac{\partial \rho^{\prime} T^{{\prime\prime}}}{\partial t} + \frac{\partial \rho^{\prime} T^{{\prime\prime}} \tilde{u}_j}{\partial x_j} + \rho^{\prime} u_j^{{\prime\prime}} \frac{\partial \tilde{T}}{\partial x_j} \right)}_{N_T^{{\prime\prime}}}. \end{align}

Likewise, $N_T^ {{\prime \prime }}$ collects all the terms related to $\rho ^ {\prime }$. There are additional terms of fluctuating turbulent heat flux, pressure dilatation and dissipation rate, as underbraced. The governing equation for $k^ {{\prime \prime }}\equiv (u_j^ {{\prime \prime }} u_j^ {{\prime \prime }}/2-k)$ (Huang et al. 1995) can be similarly derived, whose full expression is listed in Appendix A.

As seen, (3.4) is far more complicated than its incompressible counterpart for deriving a linear model. The nonlinear terms $N_\rho ^ {{\prime \prime }}$, $N_{u_i}^ {{\prime \prime }}$ and $N_T^ {{\prime \prime }}$ all result from $\rho ^ {\prime }$, not present in incompressible cases. The remaining four nonlinear terms are classic ones in the turbulence modelling theory.

3.2. Modelling issues and linearization

In general, there are two types of strategies to linearize (3.4). The first one, as adopted by Bae et al. (2020) and Dawson & McKeon (2020), is to collect all the second- and higher-order terms of ${\boldsymbol q}^ {{\prime \prime }}$, i.e. those underbraced in (3.4), into the nonlinear forcing term. The resulting form is termed the LNS equation here, written in a standard operator form as

(3.5)\begin{equation} \frac{\partial {\boldsymbol q}^{{\prime\prime}}}{\partial t} = \mathcal{L}_{LNS} {\boldsymbol q}^{{\prime\prime}} + {\boldsymbol f}_{LNS}^{{\prime\prime}}. \end{equation}

Here, the linear operator $\mathcal {L}$ is only mean flow related, and ${\boldsymbol f}^ {{\prime \prime }}$ collects all the nonlinear terms. For later use, the temporal derivative term is not included in $\mathcal {L}$. Notably, the expression of $\mathcal {L}_{LNS}$ is exactly the same as that in the linear stability theory for laminar flows (Mack 1984).

The second strategy is to utilize turbulence modelling relations for possible linearization of the nonlinear terms, which is the counterpart of the eddy-viscosity-enhanced model for incompressible flows. However, the extension to compressible flows is much more complicated with multiple combinations. More terms require linearization besides the Reynolds stress term, which, as will be shown later, can lead to difficulties in closing the equations. Our principle here is to avoid ad hoc fittings by ourselves and use only the relations from classic modelling theory, which have been widely verified.

Analogous to the incompressible case, the fluctuation of the Reynolds stress term is modelled using the linearized Boussinesq assumption,

(3.6)\begin{equation} - \bar{\rho} ({ u_i^{{\prime\prime}} u_j^{{\prime\prime}} - {\widetilde{u_i^{{\prime\prime}} u_j^{{\prime\prime}}}} }) = \mu_t \left({ \frac{\partial u_i^{{\prime\prime}}}{\partial x_j} + \frac{\partial u_j^{{\prime\prime}}}{\partial x_i} - \frac{2}{3} \frac{\partial u_m^{{\prime\prime}}}{\partial x_m} \delta_{ij} }\right) - \frac{2}{3} \bar{\rho} k^{{\prime\prime}} \delta_{ij}. \end{equation}

The turbulent heat flux is usually modelled using the SRA. The modified version by Huang et al. (1995) takes the form of

(3.7)\begin{equation} T^{{\prime\prime}} = \frac{1}{{\textit{Pr}}_t} \frac{\partial \tilde{T}}{\partial \tilde{u}} u^{{\prime\prime}}. \end{equation}

It assumes that $T^ {{\prime \prime }}$ also satisfies the mixing-length relation as $u^ {{\prime \prime }}$, and hints that $T^ {{\prime \prime }}$ is advection dominated since it is determined by $u^ {{\prime \prime }}$ and the mean-flow gradients. Based on (3.7), there are two means to linearize the turbulent heat flux term. The first one is to follow the Reynolds-averaged NS (RANS) modelling, e.g. the $k$–$\omega$ model by Wilcox (2006), which extends (3.7) to all spatial directions. The resulting form is

(3.8a)\begin{equation} - \bar{\rho} {\widetilde{u_j^{{\prime\prime}} T^{{\prime\prime}}}} = \frac{\mu_t}{{\textit{Pr}}_t} \frac{\partial \tilde{T}}{\partial x_j} \rightarrow - \bar{\rho} (u_j^{{\prime\prime}} T^{{\prime\prime}} - {\widetilde{u_j^{{\prime\prime}} T^{{\prime\prime}}}}) = \frac{\mu_t}{{\textit{Pr}}_t} \frac{\partial T^{{\prime\prime}}}{\partial x_j}. \end{equation}

The second approach is a direct application of (3.7), as

(3.8b)\begin{equation} - \bar{\rho} (u_j^{{\prime\prime}} T^{{\prime\prime}} - {\widetilde{u_j^{{\prime\prime}} T^{{\prime\prime}}}}) ={-} \frac{1}{{\textit{Pr}}_t} \frac{\partial \tilde{T}}{\partial \tilde{u}} \bar{\rho} (u^{{\prime\prime}} u_j^{{\prime\prime}} - {\widetilde{u^{{\prime\prime}} u_j^{{\prime\prime}}}} ), \end{equation}

which is evaluated using (3.6). Consequently, the heat flux fluctuation is independent of $T^ {{\prime \prime }}$ but related to $u_i^ {{\prime \prime }}$ instead. Equation (3.8a) is deployed as we prefer to treat $T^ {{\prime \prime }}$ as an independent variable. The model using (3.8b) will be discussed in § 7.

Afterward, we focus on the pressure dilatation and dissipation terms in (3.4c). These two also appear in the $k^ {{\prime \prime }}$ equation. Little can we learn from the modelling theory on the linearization of pressure dilatation. The budget analysis for $k$ shows that the pressure dilatation and dilatational dissipation in (3.3) are negligible in the supersonic channel cases (Huang et al. 1995). So intuitively, their fluctuations are also negligible. The remaining part is the solenoidal dissipation, denoted as $\bar {\rho }\epsilon$. It can be evaluated from, e.g. the $k$–$\epsilon$ model (Jones & Launder 1972), as $\bar {\rho } \epsilon =\bar {\rho }^2 C_\mu {k^2}/{\mu _t}$ with $C_\mu$ a constant, thus linearized in terms of $k^ {{\prime \prime }}$. Note that $k^ {{\prime \prime }}$ cannot be expressed as a linear function of $\hat {\boldsymbol q}^ {{\prime \prime }}$, so it needs to be treated as an independent variable, which means that the $k^ {{\prime \prime }}$ equation needs to be included in (3.4) to solve the system of six equations. We prefer not to introduce $k^ {{\prime \prime }}$ as an independent variable for three reasons. First, there are difficulties in deriving a physically reasonable energy norm for the fluctuation, including $k^ {{\prime \prime }}$ (see § 3.3). Second, more terms (e.g. the diffusion term of $k^ {{\prime \prime }}$) appear in the $k^ {{\prime \prime }}$ equation, and their linearization adds to the uncertainty and complexity of the linear model before a careful assessment. Third, including the $k^ {{\prime \prime }}$ equation requires the profile of $k$ as the input, which is harder to obtain than $\tilde {\boldsymbol q}$. As an alternative to exclude $k^ {{\prime \prime }}$, the turbulent production and dissipation terms can be combined in the spirit of algebraic models (see Appendix A), and the former can be calculated using (3.6). The effects of including $k^ {{\prime \prime }}$ will be discussed in § 7.

Finally, there are nonlinear terms $N_\rho ^ {{\prime \prime }}$, $N_{u_i}^ {{\prime \prime }}$ and $N_T^ {{\prime \prime }}$ in (3.4), with no counterparts in the modelling theory. These three terms are all related to $\rho ^ {\prime }$, so in the spirit of Morkovin's hypothesis, they are of secondary importance. Therefore, they are not linearized but included in the nonlinear forcing. Some supporting evidence will be provided in § 7.

Following the nomenclature for incompressible flows, the linear model utilizing turbulence modelling is termed the eLNS (‘e’ for eddy-viscosity-enhanced) model. By collecting the residual nonlinearity into the forcing term, the operator form is

(3.9)\begin{equation} \frac{\partial {\boldsymbol q}^{{\prime\prime}}}{\partial t} = \mathcal{L}_{eLNS} {\boldsymbol q}^{{\prime\prime}} + {\boldsymbol f}_{eLNS}^{{\prime\prime}}. \end{equation}

If using (3.6) and (3.8a) and combining the turbulent production and dissipation terms, then $\mathcal {L}_{eLNS}$ is in the same form as $\mathcal {L}_{LNS}$ except for two substitutions,

(3.10a,b)\begin{equation} \tilde{\mu} \rightarrow \tilde{\mu} + \mu_t, \quad \tilde{\kappa} \rightarrow \tilde{\kappa} + \kappa_t = \frac{c_p \tilde{\mu}}{{\textit{Pr}}} + \frac{c_p \mu_t}{{\textit{Pr}}_t}, \end{equation}

where $\kappa _t$ is the eddy diffusivity. Equation (3.10a,b) is used likewise in the algebraic models, and is one of the simplest eLNS models widely used in previous works (Alizard et al. 2015; Pickering et al. 2021; Chen et al. 2023), though the derivation and justification were not elaborated before. Detailed expressions of $\mathcal {L}_{LNS}$ and $\mathcal {L}_{eLNS}$ can be found from Chen et al. (2023). From the mathematical point of view, using (3.10a,b) means intensified damping effects due to turbulence, especially in the outer region where $\mu _t\gg \tilde {\mu }$. This point is crucial in understanding the properties of the linear operators in § 4.

Next, we address the calculation of $\mu _t$ and ${\textit {Pr}}_t$. From the Boussinesq assumption, $\mu _t$ is computed as $\mu _t=-\bar {\rho } {\widetilde {u^ {{\prime \prime }} v^ {{\prime \prime }}}}/({\partial \tilde {u}}/{\partial y})$. Also, ${\textit {Pr}}_t$ is obtained from its definition using $ {\widetilde {u^ {{\prime \prime }} v^ {{\prime \prime }}}}$ and $ {\widetilde {v^ {{\prime \prime }} T^ {{\prime \prime }}}}$. If the fluctuation statistics are not available, $\mu _t$ can be evaluated from the mean streamwise momentum (3.1b) as

(3.11)\begin{equation} 0 ={-}\frac{{\mathrm{d}} \bar{p}}{{\mathrm{d}} x} + \frac{{\mathrm{d}}}{{\mathrm{d}} y} \left[ (\tilde{\mu} + \mu_t) \frac{{\mathrm{d}} \tilde{u}}{{\mathrm{d}} y} \right] \rightarrow \frac{\mu_t}{\tilde{\mu}_w} = \frac{1 - y^+{/}Re_\tau}{{{\mathrm{d}} \tilde{u}^+}/{{\mathrm{d}} y^+}} - \frac{\tilde{\mu}}{\tilde{\mu}_w}. \end{equation}

For more general flows where (3.11) is inapplicable, $\mu _t$ can be estimated using, e.g. the algebraic models. Also, ${\textit {Pr}}_t$ can be simply assumed a constant 0.9. In this work, $\mu _t$ and ${\textit {Pr}}_t$ are determined from the DNS statistics for the benchmark case results. Equation (3.11) and the constant ${\textit {Pr}}_t$ are adopted in § 6 for the parameter study.

3.3. Response to stochastic forcing

Ensemble-averaged variable correlation $\boldsymbol {\varPhi }$ (see (2.7)) is required to compute $H_L$, $\gamma ^2$ and $A_{pm}^2$. Based on § 3.2, $\boldsymbol {\varPhi }$ can be obtained by solving a stochastically forced linear system using only mean flow input (Madhusudanan et al. 2019). As in (2.3), Fourier decomposition is applied on $\boldsymbol {f}^ {{\prime \prime }}$ for the component $\hat {\boldsymbol f}^ {{\prime \prime }}$. After substituting (2.3) into (3.5) or (3.9), the equation for a single mode with $k_x\ne 0$ or $k_z\ne 0$ is

(3.12)\begin{equation} \frac{\partial \hat{\boldsymbol q}^{{\prime\prime}}}{\partial t} = \hat{\mathcal{L}} \hat{\boldsymbol q}^{{\prime\prime}} + \hat{\boldsymbol f}^{{\prime\prime}}, \end{equation}

where $\hat {\mathcal {L}}(\tilde {\boldsymbol q},k_x,k_z)$ can be either $\hat {\mathcal {L}}_{LNS}$ or $\hat {\mathcal {L}}_{eLNS}$, $\hat {\boldsymbol f}^ {{\prime \prime }} = [\,\hat {f}_\rho ^ {{\prime \prime }},\hat {f}_u^ {{\prime \prime }},\hat {f}_v^ {{\prime \prime }},\hat {f}_w^ {{\prime \prime }},\hat {f}_T^ {{\prime \prime }}]^ {\textrm {T}}$ is the forcing term and $\hat {\boldsymbol q}^ {{\prime \prime }}$ is the response. The energy norm of $\hat {\boldsymbol q}^ {{\prime \prime }}$ is defined as (Chu 1965)

(3.13)\begin{equation} \|\hat{\boldsymbol q}^{{\prime\prime}}\|^2 = ( \hat{\boldsymbol q}^{{\prime\prime}}, \hat{\boldsymbol q}^{{\prime\prime}} )_E = \int_{0}^{2h} { \left( \bar{\rho} \hat{\boldsymbol u}^{{{\prime\prime}}{H}} \hat{\boldsymbol u}^{{\prime\prime}} + \frac{R \tilde{T}}{\bar{\rho}} \hat{\rho}^{{\prime}{\dagger}} \hat{\rho}^{\prime} + \frac{\bar{\rho} c_v}{\tilde{T}} \hat{T}^{{{\prime\prime}}{\dagger}} \hat{T}^{{\prime\prime}} \right) {\mathrm{d}} y} = \int_{0}^{2h} {\hat{\boldsymbol q}^{{{\prime\prime}}{H}} \boldsymbol{M} \hat{\boldsymbol q}^{{\prime\prime}} \, {\mathrm{d}} y}, \end{equation}

where $\boldsymbol {M}$ is the energy weight matrix.

When (3.12) is driven by a stochastic forcing, the response is also stochastic if the system is linearly stable. Following those for incompressible flows (Gupta et al. 2021), the forcing is $\hat {\boldsymbol f}=\boldsymbol {B}\hat {\boldsymbol f}_0$ where $\boldsymbol {B}(y)$ is a modelling matrix and $\hat {\boldsymbol f}_0$ is a white noise signal. Specifically, $\hat {\boldsymbol f}_0$ is assumed to be a $\delta$-correlated Gaussian white noise with zero mean,

(3.14)\begin{equation} \langle\, \hat{\boldsymbol f}_0^{{\prime\prime}} \rangle = 0, \quad \langle\, \hat{\boldsymbol f}_0^{{\prime\prime}} (y,t) \hat{\boldsymbol f}_0^{{{\prime\prime}}{H}} (y^{\prime},t^{\prime}) \rangle = \boldsymbol{I} \delta_d (y - y^{\prime}) \delta_d (t - t^{\prime}), \end{equation}

where $\boldsymbol {I}$ is the identity matrix and $\delta _d$ is the Dirac function. The introduction of $\boldsymbol {B}(y)$ allows for varying forcing amplitude for different variables at different $y$. The specific expression of $\boldsymbol {B}$ will be provided later in (3.17). Mathematically, $\boldsymbol {\varPhi }$ is obtained by solving the algebraic Lyapunov equation (see Farrell & Ioannou (1993) for more details) as

(3.15)\begin{equation} \hat{\mathcal{L}} \boldsymbol{\varPhi} + \boldsymbol{\varPhi} \hat{\mathcal{L}}^{\dagger}{+} \boldsymbol{B} \boldsymbol{B}^{H} = \boldsymbol{0}. \end{equation}

Note that $\hat {\mathcal {L}}$ contains wall-normal derivatives, so $\hat {\mathcal {L}}^{\dagger}$ is defined in an adjoint manner. Also, the vector product is with respect to (3.13). Equation (3.15) is discretized using the Chebyshev collocation point method. By default, $N_y=301$ points are used, abundant to ensure grid independence (Bae et al. 2020). The Lyapunov equation is solved using the Matlab software. For boundary conditions, a no-slip, isothermal wall is assumed on both sides. Therefore, the wall fluctuations satisfy $\hat {u}_w^ {{\prime \prime }}=\hat {v}_w^ {{\prime \prime }}=\hat {w}_w^ {{\prime \prime }}=\hat {T}_w^ {{\prime \prime }}=0$, and $\hat {\rho }_w^ {\prime }$ is solved through the continuity equation. The solver verification can be found from Chen et al. (2023).

In the incompressible model of Madhusudanan et al. (2019), $\boldsymbol {B}$ is set to $\boldsymbol {I}$ for both LNS and eLNS, meaning that the stochastic forcing has uniform amplitude along the wall-normal direction. This is also adopted in previous compressible works (Alizard et al. 2022; Chen et al. 2023). To improve the estimation below the logarithmic region, Gupta et al. (2021) design two new models for incompressible flows, termed ‘W-model’ and ‘$\lambda$-model’, where $\boldsymbol {B}$ is $y$-dependent or even scale-dependent. The simpler W-model is considered here, where the forcing amplitude for the momentum equation varies in proportion to $\nu _t=\mu _t/\bar {\rho }$. For compressible flows, $\tilde {f}_\rho$ and $\tilde {f}_T$ are present and need to be modelled. Without comprehensive knowledge of the forcing statistics, we consider the simplest form here, which is instructive and free from ad hoc implementations. First, $\boldsymbol {B}$ is set to be diagonal, i.e. $\boldsymbol {B}={diag}( [B_\rho \ \nu _t\ \nu _t\ \nu _t\ B_T]^ {\textrm {T}} )$, where the amplitude functions $B_\rho$ and $B_T$ are to be determined. In this way, the off-diagonal terms, i.e. the cospectra of the forcing components and also the anisotropy (Foysi, Sarkar & Friedrich 2004), are not modelled. Second, $\boldsymbol {B}$ is assumed to be only related to the mean flow, independent of length scales $\lambda _x$ and $\lambda _z$. The extended SRA (3.7) and the DNS data of Coleman et al. (1995) suggest that the following relations are nearly satisfied:

(3.16a,b)\begin{equation} T^{{\prime\prime}}_{rms} = \frac{1}{{\textit{Pr}}_t} \left|\frac{\partial \tilde{T}}{\partial \tilde{u}}\right| u^{{\prime\prime}}_{rms}, \quad \rho^{\prime}_{rms} = \frac{1}{{\textit{Pr}}_t} \left|\frac{\partial \bar{\rho}}{\partial \tilde{u}}\right| u^{{\prime\prime}}_{rms} \approx \frac{\bar{\rho}}{{\textit{Pr}}_t\tilde{T}} \left|\frac{\partial \tilde{T}}{\partial \tilde{u}}\right| u^{{\prime\prime}}_{rms}, \end{equation}

where the subscript ${rms}$ is for root mean square, and the approximation is due to the slightly varying mean pressure. The mean flow gradients appear in their absolute values because the fluctuation r.m.s. is positive. For channel flows, ${\partial \tilde {T}}/{\partial \tilde {u}}>0$ and ${\partial \bar {\rho }}/{\partial \tilde {u}}<0$ hold in the whole field. From (3.16a,b), the forcing amplitude is modelled as

(3.17)\begin{equation} \boldsymbol{B}={diag}\left( \left[\frac{1}{{\textit{Pr}}_t} \frac{\partial \bar{\rho}}{\partial \tilde{u}} \nu_t,\ \nu_t,\ \nu_t,\ \nu_t,\ \frac{1}{{\textit{Pr}}_t} \frac{\partial \tilde{T}}{\partial \tilde{u}} \nu_t\right]^{\text{T}} \right). \end{equation}

Since the SLSE quantities in § 2.3 contain only the ratios of cospectra, multiplying $\boldsymbol {B}$ by a non-zero constant does not affect the SLSE results. The examination for (3.16a,b) is shown in figure 2 using the present DNS data (benchmark case). Both $\rho ^ {\prime }_{rms}$ and $T^ {{\prime \prime }}_{rms}$ from (3.16a,b) are in line with DNS, though $\rho ^ {\prime }_{rms}$ is slightly underestimated. Notably, the fact that $\rho ^ {\prime }_{rms}$ and $T^ {{\prime \prime }}_{rms}$ are primarily determined by $u^ {{\prime \prime }}_{rms}$ and their mean flow gradients hints that they are advection dominated. The larger deviation in the outer region indicates more pronounced thermodynamic processes where the mean flow gradients are small and the local Mach number is high.

Figure 2. Wall-normal distributions of the (a) density and (b) temperature fluctuations (r.m.s.) from the DNS data and (3.16a,b). Note that $u^ {{\prime \prime }}_{rms}$ is also from DNS.

In summary, two linear models used in this work are (1) the LNS model (3.5) and (2) the eLNS model ((3.9) and (3.10a,b)). In both cases, (3.17) is used for variable control. There is also a notation ‘eLNS-CD’ below, to be introduced in § 5.4. In § 3.1, we use the Favre averages to derive the linear model, so we should also use $\tilde {\boldsymbol q}$ to calculate $\hat {\mathcal {L}}$ and then $\boldsymbol {\varPhi }$. As ${\textit {Ma}}_b$ here is relatively low, the difference between $\tilde {\boldsymbol q}$ and $\bar {\boldsymbol q}$ is very small (Huang et al. 1995; Cheng & Fu 2022b). We have compared $\boldsymbol {\varPhi }$ using the two mean flows, and the relative difference is less than 1 %. Thereby, for ease of later comparison with the ODE-based mean flow (§ 6), we will use $\bar {\boldsymbol q}$ to compute $\tilde {\mathcal {L}}$ and $\boldsymbol {\varPhi }$ throughout.

4. Properties of the linear operator

As introduced in § 3.2, the linear operators $\hat {\mathcal {L}}_{LNS}$ and $\hat {\mathcal {L}}_{eLNS}$ were considered in previous works, but their mathematical features were not elaborated. The framework in § 3.3 requires a globally stable system, and the response amplification depends on the non-normality of the linear operator (Trefethen 1997). Thereby, we scrutinize the eigenspectra and pseudospectra of the linear operators in this section. The categorization of different modes is also discussed. Significant differences will be noted from the incompressible case, and the discussions have important implications for later SLSE results.

4.1. Spectra of the linear operator

From (3.12), the eigenvalues of $\hat {\mathcal {L}}$ reflect the temporal modal instability of $\hat {\boldsymbol q}^ {{\prime \prime }}$, defined as $\hat {\mathcal {L}}\check {\boldsymbol q}^ {{\prime \prime }}=- {\mathrm {i}}\omega \check {\boldsymbol q}^ {{\prime \prime }}$. Here $\omega =\omega _r+ {\mathrm {i}}\omega _i$ is the eigenvalue with $\omega _r$ the frequency and $\omega _i$ the growth rate, and $\check {\boldsymbol q}^ {{\prime \prime }}$ is the right-eigenfunction. The eigenspectra of $\hat {\mathcal {L}}_{LNS}$ and $\hat {\mathcal {L}}_{eLNS}$ are displayed in figures 3(a) and 3(b) at $\lambda _x=8h$ and $\lambda _z=2h$. As will be shown in § 5, the LCS at this $(\lambda _x,\lambda _z)$ is approximately the highest within the length scales considered. All the eigenmodes in figure 3 are stable ($\omega _i<0$), meeting the requirement of a globally stable system. The distributions of $\omega _r$ are similar between the LNS and eLNS spectra. A prominent difference is that the damping rates $|\omega _i|$ of the eLNS modes are nearly two orders of magnitude higher. This is because after the substitutions in (3.10a,b), the damping second-order derivative terms are highly increased. The proof from the growth-rate decomposition is provided in Appendix C to support the explanation.

Figure 3. (a,b) Spectra, (c,d) pseudospectra ($\log _{10}\varepsilon$), and (e,f) the maximum response gain of the linear operators for (a,c,e) LNS and (b,d,f) eLNS in the benchmark case ($\lambda _x=8h$, $\lambda _z=2h$). The vertical dashed lines in (a,b,e,f) are from (4.1), and the black dots in (c,d) denote the eigenvalues.

Next, the categorization of eigenmodes is discussed, which is crucial for interpreting later linear-model results. The classic four mode branches are identified in both LNS and eLNS spectra, namely the vortical, entropy, and fast and slow acoustic branches, as denoted in different colours (Balakumar & Malik 1992). Notably, the term vortical branch hereinafter is not restricted to vortices but a general term for all kinematically dominated modes, analogous to the incompressible case. The reference frequencies of four branches, as labelled in figure 3, are obtained by solving the fluctuations on an inviscid uniform mean flow, as

(4.1)\begin{equation} \left. \begin{array}{ll@{}} \omega_{r,{vort}} = \omega_{r,{ent}} = \tilde{u}_c k_x & \text{for vortical, entropy modes,} \\ \omega_{r,{acout}} = \tilde{u}_c \left( k_x \pm \dfrac{1}{{\textit{Ma}}_c}\sqrt{k_x^2+k_z^2} \right) & \text{for fast/slow acoustic modes,} \end{array} \right\} \end{equation}

where the subscript $c$ denotes the channel centre ($y=h$), analogous to the free stream in boundary layers. The derivation for (4.1) and corresponding eigenfunctions are detailed in Appendix B. In the free stream, the acoustic modes are isentropic and dilatational, the entropy modes are isobaric and static, and the vortical modes are solenoidal with zero thermodynamic components, so they can be easily distinguished. In the channel flow, however, these modes are more complicated due to the non-uniform mean flow and viscosity. For illustration, figure 4 plots the eigenfunctions of the least stable vortical, acoustic and entropy modes in the eLNS spectra, where the linearized pressure fluctuation $\check {p}^ {\prime }/\bar {p}=\check {\rho }^ {\prime }/\bar {\rho }+\check {T}^ {{\prime \prime }}/\tilde {T}$. The modes of lower $\omega _i$ exhibit common features. The acoustic modes (figure 4a) can be recognized from others based on $\omega _r$. Also, they have large pressure fluctuations, while the vortical and entropy modes are nearly isobaric. There are a series of acoustic modes in figure 3 with different wall-normal wavenumbers ($k_y$), and $\omega _{r,{acout}}$ in (4.1) represents the lower (upper) frequency limit with zero $k_y$ for the fast (slow) branches. In comparison, the vortical and entropy modes have close $\omega _r$, so they cannot be distinguished directly in figure 3. The vortical modes (figure 4b) have large velocity components, and the non-zero $\check {\rho }^ {\prime }$ and $\check {T}^ {{\prime \prime }}$ are mainly from the passive advection due to mean-flow gradients. Their entropy fluctuations are also non-zero since the mean flow is not isentropic. In contrast, the entropy modes (figure 4c) are dominated by the thermodynamic components, while their velocity components are not zero due to the viscous and diffusive terms. The vortical modes are not isentropic and the entropy modes are kinematic, so they can be coupled with each other and the modes with smaller $\omega _i$ in the vortical/entropy branch may not be as clearly distinguished as in figure 4. The eigenfunctions of different types of modes in the LNS model are similar in shapes to figure 4. The differences are contributed by the eddy terms related to $\mu _t$ and $\kappa _t$ (see Appendix B). The effects of different modes on the response behaviour will be further discussed in § 4.2.

Figure 4. Eigenfunctions of the linear operator in the eLNS model ($\lambda _x=8h$, $\lambda _z=2h$): (a) acoustic mode ($\omega h/U_b=3.312-0.038 {\mathrm {i}}$); (b) vortical mode ($\omega h/U_b=0.870-0.059 {\mathrm {i}}$) and (c) entropy mode ($\omega h/U_b=0.866-0.075 {\mathrm {i}}$).

4.2. Pseudospectra of the linear operator

All the eigenmodes of $\hat {\mathcal {L}}$ are asymptotically stable, so its transient behaviour is determined by external forcing and the non-normality of the linear operator (Trefethen 1997). The $\varepsilon$-pseudospectra is a quantitative measure of non-normality. By definition, a smaller $\varepsilon$ means that the eigenvalues are more sensitive to perturbations on $\hat {\mathcal {L}}$, so the non-normality is stronger. Also, in response to a harmonic external forcing of frequency $\omega _f$, $\varepsilon$ measures the largest gain of the response $G(\omega _f)=max\|(- {\mathrm {i}}\omega _f\boldsymbol {I}-\hat {\mathcal {L}})^{-1}\|$, which equals $1/\varepsilon _0$ with $\varepsilon _0$ the $\varepsilon$ at the imaginary axis (see Chomaz (2005), for more details). The pseudospectra and $G(\omega _f)$ of the two linear operators are depicted in figures 3(c)–3(f) to provide insights into the response characteristics. To distinguish from the resonant effects, the gain if $\hat {\mathcal {L}}$ is normal, $G_{normal}(\omega _f)$, is also plotted. For the LNS case, strong non-normality appears for the vortical modes since $G$ can be over one order of magnitude higher than $G_{normal}$, which is known to be an essential feature of the Orr–Sommerfeld operator. The non-normality of the entropy modes cannot be directly measured because their $\omega$ are close to the vortical modes. In comparison, the acoustic modes are of little non-normality as $G$ results mostly from the resonant effects. Therefore, we can expect that the behaviour of $\hat {\boldsymbol q}^ {{\prime \prime }}$ is largely determined by the vortical modes, similar to the incompressible case. Hanifi, Schmid & Henningson (1996) show that including the acoustic modes or not negligibly affects the transient growth results for a laminar flow (LNS case). A similar examination will be performed later for the turbulent mean flow in the stochastically driven system. When the linear operator is eddy viscosity enhanced (figures 3d and 3f), the non-normality of the acoustic modes is slightly affected, but that of the vortical (and possibly entropy) modes is severely weakened, reflected from the larger $\varepsilon$ and much smaller $G/G_{normal}$. Consequently, the acoustic (and also entropy) modes, and thus the compressibility effects, are expected to play more vital roles in the fluctuation evolution.

It is worth mentioning that Symon et al. (2021) and Kuhn et al. (2022) also report in their incompressible cases that eddy viscosity decreases the non-normality of the linear operator for structures of high $\lambda _x/\lambda _z$. Symon et al. (2021) conclude that this effect of counteracting non-normality in the eLNS model mitigates the trade-off between turbulent production and nonlinear energy transfer, leading to better prediction for streamwise-elongated structures than the LNS model. In the compressible case studied here, however, additional acoustic and entropy modes are present, so the weakened non-normality of the vortical modes can lead to more complicated consequences.

To support the above observations and further quantify the effects of acoustic and entropy modes, the energy amplification factor of the response $V$ is computed. Here, $V$ is defined as the variance of $\hat {\boldsymbol q}^ {{\prime \prime }}$, obtained as the trace of $\boldsymbol {\varPhi }$ with respect to the energy norm (Hwang & Cossu 2010). Two $V$ are needed, with acoustic (or entropy) modes included and excluded, respectively, to distinguish their contributions. Equation (3.15) contains the contributions of all modes, so a framework needs to be developed first that can give the cospectra considering only partial modes. Following that in the transient growth analysis (Hanifi et al. 1996), the response is restricted to the space spanned by partial eigenfunctions,

(4.2)\begin{equation} \hat{\boldsymbol q}^{{\prime\prime}} \approx \sum_{m=1}^{N_m} { c_m(t) \check{\boldsymbol q}_m^{{\prime\prime}} } = \check{\boldsymbol Q}^{{\prime\prime}} {\boldsymbol C}(t), \end{equation}

where $m$ and $N_m$ are the indexes and total number of the modes selected, and $c_m$ is a time-dependent coefficient. In matrix form, $\check {\boldsymbol Q}^ {{\prime \prime }}=[\check {\boldsymbol q}_1^ {{\prime \prime }}, \ldots, \check {\boldsymbol q}_{N_m}^ {{\prime \prime }}]$ and ${\boldsymbol C}=[c_1, \ldots, c_{N_m}]^ {\textrm {T}}$. The resulting cospectrum is $\boldsymbol {\varPhi } = \langle \hat {\boldsymbol q}^ {{\prime \prime }} \hat {\boldsymbol q}^{ {{\prime \prime }} {H}} \rangle = \check {\boldsymbol Q}^ {{\prime \prime }} \langle {\boldsymbol C} {\boldsymbol C}^ {H} \rangle \check {\boldsymbol Q}^{ {{\prime \prime }} {H}}$. Similar to (3.15), the governing equation for $\langle {\boldsymbol C} {\boldsymbol C}^ {H} \rangle$ is derived as

(4.3)\begin{equation} \boldsymbol{\varLambda} \langle {\boldsymbol C} {\boldsymbol C}^{H} \rangle + \langle {\boldsymbol C} {\boldsymbol C}^{H} \rangle \boldsymbol{\varLambda}^{H} + [\check{\boldsymbol Q}^{{{\prime\prime}}{H}} (\boldsymbol{B} \boldsymbol{B}^{H})^{{-}1} \check{\boldsymbol Q}^{{\prime\prime}}]^{{-}1} = \boldsymbol{0}, \end{equation}

where $\boldsymbol {\varLambda }={diag}([\omega _1, \ldots, \omega _{N_m}])$. Therefore, the response when only partial eigenmodes are included can be obtained by solving (4.3). Note that $\boldsymbol {B}\boldsymbol {B}^ {H}$ can be non-invertible as the element $\nu _t$ in (3.17) tends to zero at the wall. A small positive constant is added to the diagonal of $\boldsymbol {B}$ for solving (4.3) only to avoid the singularity. Also, one drawback of (4.3) is its inapplicability to the highly non-normal case where $\check {\boldsymbol Q}^{ {{\prime \prime }} {H}} \check {\boldsymbol Q}^ {{\prime \prime }}$ has a very large condition number, hence nearly non-invertible.

Figures 5(a) and 5(b) plot the distributions of $V$ for the LNS and eLNS models at different $\lambda _z$ ($\lambda _x=8h$). The $V$ in the LNS case is over one order of magnitude higher than that in the eLNS case, consistent with the stronger non-normality in figure 3(c). Note that the higher $V$ in the LNS case does not directly influence the SLSE quantities since they contain only the ratios of cospectra. The effects of acoustic modes are studied first. In figure 5(a), the energy growth with acoustic modes excluded ($V_{nacout}$) is quite close to $V$, supporting the observation in figures 3(c) and 3(e) that the acoustic modes are in secondary roles for the linear operator in terms of non-normality. As a quantitative measure, the ratio $r_{V_{na}}=V_{nacout}/V$ is plotted in percentage. In most of the $\lambda _z$ range displayed, the non-acoustic components contribute over 90 % of the total energy growth of the response, similar to the laminar case (Hanifi et al. 1996). For the eLNS model in figure 5(b), however, the non-acoustic components contribute only approximately 50 % of the total energy growth, and the contributions of acoustic modes are prominent, consistent with the inference from figures 3(d) and 3(f).

Figure 5. (a,b) Energy amplification factors for all modes ($V$) and when acoustic modes are excluded ($V_{nacout}$) (left axis) and their ratio $r_{V_{na}}=V_{nacout}/V$ (right axis) at different spanwise wavelengths ($\lambda _x=8h$), and the contours of (c,d) $r_{V_{na}}$ and (e) $r_{V_{ne}}=V_{nentrp}/V$ at different length scales. Panels (a,c) are for the LNS model and (b,d,e) are for eLNS.

To gain further knowledge on the scale dependence of the acoustic-mode effects, the contours of $r_{V_{na}}$ are depicted in figures 5(c) and 5(d) for the two models. The $\lambda _x$ and $\lambda _z$ ranges displayed are of interest for SLSE, as will be specified in § 5. The $r_{V_{na}}$ in the LNS case is over 80 % in most scale ranges, while that in the eLNS case is within 45 % to 60 %. More importantly, $r_{V_{na}}$ tends to increase at a larger scale in the LNS model, and is over 95 % when $\lambda _x>2h$ and $\lambda _z>0.5h$. In contrast, $r_{V_{na}}$ diminishes as $\lambda _x$ and $\lambda _z$ rise in the eLNS model, indicating that the large-scale structures are especially affected by the acoustic modes. In the same way, the effects of entropy modes are quantified in figure 5(e) for the eLNS case, using $r_{V_{ne}}=V_{nentrp}/V$ with $V_{nentrp}$ the energy growth excluding entropy modes. The $r_{V_{ne}}$ in the LNS case is not displayed as the matrix excluding the entropy modes is highly non-normal (see figure 3c). In the eLNS case, $r_{V_{ne}}$ is within 60 % to 80 % in most scale ranges, so the entropy modes are also non-negligible for the response growth, though their influences are weaker than the acoustic modes since $r_{V_{ne}}>r_{V_{na}}$. Note that the contours of $r_{V_{ne}}$ are not as smooth as $r_{V_{na}}$ in figure 5(d), due to the coupling effects between the vortical and entropy modes, as discussed in § 4.1. Nevertheless, the qualitative trend in figure 5(e) is not affected. The entropy-mode effect in the LNS case is not demonstrated, so its difference from the eLNS case is not reflected. The results in § 5 suggest that the acoustic components play a larger role than the entropy ones for SLSE in the present case, so the difference of $r_{V_{na}}$ (figures 5(c) and 5(d)) is important, and the acoustic components will be more focused on below.

We provide more interpretations on the trends in figure 5. For an inviscid uniform mean flow (free stream, $\tilde {\mu }=\mu _t=0$), the eigenfunctions of different modes in (4.1) are perpendicular to each other under the energy norm in (3.13), so $\hat {\mathcal {L}}$ is normal (George & Sujith 2011). It becomes non-normal in the presence of mean flow gradients and viscosity. The latter two effects are the strongest in the near-wall region, while the acoustic and entropy modes are prevailing in the outer region, so the non-normality between the vortical and other three branches is weak (figure 3c). Meanwhile, as $\lambda _x$ and $\lambda _z$ increase from small scales, the amplified structure moves away from the wall into the outer region. The decrease of mean flow gradients leads to less coupled mode branches and thus increasing $r_{V_{na}}$ (closer to the free stream condition), as observed in figure 5(c). When the model is $\mu _t$-enhanced (eLNS), however, the trends are reversed simply because $\mu _t$ is zero at the wall and largest in the outer region with $\mu _t\gg \tilde {\mu }$ (${\sim }70$ times in the present case). This means, from the view of equations, that the damping effects due to the modelled eddy terms are highly intensified ($\tilde {\mu }\rightarrow \tilde {\mu }+\mu _t$). Consequently, the larger-scale structure is more affected by $\mu _t$, and the dissipative second-order derivative terms for $\mu _t$ result in strong coupling between the acoustic, entropy modes and others. This is exactly what we observe in figures 5(d) and 5(e).

In summary, there are four branches of eigenmodes in both LNS and eLNS cases, and the vortical branch has strong non-normality to support the transient behaviour of the fluctuations (response) subject to stochastic forcing. For the LNS model, the acoustic modes have little non-normality and increasingly small contributions to the energy growth as the length scales (both $\lambda _x$ and $\lambda _z$) increase. In contrast, the non-normality of the vortical modes diminishes in the eLNS model due to the damping effects of $\mu _t$ on the eigenmodes. As a result, the acoustic and entropy modes can account for 20 % to 55 % of the energy growth; their contribution goes higher for larger-scale structures. The above is a significant difference between the LNS and eLNS models, not present in incompressible flows.

For the eLNS model, the fact that the acoustic and entropy modes, and thus compressibility effects, can have such a high contribution to the fluctuation growth deserves serious concern, considering the great similarities reported before to incompressible flows (see § 1). We examine below whether this prominent role of acoustic and entropy components in the eLNS model is an intrinsic defect introduced by turbulence modelling, through comparison with DNS. In fact, we show in § 5 that the prominent acoustic and entropy components can be problematic, especially for the density/temperature fields, and will seek remedies.

5. General flow structures and SLSE results

5.1. General flow structures

We first present an overview of the general structures of velocity and temperature fields in the benchmark case. Comprehensive analyses have been reported in previous works based on experiments and DNS for both channel and boundary-layer flows (Smith & Smits 1995; Ringuette et al. 2008; Pirozzoli & Bernardini 2011; Cogo et al. 2022; Cheng & Fu 2023a), so here we only provide a brief description to set the grounds for later discussions.

After the velocity transformation of Griffin et al. (2021), the mean streamwise velocity $\tilde {u}_{GFM}^+(y^*)$ matches the universal incompressible profile $\bar {u}_{incp}^+(y^+)$ within and below the logarithmic region, as shown in figure 6(d). Therefore, we can follow the criterion in incompressible flows to name different regions in terms of $y^*$ and $y/h$. Three wall-normal heights are selected, namely $y^*=11$ in the buffer layer (near-wall region), $y^*=107$ ($y=0.14h$) in the logarithmic layer and $y=0.5h$ ($y^*=385$) in the outer region, to display the instantaneous fields in different regions using DNS. Hence, the overall flow organization can be recognized. Notably, Abe & Antonia (2017) show in their incompressible case that $y\approx 0.5h$ is the region where large-scale structures of the streamwise velocity and temperature (as a passive scalar) fluctuations dominate. In the near-wall region (figures 6c and 6g), both $u^{ {{\prime \prime }}+}$ and $T^{ {{\prime \prime }}+}$ exhibit streamwise-elongated streaks, connected with the ‘sweep’ and ‘ejection’ events in the near-wall self-sustaining cycles (Hamilton, Kim & Waleffe 1995; Schoppa & Hussain 2002). Meanwhile, $u^{ {{\prime \prime }}+}$ and $T^{ {{\prime \prime }}+}$ strongly resemble each other morphologically. As a quantitative measure, Cheng & Fu (2023a) show that the correlation coefficient of $u^{ {{\prime \prime }}+}$ and $T^{ {{\prime \prime }}+}$, $R_{u^ {{\prime \prime }} T^ {{\prime \prime }}} = {\langle u^ {{\prime \prime }} T^ {{\prime \prime }} \rangle }/(u_{rms}^{ {{\prime \prime }}}T_{rms}^{ {{\prime \prime }}})$, is approximately one at $y^*\lesssim 12$. Further away from the wall into the logarithmic region (figures 6b and 6f), dominant velocity streaks are in larger length scales, while the streamwise-elongated feature is less obvious for $T^ {{\prime \prime }}$. In other words, $T^ {{\prime \prime }}$ is more isotropic than its near-wall counterpart, filled with jellyfish-shaped structures. The observation position in figures 6(a) and 6(e) is out of the logarithmic region. Large scale $u^ {{\prime \prime }}$ streaks are still identifiable, but streamwise $T^ {{\prime \prime }}$ streaks are hardly seen with obvious spanwise extensions. Cheng & Fu (2023a) demonstrate that the spanwise length scales of $u^ {{\prime \prime }}$ and $T^ {{\prime \prime }}$ corresponding to their pre-multiplied spectra peaks are actually comparable to each other throughout the channel; $T^ {{\prime \prime }}$ has a shorter peak streamwise length scale within and above the logarithmic region, thus tending to be more isotropic.

Figure 6. Overview of the turbulent field for the benchmark case: (a–c) instantaneous streamwise velocity fluctuation ($u^{ {{\prime \prime }}+}$), (d) transformed mean velocity profile and (e–g) instantaneous temperature fluctuation ($T^{ {{\prime \prime }}+}$); (a,e) are at $y=0.5h$ ($y^*=385$), (b,f) are at $y=0.14h$ ($y^*=107$) and (c,g) are at $y^*=11$ (top view).

In the following, we focus on predicting the coherent velocity and temperature fluctuations using SLSE based on the measurement at outer locations. Here $y_m^*$ is fixed at the centre of the logarithmic region (Marusic et al. 2010), and the estimations are conducted at $y_p^*< y_m^*$. In analogy to the incompressible flow (Townsend 1976), the centre of the logarithmic region is evaluated to be $3.9{Re_\tau ^{*1/2}}$, which gives $y_m^*=107$.

5.2. Estimation for streamwise velocity

The two components of $H_{L,uu}$, i.e. $\gamma _{uu}^2$ and $A_{pm,uu}$, are calculated first. They describe the spectral coherence and amplitude ratio of $\hat {u}^ {{\prime \prime }}$ at different heights. Afterwards, $H_{L,uu}$ is used to estimate the instantaneous field. In general, the SLSE results for $\hat {u}^ {{\prime \prime }}$ from DNS and linear models closely resemble those from the incompressible case (Madhusudanan et al. 2019; Gupta et al. 2021), so the discussion is presented briefly, highlighting crucial points.

Panels (a $-$c) in figures 7 and 8 are the results from the DNS data. The three estimation locations are $y_p^*=70$, 40 and 10, ranging from the logarithmic layer down to the buffer layer. Limited by the computational domain, the largest length scales displayed are $\lambda _x\approx 8h$ and $\lambda _z\approx 4h$. As in incompressible flows (Baars et al. 2016; Madhusudanan et al. 2019), large-scale structures have stronger coherence than the small scale, especially in the streamwise direction. This is in line with figure 6 where the streamwise-elongated streaky motions dominate the instantaneous velocity field. Meanwhile, the coherence is weakened with the decrease of $y_p^*$ (figures 7a–7c), simply because $y_m^*$ and $y_p^*$ are farther apart. At $y^*=10$, the region with $\gamma _{uu}^2>0.5$ only appears at $\lambda _x>4h$, or $\lambda _x^*>3600$, which is much larger than the dominant length scale of the $u^ {{\prime \prime }}$ streaks ($\lambda _x^*\approx {1000}$, Yao & Hussain 2020). However, though close to the wall, the motion of the largest scale is still coherent with that at $y_m^*=107$, reminiscent of the AEM. From figures 8(b) and 8(c), the structures with $\lambda _z<0.5h$ have relatively large $A_{pm,uu}$, which means narrow (in terms of $\lambda _z$) streaky motions are continuously enhanced approaching the wall. This is consistent with figure 6 where the fluctuation amplitude rises as $y^*$ decreases.

Figure 7. Contours of 2-D LCS for $u^ {{\prime \prime }}$ ($y_m^*=107$) at (a,d,g) $y_p^*=70$, (b,e,h) $y_p^*=40$ and (c,f,i) $y_p^*=10$. Panels (a $-$c) are from the DNS data, (d $-$f) are from the LNS model and (g $-$i) are from the eLNS model. The dotted lines denote $\lambda _x=\lambda _z$.

Figure 8. Contours of the amplitude ratio for $u^ {{\prime \prime }}$ ($A_{pm,uu}$, $y_m^*=107$) at (a,d,g) $y_p^*=70$, (b,e,h) $y_p^*=40$ and (c,f,i) $y_p^*=10$. Panels (a $-$c) are from the DNS data, (d $-$f) are from the LNS model and (g $-$i) are from the eLNS model. The dotted lines denote $\lambda _x=\lambda _z$.

We next focus on the linear model results. In general, $\gamma _{uu}^2$ and $A_{pm,uu}$ from the eLNS model are closer to DNS than the LNS model. By considering the wall-normal variation of the forcing amplitude (3.17), $\gamma _{uu}^2$ and $A_{pm,uu}$ near the wall can also be well predicted by eLNS, as shown in figures 7(i) and 8(i), though $A_{pm,uu}$ at small $\lambda _z$ ($\lesssim 0.5h$) tends to be under-predicted. In comparison, the LNS results have increasingly larger deviations towards the wall, and $\gamma _{uu}^2$ and $A_{pm,uu}$ are both underestimated. In particular, $A_{pm,uu}$ is smaller than unity, i.e. $\langle |\hat {u}_p^ {{\prime \prime }}|^2\rangle < \langle |\hat {u}_m^ {{\prime \prime }}|^2\rangle$ in most regions in figures 8(e) and 8(f), which is opposite to the DNS trend. This indicates that the structures modelled by LNS are localized in the wall-normal direction. The introduction of $\mu _t$ enhances the structural coherence between different heights, and can reasonably reflect the wall-normal coherence and amplitude variation of the velocity fluctuations. Notably, in the region $\lambda _x<\lambda _z$, $\gamma _{uu}^2$ in the eLNS model has some deviations from the DNS data, not present in incompressible cases. This difference is related to compressibility effects and will be further discussed in § 6.

After obtaining $H_L$, the large-scale coherent portion of the instantaneous field $u_p^{ {{\prime \prime }}}$ at $y_p^*$ can be estimated using (2.4a) and the DNS snapshots $u^{ {{\prime \prime }}}$ at $y_m^*$, reflecting the superposition effects in IOIM. As pointed out by Baars et al. (2016), the scales at very small $\gamma ^2$ can be erroneously estimated due to the possibly large $H_L$, then one may obtain intense small-scale motions in some regions of little coherence with the fields measured. To avoid this, a threshold $\gamma _{th}^2=0.05$ is introduced, and $H_L$ with $\gamma ^2<\gamma _{th}^2$ is forced to zero (Madhusudanan et al. 2019). In this way, the small-scale motions of little coherence are artificially ruled out. The resulting snapshots of $u_p^{ {{\prime \prime }}+}$ at $y_p^*=70$ and $10$, estimated from different models, are depicted in figure 9. The corresponding 2-D velocity spectra of $u_p^{ {{\prime \prime }}+}$,

(5.1)\begin{equation} {\varPhi}_{u_pu_p}^+(y_p) = \langle \hat{u}_p^{{{\prime\prime}}+}(y_p) \hat{u}_p^{{{\prime\prime}}+{{\dagger}}}(y_p) \rangle = |H_{L,uu}(y_p,y_m)|^2 \langle \hat{u}^{{{\prime\prime}}+}(y_m) \hat{u}^{{{\prime\prime}}+{{\dagger}}}(y_m) \rangle_{DNS}, \end{equation}

are also shown for more quantitative evaluation. Compared with figures 6(b) and 6(c), small-scale structures ($\lambda _x\lesssim h$, $\lambda _z\lesssim 0.3h$) are nearly missing in figures 9(a) and 9(d) due to the weak coherence at those scales, hence not considered by the superposition effects. Also, the dominant structures in figures 9(a) and 9(d) have similar $\lambda _x$ and $\lambda _z$. These footprints posed by the motions at $y_m^*$ on $y_p^*$ are crucial in the IOIM (Baars et al. 2016). For the linear models, the eLNS results in figures 9(c) and 9(f) bear a strong resemblance to the DNS counterparts, which is expected from the relative agreement of $\gamma _{uu}^2$ and $A_{pm,uu}$ in figures 7 and 8. Nevertheless, the structures of $\lambda _z\lesssim 0.5h$ are somewhat underestimated, which may be improved by further modelling the scale dependence of $\boldsymbol {B}$ (Gupta et al. 2021). In comparison, the LNS model estimates weak $u^{ {{\prime \prime }}+}$, especially down to the buffer layer. The velocity field and 2-D spectrum in figure 9(e) are nearly uniform under the contour levels, indicative of missing footprints from the wall-attached motions in the logarithmic region.

Figure 9. Estimated instantaneous $u_p^{ {{\prime \prime }}+}$ and the pre-multiplied 2-D spectra at (a–c) $y_p^*=70$ and (d–f) $y_p^*=10$ using the snapshot at $y_m^*=107$ and the kernel function from the (a,d) DNS data and models of (b,e) LNS and (c,f) eLNS.

In summary, the SLSE results for streamwise velocity are quite similar to the incompressible case, suggesting weak compressibility effects on the linearized momentum equation under current flow conditions. By introducing $\mu _t$ and $\kappa _t$, the spectral linear coherence and amplitude ratio of $\hat {u}^ {{\prime \prime }}$ between different heights can be well captured by the linear models (eLNS), though the structures of small $\lambda _z$ can be somewhat underestimated. This agreement gives us confidence to perform a parameter study for the velocity fluctuations in § 6. Instead of collecting large amounts of time-resolved measurements or simulations, the kernel function of velocity can be obtained using the eLNS model if one only has a mean flow. Then SLSE can be performed based on one snapshot of the velocity field.

5.3. Estimation for temperature

The SLSE results for temperature fluctuations are investigated in a similar style to § 5.2. In contrast to the velocity counterparts, the features of $\gamma _{TT}^2$ and $H_{L,TT}$ were rarely reported in previous studies, so the implications from DNS are discussed first.

The $\gamma _{TT}^2$ and $A_{pm,TT}$ from DNS are plotted in panels (a $-$c) in figures 10 and 11 at $y_p^*=70$, 40 and 10. The distribution of $\gamma _{TT}^2$ is very similar to $\gamma _{uu}^2$. In particular, figure 10(c) is nearly indistinguishable from figure 7(c) at $y_p^*=10$. An important implication is that in the buffer layer, in addition to the nearly perfect correlation between $u^ {{\prime \prime }}$ and $T^ {{\prime \prime }}$ (see $R_{u^ {{\prime \prime }} T^ {{\prime \prime }}}$ and figure 6), the spectral coherence with the measurement ($\gamma _{uu}^2$ and $\gamma _{TT}^2$) are basically identical for $\hat {u}^ {{\prime \prime }}$ and $\hat {T}^ {{\prime \prime }}$, which benefits the compressible turbulence modelling. Moreover, $\gamma _{TT}^2$ is no lower than $\gamma _{uu}^2$ at the three $y_p^*$ displayed, and is even higher in some regions with $\lambda _x<\lambda _z$. Therefore, the wall-normal spectral coherence of the temperature fluctuation is at least not weaker than the streamwise velocity. Intense inner–outer interaction is anticipated, and the view is supported that $T^ {{\prime \prime }}$ is also a wall-attached quantity (Cheng & Fu 2023a). The peaks of $\gamma _{uu}^2$ and $\gamma _{TT}^2$ are all located at $\lambda _z\approx h$ at different $y_p^*$, different from the impression from figure 6 that $T^ {{\prime \prime }}$ tends to be more isotropic as $y^*$ rises. High values of $A_{pm,TT}$ appear at $\lambda _z<0.5h$, the same as that for $A_{pm,uu}$. The maximum $A_{pm,TT}$ exceeds six at $y_p^*=10$ (figure 11c), which means that $\langle |\hat {T}_p^ {{\prime \prime }}|^2\rangle$ with small $\lambda _z$ increases rapidly towards the wall with respect to $\langle |\hat {T}_m^ {{\prime \prime }}|^2\rangle$, due to local large temperature gradients.

Figure 10. Contours of 2-D LCS for $T^ {{\prime \prime }}$ ($y_m^*=107$) at (a,d,g) $y_p^*=70$, (b,e,h) $y_p^*=40$ and (c,f,i) $y_p^*=10$. Panels (a $-$c) are from the DNS data, (d $-$f) are from the LNS model and (g $-$i) are from the eLNS model.

Figure 11. Contours of the amplitude ratio for $T^ {{\prime \prime }}$ ($A_{pm,TT}$, $y_m^*=107$) at (a,d,g) $y_p^*=70$, (b,e,h) $y_p^*=40$ and (c,f,i) $y_p^*=10$. Panels (a $-$c) are from the DNS data, (d $-$f) are from the LNS model and (g $-$i) are from the eLNS model.

The linear model results for $T^ {{\prime \prime }}$ are discussed below. In general, the LNS model underestimates $\gamma _{TT}^2$ and $A_{pm,TT}$ at all $y_p^*$ selected, the same as that for $u^ {{\prime \prime }}$. The $A_{pm,TT}$ in figures 11(e) and 11(f) are less than unity in most regions, so the $\hat {T}^ {{\prime \prime }}$ predicted is more localized with insufficient wall-normal connection. Nevertheless, the feature of $\lambda _x$ preference over $\lambda _z$ for $\gamma _{TT}^2$ is captured in the LNS model, and $\gamma _{TT}^2$ peaks at $\lambda _z\approx h$. The behaviour of the eLNS model, however, is quite different from those for $u^ {{\prime \prime }}$. The $\gamma _{TT}^2$ in figures 10(g)–10(i) differs considerably from the DNS data, nearly isotropic in terms of $\lambda _x$ and $\lambda _z$. Also, $A_{pm,TT}$ at all $y_p^*$ is underestimated. Thereby, the eLNS model fails to give agreeable results with DNS for estimating temperature fluctuations. Revisiting the assumptions in § 3.2, (3.8a) treats $T^ {{\prime \prime }}$ in an isotropic way, which may be responsible for the isotropic behaviour of $\gamma _{TT}^2$ in the eLNS model. Nevertheless, our numerical tests show that using (3.8b) instead does not improve obviously the results, though some degree of anisotropy is introduced. The discussion in § 4 suggests that acoustic and entropy components are significant in the eLNS model, but no supporting signs are observed from the displayed DNS data. Therefore, we check the eLNS results below for possible explanations of the differences.

5.4. Cospectrum decomposition

To figure out the difference between the eLNS and DNS results, $\boldsymbol {\varPhi }$ is scrutinized. As discussed in § 4, multiple branches of modes are present. One way to identify different mechanisms is to conduct a mode decomposition for $\boldsymbol {\varPhi }$, known as the Karhunen–Loéve or proper orthogonal decomposition (POD). The eigenvalue problem is constructed as $\boldsymbol {\varPhi } \check {\boldsymbol q}_j^ {{\prime \prime }} = \theta _j \check {\boldsymbol q}_j^ {{\prime \prime }}$, where $\theta _j$ is the eigenvalue, and the matrix–vector multiplication is defined in terms of the energy norm. Note that the symbol $\check {\boldsymbol q}_j^ {{\prime \prime }}$ is reused to represent the eigenfunctions for clarity, though the physical meaning is not identical to that in § 4. Since $\boldsymbol {\varPhi }$ is Hermite, $\theta _j$ is real. Also, $\check {\boldsymbol q}_j^ {{\prime \prime }}$ satisfies the orthogonal relation $(\check {\boldsymbol q}_i^ {{\prime \prime }},\check {\boldsymbol q}_j^ {{\prime \prime }})_E = E_b\delta _{ij}$ ($E_b$ is the energy dimension). Different from (4.1), $\theta _j$ measures the energy of the $j$th POD mode ($\check {\boldsymbol q}_j^ {{\prime \prime }}$) and $\sum _j \theta _j=V$. If the pairs $(\theta _j, \check {\boldsymbol q}_j^ {{\prime \prime }})$ are sorted in the descending order of $\theta _j$, then $\check {\boldsymbol q}_1^ {{\prime \prime }}$ represents the leading energy-containing structure of the response.

The response with $\lambda _x=8h$, $\lambda _z=2h$ is considered first using the eLNS model, close to the maximum $\gamma _{uu}^2$ and $\gamma _{TT}^2$. Note that $\boldsymbol {\varPhi }$ discussed below is the solution of (3.15), not (5.1). After POD, at least two types of modes are identified, resulting from the different branches in § 4. The most energetic mode $\check {\boldsymbol q}_1^ {{\prime \prime }}$ is shown in figure 12(a), where $\check {u}^ {{\prime \prime }}$ is the largest component. Meanwhile, $\check {\rho }^ {\prime }$ and $\check {T}^ {{\prime \prime }}$ (not shown) can be well approximated by (3.16a,b) and $\check {u}^ {{\prime \prime }}$, analogous to figure 2. This demonstrates that $\check {\rho }^ {\prime }$ and $\check {T}^ {{\prime \prime }}$ of this mode are dominated by advection, and $\check {\boldsymbol q}_1^ {{\prime \prime }}$ is in line with the trends in DNS. Therefore, this mode is termed an advection mode, which can be considered as a subset of the vortical modes. The dilatation of the advection mode is tiny, and $\check {p}^ {\prime }/\bar {p}$ is also small, thus analogous to the incompressible counterpart. Meanwhile, the advection mode takes the form of streamwise-elongated velocity ($\check {u}^ {{\prime \prime }}$), density and temperature streaks forced by streamwise vortices ($\check {v}^ {{\prime \prime }}$, $\check {w}^ {{\prime \prime }}$), as observed in figure 6 and reported in previous experimental, DNS and linear-model results (Coleman et al. 1995; Williams et al. 2018; Chen et al. 2023). A 3-D view of this mode will be shown later. Modes $\check {\boldsymbol q}_2^ {{\prime \prime }}$, $\check {\boldsymbol q}_3^ {{\prime \prime }}$ and $\check {\boldsymbol q}_4^ {{\prime \prime }}$ have similar features to $\check {\boldsymbol q}_1^ {{\prime \prime }}$, but $\check {\boldsymbol q}_5^ {{\prime \prime }}$ is different. As shown in figure 12(b) for $\check {\boldsymbol q}_5^ {{\prime \prime }}$, $\check {p}^ {\prime }$ and $\check {\rho }^ {\prime }$ are much larger than $\check {u}^ {{\prime \prime }}$, and $\check {\rho }^ {\prime }$ is in high amplitude at the wall, so this mode experiences strong dilatation. This mode is of an acoustic nature, as discussed in § 4.1. To prove this, the analytical eigenfunction of the acoustic mode in the free stream is used (see Appendix B), leading to the following relations:

(5.2a,b)\begin{equation} \frac{\check{\rho}^{\prime}}{\bar{\rho}} = \frac{\check{p}^{\prime}}{\gamma_0\bar{p}}, \quad \frac{\check{u}^{{\prime\prime}}}{a} = \frac{k_x}{\sqrt{k_x^2+k_y^2+k_z^2}} \frac{\check{p}^{\prime}}{\gamma_0\bar{p}}, \end{equation}

where $a$ is the speed of sound. Equation (5.2a,b) is examined in figure 12(b) with $k_y=0$. Even though (5.2a,b) is derived based on a uniform mean flow, remarkable agreement is observed for both $\check {\rho }^ {\prime }$ and $\check {u}^ {{\prime \prime }}$, confirming that $\check {\boldsymbol q}_5^ {{\prime \prime }}$ is of an acoustic nature except in the near-wall region. The large $\check {\rho }^ {\prime }$ at the wall is due to the prescribed boundary condition $\check {T}^ {{\prime \prime }}_w=0$. Furthermore, a series of acoustic modes exist with different $k_y$, as shown in figure 13. These modes with $k_y\ne 0$ strongly oscillate in the interior region with (5.2a,b) also satisfied, and the maximum amplitudes of $\check {\rho }^ {\prime }$ and $\check {p}^ {\prime }$ remain nearly constants. The $\theta _j$ of these modes are labelled in figures 12 and 13, non-dimensionalized by $h/U_b$. The energy ratio of the most energetic acoustic and advection modes is 8.2 % ($\theta _5/\theta _1$), so these acoustic modes have non-negligible contributions to $\boldsymbol {\varPhi }$. In comparison, the same ratio in the LNS case is only 0.06 %. Therefore, the conclusion in § 4 is supported since the acoustic components are negligible for the response in the LNS case, but non-negligible in the eLNS case. Moreover, no counterparts to the entropy modes are found in at least the top thirty most energetic modes in both the LNS and eLNS cases, suggesting their relatively small contribution than the acoustic components at this length scale.

Figure 12. Profiles of different components of the cospectrum POD modes in the eLNS model ($\lambda _x=8h$, $\lambda _z=2h$): (a) the most energetic mode ($\theta _1$) and (b) the fifth most energetic mode ($\theta _5$). Equations (3.16a,b) and (5.2a,b) are used to approximate other components based on $\check {u}^ {{\prime \prime }}$ (for (a)) and $\check {p}^ {\prime }$ (for (b)).

Figure 13. Profiles of different components of the cospectrum POD modes in the eLNS model ($\lambda _x=8h$, $\lambda _z=2h$): the (a) seventh, (b) eleventh and (c) fifteenth energetic modes. Equation (5.2a,b) is used to approximate other components based on $\check {p}^ {\prime }$.

Considering § 4, and observing the great similarities between $u^ {{\prime \prime }}$ and $T^ {{\prime \prime }}$ in §§ 5.2 and 5.3, it is conjectured that the significant role of acoustic and entropy components to the fluctuation growth in the eLNS model is an intrinsic defect introduced by the turbulence modelling. To verify this, $\boldsymbol {\varPhi }$ is decomposed to exclude the contributions of acoustic and entropy modes, denoted as $\boldsymbol {\varPhi }_{ae}$, and the remaining part is $\boldsymbol {\varPhi }_{nae}$. The eigenvalue problem above does not allow a strict decomposition as $\boldsymbol {\varPhi }=\boldsymbol {\varPhi }_{ae}+\boldsymbol {\varPhi }_{nae}$, so the singular value decomposition is deployed instead, taking the form of

(5.4)\begin{equation} \boldsymbol{\varPhi} = \sum_j { \sigma_j \check{\boldsymbol \psi}_j \check{\boldsymbol \phi}_j^{H}} = \underbrace{\sum_{j\in {ae}} { \sigma_j \check{\boldsymbol \psi}_j \check{\boldsymbol \phi}_j^{H}}}_{\boldsymbol{\varPhi}_{ae}} + \underbrace{\sum_{j\in {nae}} { \sigma_j \check{\boldsymbol \psi}_j \check{\boldsymbol \phi}_j^{H}}}_{\boldsymbol{\varPhi}_{nae}}, \end{equation}

where $\sigma _j$ is the singular value, and $\check {\boldsymbol \psi }_j$ and $\check {\boldsymbol \phi }_j$ are the decomposition modes. As $\boldsymbol {\varPhi }$ is Hermite, $\check {\boldsymbol \psi }_j=\check {\boldsymbol \phi }_j$; $\sigma _j$ and $\check {\boldsymbol \psi }_j$ have analogous behaviours to $\theta _j$ and $\check {\boldsymbol q}_j^ {{\prime \prime }}$, hence are not elaborated for conciseness. The same criteria as in § 4 are deployed to identify the acoustic and entropy modes, based on their eigenfunctions (see figure 4). It is worth mentioning the difference of the decompositions in figure 5 and (5.4). The former is to exclude the acoustic (or entropy) eigenmodes of the linear operator before solving the response, while the latter is a posterior decomposition of the response with all eigenmodes included. The two strategies are not equivalent, and the latter is adopted here because the acoustic and entropy modes are not fully decoupled with the vortical modes (see § 4), especially in the near-wall region.

The cospectrum decomposition results using (5.4) are displayed in figure 14, still with $\lambda _x=8h$, $\lambda _z=2h$. They are normalized by the maximum of $\langle \hat {u}^ {{\prime \prime }} \hat {u}^{ {{\prime \prime }}{\dagger} }\rangle$. Interestingly, $\langle \hat {\rho }^ {\prime }\hat {\rho }^{ {\prime }{\dagger} }\rangle$, $\langle \hat {u}^ {{\prime \prime }} \hat {u}^{ {{\prime \prime }}{\dagger} }\rangle$ and $\langle \hat {u}^ {{\prime \prime }} \hat {T}^{ {{\prime \prime }}{\dagger} }\rangle$ have clear attributions, suggesting a rather ‘perfect’ decomposition. Here, ‘perfect’ means that we can clearly identify which physical process is responsible for the generation of different components. Specifically, nearly all $\langle \hat {\rho }^ {\prime }\hat {\rho }^{ {\prime }{\dagger} }\rangle$ are contributed by $\boldsymbol {\varPhi }_{ae}$, while the contribution of $\boldsymbol {\varPhi }_{nae}$ is tiny (figures 14e and 14i). From figure 12(a), the latter part can be approximated by $u^ {{\prime \prime }}$ and the mean flow gradients, representing the advection effect. In contrast, $\langle \hat {u}^ {{\prime \prime }} \hat {u}^{ {{\prime \prime }}{\dagger} }\rangle$ and $\langle \hat {u}^ {{\prime \prime }} \hat {T}^{ {{\prime \prime }}{\dagger} }\rangle$ are mostly contributed by $\boldsymbol {\varPhi }_{nae}$, nearly unaffected by the acoustic and entropy components. Thereby, the eLNS results for $u^ {{\prime \prime }}$ in § 5.2 basically agree with DNS, as in the incompressible case. Also, it is suggested that the $u^ {{\prime \prime }}$–$u^ {{\prime \prime }}$ and $u^ {{\prime \prime }}$–$T^ {{\prime \prime }}$ couplings are dominated by advection and other vortical motions, while the acoustic- and entropy-induced motions have little coherence for the couplings. For $\langle \hat {T}^ {{\prime \prime }} \hat {T}^{ {{\prime \prime }}{\dagger} }\rangle$, $\boldsymbol {\varPhi }_{nae}$ and $\boldsymbol {\varPhi }_{ae}$ have nearly equal contributions. The former is more concentrated towards the wall, represented by (3.16a,b), while the latter extends to the channel centre populating neutral acoustic and entropy waves. To further justify the decomposition, the normalized $\boldsymbol {\varPhi }$ from DNS is also computed. Since we are mainly concerned with the relative amplitudes of $\boldsymbol {\varPhi }_{nae}$ and $\boldsymbol {\varPhi }_{ae}$, only the $\boldsymbol {\varPhi }$ with $y=y^ {\prime }$ are extracted for comparison, as shown in figures 14(m)–14(p). For $\langle \hat {u}^ {{\prime \prime }} \hat {u}^{ {{\prime \prime }}{\dagger} }\rangle$ and $\langle \hat {u}^ {{\prime \prime }} \hat {T}^{ {{\prime \prime }}{\dagger} }\rangle$, $\boldsymbol {\varPhi }$ and $\boldsymbol {\varPhi }_{nae}$ are close to each other and they basically follow $\boldsymbol {\varPhi }_{DNS}$, though differences in shapes are observable. The $\boldsymbol {\varPhi }$ in the eLNS model varies more mildly than the DNS data, possibly due to using a smooth $\mu _t$. For $\langle \hat {\rho }^ {\prime }\hat {\rho }^{ {\prime }{\dagger} }\rangle$ and $\langle \hat {T}^ {{\prime \prime }} \hat {T}^{ {{\prime \prime }}{\dagger} }\rangle$, however, $\boldsymbol {\varPhi }$ is several times higher than $\boldsymbol {\varPhi }_{DNS}$, exhibiting large values around the centreline. After the decomposition, $\boldsymbol {\varPhi }_{nae}$ is much closer to $\boldsymbol {\varPhi }_{DNS}$, also tending to diminish towards the channel centre. Thereby, we can conclude that the acoustic and entropy components $\boldsymbol {\varPhi }_{ae}$ in the eLNS model barely exist in the DNS data. In short, in the eLNS model, $\langle \hat {u}^ {{\prime \prime }} \hat {u}^{ {{\prime \prime }}{\dagger} }\rangle$ and $\langle \hat {u}^ {{\prime \prime }} \hat {T}^{ {{\prime \prime }}{\dagger} }\rangle$ remain little influenced by the acoustic and entropy components, while $\langle \hat {\rho }^ {\prime }\hat {\rho }^{ {\prime }{\dagger} }\rangle$ and $\langle \hat {T}^ {{\prime \prime }} \hat {T}^{ {{\prime \prime }}{\dagger} }\rangle$ are largely affected, which can be responsible for the deviation of $H_{L,TT}$ in § 5.3. For improvement, we use $\boldsymbol {\varPhi }_{nae}$ in the following to estimate temperature and also the velocity–temperature coupling. The results are denoted as eLNS-CD (‘CD’ for cospectrum decomposition).

Figure 14. Cospectra of different components in the eLNS model ($\lambda _x=8h$, $\lambda _z=2h$): (a,e,i,m) $\langle \hat {\rho }^ {\prime } \hat {\rho }^{ {\prime }{\dagger} } \rangle$, (b,f,j,n) $\langle \hat {u}^ {{\prime \prime }} \hat {u}^{ {{\prime \prime }}{\dagger} } \rangle$, (c,g,k,o) $\langle \hat {T}^ {{\prime \prime }} \hat {T}^{ {{\prime \prime }}{\dagger} } \rangle$ and (d,h,l,p) $|\langle \hat {u}^ {{\prime \prime }} \hat {T}^{ {{\prime \prime }}{\dagger} } \rangle |$. They are non-dimensionalized by $\rho _b$, $U_b$ and $T_w$ accordingly. Panels (a $-$d) are the original results containing all modes, (e $-$h) are from $\boldsymbol {\varPhi }_{nae}$ and (i $-$l) are from $\boldsymbol {\varPhi }_{ae}$. The coordinates $y$ and $y^ {\prime }$ are for the two variables of cospectra, respectively (see (2.7)). Panels (m $-$p) are the comparisons with $\boldsymbol {\varPhi }_{DNS}$ at $y^ {\prime }=y$.

Before presenting the eLNS-CD results, we need to examine whether the dominant advection mode in figure 12 is a universal feature at different length scales. The scales $(\lambda _x,\lambda _z)=(8h,0.3h)$, $(0.3h,2h)$ and $(0.3h,0.3h)$ are considered in addition to $(8h,2h)$ in figure 12, covering the large and small streamwise/spanwise scales studied in § 5.3. The velocity and temperature components of the most energetic modes at these scales are plotted in figure 15. Their 3-D structures are depicted in figure 16, where $u_{(k_x,k_z)}^ {{\prime \prime }}=\check {u}_{(k_x,k_z)}^ {{\prime \prime }}(y)\exp [ {\mathrm {i}} (k_x x+k_z z)]$ (from (2.3)). The energetic large-scale mode $(8h,2h)$ wanders across different layers, crucial for the inner–outer interaction as discussed above. This type of mode is also characterized by geometrically self-similar hierarchies (McKeon 2019). The mode $(8h,0.3h)$ is also an advection mode as in figure 12(a), exhibiting streamwise-elongated streaky motions very close to the wall. The peaks of $\check {u}^ {{\prime \prime }}$ and $\check {T}^ {{\prime \prime }}$ are at $y^*=14$ and 11, respectively, within the buffer layer. Also, $\check {\rho }^ {\prime }$ and $\check {T}^ {{\prime \prime }}$ result primarily from the advection process, well approximated by (3.16a,b). In comparison, the mode $(0.3h,2h)$ is spanwise elongated (since $\lambda _z>\lambda _x$) with large amplitudes of $\check {v}^ {{\prime \prime }}$ and $\check {w}^ {{\prime \prime }}$. This mode is not dominated by streamwise advection, but also results from vortical modes since the thermodynamic components are small. For the mode $(0.3h,0.3h)$ in figures 15(c) and 16(d), $\check {v}^ {{\prime \prime }}$ is the largest component and $\check {u}^ {{\prime \prime }}\sim \check {w}^ {{\prime \prime }}$ as $\lambda _x\sim \lambda _z$. This mode is mostly amplified around $y\approx 0.5h$, and its oblique structure can be associated with the spanwise meandering of large-scale motions (Hutchins & Marusic 2007; Abe et al. 2018). Meanwhile, its $\check {T}^ {{\prime \prime }}$ component is in low amplitude, much lower than the approximation by (3.16a,b). In short, with the decrease of $\lambda _x$, the most energetic structure changes gradually from the motions dominated by streamwise advection to spanwise elongation. Consequently, the cospectrum decomposition at small $\lambda _x$ may not be as ‘perfect’ as in figure 14. From §§ 5.2 and 5.3, the structures with smaller $\lambda _x$ tend to have weaker coherence for the velocity and temperature fluctuations between different heights, so they can be of secondary importance for SLSE.

Figure 15. Different components of the most energetic POD modes ($\theta _1$) for the eLNS cospectrum at different length scales ($\lambda _x, \lambda _z$): (a) $(8h,0.3h)$, (b) $(0.3h,2h)$ and (c) $(0.3h,0.3h)$.

Figure 16. Isosurface $u_{(k_x,k_z)}^ {{\prime \prime }}/ max|u_{(k_x,k_z)}^ {{\prime \prime }}|=0.6$ of the most energetic mode of the eLNS cospectrum for ($\lambda _x$, $\lambda _z$) of (a) $(8h,2h)$, (b) $(8h,0.3h)$, (c) $(0.3h,2h)$ and (d) $(0.3h,0.3h)$, flooded by the mean velocity $\tilde {u}/U_b$.

The $\gamma _{TT}^2$ and $A_{pm,TT}$ from the eLNS-CD model are displayed in figures 17 and 18 at the same $y_m^*$ and $y_p^*$ as in § 5.3. The DNS results in figures 10 and 11 are re-plotted for easier comparison. The cospectrum decomposition requires categorization of the POD modes, so the eLNS-CD results are not as smooth as those by eLNS, as discussed for figures 5(e) and 15. As can be seen, after removing the acoustic and entropy components, $\gamma _{TT}^2$ and $A_{pm,TT}$ are in good agreement with DNS. To be specific, $\gamma _{TT}^2$ from eLNS-CD is not as isotropic as that in figure 10 (eLNS results), and the scale selection in terms of $\lambda _z$ is captured, so $\gamma _{TT}^2$ does not monotonically increase with $\lambda _z$. Meanwhile, the high $A_{pm,TT}$ at $\lambda _z<0.5h$ is correctly reflected, showing a prominent improvement over figures 11(g)–11(i).

Figure 17. Contours of 2-D LCS for $T^ {{\prime \prime }}$ ($y_m^*=107$) at (a,d) $y_p^*=70$, (b,e) $y_p^*=40$ and (c,f) $y_p^*=10$. Panels (a $-$c) are from the DNS data and (d $-$f) are from the eLNS model after the cospectrum decomposition (eLNS-CD). The dotted lines denote $\lambda _x=\lambda _z$.

Figure 18. Contours of the amplitude ratio for $T^ {{\prime \prime }}$ ($A_{pm,TT}$, $y_m^*=107$) at (a,d) $y_p^*=70$, (b,e) $y_p^*=40$ and (c,f) $y_p^*=10$. Panels (a $-$c) are from the DNS data and (d $-$f) are from the eLNS-CD model. The dotted lines denote $\lambda _x=\lambda _z$.

Finally, the instantaneous $T^ {{\prime \prime }}$ at $y_m^*$ is utilized to estimate the coherent temperature structures ($T_p^ {{\prime \prime }}$) at different $y_p^*$. The estimated 2-D spectra ${\varPhi }_{T_pT_p}^+(y_p)$ as in (5.1) are also computed. First, the estimations using $H_{L,TT}$ from DNS are depicted in figures 19(a) and 19(d). Compared with figure 9(a), relatively more energy resides at $\lambda _x\sim h$, contributing to the more isotropic pattern of temperature than the streamwise velocity counterpart. The estimated $T_p^ {{\prime \prime }}$ at $y_p^*=10$ is dominated by streamwise streaks and resembles a lot the velocity counterpart in figure 9(d) except for smaller amplitudes. Similar to figure 9, the LNS model severely underestimates the amplitude of $T_p^ {{\prime \prime }}$ at nearly all scales, especially close to the wall where the two plots in figure 19(e) are nearly uniform. In comparison, the eLNS-CD results exhibit significant improvements over the LNS and eLNS models, as expected from figures 17 and 18. The estimated large-scale structures ($\lambda _x\gtrsim 3h$) agree well with the DNS data, but the smaller-scale motions in figure 19(c) are obviously under-predicted due to the underestimated $\gamma _{TT}^2$ in figure 17(d) at small $\lambda _x$.

Figure 19. Estimated instantaneous $T_p^{ {{\prime \prime }}+}$ and the pre-multiplied 2-D spectra at (a,b,c) $y_p^*=70$ and (d,e,f) $y_p^*=10$ using the snapshot at $y_m^*=107$ and the kernel function from the (a,d) DNS data and models of (b,e) LNS and (c,f) eLNS-CD.

Therefore, the eLNS model can reflect reasonably the wall-normal coherence and amplitude variation of the temperature fluctuation, but it can be deteriorated by the overdamped vortical modes and thus relatively intensified acoustic and entropy components due to the changed mathematical properties of the linear operator. Through a procedure of cospectrum decomposition (eLNS-CD), the interference from the acoustic and entropy components can be removed to obtain SLSE results that basically agree with DNS data. Meanwhile, an important implication is that the coherence of temperature is dominated by the vortical components, including advection and other kinematically dominated motions, strongly analogous to the incompressible case. The compressibility effects, exerting influence through the acoustic and entropy components, are of secondary importance.

5.5. Estimation for velocity–temperature coupling

Cheng & Fu (2023a) show that SLSE can be designed to estimate the velocity–temperature coupling (§ 2.3). Here we provide the results from the linear models to examine whether $\gamma _{uT}^2$ and $H_{L,uT}$ can be reasonably predicted. First of all, $\gamma _{uT}^2$ and $A_{pm,uT}= ({\langle |\hat {T}_p^ {{\prime \prime }}|^2\rangle }/{\langle |\hat {u}_m^ {{\prime \prime }}|^2\rangle })^{1/2}$ from DNS are displayed in figures 20(a)–20(c). Two types of $y_m^*$ and $y_p^*$ combinations are considered, namely $y_m^*=y_p^*$ and $y_m^*\ne y_p^*$. In the former case, the measured $u^ {{\prime \prime }}$ is used for predicting the local $T^ {{\prime \prime }}$ at the same height, while in the latter case, $T^ {{\prime \prime }}$ at another height is predicted. In figure 20(a), $y_m^*=y_p^*=14$ are both in the buffer layer, and $u^ {{\prime \prime }}$ and $T^ {{\prime \prime }}$ have close morphological connections (see § 5.1). Therefore, $\gamma _{uT}^2$ is larger than 0.6 in most length scales displayed, and the streamwise-elongated streaky motions (larger $\lambda _x$ and smaller $\lambda _z$) are more coherent in terms of the coupling. Meanwhile, $A_{pm,uT}\approx 0.35$ (non-dimensionalized by $U_b$, $T_w$) is nearly uniform at different scales (not shown), close to the mean flow coefficient between $u_{rms}^ {{\prime \prime }}$ and $T_{rms}^ {{\prime \prime }}$ in (3.16a,b). Figures 20(b) and 20(c) are for the case $y_m^*\ne y_p^*$, where $y_m^*=107$ is at the centre of the logarithmic region and $y_p^*$ is in the buffer layer. Here, $\gamma _{uT}^2$ and $A_{pm,uT}$ have similar distributions to those for $u^ {{\prime \prime }}$–$u^ {{\prime \prime }}$ and $T^ {{\prime \prime }}$–$T^ {{\prime \prime }}$, which is reasonable considering the strong analogies between $\gamma _{uu}^2$ and $\gamma _{TT}^2$, and $A_{pm,uu}$ and $A_{pm,TT}$.

Figure 20. Contours of 2-D LCS and amplitude ratio for the $u^ {{\prime \prime }}$–$T^ {{\prime \prime }}$ coupling from the (a $-$c) DNS data, and the models of (d $-$f) LNS, (g $-$i) eLNS and (j $-$l) eLNS-CD. Panels (a,d,g,j) are at $y_m^*=y_p^*=14$, and other panels are at $y_m^*=107, y_p^*=14$. Panels (c,f,i,l) are for the amplitude ratio $({\langle |\hat {T}^ {{\prime \prime }}(y_p)|^2\rangle }/{\langle |\hat {u}^ {{\prime \prime }}(y_m)|^2\rangle })^{1/2}$, and other panels are for $\gamma _{uT}^2$. The dotted lines denote $\lambda _x=\lambda _z$. The dash-dot lines in (f,i) are the contours in the saturating region.

From figure 20, $\gamma _{uT}^2$ and $A_{pm,uT}$ from LNS and eLNS considerably deviate from the DNS data. The $\gamma _{uT}^2$ at different $y_m^*$ and $y_p^*$ are underestimated in both models. In addition, $A_{pm,uT}$ in the LNS case is high for spanwise-elongated motions ($\lambda _z>\lambda _x$) but low for the streamwise-elongated ($\lambda _x>\lambda _z$), opposite to the DNS trend. Also, $A_{pm,uT}$ in the eLNS case exhibits an overall higher value at nearly all $\lambda _x$ and $\lambda _z$. From figure 14, $\langle \hat {u}^ {{\prime \prime }} \hat {T}^{ {{\prime \prime }}{\dagger} } \rangle$ by eLNS is little affected by the acoustic and entropy components, but $\gamma _{uT}^2$ and $A_{pm,uT}$ are influenced since they are related to $\langle \hat {T}^ {{\prime \prime }} \hat {T}^{ {{\prime \prime }}{\dagger} } \rangle$. After the cospectrum decomposition, the eLNS-CD results in figures 20(j)–20(l) are much closer to DNS. The scale dependence of $\gamma _{uT}^2$ and $A_{pm,uT}$ are correctly captured, for both $y_m^*=y_p^*$ and $y_m^*\ne y_p^*$ cases. Nevertheless, an obvious deviation is that $\gamma _{uT}^2$ for the spanwise-elongated motions is underestimated in figure 20(j). As discussed in figure 15(b), wall-normal and spanwise motions dominate at these scales, affecting the performance of the cospectrum decomposition. In summary, the eLNS model can be used to estimate the velocity–temperature coupling after the cospectrum decomposition (eLNS-CD). Analogous to the coherence of temperature, the coupling of the coherent velocity and temperature fluctuations is also dominated by advection and other vortical components, while the compressibility effects are of secondary importance. Thereby, the present results support and provide more proof for the conclusions of Cheng & Fu (2023a).

6. Parameter study

Observing the consistency between the DNS data and the linear models, a parameter study is performed to investigate the variations of $\gamma ^2$ and $A_{pm}$ with different ${\textit {Ma}}_b$ and $Re_\tau$. As the compressible DNS database is very limited, the ODE solver in § 2.2 is utilized to obtain the mean flows of varying ${\textit {Ma}}_b$ and $Re_\tau$. First, we need to confirm that the SLSE results predicted by the ODE-based mean flow are as reliable as those using the DNS mean flow. The benchmark case is calculated, and the comparisons of $\gamma _{uu}^2$ and $\gamma _{TT}^2$ are shown in figure 21 using the two mean flows. The $y_m^*$ and $y_p^*$ are the same as those in § 5, with $y_p^*=40$ in the logarithmic layer. Note that the result comparisons for other $y_p^*$ are similar and thus not shown for conciseness. For an overall comparison, the ranges of $\lambda _x$ and $\lambda _z$ are much larger than those in § 5, without the limitation of the DNS computational domain. Only the eLNS model is selected, and the LNS model is not considered due to the relatively large deviation. Basically, $\gamma _{uu}^2$ and $\gamma _{TT}^2$ using the DNS and ODE mean flows agree with each other within the length scales displayed, demonstrating the reliability of the ODE solver in predicting SLSE results.

Figure 21. Contours of (a) $\gamma _{uu}^2$ and (b) $\gamma _{TT}^2$ using the mean flows from DNS and the ODE-based solver, respectively. Here, $y_m^*$ is fixed at $107$, and $y_p^*=40$. The levels displayed are 0.1, 0.3, 0.5, 0.7 and 0.9, as labelled.

Semi-local units have been extensively proven to be capable of collapsing various turbulent quantities to their incompressible counterparts, such as the mean velocity profile (Griffin et al. 2021), Reynolds stress (Huang et al. 1995), the spanwise spacing of near-wall streaks (Morinishi et al. 2004), kinetic energy budget (Duan et al. 2010) and streamwise inclination angles (Bai, Cheng & Fu 2023). One interesting point is to examine whether semi-local units collapse the 2-D LCS and amplitude ratio. This is conducive to the modelling of inner–outer interaction for compressible flows. In addition to the benchmark case, three more cases are considered, with varying ${\textit {Ma}}_b$ but the same $Re_{\tau,c}^*$. Their computational parameters are listed in table 1. The ${\textit {Ma}}_b=0.3$ case is close to the incompressible condition with nearly uniform $\bar {\rho }$ and $\tilde {T}$, and the fourth case has a high ${\textit {Ma}}_b$ of 4.

Table 1. Parameters for channel flow cases with different ${\textit {Ma}}_b$ but the same $Re_{\tau,c}^*$.

With ${\textit {Ma}}_b$ increased, $Re_\tau$ and $Re_b$ grow rapidly due to the stronger variation of the thermodynamic properties. These two Reynolds numbers are positively correlated with the required grid points for DNS, so the computational cost dramatically increases with the rise of ${\textit {Ma}}_b$, which is largely responsible for the limited DNS database for compressible channel flows.

For different cases, $y_m^*$ and $y_p^*$ are kept the same, though the corresponding $y_m^+$ and $y_p^+$ increase by over three times with ${\textit {Ma}}_b$ lifted. Meanwhile, the streamwise and spanwise wavelengths of fluctuations are all expressed in the outer scale, i.e. $\lambda _x/h$ and $\lambda _z/h$. The contours of $\gamma _{uu}^2$ and $A_{pm,uu}$ are depicted in figure 22 at $y_p^*=70$ and 10. The $\gamma _{uu}^2$ at $y_p^*=70$ varies slightly with ${\textit {Ma}}_b$ increased from 0.3 to 4.0, which means that the compressibility effects on the spectral coherence of streamwise velocity in the logarithmic region can be scaled using semi-local units, supporting Morkovin's hypothesis. When $y_p^*$ is down to 10 in the buffer layer, the differences of $\gamma _{uu}^2$ among cases are more obvious, especially at large $\lambda _z$. Nevertheless, there is still good collapse for streamwise-elongated streaky motions with $\lambda _x\gtrsim 5h$ and $\lambda _z\lesssim 0.5h$. In figures 22(b) and 22(d), the region $A_{pm,uu}>1$, i.e. $\langle |\hat {u}_p^ {{\prime \prime }}|^2\rangle > \langle |\hat {u}_m^ {{\prime \prime }}|^2\rangle$, is mainly concentrated at $\lambda _x>\lambda _z$, little influenced by ${\textit {Ma}}_b$. In contrast, $A_{pm,uu}$ at $\lambda _x<\lambda _z$ exhibits ${\textit {Ma}}_b$ dependence, where $\gamma _{uu}^2$ is low and the dominant POD mode for the cospectrum is not an advection mode (see § 5.4). In summary, $\gamma _{uu}^2$ and $A_{pm,uu}$ of the streamwise-elongated ($\lambda _x>\lambda _z$) structures with high coherence exhibit good ${\textit {Ma}}_b$ independence within and below the logarithmic region when using semi-local units. These structures have high $A_{pm,uu}$, so they are pronounced and tend to be dominant for the estimated $u^ {{\prime \prime }}(y_p^*)$, beneficial for constructing the compressible IOIM. In contrast, the spanwise-elongated structures with $\lambda _x<\lambda _z$ are more influenced by ${\textit {Ma}}_b$. They tend to be of secondary importance for SLSE because of the relatively low $\gamma _{uu}^2$ and $A_{pm,uu}$.

Figure 22. Contours of (a,c) $\gamma _{uu}^2$ and (b,d) $A_{pm,uu}$ in the cases of varying ${\textit {Ma}}_b$ but the same $Re_{\tau,c}^*$. Here, $y_m^*$ is fixed at $107$, and (a,b) $y_p^*=70$ and (c,d) $y_p^*=10$. The levels are labelled in the figures. The green dotted lines denote $\lambda _x=\lambda _z$.

The parameter study regarding temperature is not conducted for two reasons. First, heat transfer needs to be included in the incompressible case for reference (Antonia et al. 2009), but such cases are not currently available in our dataset. Second, the conclusions presented above concerning the effects of acoustic and entropy modes on temperature are limited to the case of ${\textit {Ma}}_b=1.5$, awaiting future examination at higher ${\textit {Ma}}_b$ using the DNS data.

7. Discussions

7.1. Other linear models with turbulence modelling

As introduced in § 3.2, there are a number of combinations of turbulence modelling relations to linearize the fluctuation equations, leading to different linear models. Since (3.6) is used, they can all be considered as eddy viscosity enhanced ones, so for ease of writing, the eLNS model presented above using (3.10a,b) is termed the benchmark eLNS model. We have examined other eLNS models, including

1. (i) using (3.8b) instead of (3.8a) to linearize the turbulent heat flux;

2. (ii) adding $k^ {{\prime \prime }}$ as an independent variable (so six equations with the basic variable set ${\boldsymbol q}^ {{\prime \prime }}=[\rho ^ {\prime },u^ {{\prime \prime }},v^ {{\prime \prime }},w^ {{\prime \prime }},T^ {{\prime \prime }},k^ {{\prime \prime }}]^ {\textrm {T}}$) to linearize the dissipation rate and the $k^ {{\prime \prime }}$ equation;

3. (iii) linearizing $\rho ^ {\prime }$ (in $N_\rho ^ {{\prime \prime }}$, $N_{u_i}^ {{\prime \prime }}$ and $N_T^ {{\prime \prime }}$) using the SRA in (3.16a,b) (similar to (3.8)).

They are termed the adjusted eLNS models collectively. Same as the benchmark eLNS model, non-negligible contributions of acoustic and entropy components to $\boldsymbol {\varPhi }$ are all present in these adjusted models, because the eigenmodes are damped by $\mu _t$ and the non-normality of the linear operator (in particular, the vortical mode branch) decreases as in § 4. Consequently, the SLSE results from these adjusted models are all largely different from DNS regarding the temperature fluctuations, similar to the benchmark eLNS model (see § 5.3). Thereby, detailed comparisons among these eLNS models are not provided considering the collective failure, and the main focus has been paid to the mathematical features of the linear model for possible explanations of the deviation, as has been elaborated in §§ 4 and 5.4. The procedure of cospectrum decomposition is thus introduced to obtain basically agreeable results with the DNS data.

7.2. Implications for calculating boundary layer flows

As shown in figure 13, there are a series of acoustic modes in the POD of $\boldsymbol {\varPhi }$ with different $k_y$, oscillating in the outer region. Bounded by the two walls at $y=0$ and $2h$ in the channel flow, $k_y$ is discretized so only those with $\lambda _y=2h/n$ are allowed ($n$ is an integer). In analogy, the acoustic modes exist in boundary layer flows, and their $k_y$ are determined by the height of the computational domain. Different from the channel flow, the location of the upper boundary is in the free stream and is artificially selected in the boundary layers, so the results of the acoustic modes and thus $\boldsymbol {\varPhi }$ can be dependent on the domain height, which seems not physically appropriate. When the response mode is relatively supersonic (phase speed $\omega /k_x$ larger than the free stream acoustic wave, either downstream or upstream), the acoustic components in the free stream are highly pronounced, so the responses from input–output or resolvent analyses can be tightly dependent on the computational domain, as noted by Bae et al. (2020) and Madhusudanan & McKeon (2022). The present work suggests that such a prominent role of acoustic modes in the linear models is not in line with DNS, within at least the ${\textit {Ma}}$ considered. As one solution candidate, the cospectrum decomposition procedure can be used to identify and remove the acoustic components. Detailed comparisons with experimental or DNS data for boundary layer flows are anticipated to assess further the behaviour of the linear models in different configurations.

7.3. Implications for the turbulence modelling simulations

The discussions in § 4 concerning $\hat {\mathcal {L}}_{eLNS}$ and $\mu _t$ in compressible flows may be of greater significance in some RANS-related applications. The linear operator can be encountered, e.g. in the sensitivity analysis of RANS modelling to field variables, and in some spatial numerical schemes in the computational fluid dynamics where linearization of the inviscid and viscous fluxes towards basic variables is required (see e.g. MacCormack 1985). Also, the linear operator can, at least partially, govern the transient dynamics when the turbulent flow is perturbed around its statistically equilibrium state, simulated using, e.g. unsteady RANS simulations, detached eddy simulations or wall-modelled large eddy simulations. The diminishing non-normality of the linear operator due to $\mu _t$ can overdamp the vortical modes, hence leading to possibly overemphasized roles of the acoustic and entropy modes. Nevertheless, this effect can be difficult to isolate due to the combined complexities of flow physics (shocks, separation, etc.), more sophisticated turbulence models, numerical schemes, grid resolutions and so on. More future explorations are anticipated in this direction.

8. Summary

In this work, the linear models, based on stochastically forced linearized NS equations, are deployed for SLSE of the velocity and temperature fluctuations in compressible turbulent channel flows. The benchmark case has ${\textit {Ma}}_b=1.5$ and $Re_b={20\ 020}$. The spectral coherence and amplitude ratio ($\gamma ^2$ and $A_{pm}$) of the fluctuation between different wall-normal heights are obtained using the linear models with only mean flow input. Afterwards, the coherent portion of the instantaneous field at other heights can be estimated through $H_L$ based on one snapshot at the measurement height. Three aspects are mainly focused on, namely the derivation and mathematical features of the linear models, quantitative comparison of the SLSE quantities with the DNS data, and the structural coherence of velocity and temperature fluctuations inferred by the SLSE results.

First, the derivation of the linear model is discussed carefully. The full expressions of the fluctuation equations are provided, and the assumptions for linearization are detailed. Two types of models are considered. The first is the LNS model with all the nonlinear terms collected into the forcing term, while the second is the eddy viscosity enhanced one (eLNS model), using turbulence modelling relations for possible linearization of the nonlinear terms. The mathematical properties of the linear operator are scrutinized by studying their eigenspectra and pseudospectra. There are four branches of eigenmodes for $\hat {\mathcal {L}}_{LNS}$ and $\hat {\mathcal {L}}_{eLNS}$, and the vortical branch has strong non-normality to support the transient behaviour of the fluctuations. For the LNS model, the acoustic modes have little non-normality and increasingly small contributions to the energy growth as length scales (both $\lambda _x$ and $\lambda _z$) increase. In contrast, the non-normality of the vortical modes diminishes in the eLNS model, closely related to the damping effects of $\mu _t$ (and $\kappa _t$) on the eigenmodes. Consequently, the acoustic and entropy modes can account for 20 % to 55 % of the response energy growth, which largely affects the SLSE results for temperature.

The SLSE results ($\gamma ^2$, $A_{pm}$ and $H_L$) for streamwise velocity are quite similar to the incompressible case. The eLNS model outperforms the LNS model, giving basically agreeable results with the DNS data, though the structures of small $\lambda _z$ are still somewhat underestimated. However, both models have considerable deviations from DNS regarding temperature and velocity–temperature coupling. A decomposition of the eLNS cospectrum ($\boldsymbol {\varPhi }_{nae}$ and $\boldsymbol {\varPhi }_{ae}$) shows that the acoustic and entropy components hardly affect the streamwise-velocity cospectrum, but have non-negligible contributions (especially the acoustic ones) to those of density and temperature, which is not supported by the DNS data. Thereby, the prominent role of acoustic and entropy components for the response growth is an intrinsic defect of the eLNS model introduced by the turbulence modelling related to $\mu _t$. A procedure of cospectrum decomposition is thus designed to remove the acoustic and entropy components (eLNS-CD). The resulting SLSE quantities for temperature and velocity–temperature coupling show a noticeable improvement over the eLNS model, and are basically agreeable with DNS, except in the region of small $\lambda _x$. Thereby, the coherent velocity and temperature fluctuations are dominated by advection and other vortical motions, while the compressibility effects are of secondary importance.

Finally, a parameter study is conducted for four cases with varying ${\textit {Ma}}_b$ (0.3–4) but the same $Re_{\tau,c}^*$. Semi-local units well collapse $\gamma _{uu}^2$ and $A_{pm,uu}$ to the incompressible case for the streamwise-elongated structures of high coherence, which benefits the modelling of inner–outer interaction for compressible flows. In contrast, the spanwise-elongated structures ($\lambda _x<\lambda _z$) are more influenced by ${\textit {Ma}}_b$. Moreover, the behaviours of other linear models with different turbulence modelling relations are discussed. Some implications are presented for calculating the boundary layer flows and dealing with RANS-related computations.

Future works will be focused on assessing the behaviours of the linear models for compressible cases at higher ${\textit {Ma}}_b$ and $Re$ and incompressible cases with heat transfer, using the DNS or experimental data. Improvements on the linear models are anticipated through, e.g. modified linearization of turbulent eddy terms and more realistic modelling of the forcing.