2.1. Mean flow and fluctuation equations

We consider canonical compressible turbulent channel flow under the calorically perfect gas assumption. To obtain reliable statistics from the DNS, we need to carefully deal with each term in the mean and fluctuation equations, especially those regarding high-order fluctuation terms and the difference between Reynolds and Favre averages. The equations in compressible flows are far more complicated than the incompressible version, so it is necessary to first present the complete form of the equations to be dealt with.

The Favre-averaged quantities are chosen as basic variables for the linear models, for ease of realising a thorough linearisation (Chen et al. Reference Chen, Cheng, Gan and Fu2023b ). The continuity, momentum and enthalpy equations of the mean flow are written as (e.g. Gatski & Bonnet Reference Gatski and Bonnet2013)

(2.1a) $${\frac {\partial \bar {\rho }}{\partial t} + \frac {\partial \bar {\rho } \tilde {u}_j }{\partial x_j} = 0, }$$

(2.1b) $${\bar {\rho } \left ({ \frac {\partial \tilde {u}_i}{\partial t} + \tilde {u}_j \frac {\partial \tilde {u}_i}{\partial x_j} }\right ) = - \frac {\partial }{\partial x_j} \left ({ \bar {\rho } \widetilde {u_i^{{\prime \prime }} u_j^{{\prime \prime }}} }\right ) - \frac {\partial \bar {p}}{\partial x_i} + \frac {\partial \bar {\tau }_{{i\!j}}}{\partial x_j}, \quad \text{for}\; i=1,2,3,}$$

