2. Methodology

2.1. Governing equations

For a perfect heat-conduction gas, the non-dimensional compressible Navier–Stokes equations are given as

(2.1)$$\begin{gather} \rho\left(\frac{\partial u_{i}}{\partial t}+u_{j}\,\frac{\partial u_{i}}{\partial x_{j}}\right)=-\frac{1}{\gamma M_\infty^{2}}\,\frac{\partial p}{\partial x_{i}}+\frac{1}{Re}\,\frac{\partial}{\partial x_{j}}\left[\mu\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)+\lambda\,\frac{\partial u_{k}}{\partial x_{k}}\,\delta_{i j}\right]\!, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial \rho}{\partial t}+u_{j}\,\frac{\partial \rho}{\partial x_{j}}=-\rho\,\frac{\partial u_{j}}{\partial x_{j}}, \end{gather}$$

(2.3)\begin{align} \rho\left(\frac{\partial T}{\partial t}+u_{j}\,\frac{\partial T}{\partial x_{j}}\right)&=-(\gamma-1) p\,\frac{\partial u_{i}}{\partial x_{i}}+\frac{\gamma}{Re}\, \frac{\partial}{\partial x_{j}}\left(k\,\frac{\partial T}{\partial x_{j}}\right) \nonumber\\ &\quad +\gamma(\gamma-1)\,\frac{M_\infty^{2}}{Re} \left[\mu\,\frac{\partial u_{i}}{\partial x_{j}}\,\frac{\partial u_{i}}{\partial x_{j}}+\mu\,\frac{\partial u_{i}}{\partial x_{j}}\,\frac{\partial u_{j}}{\partial x_{i}}+\lambda\left(\frac{\partial u_{k}}{\partial x_{k}}\right)^{2}\right]\!, \end{align}

where $t$ is time, $x_k$ ($k=1,2,3$) denote the streamwise, wall-normal and spanwise directions, respectively, and $u_k$ are the corresponding velocity components. Here, $\rho$ is density, $p$ is pressure, and $T$ is temperature, and they satisfy the ideal gas state equation $p=\rho T$. Also, $\mu$ and $\lambda = -2/3\mu$ are the coefficients of the molecular and second viscosity, respectively. The heat conduction coefficient $k$ is defined by $\mu /{Pr}$, where ${Pr}$ is the molecular Prandtl number. Also, $\gamma$ is the specific heat ratio, $\delta _{ij}$ is the Kronecker delta, and $M_\infty$ is the freestream Mach number. The freestream Reynolds number ${Re}$ is defined by $\rho _\infty ^* u_\infty ^* \delta _{99}^*/\mu _\infty ^*$, with $\delta _{99}$ being the 99 % boundary-layer thickness. The subscript $\infty$ denotes the freestream state, and the superscript $*$ denotes dimensional values, otherwise all variables are non-dimensionalized by the freestream flow quantities.

2.2. Resolvent analysis

The boundary-layer flow is assumed to be locally parallel (Bae et al. Reference Bae, Dawson and McKeon2020a), that is, the effects of streamwise development are ignored. Thus flow quantities $\boldsymbol {q}=[u_{1}, u_{2}, u_{3}, \rho, T]^{\rm T}$ can be Fourier transformed in the temporal, streamwise and spanwise directions, namely

(2.4)\begin{equation} \boldsymbol{q}\left(x_{1}, x_{2}, x_{3}, t\right)=\iiint_{-\infty}^{\infty} \hat{\boldsymbol{q}}\left(x_{2} ; \kappa_{1}, \kappa_{3}, \omega\right) \exp({\mathrm{i}(\kappa_{1} x_{1}+\kappa_{3} x_{3}-\omega t)})\,\mathrm{d} \kappa_{1} \,\mathrm{d} \kappa_{3} \,\mathrm{d} \omega, \end{equation}

where $\widehat {(\cdot )}$ denotes the Fourier-transformed variable, $\kappa _1$ and $\kappa _3$ are the streamwise and spanwise wavenumbers, $\omega$ is the temporal frequency, and $\mathrm {i}=\sqrt {-1}$. The speed of the streamwise wave is defined by $c = \omega /k_x$.

Substituting (2.4) into the non-dimensional compressible Navier–Stokes equations (2.1)–(2.3), we get a linearized form for each $(\kappa _1, \kappa _3, \omega )\ne\, (0,0,0)$ (Mack Reference Mack1984),

(2.5)\begin{equation} \left[-\mathrm{i} \omega \boldsymbol{I}+\boldsymbol{L}(\kappa_{1}, \kappa_{3}, \omega)\right]\hat{\boldsymbol{q}}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega)=\hat{\boldsymbol{f}}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega), \end{equation}

where $\boldsymbol {I}$ is the identity matrix. The high-order nonlinear contributions in the Navier–Stokes equations are all represented within a single forcing term $\hat {\boldsymbol {f}}$. For the detailed expression of $\boldsymbol {L}$, readers can refer to Fan et al. (Reference Fan, Li and Sandberg2023).

Equation (2.5) leads to an input–output form

(2.6)\begin{equation} \hat{\boldsymbol{q}}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega)=\boldsymbol{H}(\kappa_{1}, \kappa_{3}, \omega)\, \hat{\boldsymbol{f}}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega), \end{equation}

