Eddy-viscosity-improved resolvent analysis of compressible turbulent boundary layers

Abstract

An improved resolvent analysis is proposed in the regime of compressible turbulent boundary layers. To better model nonlinear processes in the input, the resolvent framework is augmented by adding eddy viscosity. To this end, we propose two eddy-viscosity models: a modified Cess eddy-viscosity model coupling the compressibility transformation and outer-layer correction, and a new eddy-viscosity model based on an empirical relationship and mixing-length theory. Both are incorporated into the resolvent operator to examine the performance of the eddy-viscosity-improved resolvent-based reduced-order modelling. Results of the augmented resolvent analysis are compared qualitatively and quantitatively with the first leading mode of spectral proper orthogonal decomposition, by checking the profiles and cross-spectral densities of velocities, density and temperature in two hypersonic turbulent boundary layers under different wall conditions. Higher accuracy of the turbulence prediction is achieved by adding the proposed eddy-viscosity models, particularly for the energetic cycle in the outer-layer region where strong nonlinear energy transfer exists.

1. Introduction

Resolvent analysis is known as a promising operator-theoretic modal decomposition method in fluid dynamics (McKeon & Sharma 2010; Taira et al. 2017; Jovanović 2021). Based on linear dynamical system theory, the Navier–Stokes equations are linearized around the steady base flow to form an input–output framework, governed by the resolvent operator. Resolvent analysis has been used widely to identify the prominent linear mechanisms, predict the organization of coherent structures, and design optimal flow-control algorithms, for canonical wall-bounded turbulence (Sharma & McKeon 2013; Zare, Jovanović & Georgiou 2017; Martini et al. 2020, 2022; Karban et al. 2022), mixing jets (Schmidt et al. 2018; Towne, Schmidt & Colonius 2018; Lesshafft et al. 2019; Pickering et al. 2021), bluff-body wakes (Jin, Symon & Illingworth 2021; Jin, Illingworth & Sandberg 2022) and shock buffet (He & Timme 2020; Kojima et al. 2020), to name a few.

Performing singular value decomposition of the resolvent operator, the response to the endogenous stimulus is ranked at particular wavenumber and frequency pairs. By selecting the optimal response modes, the chaotic motions in wall-bounded turbulence can be characterized by the low-rank approximation of resolvent modes. For instance, Sharma & McKeon (2013) identified the complex coherent structures in turbulent pipe flow, based on a single optimal response mode or a linear superposition of modes at representative wavenumber combinations. Moarref et al. (2013a) exploited the optimal resolvent modes as basis functions to predict streamwise velocity fluctuations in high-Reynolds-number channel flows. A convex optimization problem was solved to determine the weight functions amplifying or attenuating the energy densities of the optimal modes, and good agreements with experimental or well-resolved numerical results were achieved. Similar optimization algorithms were also applied to estimate the distribution of spanwise and wall-normal velocity fluctuations, and the Reynolds shear stress across the wall layer (Moarref et al. 2013b, 2014). Abreu et al. (2020) used the resolvent analysis and spectral proper orthogonal decomposition (SPOD) to identify the dominant near-wall coherent structures in channel flows. Quantitative comparisons suggested that the resolvent and SPOD modes are consistent with each other, especially at the frequencies and wavenumbers where the lift-up mechanism is present. In the context of compressible wall turbulence, Dawson, McKeon & Saxton-Fox (2018) and Dawson & McKeon (2019) extended the resolvent formulation to a compressible form, by incorporating the advection–diffusion equation of temperature (or any scalar), allowing the prediction of passive scalar structures. Compressible resolvent analysis has been applied to supersonic/hypersonic wall-bounded turbulence, for various purposes, including revealing the features of energy-distribution mechanisms and examining the scaling laws of resolvent modes (Bae, Dawson & McKeon 2020a,b), identifying coherent structures under the condition of wall cooling (Fan, Li & Sandberg 2023), and analysing the mechanisms responsible for the energy amplification in subsonic and supersonic response modes (Madhusudanan & McKeon 2022; Chen et al. 2023).