(2.1c) {$$\begin{aligned}& \frac {\partial }{\partial t} \left( \bar {\rho } c_p \widetilde {T} \right) + \frac {\partial }{\partial x_j} \left ({ \bar {\rho } c_p \tilde {u}_j \widetilde {T} }\right ) - \left ( \frac {\partial \bar {p}}{\partial t} + \tilde {u}_j \frac {\partial \bar {p}}{\partial x_j} + \overline {u_j^{{\prime \prime }} \frac {\partial p^{\prime }}{\partial x_j}} + \overline {u_j^{{\prime \prime }}} \frac {\partial \bar {p}}{\partial x_j} \right ) \\[-1mm] &\quad = - \frac {\partial }{\partial x_j} \left ({ \bar {\rho } c_p \widetilde {u_j^{{\prime \prime }} T^{{\prime \prime }}} }\right ) - \frac {\partial \bar {\vartheta }_j}{\partial x_j} + \bar {\tau }_{{i\!j}} \frac {\partial \tilde {u}_i}{\partial x_j} + \overline {\tau _{{i\!j}} \frac {\partial u_i^{{\prime \prime }}}{\partial x_{j_{_{_{_{\!}}}}}}}. \end{aligned}$$

Here, $\bar {\varphi }$ is the temporal (Reynolds) average and $\tilde {\varphi }=\overline {\rho \varphi }/\bar {\rho }$ is the Favre average; $\rho$ , $u_i=[u,v,w]^{T}$ and $T$ are the fluid’s density, velocities and temperature, respectively; the mean pressure is $\bar {p}=\bar {\rho }R\widetilde {T}$ , with $R$ the gas constant and $c_p$ is the isobaric specific heat. Note that the temporal derivatives are retained in (2.1) for reference for the fluctuation equations. No-slip and isothermal walls are set on both sides as the boundary condition. The mean viscous stress $\bar {\tau }_{{i\!j}}$ and heat flux $\bar {\vartheta }_j$ are calculated as

(2.2) \begin{equation} \bar {\tau }_{{i\!j}} = \underbrace{\,\bar{\!\mu} \!\left (\!{ \frac {\partial \tilde {u}_i}{\partial x_j} \!+\! \frac {\partial \tilde {u}_j}{\partial x_i} }\!\right ) - \frac {2}{3} \,\bar {\!\mu } \frac {\partial \tilde {u}_k}{\partial x_k} \delta _{{i\!j}}}_{\equiv \bar {\tau }_{{i\!j},\textit {lin}}} + \underbrace {\overline {\mu \frac {\partial u_i^{{\prime \prime }}}{\partial x_j}} \!+\! \overline {\mu \frac {\partial u_j^{{\prime \prime }}}{\partial x_i}} - \frac {2}{3} \overline {\mu \frac {\partial u_k^{{\prime \prime }}}{\partial x_k}} \delta _{{i\!j}}}_{\equiv \bar {\tau }_{{i\!j},\textit {nln}}}, \, \bar {\vartheta }_j = - \bar {\kappa } \frac {\partial \widetilde {T}}{\partial x_j} - \overline {\kappa \frac {\partial T^{{\prime \prime }}}{\partial x_j}}, \end{equation}

where $\mu (T)$ and $\kappa =\mu c_p/{\textit {Pr}}$ are the viscosity and thermal conductivity; $\mu$ is calculated by Sutherland’s law and ${\textit {Pr}}=0.72$ is the Prandtl number; $\delta _{{i\!j}}$ represents the Kronecker delta. Note that $\bar {\tau }_{{i\!j}}$ is divided into two parts, as underbraced, for later linearisation of the fluctuation equations. The first part $\bar {\tau }_{{i\!j},\textit {lin}}$ is a function of mean-flow variables only, while the second part $\bar {\tau }_{{i\!j},\textit {nln}}$ is contributed by the second-order moment of the fluctuations.

The fluctuation equations are far more complex than the incompressible counterparts. Following the derivation in Chen et al. (Reference Chen, Cheng, Gan and Fu2023b ), the fluctuation equations relative to (2.1) are written as

(2.3a) $${\qquad \qquad \qquad \qquad \qquad \underbrace {\frac {\partial \rho ^{\prime }}{\partial t} + \frac {{\partial \bar {\rho } u_j^{{\prime \prime }}}}{\partial x_j} + \frac {\partial \rho ^{\prime } \tilde {u}_j }{\partial x_j}}_{\textit{linear terms}} = \underbrace {-\frac {\partial \rho ^{\prime } u_j^{{\prime \prime }}}{\partial x_j}}_{{\textit{nonlinear term}}},}$$

(2.3b) $${\quad \begin{aligned}& \bar {\rho } \left ({ \frac {\partial u_i^{{\prime \prime }}}{\partial t} + \tilde {u}_j \frac {\partial u_i^{{\prime \prime }}}{\partial x_j} }\right ) + \big({ \rho ^{\prime } \tilde {u}_j + \bar {\rho } u_j^{{\prime \prime }}} \big) \frac {\partial \tilde {u}_i}{\partial x_j} + \frac {\partial p_{\textit {lin}}^{\prime }}{\partial x_i} - \frac {\partial \tau _{{i\!j},\textit {lin}}^{\prime }}{\partial x_j} \\[2pt] &\quad = \underbrace {- \frac {\partial }{\partial x_j} \left ({ \rho u_i^{{\prime \prime }} u_j^{{\prime \prime }} - \bar {\rho } \widetilde {u_i^{{\prime \prime }} u_j^{{\prime \prime }}} }\right )}_{{\textit{RSF}}} \underbrace {- \frac {\partial p_{\textit {nln}}^{\prime }}{\partial x_i}}_{\textit{PGNF}} + \underbrace {\frac {\partial \tau _{{i\!j},\textit {nln}}^{\prime }}{\partial x_j}}_{\textit{VSNF}} + \,\textit{MMNF}, \end{aligned}}$$

(2.3c) \begin{aligned} & \bar {\rho } \left ({ c_v \frac {\partial T^{{\prime \prime }}}{\partial t} + c_p \tilde {u}_j \frac {\partial T^{{\prime \prime }}}{\partial x_j} }\right ) - R\widetilde {T} \frac {\partial \rho ^{\prime }}{\partial t} + \big({ \rho ^{\prime } \tilde {u}_j + \bar {\rho } u_j^{{\prime \prime }} }\big ) c_p \frac {\partial \widetilde {T}}{\partial x_j} - \left ( \tilde {u}_j \frac {\partial p_{\textit {lin}}^{\prime }}{\partial x_j} + u_j^{{\prime \prime }} \frac {\partial \bar {p}}{\partial x_j} \right ) \\[2pt] &\qquad + \frac {\partial \vartheta _{j,\textit {lin}}^{\prime }}{\partial x_j} - \bar {\tau }_{{i\!j},\textit {lin}} \frac {\partial u_i^{{\prime \prime }}}{\partial x_j} - \tau _{{i\!j},\textit {lin}}^{\prime } \frac {\partial \tilde {u}_i}{\partial x_j} = \underbrace {- c_p \frac {\partial }{\partial x_j} \left ({ \rho u_j^{{\prime \prime }} T^{{\prime \prime }} - \bar {\rho } \widetilde {u_j^{{\prime \prime }} T^{{\prime \prime }}} }\right )}_{\textit{THF}} \underbrace {- \frac {\partial \vartheta _{j,\textit {nln}}^{\prime }}{\partial x_j}}_{\textit{MHNF}} \\[2pt] &\qquad + \underbrace {\bar {\tau }_{{i\!j},\textit {nln}} \frac {\partial u_i^{{\prime \prime }}}{\partial x_j} + \tau _{{i\!j},\textit {nln}}^{\prime } \frac {\partial \tilde {u}_i}{\partial x_j} + \left ({ \tau _{{i\!j}}^{\prime } \frac {\partial u_i^{{\prime \prime }}}{\partial x_j} - \overline {\tau _{{i\!j}} \frac {\partial u_i^{{\prime \prime }}}{\partial x_j}} }\right )}_{\textit{VDNF}} + \,{\textit{PCNF}} + {\textit{EMNF}}. \end{aligned}

Here, $\varphi ^{\prime }$ and $\varphi ^{{\prime \prime }}$ are the Reynolds and Favre fluctuations, respectively; the linear and nonlinear fluctuation terms are listed on the two sides of each equation. There are five independent variables in the system. We choose the basic variable set as ${\boldsymbol q}^{{\prime \prime }}=[\rho ^{\prime },\,u^{{\prime \prime }}, v^{{\prime \prime }},\,w^{{\prime \prime }},\,T^{{\prime \prime }}]^{T}$ , then the derived variables $p^{\prime }$ , $\vartheta _j^{\prime }$ and $\tau _{{i\!j}}^{\prime }$ are not strictly linear functions of ${\boldsymbol q}^{{\prime \prime }}$ . To achieve a thorough linearisation of the system, $p^{\prime }$ and $\vartheta _j^{\prime }$ are also divided into the linear and nonlinear parts in terms of ${\boldsymbol q}^{{\prime \prime }}$ , as

(2.4) \begin{equation} p^{\prime } = \underbrace {\bar {\rho } T^{{\prime \prime }} + \rho ^{\prime } \widetilde {T}}_{\equiv\; p_{lin}^{\prime }} + \underbrace {\rho ^{\prime } T^{{\prime \prime }}}_{\equiv\; p_{nln}^{\prime }}, \qquad \vartheta _j^{\prime } = \underbrace {-\bar {\kappa } \frac {\partial T^{{\prime \prime }}}{\partial x_j} - \kappa ^{\prime } \frac {\partial \widetilde {T}}{\partial x_j}}_{\equiv\; \vartheta _{j,{lin}}^{\prime }} \underbrace {- \left ( \kappa ^{\prime } \frac {\partial T^{{\prime \prime }}}{\partial x_j} - \overline {\kappa \frac {\partial T^{{\prime \prime }}}{\partial x_j}} \right )}_{\equiv\; \vartheta _{j,{nln}}^{\prime }}, \end{equation}

and so is $\tau _{{i\!j}}^{\prime }=\tau _{{i\!j},\textit {lin}}^{\prime } + \tau _{{i\!j},\textit {nln}}^{\prime }$ (not detailed for brevity). As a result, the left-hand sides in (2.3) are all linear functions of ${\boldsymbol q}^{{\prime \prime }}$ , ready for constructing the linear operator. The various nonlinear fluctuation terms reside on the right-hand sides of (2.3), classified as underbraced. Their names are collectively defined in table 1. The expressions of the unspecified terms MMNF, PCNF, EMNF are detailed in Appendix A.

Table 1. Names and abbreviations of the nonlinear terms in the fluctuation equations.

Notably, the term RSF is the only nonlinear term in incompressible flows. For compressible cases, however, many other terms appear due to the thermodynamic processes. In previous works (e.g. Alizard et al. Reference Alizard, Pirozzoli, Bernardini and Grasso2015), RSF and THF were of the primary concern and were modelled using algebraic RANS relations to formulate the eLNS model; other terms were usually neglected. But here, we provide the complete form of the fluctuation equations and will clarify in § 4.1 the contribution of each nonlinear term.

2.2. The LNS and eLNS models

Equation (2.3) constitutes a closed set of equations for ${\boldsymbol q}^{{\prime \prime }}$ , and can be rearranged into an operator form as

(2.5a) \begin{equation} \dfrac {\partial \rho ^{\prime }}{\partial t} - \mathcal{L}_\rho {\boldsymbol q}^{{\prime \prime }} = f_\rho ^{{\prime \prime }}, \quad \bar {\rho } \dfrac {\partial u_i^{{\prime \prime }}}{\partial t} - \mathcal{L}_{m_i} {\boldsymbol q}^{{\prime \prime }} = f_{m_i}^{{\prime \prime }}, \quad \bar {\rho } c_v \dfrac {\partial T^{{\prime \prime }}}{\partial t} - R\widetilde {T} \dfrac {\partial \rho ^{\prime }}{\partial t} - \mathcal{L}_{h} {\boldsymbol q}^{{\prime \prime }} = f_h^{{\prime \prime }}, \end{equation}

where $\mathcal{L}$ denotes the linear operator determined only by the mean flow $\widetilde {\boldsymbol q}$ , and $f^{{\prime \prime }}$ represents the nonlinear terms. Detailed expressions of $\mathcal{L}$ are readily contained in (2.3), also available in Dawson & McKeon (Reference Dawson and McKeon2020) and Chen et al. (Reference Chen, Cheng, Fu and Gan2023a ). Equation (2.5a ) is rewritten into a matrix form as

(2.5b) \begin{equation} \boldsymbol{A} \frac {\partial }{\partial t} \left [\!\begin{array}{*{20}{c}} \rho ^{\prime } \\ u_j^{{\prime \prime }} \\ T^{{\prime \prime }} \end{array}\!\right ] = \left [\!\begin{array}{*{20}{c}} \mathcal{L}_\rho \\ \mathcal{L}_{m_i} \\ \mathcal{L}_{h} \end{array}\!\right ] {\boldsymbol q}^{{\prime \prime }} + \left [\!\begin{array}{*{20}{c}} f_\rho ^{{\prime \prime }} \\ f_{m_i}^{{\prime \prime }} \\ f_h^{{\prime \prime }} \end{array}\!\right ], \qquad \boldsymbol{A} \equiv \left [ \begin{array}{*{20}{c}} 1&{}&{}\\ {}&\bar {\rho }\delta _{{i\!j}}&{} \\ -R\widetilde {T}&{}&\bar {\rho } c_v \end{array} \right ], \end{equation}

or equivalently

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

where $\boldsymbol{A}^{-1}$ has been absorbed into $\mathcal{L}_q$ and ${\boldsymbol f}_q^{{\prime \prime }}\equiv [f_\rho ^{{\prime \prime }},\,f_u^{{\prime \prime }},\,f_v^{{\prime \prime }},\,f_w^{{\prime \prime }},\,f_T^{{\prime \prime }}]^{\textrm{T}}$ . Equation (2.6) is the linearised NS equation and serves as a classic model problem. The linear model built upon (2.6) is often referred to as the LNS model.

In addition to the LNS model, the compressible eLNS model is frequently considered for model improvement, which is built upon algebraic RANS models (Alizard et al. Reference Alizard, Pirozzoli, Bernardini and Grasso2015; Chen et al. Reference Chen, Gan and Fu2024, Reference Chen, Gan and Fu2025). Specifically, the Boussinesq assumption and the strong Reynolds analogy (SRA) are introduced to linearise the Reynolds stress $\tau _{R,{i\!j}}$ and turbulent heat flux $\vartheta _{R,j}$

(2.7) \begin{equation} \tau _{R,{i\!j}}^{\prime } = \mu _t \left ({ \frac {\partial u_i^{{\prime \prime }}}{\partial x_j} + \frac {\partial u_j^{{\prime \prime }}}{\partial x_i} }\right ) - \frac {2}{3}\mu _t \frac {\partial u_k^{{\prime \prime }}}{\partial x_k} \delta _{{i\!j}}, \qquad \vartheta _{R,j}^{\prime } = -\kappa _t \frac {\partial T^{{\prime \prime }}}{\partial x_j}. \end{equation}

Here, $\kappa _t$ is the eddy diffusivity. In practice, $\mu _t$ and $\kappa _t$ are either assumed known from the DNS data or from semi-empirical models (e.g. Fan et al. Reference Fan, Kozul, Li and Sandberg2024). The terms $\tau _R^{\prime }$ and $\vartheta _R^{\prime }$ were originally designed to partially model the colour of the nonlinear terms. Taking the RSF term as an example, it is modelled as

(2.8) \begin{equation} -\frac {\partial }{\partial x_j} \left ({ \rho u_i^{{\prime \prime }} u_j^{{\prime \prime }} - \bar {\rho } \widetilde {u_i^{{\prime \prime }} u_j^{{\prime \prime }}} }\right ) \equiv {\textit{RSF}}_i = \frac {\partial \tau _{R,{i\!j}}^{\prime }}{\partial x_j} + {\textit{eRSF}}_i, \end{equation}

where ${\rm eRSF}$ represents the residual nonlinear term of RSF, after subtracting the modelled stress flux ${\partial \tau _{R,{i\!j}}^{\prime }}/{\partial x_j}$ . The other two residual terms $\rm{eTHF}$ and ${\rm eVDNF}$ , after modelling the turbulent heat flux and viscous dissipation, are similarly defined (e.g. Chen et al. Reference Chen, Cheng, Gan and Fu2023b ).

It is worth emphasising that (2.8) is the core assumption for the eLNS model. However, we will show in § 4 that (2.8) can deviate from its original intention. The supposedly small residual terms, ${\rm{eRSF}}$ and $\rm{eTHF}$ , turn out to essentially rebuild the forcing distributions.

Since $\mu _t$ and $\kappa _t$ profiles are assumed fixed, $\tau _R^{\prime }$ and $\vartheta _R^{\prime }$ are linear in terms of ${\boldsymbol q}^{{\prime \prime }}$ , so these eddy-viscosity-related terms in the momentum and enthalpy equations are linearised as

(2.9) \begin{equation} \frac {\partial \tau _{R,{i\!j}}^{\prime }}{\partial x_j} \equiv \mathcal{M}_{m_i} {\boldsymbol q}^{{\prime \prime }}, \quad - \frac {\partial \vartheta _{R,j}^{\prime }}{\partial x_j} + \bar {\tau }_{R,{i\!j}} \frac {\partial u_i^{{\prime \prime }}}{\partial x_j} + \tau _{R,{i\!j}}^{\prime } \frac {\partial \tilde {u}_i}{\partial x_j} \equiv \mathcal{M}_{h} {\boldsymbol q}^{{\prime \prime }}. \end{equation}

Consequently, parts of the LNS nonlinear terms can be linearised and moved into the linear operator, so that the model is eddy-viscosity enhanced, following the term for incompressible flows (Madhusudanan, Illingworth & Marusic Reference Madhusudanan, Illingworth and Marusic2019). The resulting operator form of the eLNS model reads

(2.10) \begin{equation} \left \{ \!\!\begin{array}{l} \dfrac {\partial \rho ^{\prime }}{\partial t} = \mathcal{E}_\rho {\boldsymbol q}^{{\prime \prime }} + e_\rho ^{{\prime \prime }}, \quad \text{with $\mathcal{E}_{\rho }\equiv \mathcal{L}_{\rho }$ and $e_\rho ^{{\prime \prime }}\equiv f_\rho ^{{\prime \prime }}$},\\[10pt] \displaystyle \bar {\rho } \dfrac {\partial u_i^{{\prime \prime }}}{\partial t} = \underbrace {\left ( \mathcal{L}_{m_i} + \mathcal{M}_{m_i} \right )}_{\equiv\; \mathcal{E}_{m_i}} {\boldsymbol q}^{{\prime \prime }} + \underbrace {\left ( f_{m_i}^{{\prime \prime }} - \mathcal{M}_{m_i} {\boldsymbol q}^{{\prime \prime }} \right )}_{\equiv\; e_{m_i}^{{\prime \prime }}},\\[10pt] \displaystyle \bar {\rho } c_v \dfrac {\partial T^{{\prime \prime }}}{\partial t} - R\widetilde {T} \dfrac {\partial \rho ^{\prime }}{\partial t} = \underbrace {\left ( \mathcal{L}_{h} + \mathcal{M}_{h} \right )}_{\equiv\; \mathcal{E}_{h}} {\boldsymbol q}^{{\prime \prime }} + \underbrace {\left ( f_{h}^{{\prime \prime }} - \mathcal{M}_{h} {\boldsymbol q}^{{\prime \prime }} \right )}_{\equiv\; e_{h}^{{\prime \prime }}}, \end{array} \right . \end{equation}

where the linear operator $\mathcal{E}$ and the nonlinear term ${\boldsymbol e}^{{\prime \prime }}$ for eLNS are defined as underbraced. The final standard form similar to (2.6) is

(2.11) \begin{equation} \frac {\partial {\boldsymbol q}^{{\prime \prime }}}{\partial t} = \mathcal{E}_q {\boldsymbol q}^{{\prime \prime }} + {\boldsymbol e}_q^{{\prime \prime }}. \end{equation}

Note that $\mathcal{E}_q$ is the same as $\mathcal{L}_q$ except for replacing $\,\bar {\!\mu }$ with $\,\bar {\!\mu }+\mu _t$ and $\bar {\kappa }$ with $\bar {\kappa }+\kappa _t$ . Equations (2.6) and (2.11) constitute the model equations for LNS and eLNS, respectively.

It is worth mentioning that the eLNS model was originally derived based on the triple decomposition of $\boldsymbol q$ (Reynolds & Hussain Reference Reynolds and Hussain1972; Alizard et al. Reference Alizard, Pirozzoli, Bernardini and Grasso2015), where $\boldsymbol q$ is decomposed into the mean part, the turbulent fluctuation part and a small-amplitude organised perturbation. This approach is more physically interpretable, whereas the final linearised equation is mathematically equivalent to (2.11). Therefore, we do not present the triply decomposed fluctuation equations, considering the already complicated (2.3).

2.3. The DNS datasets

A series of DNSs for compressible turbulent channel flows have been conducted by the present authors’ group, as described at length in Cheng & Fu (Reference Cheng and Fu2022, Reference Cheng and Fu2023). The simulations adopt three bulk Mach numbers ${\textit {Ma}}_b=U_b/a_w=0.8, 1.5, 3.0$ and a series of bulk Reynolds numbers $\textit{Re}_b=\rho _b U_b h/\mu _w$ . Here, $\rho _b$ and $U_b$ are the bulk density and streamwise velocity, $a=(\gamma R T)^{1/2}$ is the speed of sound with the specific ratio $\gamma =1.4$ , $h$ is the channel half-width and the subscript $w$ denotes wall quantities. Our focus is the compressibility effects on the forcing models, so we first select two DNS cases at ${\textit {Ma}}_b=1.5$ and 3.0 with nearly equal $\textit{Re}_\tau ^*$ . A third case with a higher $\textit{Re}$ at ${\textit {Ma}}_b=1.5$ is also included, to examine potential scalings. Additionally, two incompressible DNS cases at different $\textit{Re}_\tau$ are adopted, with data from Ying et al. (Reference Ying, Liang, Li and Fu2023), to facilitate quantitative measure of the compressibility effects. Detailed parameters of the five cases are listed in table 2. Here, $\textit{Re}_\tau =\rho _w u_\tau h/\mu _w$ and $\textit{Re}_\tau ^*=\bar {\rho }u_\tau ^* h/\,\bar {\!\mu }$ are the friction and semi-local Reynolds numbers; $u_\tau$ is the friction velocity and $u_\tau ^*$ is the semi-local counterpart; $\widetilde {T}_c$ is the temperature at the channel centre; $L_x$ and $L_z$ are the domain sizes. For later use, the wall viscous and semi-local length scales are defined as $\delta _\nu =\mu _w/(\rho _wu_\tau )$ and $\delta _\nu ^*=\,\bar {\!\mu }/(\bar {\rho }u_\tau ^*)$ , based on which the two sets of non-dimensional lengths are denoted by superscripts $+$ and $*$ , respectively.

Table 2. Parameters of the DNS cases for turbulent channel flows, where $t_{\textit {total}} u_{\tau }/h$ is the total eddy turnover time to accumulate statistics. Two incompressible cases are also included for later reference.

To obtain the real forcing statistics, the nonlinear terms ${\boldsymbol f}_q^{{\prime \prime }}$ and ${\boldsymbol e}_q^{{\prime \prime }}$ are calculated by their definitions in §§ 2.1, 2.2 using the DNS data. The linear models in this work are solved within the framework of the stochastic Lyapunov equation (see § 2.4), so we do not need time-resolved DNS data of massive memory requirement as in Morra et al. (Reference Morra, Nogueira, Cavalieri and Henningson2021) and Ying et al. (Reference Ying, Chen, Li and Fu2024a ) for incompressible cases. Instead, we focus on the ensemble (temporal) average $\langle \cdot \rangle$ of variables, where time-not-resolved data are adequate. Therefore, the sampling number of the instantaneous fields $N_t$ can be orders of magnitude lower than that required by the time-resolved data; for the latter case, $N_t$ can reach ${\sim}10^4$ at current $\textit{Re}_\tau$ . The price is that the frequency spectrum is not resolved, but we can still reveal the forcing statistics in the spatial spectra, which provides adequate information to improve the stochastic linear models. The Fourier decomposition is applied on ${\boldsymbol q}^{{\prime \prime }}$ in two homogenous directions as

(2.12) \begin{equation} {\boldsymbol q}^{{\prime \prime }}\! \left (x,y,z,t\right ) = \iint \nolimits _{-\infty }^{\infty } {\hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \left (y, t; k_x, k_z\right ) \exp \left [{\mathrm{i}} \left (k_x x + k_z z\right )\right ] {\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}^{{\kern0.3pt}{\prime \prime }}$ is the shape function. The Fourier components of $\,\hat {\!\boldsymbol f}^{{\prime \prime }}$ and $\,\hat {\!\boldsymbol e}^{{\prime \prime }}$ are similarly defined. Thevalue of $N_t$ used for different cases is listed in table 2, which will be shown in § 3 to be large enough to provide converged spectral statistics.

2.4. Fluctuations in response to modelled and DNS forcing

In this subsection, we first present how the spectral correlation tensor $\widehat {\boldsymbol{\Phi }}(y,y^{\prime };k_x,k_z)= \langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }}(y,t;k_x,k_z) \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}}(y^{\prime },t;k_x,k_z) \rangle$ in response to a modelled forcing is obtained in the linear model, which requires only the mean flow (including $\mu _t$ , $\kappa _t$ ) as the input. Afterwards, we present how the real forcing can be computed from the DNS data, hence enabling a direct interpretation of the model errors.