where $\boldsymbol {H}(\kappa _{1}, \kappa _{3}, \omega )=(-\mathrm {i} \omega \boldsymbol {I}+\boldsymbol {L}(\kappa _{1}, \kappa _{3}, \omega ))^{-1}$ is termed the resolvent operator, which linearly relates the input forcing to the output state. To characterize the resolvent operator $\boldsymbol {H}$, a singular value decomposition is used, with a weight matrix taken into account for a physically appropriate norm for compressible turbulent boundary layers, that is,

(2.7a–c)\begin{equation} \boldsymbol{W}^{1/2}\boldsymbol{H}\boldsymbol{W}^{-1/2} = \sum_{p=1}^{N}\tilde{\boldsymbol{\psi}}_p \sigma_p \tilde{\boldsymbol{\phi}}_p^{\dagger},\quad \boldsymbol{\phi}_p =\boldsymbol{W}^{-1/2} \tilde{\boldsymbol{\phi}}_p,\quad \boldsymbol{\psi}_p =\boldsymbol{W}^{-1/2} \tilde{\boldsymbol{\psi}}_p, \end{equation}

where $N$ is the number of resolvent modes, $\boldsymbol {\psi }_p$ and $\boldsymbol {\phi }_p$ are the $p$th orthogonal basis functions of the response and forcing modes, and the real singular value $\sigma _p$ represents the linear amplification rate ranked by $\sigma _{p}>\sigma _{p+1}$. The weight matrix selected here is defined by Chu (Reference Chu1965), i.e. $\boldsymbol {W}={\rm diag}(\bar {\rho }, \bar {\rho }, \bar {\rho }, {\bar {T}}/{\gamma M^2 \bar {\rho }}, {\bar {\rho }}/{\gamma (\gamma -1)M^2\bar {T}}) \,\mathrm {d}\kern0.06em x_2$, where the overline represents Reynolds averaging. In such a way, the corresponding norm $\left \langle \hat {\boldsymbol {q}},\hat {\boldsymbol {q}}\right \rangle =\hat {\boldsymbol {q}}^{\dagger} \boldsymbol {W}\hat {\boldsymbol {q}}$ (where the superscript ${\dagger}$ denotes the complex conjugate transpose) eliminates the pressure-related energy transfer terms (compression work) and has been used widely in studies of compressible flows.

2.3. Improved resolvent analysis with addition of eddy viscosity

To model part of the pivotal nonlinear contribution, a widely acknowledged method is to add an eddy-viscosity term in the linearized Navier–Stokes equations. The role of the eddy viscosity is to introduce a closure, establishing a constitutive relation between the oscillation of Reynolds stresses associated with the passage of the disturbance and the strain rate of the turbulent motions (Illingworth et al. Reference Illingworth, Monty and Marusic2018; Symon et al. Reference Symon, Illingworth and Marusic2021; Holford, Lee & Hwang Reference Holford, Lee and Hwang2023). The same approach is applied to the relation between the fluctuations of turbulent heat flux and the temperature gradient in the energy equation. In this way, a new set of non-dimensional linearized Navier–Stokes equations is obtained:

(2.8)\begin{align} \hat{\boldsymbol{q}}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega)&= \boldsymbol{H}_e(\kappa_{1}, \kappa_{3}, \omega)\,\hat{\boldsymbol{f}}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega)\nonumber\\ &=\left(-\mathrm{i} \omega \boldsymbol{I}+\boldsymbol{L}_e(\kappa_{1}, \kappa_{3}, \omega)\right)^{-1} \hat{\boldsymbol{f}}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega). \end{align}

In the eddy-embedded resolvent operator $\boldsymbol {H}_e$, substitutions are made for the momentum and energy equations, that is,

(2.9a–c)\begin{equation} \mu \rightarrow \mu + \mu_t, \quad \lambda \rightarrow \lambda+ \lambda_t \quad {\rm and}\quad k \rightarrow k + k_t = \frac{\mu}{Pr} + \frac{\mu_t}{{Pr}_t},\end{equation}

where $\mu _t$ and $\lambda _t$ are the eddy molecular and second viscosity coefficients, respectively. Here, ${Pr}_t$ is the turbulent Prandtl number, which is nearly constant in most regions of the boundary layer and is mostly insensitive to the freestream Mach number and wall temperature conditions (Zhang, Duan & Choudhari Reference Zhang, Duan and Choudhari2018). In this study, we set ${Pr}_t$ to 0.85, which is informed by previous studies (Zhang et al. Reference Zhang, Bi, Hussain and She2014), and our preliminary tests showed that this value does not influence the prediction results, which was also reported by Chen et al. (Reference Chen, Cheng, Fu and Gan2023). As the additional part compensates for the inter-scale energy transfer in the physical process, the resolvent operator $\boldsymbol {H}_e$ is able to provide a more promising tool for prediction than $\boldsymbol {H}$. Accordingly, this eddy-viscosity-improved resolvent operator can also be characterized by singular value decomposition, leading to the response and forcing modes as $\boldsymbol {\psi }_{e,p}$ and $\boldsymbol {\phi }_{e,p}$. The definition of the norm $\left \langle \hat {\boldsymbol {q}},\hat {\boldsymbol {q}}\right \rangle =\hat {\boldsymbol {q}}^{\dagger} \boldsymbol {W}\hat {\boldsymbol {q}}$ given by Chu (Reference Chu1965) is used here, yielding an energy eliminating the pressure-related work.