Recently, various modifications to the standard resolvent formulation have been proposed to improve the accuracy of resolvent-based modelling and prediction (McKeon 2017; Gupta et al. 2021; Wu & He 2023). Among the modifications, a critical point is that the nonlinear process in the resolvent formulation needs to be modelled appropriately. Generally, a stochastic excitation is imposed on the linear system as an input, to avoid the calculation or modelling of the complex nonlinear terms. However, this strategy may lead to a poor performance of prediction. For instance, Symon, Illingworth & Marusic (2021) reported that given white noise inputs, the predicted production term in the turbulent kinetic energy budgets is much larger than the inter-scale transfer term, which conflicts with the real physics. They found that this unphysical phenomenon is associated directly with the non-normality in the resolvent operator, and it can be alleviated by adding appropriate eddy viscosity to the resolvent formulation, so as to account for part of the nonlinear process. In the resolvent analysis of incompressible turbulent channel flows, the Cess eddy-viscosity model (Cess 1958; Reynolds & Hussain 1972) is used widely. Significant improvements of accuracy were achieved, for the prediction of instantaneous velocity fields (Illingworth, Monty & Marusic 2018; Madhusudanan, Illingworth & Marusic 2019), energy transfer (Symon et al. 2021) and spatio-temporal power spectral densities (Morra et al. 2019, 2021; Symon, Illingworth & Marusic 2020), especially for the high-energy-containing scales. The turbulent dissipation introduced by the eddy viscosity plays an important role in modelling the nonlinear transfer from large- to small-scale motions (Symon et al. 2021, 2023). Aiming to model the remaining part of the nonlinear process, Gupta et al. (2021) decomposed the nonlinear terms into two parts, the eddy-viscosity term and the improved stochastic-forcing term. The two terms were modelled separately as functions of the wall distance and turbulence scales, to estimate the large-scale fluctuations at different wall-normal locations in channel flows.

So far, no eddy-viscosity model has been developed for the resolvent analysis in the case of compressible turbulent boundary layers, to the best of the authors’ knowledge. Due to the different geometries between the channel and boundary layer, the classical Cess model, without modifications, may not be suitable for turbulent boundary layers, especially for the dynamics in the outer region. On the other hand, in compressible turbulent boundary layers, the mean velocity profiles differ from their incompressible counterparts due to the variations of thermodynamic properties. The Cess eddy-viscosity model (Cess 1958; Reynolds & Hussain 1972) is defined such that the universal mean velocity profile in incompressible channel flows is obtained by simply integrating the expression ${Re}_\tau (1-x_2)\mu _w/\mu _T$ in the wall-normal direction. (In this expression, ${Re}_\tau$ is the friction Reynolds number, $\mu _w$ and $\mu _T$ are the wall dynamic and total viscosities, and $x_2$ is the normalized wall-normal distance to the surface.) However, with rapid density variations in the compressible cases, the velocity profiles deviate from those resulting from the classical Cess model. Therefore, compressibility effects need to be taken into account in the development of eddy-viscosity models for compressible flows.

In the present study, two eddy-viscosity models are proposed to improve the resolvent analysis for compressible turbulent boundary layers: (i) the classical Cess eddy-viscosity model is modified by coupling a compressibility transformation with an outer-layer correction; (ii) we develop a new eddy-viscosity model based on an empirical expression and the mixing-length theory. To examine their capability in resolvent-based modelling, we apply them to the resolvent formulation for two hypersonic turbulent boundary layers under different wall conditions. The remainder of the paper is outlined as follows. Section 2 introduces the resolvent formulation for the compressible turbulent boundary layers, with/without the addition of eddy viscosity. Two eddy-viscosity models are proposed herein, with a priori examinations conducted in contrast to the direct numerical simulations (DNS) data. In § 3, to evaluate the performance of eddy-viscosity-improved resolvent-based reduced modelling, we compare the predicted perturbations in velocities, density and temperature to those issued from the SPOD. Finally, concluding remarks are given in § 4.

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. 2020a), 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 1984),

(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. (2023).

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 (1965), 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. 2018; Symon et al. 2021; Holford, Lee & Hwang 2023). 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 2018). In this study, we set ${Pr}_t$ to 0.85, which is informed by previous studies (Zhang et al. 2014), and our preliminary tests showed that this value does not influence the prediction results, which was also reported by Chen et al. (2023). 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 (1965) 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 1958) 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. 2009; Amaral et al. 2021). 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 1971). 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 2000), 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 2000).

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 1951) (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 2004). For example, Escudier (1966) 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 (2006) 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. 2006). 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. 2006; Xia, Zhang & Yang 2021).

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 (2011) and Bernardini & Pirozzoli (2011). 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 (2022) and Fan & Li (2023). 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 2011); $M=3$, $Re_\tau \approx 500$ and $M=4$, $Re_\tau \approx 500$ (Bernardini & Pirozzoli 2011); and $M =5.86$, $Re_\tau \approx 420$ under the adiabatic and cold-wall conditions (Fan et al. 2022; Fan & Li 2023). 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 1967; Berkooz, Holmes & Lumley 1993; Towne et al. 2018) 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. 2017; Towne et al. 2018). 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. 2020; Pickering et al. 2021) 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 2023). 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. 2023), 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. (2023).

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 2010). 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 2009).

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. 2021). 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. 2023). 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 2022) 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 (1965), 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. 2020; Pickering et al. 2021). 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. 2018). 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 2019; Fan et al. 2022). Moreover, the additional eddy diffusion also plays a critical role in the positive energy receipt in the near-wall region (Symon et al. 2023). 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.