In the Fourier space, the model equations (2.6) and (2.11) are expressed as

(2.13a,b) \begin{equation} \frac {\partial \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }}}{\partial t} = \widehat {\mathcal{L}}_q \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} + \,\hat {\!\boldsymbol f}_q^{{\prime \prime }}, \qquad \frac {\partial \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }}}{\partial t} = \widehat {\mathcal{E}}_q \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} + \,\hat {\!\boldsymbol e}_q^{{\prime \prime }}. \end{equation}

When the real forcings $\,\hat {\!\boldsymbol f}_q^{{\prime \prime }}$ and $\,\hat {\!\boldsymbol e}_q^{{\prime \prime }}$ are unknown, the stochastic linear models provide a mean to estimate $\widehat {\boldsymbol{\Phi }}$ . An idealised stochastic forcing is assumed as $\,\hat {\!\boldsymbol f}_q=\boldsymbol{F}\,\hat {\!\boldsymbol f}_0$ (so is $\,\hat {\!\boldsymbol e}_q=\boldsymbol{E}\,\hat {\!\boldsymbol e}_0$ ), where the matrix $\boldsymbol{F}(y)$ models the wall-normally varied forcing amplitude, and $\,\hat {\!\boldsymbol f}_0$ (and $\,\hat {\!\boldsymbol e}_0$ ) is a delta-correlated Gaussian white noise with zero mean. Consequently, the correlation tensor $\widehat {\boldsymbol{\Phi }}$ is obtained by solving the Lyapunov equation (Farrell & Ioannou Reference Farrell and Ioannou1993), as

(2.14a) $${\widehat {\mathcal{L}}\, \widehat {\boldsymbol{\Phi }} + \widehat {\boldsymbol{\Phi }} \widehat {\mathcal{L}}^{{H}} + (\boldsymbol{{FF}}^{{H}}) = \boldsymbol{0}, \qquad \text{for the LNS model},}$$

(2.14b) $${\widehat {\mathcal{E}}\, \widehat {\boldsymbol{\Phi }} + \widehat {\boldsymbol{\Phi }} \widehat {\mathcal{E}}^{{H}} + (\boldsymbol{{EE}}^{{H}}) = \boldsymbol{0}, \qquad \text{for the eLNS model},}$$

where the Hermitian transpose of $\widehat {\mathcal{L}}$ and $\widehat {\mathcal{E}}$ is defined as the discrete adjoint. In the simplest models, $\boldsymbol{F}$ and $\boldsymbol{E}$ are assumed to be diagonal, suggesting perfect zero correlation of the forcing between different variables and between different wall-normal heights. The simple yet well-behaved W-model (Gupta et al. Reference Gupta, Madhusudanan, Wan, Illingworth and Juniper2021; Chen et al. Reference Chen, Cheng, Gan and Fu2023b ) is considered here, where the forcing amplitude varies in proportion to the kinematic eddy viscosity $\nu _t=\mu _t/\bar {\rho }$ , as

(2.15) \begin{equation} \boldsymbol{F},\boldsymbol{E} = \text{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 \widetilde {T}}{\partial \tilde {u}} \nu _t\right ]^{\textrm{T}} \right ) \,\boldsymbol{M}^{-1/2}. \end{equation}

The weight matrix $\boldsymbol{M}$ will be defined later, and the extended SRA (Huang, Coleman & Bradshaw Reference Huang, Coleman and Bradshaw1995) has been utilised in deriving (2.15) as

(2.16) \begin{equation} T^{{\prime \prime }}_{\textit {rms}} \approx \frac {1}{{\textit {Pr}}_t} \left |\frac {\partial \widetilde {T}}{\partial \tilde {u}} \right | u^{{\prime \prime }}_{\textit {rms}}, \qquad \rho ^{\prime }_{\textit {rms}} \approx \frac {1}{{\textit {Pr}}_t} \left |\frac {\partial \bar {\rho }}{\partial \tilde {u}}\right | u^{{\prime \prime }}_{\textit {rms}}, \end{equation}

where $rms$ denotes the root mean square. Equation (2.16) also suggests two wall viscous units for the temperature and density, $T_\tau =(\partial \widetilde {T}/\partial \tilde {u})_w u_\tau /{\textit {Pr}}$ and $\rho _\tau =(\partial \bar {\rho }/\partial \tilde {u})_w u_\tau$ .

The eigenmodes of $\widehat {\boldsymbol{\Phi }}$ , known as the proper orthogonal decomposition (POD) or Karhunen–Loève modes, are of interest in various analyses on turbulence. They are defined as

(2.17) \begin{equation} \left(\widehat {\boldsymbol{\Phi }}, \,\check {\!\boldsymbol q}_j^{{\prime \prime }}\right)_E = \theta _j \,\check {\!\boldsymbol q}_j^{{\prime \prime }}, \qquad j =1,2,\ldots , \end{equation}

where the eigenvalues and eigenfunctions $\theta _j$ and $\,\check {\!\boldsymbol q}_j^{{\prime \prime }}$ are sorted in the descending order of the energy of the $j$ th mode ( $\theta _j$ ). These POD modes are orthogonal to each other under the energy norm $(\cdot ,\cdot )_E$ , which is defined following common usage (Chu Reference Chu1965) as

(2.18) \begin{equation} \|\,\check {\!\boldsymbol q}^{{\prime \prime }}\|^2 = \left( \,\check {\!\boldsymbol q}^{{\prime \prime }}, \,\check {\!\boldsymbol q}^{{\prime \prime }} \right)_E = \int _{0}^{2h} { \left ( \bar {\rho } \,\check {\!\boldsymbol u}^{{{\prime \prime }}{{H}}} \,\check {\!\boldsymbol u}^{{\prime \prime }} + \frac {R \widetilde {T}}{\bar {\rho }} \check {\rho }^{{\prime }\dagger } \check {\rho }^{\prime } + \frac {\bar {\rho } c_v}{\widetilde {T}} \check {T}^{{{\prime \prime }}\dagger } \widehat {T}^{{\prime \prime }} \right ) {\mathrm{d}} y} = \,\check {\!\boldsymbol Q}^{{{\prime \prime }}{{H}}} \boldsymbol{M} \,\check {\!\boldsymbol Q}^{{\prime \prime }}. \end{equation}

Here, $\dagger$ denotes complex conjugate, the global vector $\,\check {\!\boldsymbol Q}^{{\prime \prime }}$ contains $\,\check {\!\boldsymbol q}^{{\prime \prime }}$ at all wall-normal grids and $\boldsymbol{M}$ is the weight matrix. Regarding the numerics, (2.14) is discretised using the Chebyshev collocation point method; 241 points are used in default, which is abundant to ensure grid independence. The solver verification can be found in Chen et al. (Reference Chen, Cheng, Fu and Gan2023a ).

Instead of the ideal modelling in (2.15), the forcings can be realistically computed by their definitions from the DNS data. Equation (2.13) indicates that the correlation tensor satisfies the following two Lyapunov equations for the LNS and eLNS models, respectively (derivations in Appendix A):