2.3.1. A modified Cess eddy-viscosity model

An issue here is how to define the profile of eddy viscosity ($\mu _t$) in the compressible boundary layer. In resolvent analysis of incompressible channel flows, the Cess approximation has been used commonly and proved to be efficient in recovering the energetic modes. As is derived from the mean momentum equation for incompressible channel flows, this algebraic total-viscosity model (Cess Reference Cess1958) is given as

(2.10)\begin{align} \frac{\mu_{T,Cess}}{\mu_w} &= \frac{1}{2}\left( 1+ \frac{\kappa^2 \,{Re}_\tau^2}{9} \left(2x_2-x_2^2\right)^2 \left(3-4x_2+2x_2^2 \right)^2 \left[1-\exp\left( \frac{-{Re}_\tau\,x_2}{A} \right) \right]^2\right)^{1/2} \nonumber\\ &\quad + \frac{1}{2}, \end{align}

where $\kappa =0.426$ and $A=25.4$, following previous studies (Pujals et al. Reference Pujals, García-Villalba, Cossu and Depardon2009; Amaral et al. Reference Amaral, Cavalieri, Martini, Jordan and Towne2021). The subscript $w$ denotes quantities at the wall.

However, in the case of compressible turbulent boundary layers, the situation differs (see the Introduction). The turbulent boundary layer is a composite layer, which can be characterized by two regions associated with different fluid responses to shear and pressure gradient (Cebeci Reference Cebeci1971). Therefore, in the present study, we decompose the boundary layer into two separate regions, i.e. the inner and outer regions, where different strategies of modelling are developed. Note that the so-called ‘inner’ and ‘outer’ regions in the remainder of the paper are here used as broad concepts, distinguishing the regions nearer the wall ($x_2\leq \delta _c$) and further away from the wall ($x_2 > \delta _c$), where $\delta _c$ is the crossover wall-normal position. This differs from the classical identification of layers related to the mean-velocity behaviour.

In the inner region, we adopt the Cess approximation, considering that the dynamics there is very similar to that in channel flows (Pope Reference Pope2000), whereas in the outer region, since the boundary-layer flow behaviour departs from that in channels, the Cess model (2.10), which is specifically derived for incompressible channel flows, is not applicable. Moreover, the Cess formulation does not consider the region outside the boundary-layer thickness. In this sense, we propose to define the eddy viscosity based on Prandtl's mixing-length hypothesis in the outer region. Consequently, a two-layer algebraic eddy-viscosity model is proposed, generalizing the Cess model to the regime of compressible turbulent boundary layers:

(2.11)\begin{equation} \mu_t = \begin{cases} \sqrt{\dfrac{\bar{\rho}}{\rho_w}}\,\mu_{T,Cess}-\bar{\mu}, & {\rm for} \ x_2\leq\delta_c, \\ {Re}\,\bar{\rho} l_m^2 \left|\dfrac{{\rm d}\bar{u}_1}{{\rm d}x_2}\right|, & {\rm for} \ x_2>\delta_c, \end{cases} \end{equation}

where $l_m$ is the mixing length (Pope Reference Pope2000).

In the modified model (2.11), when $x_2 \leq \delta _c$, the original Cess total-viscosity definition is premultiplied by $\sqrt {\bar {\rho }/\rho _w}$, so that the velocity profile obtained from the integral of ${Re}_\tau \,(1-x_2)\mu _w/\mu _T$ still follows the classical velocity laws, upon van Driest transformation (van Driest Reference van Driest1951) (especially suitable for the adiabatic cases). The subscript $T$ denotes the total quantity, which means $\mu _T=\mu _t+\bar {\mu }$, with the first term representing the turbulent contribution, and the second the mean molecular viscosity. When $x_2 > \delta _c$, the eddy viscosity is determined by Prandtl's mixing-length hypothesis. It is necessary to specify the expression of $l_m$ that we are concerned about. Normally, it can be approximated as a constant in the outer region (core region) (Cebeci Reference Cebeci2004). For example, Escudier (Reference Escudier1966) set $l_m$ to be the minimum value of $\kappa x_2$ (which is assumed to be the $l_m$ expression in the logarithmic region) in the core region and 0.09. In the present study, to ensure the continuity of the $\mu _t$ profile across the boundary layer, $l_m$ is determined so that the $\mu _t$ at $x_2=\delta _c$ calculated from the first and second parts of (2.11) are equal. This treatment is believed to be reasonable, and the following results suggest that it yields a good approximation of turbulent eddy viscosity in the outer region.

2.3.2. A new semi-empirical eddy-viscosity model

We propose another approach to estimate the eddy-viscosity profile, based on an empirical law and the mixing-length theory, which is especially suitable for zero-pressure-gradient turbulent boundary layers. It is still a two-layer eddy-viscosity model, with major differences from (2.11) in the inner region. Instead of a Cess-based correction originated from the channel-flow models, here we start from the universal scaling of the total shear stress, which is exclusive for turbulent boundary layers and hence will potentially lead to a more accurate eddy-viscosity profile.

