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 Reference Cheng and Fu2022b). 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 Reference Agostini and Leschziner2014). Readers can refer to Cheng & Fu (Reference Cheng and Fu2022b, Reference Cheng and Fu2023a) 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 Reference Griffin, Fu and Moin2022, Reference Griffin, Fu and Moin2023; Chen et al. Reference Chen, Cheng, Fu and Gan2023; Song, Zhang & Xia Reference Song, Zhang and Xia2023). The ODE-based method of Chen et al. (Reference Chen, Cheng, Fu and Gan2023) 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. Reference Tinney, Coiffet, Delville, Hall, Jordan and Glauser2006),

(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. Reference Baars, Hutchins and Marusic2016). 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 (Reference Cheng and Fu2022b) 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 Reference Cheng and Fu2023a). 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. (Reference Tinney, Coiffet, Delville, Hall, Jordan and Glauser2006), 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. Reference Baars, Hutchins and Marusic2016). 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.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 (Reference Gatski and Bonnet2013), 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. Reference Huang, Coleman and Bradshaw1995). 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 Reference Lele1994). 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 Reference Yao and Hussain2020). 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. Reference Huang, Coleman and Bradshaw1995) 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. (Reference Bae, Dawson and McKeon2020) and Dawson & McKeon (Reference Dawson and McKeon2020), 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 Reference Mack1984).

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. (Reference Huang, Coleman and Bradshaw1995) 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 (Reference Wilcox2006), 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. Reference Huang, Coleman and Bradshaw1995). 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 Reference Jones and Launder1972), 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. Reference Alizard, Pirozzoli, Bernardini and Grasso2015; Pickering et al. Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021; Chen et al. Reference Chen, Cheng, Fu and Gan2023), 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. (Reference Chen, Cheng, Fu and Gan2023). 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. Reference Madhusudanan, Illingworth and Marusic2019). 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 Reference Chu1965)

(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. Reference Gupta, Madhusudanan, Wan, Illingworth and Juniper2021), 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 (Reference Farrell and Ioannou1993) 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. Reference Bae, Dawson and McKeon2020). 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. (Reference Chen, Cheng, Fu and Gan2023).

In the incompressible model of Madhusudanan et al. (Reference Madhusudanan, Illingworth and Marusic2019), $\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. Reference Alizard, Gibis, Selent, Rist and Wenzel2022; Chen et al. Reference Chen, Cheng, Fu and Gan2023). To improve the estimation below the logarithmic region, Gupta et al. (Reference Gupta, Madhusudanan, Wan, Illingworth and Juniper2021) 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 Reference Foysi, Sarkar and Friedrich2004), 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. (Reference Coleman, Kim and Moser1995) 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. Reference Huang, Coleman and Bradshaw1995; Cheng & Fu Reference Cheng and Fu2022b). 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.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 Reference Balakumar and Malik1992). 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 Reference Trefethen1997). 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 (Reference Chomaz2005), 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 (Reference Hanifi, Schmid and Henningson1996) 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. (Reference Symon, Illingworth and Marusic2021) and Kuhn et al. (Reference Kuhn, Müller, Knechtel, Soria and Oberleithner2022) 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. (Reference Symon, Illingworth and Marusic2021) 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 Reference Hwang and Cossu2010). 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. Reference Hanifi, Schmid and Henningson1996), 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. Reference Hanifi, Schmid and Henningson1996). 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 Reference George and Sujith2011). 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.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 Reference Smith and Smits1995; Ringuette et al. Reference Ringuette, Wu and Martín2008; Pirozzoli & Bernardini Reference Pirozzoli and Bernardini2011; Cogo et al. Reference Cogo, Salvadore, Picano and Bernardini2022; Cheng & Fu Reference Cheng and Fu2023a), so here we only provide a brief description to set the grounds for later discussions.

After the velocity transformation of Griffin et al. (Reference Griffin, Fu and Moin2021), 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 (Reference Abe and Antonia2017) 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 Reference Hamilton, Kim and Waleffe1995; Schoppa & Hussain Reference Schoppa and Hussain2002). Meanwhile, $u^{ {{\prime \prime }}+}$ and $T^{ {{\prime \prime }}+}$ strongly resemble each other morphologically. As a quantitative measure, Cheng & Fu (Reference Cheng and Fu2023a) 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 (Reference Cheng and Fu2023a) 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. Reference Marusic, Mathis and Hutchins2010), and the estimations are conducted at $y_p^*< y_m^*$. In analogy to the incompressible flow (Townsend Reference Townsend1976), 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. Reference Madhusudanan, Illingworth and Marusic2019; Gupta et al. Reference Gupta, Madhusudanan, Wan, Illingworth and Juniper2021), 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. Reference Baars, Hutchins and Marusic2016; Madhusudanan et al. Reference Madhusudanan, Illingworth and Marusic2019), 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 Reference Yao and Hussain2020). 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. (Reference Baars, Hutchins and Marusic2016), 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. Reference Madhusudanan, Illingworth and Marusic2019). 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. Reference Baars, Hutchins and Marusic2016). 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. Reference Gupta, Madhusudanan, Wan, Illingworth and Juniper2021). 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 Reference Cheng and Fu2023a). 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. Reference Coleman, Kim and Moser1995; Williams et al. Reference Williams, Sahoo, Baumgartner and Smits2018; Chen et al. Reference Chen, Cheng, Fu and Gan2023). 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 Reference McKeon2019). 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 Reference Hutchins and Marusic2007; Abe et al. Reference Abe, Antonia and Toh2018). 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 (Reference Cheng and Fu2023a) 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 (Reference Cheng and Fu2023a).