(2.19a) $${{\boldsymbol 0} \approx \dfrac {\partial \widehat {\boldsymbol{\Phi }}}{\partial t} = \widehat {\mathcal{L}}_q \widehat {\boldsymbol{\Phi }} + \widehat {\boldsymbol{\Phi }} \widehat {\mathcal{L}}_q^{{H}} + \underbrace {\left\langle \,\hat {\!\boldsymbol f}_q^{{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} + \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \,\hat {\!\boldsymbol f}_q^{{{\prime \prime }}{{H}}} \right\rangle }_{\equiv (\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}},\qquad \text{for the LNS model},}$$

(2.19b) $${{\boldsymbol 0} \approx \dfrac {\partial \widehat {\boldsymbol{\Phi }}}{\partial t} = \widehat {\mathcal{E}}_q \widehat {\boldsymbol{\Phi }} + \widehat {\boldsymbol{\Phi }} \widehat {\mathcal{E}}_q^{{H}} + \underbrace {\left\langle \,\hat {\!\boldsymbol e}_q^{{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} + \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \,\hat {\!\boldsymbol e}_q^{{{\prime \prime }}{{H}}} \right\rangle }_{\equiv (\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}},\qquad \text{for the eLNS model}.}$$

Here, the approximation on the left-hand side holds subject to a sufficiently large $N_t$ . Comparing between (2.14), (2.19), we can define the real forcing matrix from DNS as $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ and $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ , as underbraced in (2.19). A reasonable consequence is that when $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ and $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ are input into the stochastic linear models (2.14), the output $\widehat {\boldsymbol{\Phi }}$ should be identical to the DNS counterpart $\widehat {\boldsymbol{\Phi }}_{\textit{DNS}}$ , which is required by the models’ mathematical self-consistency. On the other hand, $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ and $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ naturally provide the ‘standard answers’ to how to model $\boldsymbol{{FF}}^{{H}}$ and $\boldsymbol{{EE}}^{{H}}$ , and can directly tell why the current LNS and eLNS models succeed or fail in specific cases.

It is worth noting that, when deriving (2.19), we do not and need not assume the DNS signal to be a white noise in the temporal domain; it is indeed not a white noise. The consequence is that $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ (and also $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ ) may be negative definite for some scales, so the matrix $(\boldsymbol{F})_{\textit{DNS}}$ may not really exist to satisfy this product. That is also the case when the white noise assumption and (2.15) fail. Nonetheless, this consequence does not affect the proof of the models’ self-consistency because we simply treat $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ as a whole and do not pursue further matrix decompositions.

4.1. Distributions of the nonlinear terms

Given the significance of accurately computing the forcing (§ 3), we scrutinise each nonlinear fluctuation term defined in (2.3). Their root mean squares in the LNS and eLNS models are computed from the DNS. The wall-normal distributions of these components are plotted in figures 2 and 3 for cases Ma15Re3k and Ma30Re5k, respectively. The LNS model is discussed first. For $f_{u,\textit {rms}}^{{\prime \prime }}$ and $f_{T,\textit {rms}}^{{\prime \prime }}$ , RSF and THF are respectively the dominant terms in nearly all regions. The terms related to molecular viscosity and conductivity (VSNF, VDNF and MHNF) are negligible except in the viscous sublayer, and the pressure-related terms PGNF and PCNF remain small throughout. For $f_{v,\textit {rms}}^{{\prime \prime }}$ and $f_{w,\textit {rms}}^{{\prime \prime }}$ , RSF still dominates away from the wall, but within and below the buffer layer, PGNF is significant and leads to an additional peak of $f_{\textit {rms}}^{{\prime \prime }}$ , especially for case Ma30Re5k (figures 3 b and 3 c). Therefore, the compressibility effect is prominent near the wall for the nonlinear terms in the wall-normal and spanwise momentum equations. The physics behind this additional peak will be clarified later using its spectra. Moreover, the convection-related terms MMNF and EMNF remain moderate throughout the field. They become increasingly important as ${\textit {Ma}}_b$ rises, suggesting notable convection of momentum and energy caused by $\rho ^{\prime }$ . The results of case Ma15Re9k are qualitatively the same as case Ma15Re3k, and are hence not displayed.

Next, we turn to the eLNS model. The original intention of introducing $\tau _R^{\prime }$ and $\vartheta _R^{\prime }$ is to partially model RSF and THF, so the residual terms eRSF and eTHF should be smaller than RSF and THF (see (2.8)). In figures 2 and 3, however, the results are the opposite. For all the equations shown, eRSF and eTHF are close to RSF and THF at $y^*\lesssim 15$ , where $\mu _t/\,\bar {\!\mu }$ (and $\kappa _t/\bar {\kappa }$ ) is small. In the outer region with $\mu _t\gg \,\bar {\!\mu }$ , however, eRSF and eTHF can be over two times larger than RSF and THF, so $e_{\textit {rms}}^{{\prime \prime }}$ are much higher than $f_{\textit {rms}}^{{\prime \prime }}$ in all the momentum and enthalpy equations. The $\textit{Re}_\tau$ in the present cases are relatively low, so the ratio $e_{\textit {rms}}^{{\prime \prime }}/f_{\textit {rms}}^{{\prime \prime }}$ will become even larger under higher $\textit{Re}_\tau$ , due to the increased $\mu _t/\,\bar {\!\mu }$ .

Figure 4. Variance of the three terms in (2.8) and the correlation between RSF and the modelled stress flux (all normalised by $u_\tau ^2$ ) for case Ma30Re5k: (a) streamwise and (b) wall-normal momentum equations.

Considering the unusual trend of eRSF and eTHF, the modelling assumption (2.8) deserves more discussion. The variances of the three components in (2.8), i.e. ${\textit{RSF}}_i$ , ${\partial \tau _{R,{i\!j}}^{\prime }}/{\partial x_j}$ and ${\textit{eRSF}}_i$ , are plotted in figure 4 for case Ma30Re5k, along with the correlation between ${\textit{RSF}}_i$ and ${\partial \tau _{R,{i\!j}}^{\prime }}/{\partial x_j}$ . For both the streamwise and wall-normal momentum equations (and also the spanwise one not shown), the modelled stress flux ${\partial \tau _{R,{i\!j}}^{\prime }}/{\partial x_j}$ does not follow ${\textit{RSF}}_i$ throughout, regarding both shapes and amplitudes. The modelled flux highly overestimates, and has moderately negative correction with RSF in the outer layer, which leads to the much higher amplitudes of eRSF and RSF. The same conclusion applies to THF and eTHF. Also, we have examined that the trends in figure 4 are analogous for other cases, insensitive to $\textit{Re}_\tau$ and ${\textit {Ma}}_b$ .

Consequently, we arrive at a conclusion that, the $\tau _R^{\prime }$ and $\vartheta _R^{\prime }$ fluxes are not partial modellings of RSF and THF, as previously believed (Hwang & Cossu Reference Hwang and Cossu2010; Pickering et al. Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021; Chen et al. Reference Chen, Cheng, Fu and Gan2023a ; Ying et al. Reference Ying, Liang, Li and Fu2023). On the contrary, the ‘residual’ terms eRSF and eTHF can be several times larger than RSF and THF, so $e_{\textit {rms}}^{{\prime \prime }}$ is more dominant by eRSF and eTHF, than the dominance of RSF and THF over $f_{\textit {rms}}^{{\prime \prime }}$ . As a result, the effects of other unmodelled nonlinear fluctuations (VSNF, MMNF, etc.) become minor, which increases the robustness of the shapes and amplitudes of the forcing. More discussions on this point will be provided in § 4.2. This robustness of the eLNS forcing also explains the trend in figure 1 that its ensemble average can reach convergence using a lower $N_t$ than the LNS counterpart.

Figure 5. Pre-multiplied one-dimensional spectra of the nonlinear fluctuation terms in the (a,b) streamwise and (c,d) spanwise directions in semi-local units, for the (a,c) LNS model ( $\,\hat {\!\boldsymbol f}_q^{{\prime \prime }}$ ) and (b, d) eLNS model ( $\,\hat {\!\boldsymbol e}_q^{{\prime \prime }}$ ), respectively, for case Ma30Re5k. The contours in each panel are normalised by their extreme values labelled on the top (in wall units $\rho _\tau$ , $u_\tau$ and $T_\tau$ ). The blue dotted lines denote the peak location of the $u$ -spectrum, and the blue dashed lines denote the $u$ -contour of 0.4.

The spectra of these nonlinear fluctuations are calculated, to identify the active wavenumbers of different components. Case Ma30Re5k is studied first, as it features the strongest density and temperature fluctuations. The pre-multiplied one-dimensional (1-D) spectra of $\,\hat {\!\boldsymbol f}_q^{{\prime \prime }}$ and $\,\hat {\!\boldsymbol e}_q^{{\prime \prime }}$ in the streamwise and spanwise directions are displayed in figure 5. The semi-local coordinates $y^*$ , $\lambda _x^*$ and $\lambda _z^*$ are adopted, and the spectra of the $u$ component are used as references in each row to examine the spectral similarity among different variables. A notable observation is that, compared with the LNS model, the spectra of $\,\hat {\!\boldsymbol e}_u^{{\prime \prime }}$ , $\,\hat {\!\boldsymbol e}_v^{{\prime \prime }}$ , $\,\hat {\!\boldsymbol e}_w^{{\prime \prime }}$ and $\,\hat {\!\boldsymbol e}_T^{{\prime \prime }}$ in eLNS have greater resemblance to each other in shapes and peak locations for both the streamwise and spanwise ones; see figures 5(b) and 5(d). This similarity among variables can be explained by the dominance of eRSF and eTHF over other terms, as observed in figure 3. A further decomposition of $\partial \tau _{R,{i\!j}}^{\prime }/\partial x_j$ and $\partial \vartheta _{R,j}^{\prime }/\partial x_j$ (difference between ${\boldsymbol f}_q^{{\prime \prime }}$ and ${\boldsymbol e}_q^{{\prime \prime }}$ ) shows that the dominant terms for the four terms are $\mu _t \Delta \hat {u}^{{\prime \prime }}$ , $\mu _t \Delta \hat {v}^{{\prime \prime }}$ , $\mu _t \Delta \hat {w}^{{\prime \prime }}$ and $\kappa _t \Delta \hat {T}^{{\prime \prime }}$ , respectively, with $\Delta$ the Laplace operator. Therefore, $\,\hat {\!\boldsymbol e}_q^{{\prime \prime }}$ resides in smaller scales than $\hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }}$ and $\,\hat {\!\boldsymbol f}_q^{{\prime \prime }}$ , and the similarity among $\,\hat {\!\boldsymbol e}_u^{{\prime \prime }}$ , $\,\hat {\!\boldsymbol e}_v^{{\prime \prime }}$ , $\,\hat {\!\boldsymbol e}_w^{{\prime \prime }}$ and $\,\hat {\!\boldsymbol e}_T^{{\prime \prime }}$ is more easily established than $\,\hat {\!\boldsymbol f}^{{\prime \prime }}$ .

An additional peak appears in the $\hat {f}_v^{{\prime \prime }}$ and $\hat {e}_v^{{\prime \prime }}$ spectra in the viscous sublayer, which is contributed by PGNF (equal to $-\partial (\rho ^{\prime } T^{{\prime \prime }})/\partial y$ ) from figure 3. Compared with the outer peak, this inner peak resides in a larger scale of $\lambda _z^*\approx 180$ . The value of $\lambda _x^*$ of this inner peak exceeds 1500, but due to the relatively low $\textit{Re}_\tau$ in this case, its streamwise spectrum is not fully resolved. The $\lambda _x^*\ll \lambda _z^*$ feature of this inner peak indicates that it is not caused by the near-wall acoustic motions (where $\lambda _x\lt \lambda _z$ ; see § 5), but is attributed to the near-wall density and temperature streaks which are most active at $\lambda _x^+\sim 1000$ and $\lambda _z^+\sim 100$ ; these two streaks are also in nearly perfect negative correlation (Coleman, Kim & Moser Reference Coleman, Kim and Moser1995; Patel, Boersma & Pecnik Reference Patel, Boersma and Pecnik2016; Cheng & Fu Reference Cheng and Fu2024).

Figure 6. Same as figure 5, but for (a,b) case Ma15Re3k and (c,d) case Ma15Re9k and only the results of the eLNS model are shown. Panels (a,c) are the streamwise spectra and panels (b,d) are the spanwise ones.