In zero-pressure-gradient turbulent boundary layers, Hou, Somandepalli & Mungal (Reference Hou, Somandepalli and Mungal2006) found that the $(1-x_2)$ weighted total shear stress ($\tau$) can be fitted by a linear expression in the region $x_2\leq 0.5$. This $(1-x_2)$ fit is expressed as

(2.12)\begin{equation} (1-x_2)\,\frac{\tau}{\tau_w} = ax_2+b, \quad {\rm for} \ x_2\le0.5, \end{equation}

where ${\tau _w}$ is the wall shear stress. The intercept $b$ is fixed as $b=1$ to ensure $\tau =\tau _w$ at $x_2=0$. The slope $-a$ lies in the range from $1.2$ to $1.6$ (Hou et al. Reference Hou, Somandepalli and Mungal2006). The databases of numerical and experimental results for a wide range of Reynolds numbers suggested that there is no clear dependence of $-a$ on the Reynolds number (Hou et al. Reference Hou, Somandepalli and Mungal2006; Xia, Zhang & Yang Reference Xia, Zhang and Yang2021).

To examine the feasibility of (2.12) for compressible turbulent boundary layers, the profiles of the weighted total shear stress at various Mach numbers and Reynolds numbers are depicted in figure 1. Details of the cases at the freestream Mach number $M=2,3,4$ with the friction Reynolds number $Re_\tau$ ranging from 250 to 1110 can be retrieved from Pirozzoli & Bernardini (Reference Pirozzoli and Bernardini2011) and Bernardini & Pirozzoli (Reference Bernardini and Pirozzoli2011). The hypersonic databases at $M=5.86$, $Re_\tau \approx 420$ and the wall-to-recovery temperature ratios $T_w/T_r=1.0$ and $0.25$, representing the adiabatic and cold-wall conditions, are obtained from DNS in Fan, Li & Pirozzoli (Reference Fan, Li and Pirozzoli2022) and Fan & Li (Reference Fan and Li2023). In figure 1, good alignment is observed between the weighted total shear stress and the fitting law $-1.3x_2+1$ in $x_2 \le 0.5$, for all of the compressible cases considered, confirming its suitability for the empirical modelling.

Figure 1. The DNS profiles of $(1-x_2){\tau }/{\tau _w}$ in compressible zero-pressure-gradient turbulent boundary layers, at: $M=2$, $Re_\tau \approx 250\unicode{x2013}1110$ (Pirozzoli & Bernardini Reference Pirozzoli and Bernardini2011); $M=3$, $Re_\tau \approx 500$ and $M=4$, $Re_\tau \approx 500$ (Bernardini & Pirozzoli Reference Bernardini and Pirozzoli2011); and $M =5.86$, $Re_\tau \approx 420$ under the adiabatic and cold-wall conditions (Fan et al. Reference Fan, Li and Pirozzoli2022; Fan & Li Reference Fan and Li2023). The black dashed line denotes the linear expression $-1.3 x_2+1$.

In the region further away from the wall ($x_2>\delta _c$), mixing-length theory is also exploited, leading to the analytical description of the new eddy-viscosity profile across the boundary layer:

(2.13)\begin{equation} \mu_t = \begin{cases} \dfrac{\mu_w}{{\rm d} \bar{u}_1^+{/}{\rm d} \bar{x}_2^+}\,\dfrac{-1.3x_2+1}{1-x_2} -\bar{\mu} , & {\rm for} \ x_2\leq\delta_c, \\ {Re}\,\bar{\rho} l_m^2 \left|\dfrac{{\rm d}\bar{u}_1}{{\rm d}\kern 0.06em x_2}\right| , & {\rm for} \ x_2>\delta_c, \end{cases} \end{equation}

where the superscript $+$ denotes normalization by wall units. To ensure the applicability of the model in $x_2\leq \delta _c$, $\delta _c$ should not be larger than 0.5. On the other hand, a crossover position too close to the wall is also avoided under the constant-$l_m$ assumption. Conservatively, we make $\delta _c= 0.5$ for the cases under scrutiny in the following study.

2.3.3. A priori comparison with DNS results

Hereafter, we focus on hypersonic turbulent boundary layers at $M=5.86$ with adiabatic and cold-wall conditions, since the compressibility effects in these cases are strong and, more importantly, we have detailed DNS solutions to evaluate the performance of the resolvent-based modelling. With the models (2.11) and (2.13), the distribution of the total shear stress across the boundary layer can be obtained on the basis of known mean flow quantities (a requirement for resolvent analysis), including the velocities, density and temperature. Figure 2 compares the modelled weighted total shear stress ($(1-x_2)\tau /\tau _w = (1-x_2) (\mu _t +\bar {\mu })({\rm d}\bar {u}_1/{\rm d}\kern 0.06em x_2)/\tau _w$) with the DNS results. The total stresses predicted with the semi-empirical model are found to collapse perfectly onto the DNS profiles, whereas with the Cess-correction model, there are obvious deviations, especially in the inner region.

Figure 2. Modelled profiles of $(1-x_2){\tau }/{\tau _w}$ in the adiabatic and cold-wall hypersonic turbulent boundary layer at $M =5.86$.

Figure 3 quantifies the total viscosity ($\mu _T/\mu _w$) profile across the boundary layer. It can be seen that both models proposed here capture the trends of the DNS profiles. Note that the eddy viscosity in DNS is solved from the Reynolds stress expressions based on the turbulent viscosity hypothesis. Nonetheless, the semi-empirical eddy-viscosity model is more accurate, with the profiles almost collapsed onto those of the DNS, according to the a priori examinations. In the next section, both eddy-viscosity models will be utilized in the resolvent operator, to examine their roles in the performance of resolvent-based flow modelling.

Figure 3. Profiles of the total viscosity $\mu _T/\mu _w$ in the (a) adiabatic and (b) cold-wall hypersonic turbulent boundary layer at $M =5.86$.

3. Consistency between the resolvent modes and SPOD modes

Spectral proper orthogonal decomposition (Lumley Reference Lumley1967; Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993; Towne et al. Reference Towne, Schmidt and Colonius2018) is known as a promising data-driven tool to construct a reduced-order model of fluctuating flow fields, from the perspective of energy. It extracts the orthogonal basis functions at discrete frequencies from a series of snapshots using Welch's method, with the leading SPOD modes to be optimally representative of the turbulent energy of the system (Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017; Towne et al. Reference Towne, Schmidt and Colonius2018). This section aims to examine the similarity between the SPOD modes issued from the DNS data and the response modes from the resolvent analysis with/without eddy viscosity, to evaluate the feasibility of resolvent-based reduced-order modelling (Abreu et al. Reference Abreu, Cavalieri, Schlatter, Vinuesa and Henningson2020; Pickering et al. Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021) in the compressible regime.

Figure 4 quantifies the energy ratios of the first eight leading resolvent and SPOD modes, at two representative scales in the adiabatic and cold-wall turbulent boundary layers. The scales are selected at positions of local energy peaks characterizing the inner- and outer-layer dynamics, referring to our previous study (Fan & Li Reference Fan and Li2023). It is observed that the first mode in the standard resolvent and SPOD analysis captures near or more than 50 % of the total energy. The first six modes are capable of recovering more than 85 % of the small-scale energy, and more than 90 % of the large-scale energy. This feature is not evident in the eddy-viscosity-improved resolvent analysis, with the first mode capturing only 2 %–20 % of the total energy, and the higher-order modes still holding significant energy. Nevertheless, such results do not reduce the accuracy of the response modes in the prediction of flow structures (Symon et al. Reference Symon, Madhusudanan, Illingworth and Marusic2023), since the singular value does not strictly convey the weight of the corresponding mode in recovering the flow field in resolvent analysis. The weight function is determined by both the singular value and the projection of the forcing mode on the real nonlinear process; see (14) in Fan et al. (Reference Fan, Li and Sandberg2023).

Figure 4. Energy contained in the $i$th resolvent and SPOD modes relative to the total, defined by $\sigma _i^2/\sum _j \sigma _j^2$ and $\varLambda _i/\sum _j \varLambda _j$ (where $\varLambda$ denotes the eigenvalue in the SPOD analysis), respectively, at (a) $(\lambda _1^+,\lambda _3^+,c)=(465,96,0.52)$ and (b) $(\lambda _1^+,\lambda _3^+,c)=(1210,450,0.89)$ in the adiabatic turbulent boundary layer, and (c) $(\lambda _1^+,\lambda _3^+,c)=(929,208,0.66)$ and (d) $(\lambda _1^+,\lambda _3^+,c)=(1239,459,0.88)$ in the cold-wall turbulent boundary layer.

Considering the representativeness of the first SPOD mode in terms of energy, we first examine its alignment with the optimal resolvent response mode with/without eddy viscosity. Figure 5 presents the wall-normal profiles of the streamwise velocity ($\hat {u}_1$), wall-normal velocity ($\hat {u}_2$), spanwise velocity ($\hat {u}_3$), density ($\hat {\rho }$) and temperature ($\hat {T}$) components of the optimal resolvent and the first SPOD mode at the representative scales in the adiabatic turbulent boundary layer. Generally, in figures 5(a–e), all three types of resolvent-based results are able to capture approximately the features in the distribution of perturbations, especially with respect to the correctly identified most energetic location following the critical-layer arguments (McKeon & Sharma Reference McKeon and Sharma2010). More specifically, better agreements can be observed at the tails between the eddy-viscosity-improved resolvent and SPOD modes for all components, in contrast to the performance of the standard resolvent mode. It thus suggests a higher prediction accuracy of resolvent analysis with the addition of eddy viscosity in the near-wall cycle. The profiles with the addition of the Cess-correction and semi-empirical models are almost collapsed onto each other. A further comparison of the two eddy-viscosity models is performed, in terms of the mode shapes of the streamwise velocity. Details are shown in the Appendix. The remarkable consistency between the two models confirms that the optimal resolvent mode is not sensitive to the small differences in the eddy-viscosity quantities, as long as the $\mu _t$ model can resemble its real distribution appropriately (Cossu, Pujals & Depardon Reference Cossu, Pujals and Depardon2009).

Figure 5. The streamwise velocity, wall-normal velocity, spanwise velocity, density and temperature components of the optimal resolvent with/without adding the semi-empirical eddy viscosity model and the leading SPOD mode at (a–e) $(\lambda _1^+,\lambda _3^+,c)=(465,96,0.52)$ and (f–j) $(\lambda _1^+,\lambda _3^+,c)=(1210,450,0.89)$, in the adiabatic turbulent boundary layer.