To examine the robustness of the spectral features of the eLNS nonlinear terms at different ${\textit {Ma}}_b$ and $\textit{Re}_\tau$ , the 1-D spectra of $\,\hat {\!\boldsymbol e}_q^{{\prime \prime }}$ for cases Ma15Re3k and Ma15Re9k are depicted in figure 6. As in case Ma30Re5k, similarity between the spectra of $\,\hat {\!\boldsymbol e}_u^{{\prime \prime }}$ , $\,\hat {\!\boldsymbol e}_v^{{\prime \prime }}$ , $\,\hat {\!\boldsymbol e}_w^{{\prime \prime }}$ and $\,\hat {\!\boldsymbol e}_T^{{\prime \prime }}$ is observed in both the streamwise and spanwise directions. Meanwhile, the inner peak of $\,\hat {\!\boldsymbol e}_v^{{\prime \prime }}$ due to the PGNF is less evident than in figure 5 owing to the reduced ${\textit {Ma}}_b$ . Comparing between figures 6 and 5, we note that the (outer) peak location of these 1-D spectra is not sensitive to ${\textit {Ma}}_b$ and $\textit{Re}_\tau$ after using the semi-local units, which is measured to be $\lambda _x^*\approx 90$ , $\lambda _x^*\approx 30$ and $y^*\approx 40$ .

In short, we emphasise that the introduction of $\tau _R^{\prime }$ and $\vartheta _R^{\prime }$ is not a partial modelling of RSF and THF. On the contrary, it highly lifts the amplitudes of eRSF and eTHF, and hence the dominance of $\,\hat {\!\boldsymbol e}_q^{{\prime \prime }}$ by eRSF and eTHF. Therefore, $\,\hat {\!\boldsymbol e}_q^{{\prime \prime }}$ is more robust than $\,\hat {\!\boldsymbol f}_q^{{\prime \prime }}$ for different variables in terms of the peaks and shapes of their physical and spectral distributions. Consequently, $\,\hat {\!\boldsymbol e}_q^{{\prime \prime }}$ can be more easily modelled than $\,\hat {\!\boldsymbol f}_q^{{\prime \prime }}$ by featuring more robust composition and structures, although this robustness may not always improve the model results; see § 4.2.

4.2. Distributions of the forcing matrix

As mentioned in § 3, $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ and $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ provide standard answers to how to model $(\boldsymbol{{FF}}^{{H}})$ and $(\boldsymbol{{EE}}^{{H}})$ in the stochastic linear models, so the distributions of these forcing are discussed at different $k_x$ and $k_z$ . The spectral fluctuation is usually classified into two branches according to its aspect ratio $\lambda _x/\lambda _z$ . The structure with $\lambda _x\gt \lambda _z$ is considered as streamwise elongated, which is energetic and has strong linear coherence in the wall-normal direction (e.g. del Álamo & Jiménez Reference del Álamo and Jiménez2006). In comparison, the structure with $\lambda _x\lt \lambda _z$ is regarded as spanwise elongated, which receives energy from the nonlinear terms and can be worse predicted by the linear models than the streamwise-elongated one (Gupta et al. Reference Gupta, Madhusudanan, Wan, Illingworth and Juniper2021; Symon et al. Reference Symon, Illingworth and Marusic2021). These two types of fluctuations are discussed separately below.

Figure 7. Contours of the forcing matrix (a) $|(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}|$ and (b) $|(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}|$ , at $(k_x,k_z)h=(0.5,15)$ for case Ma30Re5k. Their diagonal terms, as labelled in blue dotted lines, are plotted in panels (c,d) in normalised values; the purple reference lines are from (2.15).

The streamwise-elongated fluctuations are considered first. Figure 7 vividly demon- strates what the real forcing matrices $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ and $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ look like for case Ma30Re5k at $(k_x,k_z)h=(0.5,15)$ ( $\lambda _x=12.6h$ , $\lambda _z=0.4h$ ). A direct observation is that these two matrices are far more complicated than the simple diagonal modelling in (2.15). Recall that two levels of simplification are made in (2.15). First, only the self-correlation matrices of the forcing variables ( $\hat {f}_\rho ^{{\prime \prime }}\hat {f}_\rho ^{{{\prime \prime }}{{H}}}$ , $\hat {f}_u^{{\prime \prime }}\hat {f}_u^{{{\prime \prime }}{{H}}}$ , etc.) are retained, while all the inter-variable correlations ( $\hat {f}_\rho ^{{\prime \prime }}\hat {f}_u^{{{\prime \prime }}{{H}}}$ , etc.) are neglected. Second, each self-correlation matrix is set to be diagonal, signifying perfect zero correlation between different $y$ . In figure 7, however, there are significant off-diagonal terms for both $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ and $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ , suggesting non-negligible inter-variable and inter-height correlations. Regarding the difference between the LNS and eLNS forcing, a visual observation seems to suggest that $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ is more diagonal than $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ , in terms of wall-normal correlations. This is reasonable because the dominant terms $\tau _R^{\prime }$ and $\vartheta _R^{\prime }$ in eLNS are linear functions of the fluctuations and are related to local mean-flow gradients, so they are more localised than the other nonlinear terms which contain high-order fluctuation moments. More quantitative evidence will be provided later. Consequently, the diagonal modelling in (2.15) benefits more the forcing in the eLNS model, than the LNS one. The diagonal elements of $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ and $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ are extracted in figures 7(c) and 7(d). The modelled distributions in (2.15) are also displayed for reference. As discussed in § 4.1, the diagonal terms for eLNS reflect the large contributions of eRSF and eTHF, so the diagonal terms of the $u$ , $v$ , $w$ and $T$ components exhibit stronger resemblance to each other, hence they are more robust for modelling compared with the LNS model. Nevertheless, the prediction of (2.15) is not very consistent with the DNS results; the real forcing is more concentrated near the wall under this small $\lambda _z$ .

Figure 8. (a,c) Energy ratios occupied by the leading 10 singular values of the forcing matrix in figure 7, and (b,d) shape functions of the principal forcing mode (case Ma30Re5k, $(k_x,k_z)h=(0.5,15)$ ). Panels (a,b) are for the LNS model, and (c,d) for the eLNS one.

From the modelling perspective, it is interesting to show whether the DNS forcing is of low rank. The diagonal modelling in (2.15) means that the forcing is assumed to be of full rank, while Morra et al. (Reference Morra, Nogueira, Cavalieri and Henningson2021) demonstrate using the incompressible DNS data that the forcing has clear spatial-temporal coherence. The rank characteristics of the forcing matrix can be measured using the singular value decomposition (SVD), as

(4.1) \begin{equation} (\boldsymbol{{FF}}^{{H}})_{\textit{DNS}} = \sum _j {\sigma _j \,\check {\!\boldsymbol f}_j^{{\prime \prime }} \,\check {\!\boldsymbol f}_j^{{{\prime \prime }}{{H}}}} = \sum _j { \sigma _j \,\check {\!\boldsymbol{\Psi }}_j }, \end{equation}

where $\sigma _j$ is the singular value sorted in the descending order; $\,\check {\!\boldsymbol f}_j^{{\prime \prime }}$ is the singular vector satisfying the orthogonal condition $(\,\check {\!\boldsymbol f}_i^{{\prime \prime }},\,\check {\!\boldsymbol f}_j^{{\prime \prime }})_E=\delta _{{i\!j}}$ , so each component $\,\check {\!\boldsymbol{\Psi }}_j$ is of rank one. The left and right singular vectors are the same because the forcing matrix is Hermitian. Besides, the sum of the singular values, $\sum _j{\sigma _j}\equiv V_\sigma$ , represents the total energy input by the forcing. The decomposition for $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ is defined likewise. Based on (4.1), we would regard a forcing matrix as of low rank if the first or first few SVD modes occupy a large portion of the total forcing energy.

The leading ten $\sigma _{F,j}$ and $\sigma _{E,j}$ of the matrices in figure 7 are plotted in figures 8(a) and 8(c), expressed in relative ratios $r_{\sigma ,j}=\sigma _j/V_\sigma$ . Note that $\sigma _j$ appears in pairs due to the flow symmetry or antisymmetry on the two sides of the channel regarding the centreline, so every two $\sigma _j$ are combined to reflect their joint contributions. The forcing $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ for the LNS model is of relatively low rank, with $r_{\sigma F,1}$ over 57 %. This low-rank feature is very different from the full-rank $\boldsymbol{{FF}}^{{H}}$ (fully diagonal) assumed in the linear model. The shape of the principal mode $\,\check {\!\boldsymbol f}_1^{{\prime \prime }}$ is plotted in figure 8(b). The $\check {f}_{u,1}^{{\prime \prime }}$ and $\check {f}_{T,1}^{{\prime \prime }}$ components turn out to be the two largest, which is different from a previous suggestion of the linear models that, for the principal forcing mode, the $\check {f}_{v,1}^{{\prime \prime }}$ and $\check {f}_{w,1}^{{\prime \prime }}$ components are dominant for streamwise-elongated modes in the form of streamwise vortices to induce streaks (Hwang & Cossu Reference Hwang and Cossu2010; Chen et al. Reference Chen, Cheng, Fu and Gan2023a ). In comparison, $r_{\sigma E,1}={26}\,\%$ for the eLNS model is lower than $r_{\sigma F,1}$ . In other words, the modelled terms eRSF and eTHF lessen the spatial coherence of the forcing, and add to its diagonal dominance. Meanwhile, eRSF and eTHF markedly enhance the outer-layer forcing (figure 3), so the shapes of the principal mode in figure 8(d) is also lifted away from the wall, compared with its LNS counterpart. This amplified outer-layer forcing energy can be rapidly dissipated within the linear operator, where the stronger damping in eLNS is due to the replacement of $\,\bar {\!\mu }$ with $\,\bar {\!\mu }+\mu _t$ in the linear operator (see § 2.2). Consequently, the response mode can still take the form of near-wall streaky motions (c.f. figure 1 b also with $\lambda _x\gt \lambda _z$ ).

Figure 9. Same as figure 7 except for scale $(k_x,k_z)h=(20,3)$ .

The forcing pattern for spanwise-elongated modes is demonstrated in figure 9, where the wavenumbers are selected to be $(k_x,k_z)h=(20,3)$ ( $\lambda _x=0.6h$ , $\lambda _z=4.2h$ ). Compared with figure 7, the matrices $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ and $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ become visually more diagonal in the wall-normal direction. Equivalently, the two forcings are both of higher rank, compared with those in figures 7, 8, which is inferred from their more uniform $\sigma _j$ (not shown; $r_{\sigma F,1},r_{\sigma E,1}\lt {12}\,\%$ ). This stronger diagonality for spanwise-elongated structures is reasonable, because they tend to be more localised and exhibit weak inter-layer linear coherence (e.g. Marusic, Baars & Hutchins Reference Marusic, Baars and Hutchins2017). The diagonal elements of the two forcing matrices are plotted in figures 9(c) and 9(d). It is interesting to show that the prediction of (2.15) is in good agreement with the DNS results, especially for the eLNS model. Nevertheless, both incompressible and compressible works suggest that the linear models perform worse for spanwise-elongated fluctuations (Gupta et al. Reference Gupta, Madhusudanan, Wan, Illingworth and Juniper2021; Symon et al. Reference Symon, Illingworth and Marusic2021; Chen et al. Reference Chen, Cheng, Gan and Fu2023b ). From the modelling aspect, a prominent disrupting factor is the unmodelled inter-variable correlations, particularly $\hat {f}_\rho ^{{\prime \prime }}\hat {f}_u^{{{\prime \prime }}{{H}}}$ , $\hat {f}_\rho ^{{\prime \prime }}\hat {f}_T^{{{\prime \prime }}{{H}}}$ , $\hat {f}_u^{{\prime \prime }}\hat {f}_v^{{{\prime \prime }}{{H}}}$ and $\hat {f}_u^{{\prime \prime }}\hat {f}_T^{{{\prime \prime }}{{H}}}$ (and also the $\hat {e}^{{\prime \prime }}$ counterparts). They are all assumed zero in (2.15). From the physical perspective, previous studies showed that the linear dynamics (lift-up mechanism, transient growth, etc.) is important for generating streamwise-elongated structures, but contrarily, the travelling-wave-like spanwise-elongated fluctuations are mainly caused by nonlinear effects related to turbulent bursting, vortex stretching and streak breakdown (Schoppa & Hussain Reference Schoppa and Hussain2002; Yu et al. Reference Yu, Zhou, Dong, Yuan and Xu2024).