For the larger energetic scales in the outer region, e.g. $(\lambda _1^+,\lambda _3^+,c)=(1210,450,0.89)$ in figures 5(f–j), the resolvent modes without eddy viscosity are localized narrowly, since they neglect the true nonlinear transfer of energy, hence the local production is dissipated rapidly by molecular viscosity. On the other hand, the predictive performances of the proposed eddy-viscosity-improved resolvent analysis are nearly perfect for the components $\hat {u}_1$, $\hat {\rho }$ and $\hat {T}$, by modulating the energy transfer mechanism upon adding the eddy viscous work (Symon et al. Reference Symon, Illingworth and Marusic2021). As for the components $\hat {u}_2$ and $\hat {u}_3$ (in figures 5g,h), the resolvent-based predictions are degraded when compared with the SPOD profiles. Underestimations of these components are observed without eddy viscosity, which is associated with the nature of high non-normality in the shear flow operator (Symon et al. Reference Symon, Madhusudanan, Illingworth and Marusic2023). The addition of either eddy-viscosity profile is able to reduce the deviation between the resolvent and SPOD modes to some extent, exhibiting larger values of $|\hat {u}_2|$, $|\hat {u}_3|$ and $|\hat {\rho }|$ than the mode without eddy viscosity.

In the cold-wall case, a similar phenomenon is seen in figure 6. The addition of eddy viscosity makes the response modes richer in light of the wide spanning in the wall-normal direction, as mentioned earlier. Hence it provides a more efficient basis for turbulence prediction. In particular for the prediction of density and temperature perturbations close to the wall (in figures 6d,e), the inherent dual-peak feature is identified reasonably well in all types of resolvent modes, though no significant improvement is observed in the resolvent analysis with eddy viscosity. The dual-peak feature is ascribed to the non-monotonous profile of the mean temperature in the cold-wall turbulent boundary layer. The sign change in the mean temperature gradient leads to a reversed conductive and turbulent heat transfer in the wall-normal direction, which challenges the basic assumptions in the streamwise-velocity–temperature correlation (Wenzel, Gibis & Kloker Reference Wenzel, Gibis and Kloker2022) and makes the prediction of the thermal statistics more difficult.

Figure 6. The streamwise velocity, wall-normal velocity, spanwise velocity, density and temperature components of the optimal resolvent with/without adding the semi-empirical eddy viscosity model and the leading SPOD mode at (a–e) $(\lambda _1^+,\lambda _3^+,c)=(929,208,0.66)$ and (f–j) $(\lambda _1^+,\lambda _3^+,c)=(1239,459,0.88)$, in the cold-wall turbulent boundary layer.

To assess quantitatively the alignment between the resolvent and SPOD modes, a coefficient $\beta _{11}$, which is defined by the projection of the optimal resolvent mode ($\boldsymbol {\psi }_{1}$ or $\boldsymbol {\psi }_{e,1}$) onto the first SPOD mode ($\hat {\boldsymbol {q}}_{SPOD1}$) at each wavenumber combination, is introduced. It is written as

(3.1)\begin{equation} \beta_{11}(\kappa_{1}, \kappa_{3}, \omega) = \frac{\left|\left\langle\hat{\boldsymbol{q}}_{SPOD1}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega),\boldsymbol{\psi}_{1}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega)\right\rangle\right|}{ \|\hat{\boldsymbol{q}}_{SPOD1}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega)\| \times \|\boldsymbol{\psi}_{1}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega)\|}, \end{equation}

where $\langle \cdot,\cdot \rangle$ denotes the inner product defined by Chu (Reference Chu1965), and $\|\cdot \|$ represents the weighted $L_2$ norm. Note that $\boldsymbol {\psi }_{1}$ is replaced by $\boldsymbol {\psi }_{e,1}$ when taking into account the eddy viscosity. Having $\beta _{11} = 1$ indicates perfect agreement between the two, giving identical identifications of turbulent structures and energy distribution, whereas $\beta _{11} = 0$ suggests that the SPOD mode and the resolvent mode are orthogonal to each other (Abreu et al. Reference Abreu, Cavalieri, Schlatter, Vinuesa and Henningson2020; Pickering et al. Reference Pickering, Rigas, Schmidt, Sipp and Colonius2021). For the scales under scrutiny in figures 5 and 6, the projection coefficient $\beta _{11}$ is calculated and listed in table 1. Remarkable improvement of the representativeness of the optimal resolvent mode is seen clearly, especially for the larger scale where the added eddy viscosity plays a significant role in modelling the nonlinear energy transfer.

Table 1. The alignment between the optimal resolvent and the first SPOD mode at two selected scales.