Figure 10. Projection coefficients of the DNS forcing to the space spanned by the POD modes for case Ma30Re5k: (a,b,c,d) LNS part $\alpha _{{i\!j}}$ , (e,f,g,h) eddy-viscosity part $\xi _{{i\!j}}$ and (i,j,k,l) eLNS part $\beta _{{i\!j}}$ ; see (4.2). Panels show (a,e,i) $\lambda _x=6.3h$ , $\lambda _z=3.1h$ ; (b,f,j) $\lambda _x=6.3h$ , $\lambda _z^+=179$ ; (c,g,k) $\lambda _x^+=143$ , $\lambda _z=3.1h$ ; (d,h,l) $\lambda _x^+=143$ , $\lambda _z^+=179$ . The input energy $V_\sigma$ shown is amplified by 10⁵ for all panels for convenience.

The diagonality of the forcing matrices in figures 7 and 9 relies on visual measurement. For more quantitative estimation, we project the forcing into the space spanned by the POD modes of $\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }}\hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle$ . The projection coefficients for models LNS and eLNS are then computed as

(4.2) \begin{equation} \alpha _{{i\!j}} = \,\check {\!\boldsymbol q}_i^{{{\prime \prime }}{{H}}}\boldsymbol{M}\big(\boldsymbol{{FF}}^{{H}}\big)_{\textit{DNS}} \boldsymbol{M} \,\check {\!\boldsymbol q}_j^{{\prime \prime }}, \qquad \beta _{{i\!j}} = \,\check {\!\boldsymbol q}_i^{{{\prime \prime }}{{H}}}\boldsymbol{M}\big(\boldsymbol{{EE}}^{{H}}\big)_{\textit{DNS}} \boldsymbol{M} \,\check {\!\boldsymbol q}_j^{{\prime \prime }}. \end{equation}

This projection reduces the forcing in the physical space into discrete coefficients, and can jointly reflect inter-variable and inter-height correlations. If the forcing is ideally a white noise modelled by (2.15), then $\alpha _{{i\!j}}$ (and $\beta _{{i\!j}}$ ) is also a diagonally dominant matrix for nearly all spatial scales. The difference between $\alpha _{{i\!j}}$ and $\beta _{{i\!j}}$ is also of interest, denoted as $\xi _{{i\!j}}=\beta _{{i\!j}}-\alpha _{{i\!j}}$ , which reflects the contribution of the eddy-viscosity-related terms.

The Ma30Re5k case is discussed first using four representative scales. The first two are streamwise- and spanwise-elongated fluctuations discussed above; the third type is a large-scale fluctuation in the outer region with $(\lambda _x,\lambda _z)/h=(6.3,3.1)$ , and the fourth is a small-scale fluctuation with $\lambda _x^+,\lambda _z^+\sim 100$ . Their $\alpha _{{i\!j}}$ , $\xi _{{i\!j}}$ and $\beta _{{i\!j}}$ values are depicted in figure 10 for the leading indices $i,j\lt 10$ . These coefficients are normalised by their input energy $V_{\sigma }$ , as labelled, such that their trace is unity: $\Sigma _i \alpha _{ii}/V_{\sigma ,\alpha }=1$ (so are $\xi _{{i\!j}}$ and $\beta _{{i\!j}}$ ). For the large-scale fluctuation in figure 10(a), $\alpha _{{i\!j}}$ for the LNS model has pronounced off-diagonal elements, indicating strong inter-variable and/or inter-height correlations. In comparison, the eddy-viscosity term $\xi _{{i\!j}}$ has much higher energy ( $V_\sigma$ ) than $\alpha _{{i\!j}}$ , as discussed for figure 3, and is highly diagonally dominant. As a result, their summation $\beta _{{i\!j}}=\alpha _{{i\!j}}+\xi _{{i\!j}}$ for the eLNS model has a similar distribution to $\xi _{{i\!j}}$ and features high diagonal dominance. The above trend is analogous to streamwise-elongated structures in the second column of figure 10, which has been visually shown in figure 9. In comparison, for the structures of small $\lambda _x$ (irrespective of $\lambda _z$ ), $\alpha _{{i\!j}}$ , $\xi _{{i\!j}}$ and $\beta _{{i\!j}}$ are all diagonally dominant, as shown in the last two columns of figure 10, which is in line with the observations in figure 9.

Figure 11. Same as figure 10 but for (a,b,e,f) case Ma15Re3k and (c,d,g,h) case Ma15Re9k. Panels (a,c,e,g) are fluctuations of large $\lambda _x$ , $\lambda _z$ and panels (b,d,f,h) are of small $\lambda _x$ , $\lambda _z$ .

The projection results for cases Ma15Re3k and Ma15Re9k are shown in figure 11. The distributions of $\alpha _{{i\!j}}$ and $\beta _{{i\!j}}$ for different scales resemble a lot those in figure 10, suggesting weak dependence on ${\textit {Ma}}_b$ and $\textit{Re}_\tau$ . Therefore, the role of the eddy-viscosity terms is clarified again (as in § 4.1) that it adds a large-amplitude diagonally dominant term to the LNS forcing matrix (specifically, the self-correlation part of the forcing variables), so that the eLNS forcing is more diagonally dominant and has similar shapes at different scales. From the physical point of view, this enhanced diagonal dominance was shown to build a better balance between local turbulent energy dissipation and transport (Symon et al. Reference Symon, Madhusudanan, Illingworth and Marusic2023; Ying et al. Reference Ying, Chen, Li and Fu2024a ). Consequently, the eLNS model is more robust among different scales in predicting variable correlations than LNS.

Based on the above conclusion, a natural thought to improve the model arises that the robustness of the eLNS model can be further strengthened by artificially adding very large diagonal elements, for example, 10 times ${\rm{eRSF}}$ , although unphysical, to the forcing matrix $(\boldsymbol{{EE}}^{{H}})$ . Nevertheless, we note that the model behaviours will not be necessarily improved for two reasons. First, although the inter-height correlation is more robustly modelled, the off-diagonal inter-variable correlation for the forcing also matters (see figure 9), so the relative amplitudes of different variables still need to be carefully determined. Second, artificial enhancement of the diagonal terms can destroy the non-normality of the NS operator, as already reflected in the eLNS model to some extent (Symon et al. Reference Symon, Illingworth and Marusic2021; Chen et al. Reference Chen, Cheng, Gan and Fu2023b ). As a result, not only can the linear amplification mechanism be disrupted, but also the energies of different branches of modes can be disorderly sorted, possibly leading to excessive amplification of the acoustic and other components for compressible flows.

Finally, we provide a global view of the coherence of the LNS and eLNS forcing at all scales. The contours of $r_{\sigma ,1}$ , which measures the energy ratio of the most energetic forcing mode to all modes, are plotted in figure 12 for all three cases. For the LNS model, the forcing for streamwise-elongated fluctuations ( $\lambda _x\gt \lambda _z$ ) has relatively high coherence and low rank, compared with the spanwise-elongated parts. The maximum $r_{\sigma ,1}$ for case Ma15Re3k reaches over 70 % at $(\lambda _x,\lambda _z)/h\approx (12.6,0.7)$ . The maxima in other two cases (figures 12 b and 12 c) also appear at approximately the same scale, suggesting a robust mechanism of scale selection for the forcing at different ${\textit {Ma}}_b$ and $\textit{Re}_\tau$ . As discussed above, this low-rank feature of the LNS forcing at $\lambda _x\gt \lambda _z$ explains its poor behaviours in predicting these energetic fluctuations (e.g. Madhusudanan et al. Reference Madhusudanan, Illingworth and Marusic2019). Meanwhile, we notice that the distribution of $r_{\sigma ,1}$ resembles a lot the linear coherence spectra of $u$ and $T$ between different heights (Marusic et al. Reference Marusic, Baars and Hutchins2017; Gupta et al. Reference Gupta, Madhusudanan, Wan, Illingworth and Juniper2021; Chen et al. Reference Chen, Cheng, Gan and Fu2023b ) within the inner–outer interaction model. Therefore, the LNS forcing for the scales with strong inter-height linear coherence of velocity and temperature is of particularly low rank.

Figure 12. Energy ratio occupied by the leading POD mode ( $r_{\sigma ,1}=\sigma _1/\sum _j\sigma _j$ ) for the forcing in the (a,b,c) LNS and (d,e,f) eLNS models at different $\lambda _x$ , $\lambda _z$ for cases (a,d) Ma15Re3k, (b,e) Ma30Re5k and (c,f) Ma15Re9k. The black dashed line denotes $\lambda _x=\lambda _z$ .

For a more quantitative comparison of $r_{\sigma ,1}$ between different cases, its variation with $\lambda _z/h$ is plotted in figure 13 for all three compressible cases, at the largest streamwise scale $\lambda _x=4\pi h$ . The two incompressible cases are also included for reference. Cases Ma00Re3k, Ma15Re3k and Ma30Re5k have comparable $\textit{Re}_\tau ^*=141\!\sim \!186$ , and their $r_{\sigma ,1}$ values are also close to each other, except at the very small scale. For the two higher $\textit{Re}_\tau ^*$ cases Ma15Re9k and Ma00Re10k, their $r_{\sigma ,1}$ also have analogous distributions. Therefore, it is concluded that the low-rank feature of the LNS forcing is not sensitive to ${\textit {Ma}}_b$ , although a higher ${\textit {Ma}}_b$ can slightly decrease $r_{\sigma ,1}$ due to the increasing acoustic components (see § 5). Meanwhile, the low-rank feature is weakened by a rising $\textit{Re}_\tau$ , suggesting reduced spatial-temporal coherence as turbulence intensifies.

Figure 13. Energy ratio of the leading POD mode ( $r_{\sigma ,1}$ ) for the LNS forcing with different $\lambda _z$ for different cases (at largest $\lambda _x=4\pi h$ ). Cases Ma00Re3k, Ma15Re3k and Ma30Re5k have comparable $\textit{Re}_\tau ^*=141\!\sim \!186$ ; cases Ma15Re9k and Ma00Re10k have higher $\textit{Re}_\tau ^*=393$ and 547 (see table 2).

Regarding the eLNS results, figure 12 shows that the eLNS model largely disrupts the low-rank feature of the LNS forcing for scales $\lambda _x\gt \lambda _z$ , by including the eddy-viscosity-enhanced terms. The large-scale fluctuations of high linear coherence are particularly affected, with $r_{\sigma ,1}$ down to $\lesssim {30}\,\%$ , so these scales are the region where the prediction of the linear coherence is mostly improved by the eLNS model over LNS. For other scales ( $\lambda _x\lt \lambda _z$ , or small $\lambda _z$ ), the improvement by eLNS is not obvious. Further improvements are anticipated by optimising the $\mu _t$ profile for different scales, as realised by Symon et al. (Reference Symon, Madhusudanan, Illingworth and Marusic2023) and Ying et al. (Reference Ying, Chen, Li and Fu2024a ) for incompressible flows, and by considering inter-variable correlations.

5.1. Pressure decomposition for acoustic fluctuations

One viable path to extract the acoustic fluctuations is to first decompose the pressure fluctuation $p^\prime$ based on the pressure Poisson equation, and then formulate other acoustic variables based on the compressible part of $p^\prime$ . The first step is conducted in this subsection.