To examine the alignment at other scales, we sample at varying spanwise wavelengths and a fixed streamwise wavelength and wave speed, which is cross-sectional in the near-wall and outer-layer cycle, i.e. $(\lambda _1^+,c)=(430,0.57)$ and $(\lambda _1^+,c)=(1291,0.8)$ for the adiabatic case, and $(\lambda _1^+,c)=(929,0.66)$ and $(\lambda _1^+,c)=(1239,0.88)$ for the cold-wall case. The variation of $\beta _{11}$ as a function of $\lambda _3^+$ is shown in figures 7 and 8. The background colour represents the turbulent energy spectra at the corresponding wall-normal location in the sense of the critical layer. As anticipated, for the near-wall scales shown in figures 7(a) and 8(a), $\beta _{11}$ is larger than or close to 0.8 in the high-energy-containing region (e.g. $20<\lambda _3^+<175$ in the adiabatic case, and $60<\lambda _3^+<400$ in the cold-wall case), between either type of resolvent and SPOD modes. Only small variations are induced by the addition of eddy viscosity. As highlighted by the red dashed box, the so-called ‘high-energy-containing region’ covers more than 70 % of the total energy at respective heights. Beyond this region, the value of $\beta _{11}$ reduces dramatically, thus indicating a poor prediction of the low-energy-containing structures. Nonetheless, this biased predictive capability is believed to have limited impacts on constructing the full flow field, as the contribution of low-energy-containing structures is minor with low energy amplification.

Figure 7. Variation of $\beta _{11}$ with regard to the spanwise wavelength, in the adiabatic turbulent boundary layer. The background colour represents the normalized energy spectra at matching streamwise wavelengths and wall-normal locations.

Figure 8. Variation of $\beta _{11}$ with regard to the spanwise wavelength, in the cold-wall turbulent boundary layer.

To be more specific, we quantify the energy captured by the resolvent mode in contrast to that contained in the first SPOD mode ($E_{qq,SPOD1}$) at the characteristic wall-normal location, as shown in figure 9. The $E_{qq,SPOD1}$ value is calculated based on the energy spectra at matching streamwise wavelengths and wall-normal locations, which are obtained from DNS, premultiplied by the rate of the first eigenvalue to the total (i.e. $\varLambda _1/\sum _j \varLambda _j$). In the adiabatic case, the energy distribution exhibited by the eddy-viscosity-improved resolvent mode aligns quite well with that of the first SPOD mode, especially in the highlighted ‘high-energy-containing region’. In general, both types of eddy-viscosity-improved modes are able to capture up to 53 % of the SPOD energy when integrated in the spanwise wavenumber direction, while the standard resolvent mode captures approximately 38 %. The situation is similar in the cold-wall turbulent boundary layer, with approximately 50 % of the energy captured in the resolvent mode with eddy viscosity, while 45 % is captured in the standard resolvent mode.

Figure 9. Energy distribution in the optimal resolvent modes and the first SPOD mode, for the near-wall cycle in the (a) adiabatic and (b) cold-wall turbulent boundary layers. The grey background covers the high-energy-containing region.

As for the energetic large-scale motions characterizing the outer-layer dynamics in figures 7(b) and 8(b), although the general performance is poorer than that for the small scales, the addition of the proposed eddy-viscosity model indeed enhances the projection coefficients. Hence it confirms the efficient role of the rational eddy-viscosity models in resolvent-based reduced-order modelling.

In addition to the wall-normal profiles of the state variables, the distribution of the cross-spectral density (CSD) is of great concern in turbulence prediction, since it can properly reconstruct the second-order statistics of the flow (Towne et al. Reference Towne, Schmidt and Colonius2018). The low-rank approximation of the CSD tensor is calculated by the optimal response mode at a particular wavenumber pair:

(3.2)\begin{equation} S_{11}(x_{2},x_{2}' ; \kappa_{1}, \kappa_{3}, \omega) = \hat{\boldsymbol{\psi}}_1(x_{2} ; \kappa_{1}, \kappa_{3}, \omega)\,\hat{\boldsymbol{\psi}}_1^{\dagger} (x_{2}'; \kappa_{1}, \kappa_{3}, \omega). \end{equation}

The absolute value of the modelled rank-1 CSD for each component (including $\hat {u}_1$, $\hat {u}_2$, $\hat {u}_3$, $\hat {\rho }$ and $\hat {T}$) is plotted as a function of $x_2$ and $x_2'$ in figure 10, at the two selected scales in the adiabatic turbulent boundary layer. At small scales (figure 10a), the three types of resolvent-based models are expected to reproduce faithfully the CSD from SPOD results. A slight exception occurs in the autocorrelation and cross-correlations of $\hat {u}_3$. As we see in the third row and the third column, the $\mu _t$-modelled energy densities appear at a lower height. This indicates that the addition of eddy viscosity introduces a strong energy transfer among components near the wall, since the spanwise velocity receives energy from the streamwise and wall-normal components, due to the effect of pressure strain and splat (Lee & Moser Reference Lee and Moser2019; Fan et al. Reference Fan, Li and Pirozzoli2022). Moreover, the additional eddy diffusion also plays a critical role in the positive energy receipt in the near-wall region (Symon et al. Reference Symon, Madhusudanan, Illingworth and Marusic2023). At a larger energetic scale (figure 10b), the CSD results using the Cess-correction and semi-empirical models are also observed to be very similar to each other, and they exhibit much better alignment with the DNS data in contrast to those without an eddy-viscosity model. A very narrow distribution is observed in the CSD without eddy viscosity, and the turbulent energy decays from its peak dramatically, consistent with the observations in the profiles of the response mode. This issue can be alleviated efficiently with the addition of eddy viscosity, suggesting that the eddy-viscosity models help to improve the mode representation of the second-order turbulent statistics. Hence the influence of the dominant turbulent structures can be identified more accurately in the eddy-viscosity-improved resolvent analysis.