Following Sarkar (Reference Sarkar1992) and Zhang et al. (Reference Zhang, Wan, Sun and Lu2024), $p^{\prime }$ is decomposed into the quasi-incompressible part $p_{ic}^{\prime }$ and the compressible part $p_c^{\prime }$ . The governing equations for the two are obtained by classifying the right-hand-side terms of the pressure Poisson equation, as

(5.1a) $${\frac {\partial ^2 p_{ic}^{\prime }}{\partial x_i \partial x_i} = - 2\frac {\partial \tilde {u}_i}{\partial x_j} \frac {\partial \rho u_j^{{\prime \prime }}}{\partial x_i} - \frac {\partial ^2}{\partial x_i \partial x_j} \left ( \rho u_i^{{\prime \prime }} u_j^{{\prime \prime }} - \bar {\rho }\widetilde {u_i^{{\prime \prime }} u_j^{{\prime \prime }}} \right ) + \frac {\partial ^2 \tau _{ij}^{\prime }}{\partial x_i \partial x_j},}$$

(5.1b) $${\frac {\partial p_c^{\prime }}{\partial x_i \partial x_i} = \frac {\partial ^2 \rho ^{\prime }}{\partial t^2} - \frac {\partial ^2 (\rho ^{\prime } \tilde {u}_i \tilde {u}_j )}{\partial x_i \partial x_j} - \tilde {u}_i \frac {\partial ^2 (2\rho u_j^{{\prime \prime }})}{\partial x_i \partial x_j} - \frac {\partial }{\partial x_i} \left (2 \rho u_i^{{\prime \prime }} \frac {\partial \tilde {u}_j}{\partial x_j}\right ).}$$

The boundary conditions are $\partial p_{ic}^{\prime }/\partial y=\partial \tau _{i2}^{\prime }/\partial x_i$ and $\partial p_c^{\prime }/\partial y=0$ at both wall sides (Yu, Xu & Pirozzoli Reference Yu, Xu and Pirozzoli2020). For incompressible flows, $p_c^{\prime }$ naturally reduces to zero since $\rho ^{\prime }=\partial \tilde {u}_j/\partial x_j=\partial {u}_j^{{\prime \prime }}/\partial x_j=0$ . The part $p_{ic}^{\prime }$ can be further decomposed (into a rapid term, etc.), but here, we simply treat it as a whole, as we focus on the compressibility effects. In this sense, $p_{ic}^{\prime }$ reflects the collective contribution of vortical modes, while $p_c^{\prime }$ is mainly contributed by acoustic modes. Equation (5.1) is solved using the Fourier–Galerkin scheme. Since the temporal derivative in (5.1b ) is not available in our time-not-resolved data, we obtain $p_c^{\prime }$ simply from $p_c^{\prime }=p^{\prime }-p_{ic}^{\prime }$ after solving $p_{ic}^{\prime }$ from (5.1a ). The accuracy of $p_c^{\prime }$ is reliable because the realistic $p^{\prime }$ and that solved from the Poisson equation were shown to differ by less than 1 % for similar channel cases (Yu et al. Reference Yu, Xu and Pirozzoli2020).

Figure 14. (a) Variance and correlation of pressure fluctuation components for cases Ma15Re3k and Ma30Re5k. (b $-$ e) Pre-multiplied 2-D spectra for the (b,d) incompressible and (c,e) compressible parts of the wall pressure fluctuations for cases (b,c) Ma15Re3k and (d,e) Ma30Re5k. The black dashed line in the right four panels denotes $\lambda _x=\lambda _z$ .

The variances of the pressure and its components for cases Ma15Re3k and Ma30Re5k are plotted in figure 14(a), where the normalisation is $p^{{\prime }+}=p^{\prime }/\tau _w$ . After the decomposition, the $p_{ic,\textit {rms}}^{{\prime }+}$ in the two cases are close to each other under the same $\textit{Re}_\tau ^*$ , suggesting robust contributions from the quasi-incompressible part. If $p_{ic}^{{\prime }+}$ is regarded as the result at ${\textit {Ma}}_b=0$ , then the wall pressure variance roughly follows a linear relation

(5.2) \begin{equation} \overline {p_w^{{\prime }+2}}({\textit {Ma}}_b) \approx \overline {p_w^{{\prime }+2}}(0) + 0.216{\textit {Ma}}_b^2, \quad \text{for}\; \textit{Re}_\tau ^*\approx 140. \end{equation}

This scaling is consistent with Yu et al. (Reference Yu, Xu and Pirozzoli2020) although the slope differs due to a different $\textit{Re}$ . With ${\textit {Ma}}_b$ increased, $p_{c,\textit {rms}}^{{\prime }+}$ quickly rises and surpasses $p^{{\prime }+}_{ic,\textit {rms}}$ at ${\textit {Ma}}_b=3$ at the wall. Moreover, $p_{c,\textit {rms}}^{{\prime }+}$ is mostly amplified near the wall, suggesting locally intensified fluctuations of mass flux. The correlation $\overline {p_{ic}^{{\prime }+}p_c^{{\prime }+}}$ is also plotted, which remains small ( ${\lt } 0.2$ ) throughout the field for both cases.

Furthermore, the pre-multiplied 2-D spectra of $\hat {p}_{ic}^{\prime }$ and $\hat {p}_{c}^{\prime }$ at the wall are displayed in figure 14. The pressure spectra in compressible flows have been comprehensively studied before (Bernardini & Pirozzoli Reference Bernardini and Pirozzoli2011; Duan, Choudhari & Zhang Reference Duan, Choudhari and Zhang2016). Our intention here is to provide further support to the decomposition (5.1), by demonstrating the resemblance of the $\hat {p}_{ic}^{\prime }$ spectra under different ${\textit {Ma}}_b$ and the clear scale separation between $\hat {p}_{ic}^{\prime }$ and $\hat {p}_{c}^{\prime }$ . The part $\hat {p}_{c}^{\prime }$ resides at much smaller streamwise scales than $\hat {p}_{ic}^{\prime }$ , taking the form of travelling waves, so $\hat {p}_{c}^{\prime }$ turns out to be spanwise elongated with the peak locating at $(\lambda _x,\lambda _z)/h\approx (0.5,1)$ for case Ma30Re5k.

5.2. Variance and structure of acoustic components

The $\hat {p}_c^{\prime }$ obtained in § 5.1 can be utilised to recover the other variables of acoustic modes. From 1-D eigenmode analysis, the $\hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }}$ of acoustic modes, denoted as $\hat {\boldsymbol q}_{\textit {ac}}^{{\prime \prime }}$ , is related to $\hat {p}_c^{\prime }$ as

(5.3) \begin{equation} \hat {\boldsymbol q}_{\textit {ac}}^{{\prime \prime }} = \big[\rho _{\textit {ac}}^{\prime },\,u_{\textit {ac}}^{{\prime \prime }},\,v_{\textit {ac}}^{{\prime \prime }},\,w_{\textit {ac}}^{{\prime \prime }},\,T_{\textit {ac}}^{{\prime \prime }}\big]^{T} = \frac {\boldsymbol{C}_{\textit {vis}}}{\gamma \bar {p}} \left [ \bar {\rho },\, \frac {k_xa}{k}, -\frac {{\mathrm{i}} a}{k}\frac {{{d}}}{{\mathrm{d}} y},\, \frac {k_za}{k},\, (\gamma -1)\widetilde {T} \right ]^{T} \hat {p}_c^{\prime }, \end{equation}

where $k^2=k_x^2+k_z^2$ . Derivation of (5.3) is detailed in Appendix C. Consequently, the acoustic part $\langle \hat {\boldsymbol q}_{\textit {ac}}^{{\prime \prime }} \hat {\boldsymbol q}_{\textit {ac}}^{{{\prime \prime }}{{H}}} \rangle \equiv \langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle _{\textit {ac}}$ can be extracted out of $\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle$ , based on $\langle \hat {p}_c^{\prime }\hat {p}_c^{{\prime }{{H}}} \rangle$ and (5.3). The residual $\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle _{\textit {nac}}=\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle -\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle _{\textit {ac}}$ is regarded as the non-acoustic part, which mostly comes from vortical modes although it also includes weak vortical–acoustic coupling. In the following, the relative contributions of $\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle _{\textit {ac}}$ and $\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle _{\textit {nac}}$ are assessed, and the two are further utilised to decompose the forcing, which will uncover the quantitative contribution of the acoustic components in the stochastic linear model.

Equation (5.3) explicitly indicates the scale dependence of $\hat {\boldsymbol q}_{\textit {ac}}^{{\prime \prime }}$ . The ratios $\hat {\rho }_{\textit {ac}}^{\prime }/\hat {p}_c^{\prime }$ and $\hat {T}_{\textit {ac}}^{{\prime \prime }}/\hat {p}_c^{\prime }$ turn out to be nearly independent of $k_x$ and $k_z$ , so the spectral distributions of $\hat {\rho }_{\textit {ac}}^{\prime }$ and $\hat {T}_{\textit {ac}}^{{\prime \prime }}$ are primarily determined by $\hat {p}_c^{\prime }$ and the mean flow. Regarding the velocity components, the relation $\hat {u}_{\textit {ac}}^{{\prime \prime }}\ll \hat {w}_{\textit {ac}}^{{\prime \prime }}$ holds if $k_x\ll k_z$ , so the acoustic motion for streamwise-elongated structures is mainly active in the $y$ – $z$ plane with velocities $\hat {v}_{\textit {ac}}^{{\prime \prime }}$ and $\hat {w}_{\textit {ac}}^{{\prime \prime }}$ . Conversely, the acoustic motion is constrained in the $x$ – $y$ plane for spanwise-elongated structures when $k_x\gg k_z$ .

The $\hat {p}_c^{\prime }$ value in the two ${\textit {Ma}}_b=1.5$ cases is relatively low, so we focus on case Ma30Re5k to see the potential effects of acoustic modes. Before inspecting $\hat {\boldsymbol q}_{\textit {ac}}^{{\prime \prime }}$ , we first examine the spectral distribution of the relative strength of $\hat {p}_c^{\prime }$ and $\hat {p}_{ic}^{\prime }$ . To reflect their amplitudes at different scales, the energy norm for pressure is defined as (Chu Reference Chu1965)

(5.4) \begin{equation} V_p = ( \hat {p}^{\prime }, \hat {p}^{\prime } )_E = \int _{0}^{2h} { \left [ \bar {\rho } a^2 \frac {\hat {p}^{{\prime }\dagger }\hat {p}^{{\prime }}}{(\gamma \bar {p})^2} \right ] {\mathrm{d}} y}. \end{equation}

Figure 15 provides the distributions of three norm ratios of the pressure components among all scales: the compressible part $r_{p_c}=(\hat {p}_c^{\prime },\hat {p}_c^{\prime })_E/V_p$ , the incompressible part $r_{p_i}=(\hat {p}_{ic}^{\prime },\hat {p}_{ic}^{\prime })_E/V_p$ and their coupling $r_{p_cp_i}=(\hat {p}_{c}^{\prime },\hat {p}_{ic}^{\prime })_E/V_p$ . Based on these ratios, three distinct regions in terms of the horizontal scales $(\lambda _x,\lambda _z)$ can be divided, as denoted by the dashed lines (see figure 15 c). The first region (I) is of either small $\lambda _x^+$ or $\lambda _z^+$ ( $\lesssim 30$ ), where $r_{p_c}$ exceeds 40 % and $r_{p_i}$ falls below 50 % and can down to 10 %. These fluctuations reside near the wall where $p_{c}^{{\prime }}$ is the strongest, so they are most affected by acoustic modes, resulting in deviations from their incompressible counterpart. The second region (II) comprises the streamwise-elongated fluctuations of high linear coherence ( $\lambda _z^+\gt 30$ ). They are nearly unaffected by $p_c^{\prime }$ with $r_{p_i}$ over 90 %, so they are dominated by vortical modes. The slightly higher $r_{p_c}$ for the largest scale in the domain $(\lambda _x,\lambda _z)/h=(4\pi ,2\pi )$ may be attributed to the artificial scale truncation by the domain size. In comparison, the spanwise-elongated ones in region III feature higher $r_{p_c}\gtrsim {20}\,\%$ and stronger $\hat {p}_c^{\prime }$ – $\hat {p}_i^{\prime }$ coupling with $|r_{p_cp_i}|$ over 10 %, which is in line with the active zone of $p_{c,w}^{\prime }$ in figure 14(e).