Figure 10. Absolute value of the CSD for each component, at (a) $(\lambda _1^+,\lambda _3^+,c)=(465,96,0.52)$ and (b) $(\lambda _1^+,\lambda _3^+,c)=(1210,450,0.89)$, in the adiabatic turbulent boundary layer. The background colour represents the CSD of the first SPOD results. The contour lines indicate 0.3, 0.6 and 0.9 of their respective maximum, moving inwards, and the legend for the coloured lines is as in figure 7.

For the CSD prediction in the cold-wall cases, as shown in figure 11, the situation is similar. The predictive performance is greatly improved by using the eddy-viscosity-improved resolvent analysis. Therefore, the reliability and improvement of the proposed models in comparison with the standard resolvent-based modelling, for the prediction of turbulent fluctuations in the compressible turbulent boundary layers, are believed to be achieved.

Figure 11. Absolute value of the CSD for each component, at (a) $(\lambda _1^+,\lambda _3^+,c)=(929,208,0.66)$ and (b) $(\lambda _1^+,\lambda _3^+,c)=(1239,459,0.88)$, in the cold-wall turbulent boundary layer.

At last, although the first (optimal) mode is believed to be of the greatest importance in terms of energy occupation, the alignment of the other resolvent and SPOD modes is also of interest. We project the first four response modes onto the corresponding SPOD modes, and the projection coefficients are shown in figures 12 and 13 for the adiabatic and cold-wall turbulent boundary layers, respectively. The projection coefficient is formulated as

(3.3)\begin{equation} \beta_{ij}(\kappa_{1}, \kappa_{3}, \omega) = \frac{\left|\left\langle \hat{\boldsymbol{q}}_{SPODi}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega),\boldsymbol{\psi}_{j} (x_{2} ; \kappa_{1}, \kappa_{3}, \omega)\right\rangle\right|}{ \| \hat{\boldsymbol{q}}_{SPODi}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega)\| \times \|\boldsymbol{\psi}_{j}(x_{2} ; \kappa_{1}, \kappa_{3}, \omega)\|}, \end{equation}

where $\hat {\boldsymbol {q}}_{SPODi}$ denotes the $i$th SPOD mode. For the near-wall cycle in both cases, evident promotion of the eddy-viscosity models is observed only in some of the modes, e.g. the first and third resolvent modes at the scale under scrutiny. Similarities between the two proposed eddy-viscosity-improved resolvent analyses are also observed in higher-order modes. On the other hand, for the larger scales in the outer region, the plots in figures 12(d–f) and 13(d–f) show that the projection coefficients of the eddy-viscosity-improved modes at the diagonal line (i.e. $\beta _{ii}$) are all superior to those of the standard modes. Hence in these cases, higher-order resolvent modes can help to improve the capability of turbulence modelling/prediction, in particular in that the second SPOD modes capture more than 20 % of the total energy for large-scale motions (see figures 4b,d).

Figure 12. Projections of the first four response modes onto the corresponding SPOD modes, at (a–c) $(\lambda _1^+,\lambda _3^+,c)=(465,96,0.52)$ and (d–f) $(\lambda _1^+,\lambda _3^+,c)=(1210,450,0.89)$, in the adiabatic turbulent boundary layer.

Figure 13. Projections of the first four response modes onto the corresponding SPOD modes, at (a–c) $(\lambda _1^+,\lambda _3^+,c)=(929,208,0.66)$ and (d–f) $(\lambda _1^+,\lambda _3^+,c)=(1239,459,0.88)$, in the cold-wall turbulent boundary layer.

4. Summary

This paper develops two eddy-viscosity models to improve the low-rank approximation of resolvent analysis in compressible turbulent boundary layers, i.e. a modified Cess eddy-viscosity model and a new semi-empirical eddy-viscosity model. The models are represented with a two-layer algebraic expression, considering the composite-layer nature of the turbulent boundary layer. The first model is a correction of the classical Cess model, based on the compressibility transformation in the inner region and the mixing-length hypothesis in the outer region, such that the density variation and the outer-layer difference between the boundary layer and the channel flow are properly modelled. Considering the exclusive features of zero-pressure-gradient turbulent boundary layers, the second model is proposed based on an empirical fitting law of the total shear stress and the mixing-length theory.

In the a priori tests, both models are able to capture the tendency of the eddy-viscosity profile reasonably well; nonetheless, the semi-empirical model seems to outperform the corrected Cess model. With the DNS data of two hypersonic turbulent boundary layers, comparisons between the resolvent response modes with/without eddy viscosity and the SPOD modes are conducted in terms of the perturbations of velocities, density and temperature. Results show that the predictive performance of the optimal resolvent mode in turbulent motions is improved significantly by adding eddy viscosity in the resolvent operator, especially at the high-energy-containing scales, confirming the promising low-rank approximations of resolvent analysis with the proposed eddy-viscosity models. It is also noted that there is no significant difference observed between the predictions with the two models, suggesting that the optimal response mode is not sensitive to the relatively small differences in the eddy-viscosity quantities, as long as its distribution can be represented properly.