Figure 15. Energy norm ratios of (a) the compressible part $\langle \hat {p}_c^{\prime }\hat {p}_c^{{\prime }{{H}}} \rangle$ , (b) the incompressible part $\langle \hat {p}_{ic}^{\prime }\hat {p}_{ic}^{{\prime }{{H}}} \rangle$ and (c) their coupling $\langle \hat {p}_c^{\prime }\hat {p}_{ic}^{{\prime }{{H}}} \rangle$ relative to $\langle \hat {p}^{\prime }\hat {p}^{{\prime }{{H}}} \rangle$ for case Ma30Re5k. Three featured regions (I, II, III) are divided by the three black dashed lines $\lambda _x=\lambda _z$ , $\lambda _x^+=30$ and $\lambda _z^+=30$ . Extra contours outside the shaded levels are shown in grey dashed lines with labelled levels.

Figure 16. (a) Pre-multiplied 2-D spectrum of the acoustic energy $V_{q_{\textit {ac}}}$ (normalised by $\rho _bU_b^2/h$ ) and the energy-norm ratios of (b) $\langle \,\hat {\!\boldsymbol q}^{{\prime \prime }}\,\hat {\!\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle _{\textit {ac}}$ and (c) $\langle \hat {\rho }^{\prime } \hat {\rho }^{{\prime }{{H}}} \rangle _{\textit {ac}}$ relative to $\langle \,\hat {\!\boldsymbol q}^{{\prime \prime }}\,\hat {\!\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle$ and $\langle \hat {\rho }^{\prime } \hat {\rho }^{{\prime }{{H}}} \rangle$ , respectively, for case Ma30Re5k. Three regions (I, II, III) are divided as in figure 15. The dotted lines in panel (a) are the contours (0.2, 0.4, 0.6, 0.8) of the normalised pre-multiplied spectrum of the total energy $V_q$ .

Next, the variance of $\hat {\boldsymbol q}_{\textit {ac}}^{{\prime \prime }}$ , i.e. $V_{q_{\textit {ac}}}=\|\hat {\boldsymbol q}_{\textit {ac}}^{{\prime \prime }}\|^2$ , is computed, which reflects how much energy is contributed by acoustic modes to the total fluctuation energy in the linear models. The pre-multiplied 2-D spectrum of $V_{q_{\textit {ac}}}$ is first plotted in figure 16(a), along with the $V_q$ spectrum to show the region of energetic total fluctuations. The energy $V_{q_{\textit {ac}}}$ is more pronounced at scales $\lambda _x\lt \lambda _z$ , consistent with figure 14(e). The ratio $r_{q_{\textit {ac}}}=V_{q_{\textit {ac}}}/V_{q}$ for all scales is displayed in figure 16(b), whose pattern is qualitatively the same as in figure 15(a). The three distinct regions (I, II, III) analogous to figure 15 can be identified. The near-wall fluctuations of either small $\lambda _x^+$ or $\lambda _z^+$ ( ${\lt } 30$ ) are most susceptible to acoustic modes, with $r_{q_{\textit {ac}}}\gt {10}\,\%$ . This is consistent with Chen et al. (Reference Chen, Cheng, Fu and Gan2023a ) that the near-wall inner peak of the energy amplification curves, representing the optimal near-wall streak in the buffer layer, becomes non-existent with the rise of ${\textit {Ma}}_b$ . In other words, the travelling-wave-like acoustic components disrupt the near-wall scaling of streamwise-elongated streaks and tend to dominate the near-wall turbulent motions.

When $\lambda _x^+,\lambda _z^+\gt 30$ , $r_{q_{\textit {ac}}}$ is generally less than 2 % for this ${\textit {Ma}}_b=3$ case, indicative of minor contributions from acoustic modes and dominance of vortical modes. This fact justifies the post-processing procedure of the cospectrum decomposition designed by Chen et al. (Reference Chen, Cheng, Gan and Fu2023b ), where the acoustic components are artificially removed from $\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }}\hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle$ to improve the prediction of the linear coherence spectra of velocity and temperature in the eLNS model. The overly amplified acoustic components before the processing ( $r_{q_{\textit {ac}}}$ can reach 20%–55 % in their eLNS model under (2.15)) is due to misordered energies of the POD modes subject to an inaccurate forcing model.

Regarding different components of $\hat {\boldsymbol q}_{\textit {ac}}^{{\prime \prime }}$ , the density fluctuation is most affected by the acoustic components. As shown in figure 16(c), $r_{\rho _{\textit {ac}}}=\|\hat {\rho }_{\textit {ac}}^{\prime }\|^2/\|\hat {\rho }^{\prime }\|^2$ can surpass 10 % for spanwise-elongated fluctuations (region III). In comparison, the corresponding maximum ratios of the velocity components are all less than 3 %, so the velocities are dominated by vortical modes.

Figure 17. Profiles of the total density fluctuation variance and those contributed by the acoustic and non-acoustic parts for scales of (a) $(\lambda _x,\lambda _z)/h=(0.6,6.3)$ , (b) $(\lambda _x,\lambda _z)/h=(12.6,6.3)$ and (c) $(\lambda _x,\lambda _z)/h=(12.6,2.1)$ for case Ma30Re5k. Those predicted by the SRA relation (2.16) are also shown.

Considering the sensitivity of $r_{\rho _{\textit {ac}}}$ to acoustic modes, the wall-normal structure of $\hat {\rho }_{\textit {ac}}^{\prime }$ is presented in figure 17 using three selective scales. The first scale is $(\lambda _x,\lambda _z)/h=(0.6,6.3)$ , corresponding to the peak of $r_{q_{\textit {ac}}}$ (and $r_{\rho _{\textit {ac}}}$ ) in figure 16(b) in region III. The other two scales are both streamwise-elongated ones with $\lambda _x=12.6h$ and in different $\lambda _z$ . In figures 17(b) and 17(c), $\hat {\rho }_{\textit {ac}}^{\prime }$ is much smaller than $\hat {\rho }_{\textit {nac}}^{\prime }$ , and the latter can be well approximated by the SRA relation (2.16). Therefore, $\hat {\rho }^{\prime }$ is transported like a passive scalar controlled by the mean-flow gradients for these streamwise-elongated motions (Coleman et al. Reference Coleman, Kim and Moser1995). Compared with the temperature part, the SRA for density has received less attention. Its effectiveness for these energetic motions also supports the usage of (2.15) for the density and temperature components. In figure 17(a) where $r_{q_{\textit {ac}}}$ peaks, however, $\hat {\rho }_{\textit {ac}}^{\prime }$ is more prominent and is higher than $\hat {\rho }_{\textit {nac}}^{\prime }$ near the wall. The $\hat {\rho }^{\prime }$ at the wall is almost solely contributed by $\hat {\rho }_{\textit {ac}}^{\prime }$ , leaving the residual $\hat {\rho }_{\textit {nac},w}^{\prime }$ nearly zero. This is reasonable because $\hat {\rho }_{\textit {nac}}^{\prime }$ is induced by vortical modes, then $\hat {\rho }_{\textit {nac}}^{\prime }$ should tend to zero near the no-slip wall. The SRA relation (2.16) severely underestimates $\hat {\rho }_{\textit {nac}}^{\prime }$ at this scale because the localised spanwise motion ( $\hat {w}^{{\prime \prime }}$ ) is active for these spanwise-elongated fluctuations.

The above decomposition of $\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle$ into the acoustic and non-acoustic elements enables a decomposition of $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ and $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ , which aids the respective modelling of the acoustic and non-acoustic parts of the forcing. From (2.19), the LNS forcing can be decomposed into $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}=(\boldsymbol{{FF}}^{{H}})_{\textit{DNS},\textit {ac}}+(\boldsymbol{{FF}}^{{H}})_{\textit{DNS},\textit {nac}}$ , with the two parts satisfying

(5.5a) $${\big(\boldsymbol{{FF}}^{{H}}\big)_{\textit{DNS},\textit {ac}} = - \big ( \widehat {\mathcal{L}}_q \big\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \big\rangle _{\textit {ac}} + \big\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \big\rangle _{\textit {ac}} \widehat {\mathcal{L}}_q^{{H}} \big ),}$$

(5.5b) $${\big(\boldsymbol{{FF}}^{{H}}\big)_{\textit{DNS},\textit {nac}} = - \big ( \widehat {\mathcal{L}}_q \big\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \big\rangle _{\textit {nac}} + \big\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \big\rangle _{\textit {nac}} \widehat {\mathcal{L}}_q^{{H}} \big ),}$$

respectively. The decomposition for the eLNS forcing is defined likewise. Different from vortical modes featuring strong non-normality, acoustic modes are nearly perpendicular to each other. Taking the leading acoustic modes (lowest frequency $|\omega |$ ) as a representative, $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS},\textit {ac}}$ can be approximated as $-2\omega _{\textit {ac}} \langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle _{\textit {ac}}$ since $\widehat {\mathcal{L}}_q\,\check {\!\boldsymbol q}^{{\prime \prime }}=\omega _{\textit {ac}} \,\check {\!\boldsymbol q}^{{\prime \prime }}$ (details in Appendix C). Thereby, it can be conjectured that $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS},\textit {ac}}$ resembles the distribution of $\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle _{\textit {ac}}$ except in different amplitudes.

Figure 18. Contours of the decomposed eLNS forcing $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ into the (a) non-acoustic part and (b) acoustic part for the scale $(\lambda _x,\lambda _z)/h=(0.6,6.3)$ for case Ma30Re5k.

This resemblance also applies to $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS},\textit {ac}}$ because the eddy-viscosity-enhanced terms reside in the outer layer, while the acoustic components are mostly amplified near the wall (see figure 14). As a demonstration, $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS},\textit {nac}}$ and $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS},\textit {ac}}$ are depicted in figure 18 for the scale in figure 17(a), at which acoustic modes are prominent. It is observed that $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS},\textit {ac}}$ is indeed confined in the very near-wall region like $\langle \hat {\boldsymbol q}^{{\kern0.3pt}{\prime \prime }} \hat {\boldsymbol q}^{{{\prime \prime }}{{H}}} \rangle _{\textit {ac}}$ . The relative amplitudes of different components ( $\langle \hat {e}_\rho ^{\prime } \hat {e}_\rho ^{{\prime }{{H}}} \rangle _{\textit {ac}}$ , $\langle \hat {e}_\rho ^{{\prime \prime }} \hat {e}_u^{{{\prime \prime }}{{H}}} \rangle _{\textit {ac}}$ , etc.) are thus readily obtained in (5.3). Most of the eLNS forcing is produced by $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS},\textit {nac}}$ resulting from vortical modes, and so is the LNS forcing. The modelling of these non-acoustic parts can learn a lot from the various forcing models for incompressible flows, as introduced in § 1. Systematic improvements on the modelling of $(\boldsymbol{{FF}}^{{H}})_{\textit{DNS}}$ and $(\boldsymbol{{EE}}^{{H}})_{\textit{DNS}}$ for compressible flows will be explored in future works.