2.1 Flow conservation equations

The laws of conservation of mass, momentum and energy (known as the Navier–Stokes equations), in differential and dimensionless form, are

(2.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}u_{j})}{\unicode[STIX]{x2202}x_{j}}=0,\\ \displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}u_{i})}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}u_{i}u_{j}+p\unicode[STIX]{x1D6FF}_{ij}-\unicode[STIX]{x1D70F}_{ij})}{\unicode[STIX]{x2202}x_{j}}=0,\\ \displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}E)}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}Eu_{j}+pu_{j}+q_{j}-u_{i}\unicode[STIX]{x1D70F}_{ij})}{\unicode[STIX]{x2202}x_{j}}=0,\end{array}\right\}\end{eqnarray}$$

where $x_{i}=(x,y,z)$ are the coordinates in the streamwise, wall-normal and spanwise directions, $u_{i}=(u,v,w)$ are the corresponding velocity components, $t$ the time, $\unicode[STIX]{x1D70C}$ the fluid density, $E=e+u_{i}u_{i}/2$ the total energy, $e$ the internal energy and $p$ is the pressure. The viscous stress tensor, $\unicode[STIX]{x1D70F}_{ij}$ , and the heat flux vector, $q_{j}$ , are given by

(2.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D70F}_{ij}=\frac{\unicode[STIX]{x1D707}}{\mathit{Re }_{\infty }}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}\right)+\frac{\unicode[STIX]{x1D706}}{\mathit{Re }_{\infty }}\unicode[STIX]{x1D6FF}_{ij}\frac{\unicode[STIX]{x2202}u_{k}}{\unicode[STIX]{x2202}x_{k}},\\ \displaystyle q_{j}=-\frac{\unicode[STIX]{x1D705}}{\mathit{Re }_{\infty }\mathit{Pr}_{\infty }\mathit{Ec}_{\infty }}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}x_{j}}.\end{array}\right\}\end{eqnarray}$$

Here $\unicode[STIX]{x1D707}$ is the dynamic viscosity, $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D707}_{b}-2/3\unicode[STIX]{x1D707}$ the second viscosity, $\unicode[STIX]{x1D707}_{b}$ the bulk viscosity and $\unicode[STIX]{x1D705}$ is the thermal conductivity. It is known that $\unicode[STIX]{x1D707}_{b}$ has a very limited effect on the linear stability of channel flows (Ren et al. Reference Ren, Fu and Pecnik2019). Results presented in the following sections are subjected to $\unicode[STIX]{x1D707}_{b}=0$ . An additional assumption is that buoyancy effects are not considered. The equations above have been non-dimensionalized by reference values, as follows

(2.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}u={\displaystyle \frac{u^{\ast }}{u_{\infty }^{\ast }}},\quad x_{i}={\displaystyle \frac{x_{i}^{\ast }}{l_{0}^{\ast }}},\quad t={\displaystyle \frac{t^{\ast }u_{\infty }^{\ast }}{l_{0}^{\ast }}},\quad p={\displaystyle \frac{p^{\ast }}{\unicode[STIX]{x1D70C}_{\infty }^{\ast }u_{\infty }^{\ast 2}}},\quad \unicode[STIX]{x1D70C}={\displaystyle \frac{\unicode[STIX]{x1D70C}^{\ast }}{\unicode[STIX]{x1D70C}_{\infty }^{\ast }}},\\ T={\displaystyle \frac{T^{\ast }}{T_{\infty }^{\ast }}},\quad E={\displaystyle \frac{E^{\ast }}{u_{\infty }^{\ast 2}}},\quad \unicode[STIX]{x1D707}={\displaystyle \frac{\unicode[STIX]{x1D707}^{\ast }}{\unicode[STIX]{x1D707}_{\infty }^{\ast }}},\quad \unicode[STIX]{x1D705}={\displaystyle \frac{\unicode[STIX]{x1D705}^{\ast }}{\unicode[STIX]{x1D705}_{\infty }^{\ast }}},\end{array}\right\}\end{eqnarray}$$

which leads to the definition of the Reynolds number, $\mathit{Re}_{\infty }$ , Prandtl number, $\mathit{Pr}_{\infty }$ , Eckert number, $\mathit{Ec}_{\infty }$ , and the Mach number, $\mathit{Ma}_{\infty }$ (all based on free-stream parameters)

(2.4a-d ) $$\begin{eqnarray}\mathit{Re}_{\infty }={\displaystyle \frac{\unicode[STIX]{x1D70C}_{\infty }^{\ast }u_{\infty }^{\ast }l_{0}^{\ast }}{\unicode[STIX]{x1D707}_{\infty }^{\ast }}},\quad \mathit{Pr}_{\infty }={\displaystyle \frac{\unicode[STIX]{x1D707}_{\infty }^{\ast }C_{p\infty }^{\ast }}{\unicode[STIX]{x1D705}_{\infty }^{\ast }}},\quad \mathit{Ec}_{\infty }={\displaystyle \frac{u_{\infty }^{\ast 2}}{C_{p\infty }^{\ast }T_{\infty }^{\ast }}},\quad \mathit{Ma}_{\infty }={\displaystyle \frac{u_{\infty }^{\ast }}{a_{\infty }^{\ast }}}.\end{eqnarray}$$

The subscript $\infty$ denotes free-stream values, superscript $\ast$ stands for dimensional variables, $l_{0}^{\ast }$ is a chosen length scale, $a_{\infty }^{\ast }$ is the speed of sound in the free stream. Note that for an ideal gas $\mathit{Ec}_{\infty }=(\unicode[STIX]{x1D6FE}-1)\mathit{Ma}_{\infty }^{2}$ , where $\unicode[STIX]{x1D6FE}$ is the heat capacity ratio. In linear stability theory, $l_{0}^{\ast }$ is chosen to be the local boundary-layer thickness scale $\unicode[STIX]{x1D6FF}^{\ast }$ , which results in the definition of $\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ ,

(2.5a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}^{\ast }=\left({\displaystyle \frac{\unicode[STIX]{x1D707}_{\infty }^{\ast }x^{\ast }}{\unicode[STIX]{x1D70C}_{\infty }^{\ast }u_{\infty }^{\ast }}}\right)^{1/2},\quad \mathit{Re}_{\unicode[STIX]{x1D6FF}}={\displaystyle \frac{\unicode[STIX]{x1D70C}_{\infty }^{\ast }u_{\infty }^{\ast }\unicode[STIX]{x1D6FF}^{\ast }}{\unicode[STIX]{x1D707}_{\infty }^{\ast }}}=\left({\displaystyle \frac{\unicode[STIX]{x1D70C}_{\infty }^{\ast }u_{\infty }^{\ast }x^{\ast }}{\unicode[STIX]{x1D707}_{\infty }^{\ast }}}\right)^{1/2}.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D6FF}^{\ast }$ measures the order of boundary-layer thickness over a flat plate. For example, applying the Blasius solution, the displacement thickness $\unicode[STIX]{x1D6FF}_{1}^{\ast }\approx 1.721\unicode[STIX]{x1D6FF}^{\ast }$ , momentum thickness $\unicode[STIX]{x1D6FF}_{2}^{\ast }\approx 0.664\unicode[STIX]{x1D6FF}^{\ast }$ (Schlichting & Gersten Reference Schlichting and Gersten2017).

2.2 The laminar base flow

The self-similar solution to the boundary-layer equation over a flat plate is employed in this work. It serves as the base flow for the stability analysis as well as for the initial state of the DNS. Applying the boundary-layer assumption and a routine order-of-magnitude analysis of the dimensional Navier–Stokes (N–S) equations, the boundary-layer equations read,

(2.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}^{\ast }u^{\ast })}{\unicode[STIX]{x2202}x^{\ast }}}+{\displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}^{\ast }v^{\ast })}{\unicode[STIX]{x2202}y^{\ast }}}=0,\\ \unicode[STIX]{x1D70C}^{\ast }u^{\ast }{\displaystyle \frac{\unicode[STIX]{x2202}u^{\ast }}{\unicode[STIX]{x2202}x^{\ast }}}+\unicode[STIX]{x1D70C}^{\ast }v^{\ast }{\displaystyle \frac{\unicode[STIX]{x2202}u^{\ast }}{\unicode[STIX]{x2202}y^{\ast }}}+{\displaystyle \frac{\text{d}p_{\infty }^{\ast }}{\text{d}x^{\ast }}}-{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y^{\ast }}}\left(\unicode[STIX]{x1D707}^{\ast }{\displaystyle \frac{\unicode[STIX]{x2202}u^{\ast }}{\unicode[STIX]{x2202}y^{\ast }}}\right)=0,\\ \unicode[STIX]{x1D70C}^{\ast }u^{\ast }{\displaystyle \frac{\unicode[STIX]{x2202}h^{\ast }}{\unicode[STIX]{x2202}x^{\ast }}}+\unicode[STIX]{x1D70C}^{\ast }v^{\ast }{\displaystyle \frac{\unicode[STIX]{x2202}h^{\ast }}{\unicode[STIX]{x2202}y^{\ast }}}-u^{\ast }{\displaystyle \frac{\text{d}p_{\infty }^{\ast }}{\text{d}x^{\ast }}}-{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y^{\ast }}}\left(\unicode[STIX]{x1D705}^{\ast }{\displaystyle \frac{\unicode[STIX]{x2202}T^{\ast }}{\unicode[STIX]{x2202}y^{\ast }}}\right)-\unicode[STIX]{x1D707}^{\ast }\left({\displaystyle \frac{\unicode[STIX]{x2202}u^{\ast }}{\unicode[STIX]{x2202}y^{\ast }}}\right)^{2}=0.\end{array}\right\}\end{eqnarray}$$

Introducing the Lees–Dorodnitsyn transformation (see introduction in Anderson Jr Reference Anderson2000; Schlichting & Gersten Reference Schlichting and Gersten2017)

(2.7) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\text{d}\unicode[STIX]{x1D709}=\unicode[STIX]{x1D70C}_{\infty }^{\ast }\unicode[STIX]{x1D707}_{\infty }^{\ast }u_{\infty }^{\ast }\,\text{d}x^{\ast },\\ \text{d}\unicode[STIX]{x1D702}={\displaystyle \frac{\unicode[STIX]{x1D70C}^{\ast }u_{\infty }^{\ast }}{\sqrt{2\unicode[STIX]{x1D709}}}}\,\text{d}y^{\ast },\end{array}\right\}\end{eqnarray}$$

for the boundary-layer equations, yields the transformed ordinary differential equations (ODE) for $f$ and $g$ given as,

(2.8) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\displaystyle \frac{\text{d}}{\text{d}\unicode[STIX]{x1D702}}}\left(C_{l}{\displaystyle \frac{\text{d}^{2}f}{\text{d}\unicode[STIX]{x1D702}^{2}}}\right)+f{\displaystyle \frac{\text{d}^{2}f}{\text{d}\unicode[STIX]{x1D702}^{2}}}=0,\\ {\displaystyle \frac{\text{d}}{\text{d}\unicode[STIX]{x1D702}}}\left({\displaystyle \frac{C_{l}}{\mathit{Pr}_{l}}}{\displaystyle \frac{\text{d}g}{\text{d}\unicode[STIX]{x1D702}}}\right)+f{\displaystyle \frac{\text{d}g}{\text{d}\unicode[STIX]{x1D702}}}+C_{l}{\displaystyle \frac{u_{\infty }^{\ast 2}}{h_{t\infty }^{\ast }}}\left({\displaystyle \frac{\text{d}^{2}f}{\text{d}\unicode[STIX]{x1D702}^{2}}}\right)^{2}=0,\end{array}\right\}\end{eqnarray}$$

where

(2.9a-d ) $$\begin{eqnarray}{\displaystyle \frac{\text{d}f}{\text{d}\unicode[STIX]{x1D702}}}={\displaystyle \frac{u^{\ast }}{u_{\infty }^{\ast }}},\quad g={\displaystyle \frac{h_{t}^{\ast }}{h_{t\infty }^{\ast }}},\quad C_{l}={\displaystyle \frac{\unicode[STIX]{x1D70C}^{\ast }\unicode[STIX]{x1D707}^{\ast }}{\unicode[STIX]{x1D70C}_{\infty }^{\ast }\unicode[STIX]{x1D707}_{\infty }^{\ast }}},\quad \mathit{Pr}_{l}={\displaystyle \frac{\unicode[STIX]{x1D707}^{\ast }C_{p}^{\ast }}{\unicode[STIX]{x1D705}^{\ast }}}.\end{eqnarray}$$

Here, $f$ and $g$ are unary functions of the transformed coordinate $\unicode[STIX]{x1D702}$ , while $h_{t}$ denotes the total enthalpy. The above ODEs are numerically integrated using the fourth-order Runge–Kutta scheme subjected to adiabatic wall boundary conditions. During the integration process, a one-dimensional (1-D) look-up table (see appendix A) is used to calculate $C_{l}$ and $\mathit{Pr}_{l}$ from the static temperature that can be obtained from the total enthalpy  $g$ .

2.3 Linear stability theory

The linear stability equations for the non-ideal gas have been derived and documented in Ren et al. (Reference Ren, Fu and Pecnik2019). It is known that for simple compressible systems (e.g. single phase, pure substances and uniform mixtures of non-reacting gases), the thermodynamic state is defined by two independent thermodynamic properties. We choose $\unicode[STIX]{x1D70C}$ and $T$ as the two independent thermodynamic quantities, while the remaining thermodynamic and transport properties (e.g. $E$ , $p$ , $\unicode[STIX]{x1D707}$ , $\unicode[STIX]{x1D705}$ ) are determined as functions of $\unicode[STIX]{x1D70C}$ and $T$ . For example, the viscosity perturbation is given by the two-dimensional Taylor expansion,

(2.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}^{\prime } & = & \displaystyle \left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}}\right|_{T}\unicode[STIX]{x1D70C}^{\prime }+\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}T_{0}}}\right|_{\unicode[STIX]{x1D70C}}T^{\prime }+{\displaystyle \frac{1}{2}}\left(\!\left.{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}^{2}}}\right|_{T}\unicode[STIX]{x1D70C}^{\prime }\unicode[STIX]{x1D70C}^{\prime }+2\left.{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}T_{0}}}\right|_{\unicode[STIX]{x1D70C}}\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}}\right|_{T}\unicode[STIX]{x1D70C}^{\prime }T^{\prime }+\left.{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}T_{0}^{2}}}\right|_{\unicode[STIX]{x1D70C}}T^{\prime }T^{\prime }\right)\nonumber\\ \displaystyle & & \displaystyle +\,\cdots \,.\end{eqnarray}$$

It can be seen that in non-ideal-gas flows, the viscosity perturbation is dependent on $\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}$ (at constant $T$ ), which is not accounted for in the conventional empirical viscosity laws of an ideal gas (e.g. Sutherland’s law and power law). The full stability equation is derived by introducing small perturbations into the N–S equations (2.1) and subtracting the governing equations of the base flow. With the nonlinear terms neglected, the linear stability equations are formulated as

(2.11) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D647}_{\unicode[STIX]{x1D669}}{\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{q}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D647}_{\unicode[STIX]{x1D66D}}{\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{q}}{\unicode[STIX]{x2202}x}}+\unicode[STIX]{x1D647}_{\unicode[STIX]{x1D66E}}{\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{q}}{\unicode[STIX]{x2202}y}}+\unicode[STIX]{x1D647}_{\unicode[STIX]{x1D66F}}{\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{q}}{\unicode[STIX]{x2202}z}}+\unicode[STIX]{x1D647}_{\unicode[STIX]{x1D666}}\boldsymbol{q}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D66D}\unicode[STIX]{x1D66D}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\boldsymbol{q}}{\unicode[STIX]{x2202}x^{2}}}+\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D66D}\unicode[STIX]{x1D66E}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\boldsymbol{q}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}y}}+\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D66D}\unicode[STIX]{x1D66F}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\boldsymbol{q}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z}}+\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D66E}\unicode[STIX]{x1D66E}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\boldsymbol{q}}{\unicode[STIX]{x2202}y^{2}}}+\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D66E}\unicode[STIX]{x1D66F}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\boldsymbol{q}}{\unicode[STIX]{x2202}y\unicode[STIX]{x2202}z}}+\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D66F}\unicode[STIX]{x1D66F}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\boldsymbol{q}}{\unicode[STIX]{x2202}z^{2}}}=0.\qquad\end{eqnarray}$$

Here $\boldsymbol{q}=(\unicode[STIX]{x1D70C}^{\prime },u^{\prime },v^{\prime },w^{\prime },T^{\prime })^{\text{T}}$ is the perturbation vector and the detailed expressions for the matrices $\unicode[STIX]{x1D647}_{\unicode[STIX]{x1D669}}$ , $\unicode[STIX]{x1D647}_{\unicode[STIX]{x1D66D}}$ , $\unicode[STIX]{x1D647}_{\unicode[STIX]{x1D66E}}$ , $\unicode[STIX]{x1D647}_{\unicode[STIX]{x1D66F}}$ , $\unicode[STIX]{x1D647}_{\unicode[STIX]{x1D666}}$ , $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D66D}\unicode[STIX]{x1D66D}}$ , $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D66E}\unicode[STIX]{x1D66E}}$ , $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D66F}\unicode[STIX]{x1D66F}}$ , $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D66D}\unicode[STIX]{x1D66E}}$ , $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D66E}\unicode[STIX]{x1D66F}}$ and $\unicode[STIX]{x1D651}_{\unicode[STIX]{x1D66D}\unicode[STIX]{x1D66F}}$ are functions of the dimensionless parameter (2.4), the base flow and thermodynamic and transport properties, and their detailed expressions can be found in Ren et al. (Reference Ren, Fu and Pecnik2019). The perturbation is assumed to have the normal-mode form,

(2.12) $$\begin{eqnarray}\boldsymbol{q}(x,y,z,t)=\hat{\boldsymbol{q}}(y)\exp (\text{i}\unicode[STIX]{x1D6FC}x+\text{i}\unicode[STIX]{x1D6FD}z-\text{i}\unicode[STIX]{x1D714}t)+\text{c.c.}\end{eqnarray}$$

where c.c. stands for the complex conjugate. Substituting (2.12) into (2.11) results in an eigenvalue problem. For boundary layers, we consider the spatial problem, where $\unicode[STIX]{x1D714}$ and $\unicode[STIX]{x1D6FD}$ are prescribed frequency and spanwise wavenumber. The real and imaginary parts of the eigenvalue $\unicode[STIX]{x1D6FC}$ give the streamwise wavenumber and its local growth rate, respectively.

To solve the eigenvalue problem, Chebyshev collocation points and Chebyshev differentiation matrices are introduced to discretize the equations. The perturbations are subjected to the following boundary conditions: $u^{\prime }=v^{\prime }=w^{\prime }=0$ and $\unicode[STIX]{x2202}T^{\prime }/\unicode[STIX]{x2202}y=0$ at the wall ( $y=0$ ), while in the free stream at $y=y_{max}$ , $u^{\prime }=v^{\prime }=w^{\prime }=0$ and $T^{\prime }=0$ .

2.4 Direct numerical simulation

In the context of direct numerical simulations, the N–S equations (2.1) are numerically integrated. Wall blowing and suction are introduced to excite the Tollmien–Schlichting (T–S) waves, which are the normal-mode solutions of the linearized N–S equations (Schmid & Henningson Reference Schmid and Henningson2001). The amplitude of the forcing has been properly chosen (between $10^{-5}$ and $10^{-3}$ of the free-stream Mach number) such that the excited perturbations stay in the linear regime. At the wall, no-slip and adiabatic boundary conditions are applied. At the inflow the self-similar solution obtained in § 2.2 is prescribed. Close to the outflow and the free-stream boundaries, a sponge region forces the solution towards the corresponding laminar state. The algorithm is based on a sixth-order compact finite-difference method with a staggered arrangement of flow variables, while the time stepping scheme is based on an explicit third-order Runge–Kutta time method. The DNS has been performed for 2-D perturbations in this study, a $N_{x}\times N_{y}=1201\times 201$ mesh has been used for most cases and is sufficient to give a grid-independent result. The non-ideal fluid properties are incorporated using 2-D look-up tables (see appendix A). Readers may refer to Nagarajan, Lele & Ferziger (Reference Nagarajan, Lele and Ferziger2003) and Marxen et al. (Reference Marxen, Magin, Iaccarino and Shaqfeh2011, Reference Marxen, Magin, Shaqfeh and Iaccarino2013) for numerical details.

To analyse the results from DNS and to compare them with LST, the results are Fourier transformed in time with a fundamental angular frequency of $\unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D714}/2$ ,

(2.13) $$\begin{eqnarray}\hat{q}_{k}(x,y,z)={\displaystyle \frac{1}{N}}\mathop{\sum }_{l=1}^{N}q(x,y,z,t_{l})\exp (2\text{i}k\unicode[STIX]{x1D6FA}t_{l}).\end{eqnarray}$$

We take $N=50$ samples within two forcing periods, $l$ is the sampling index and the discrete time is $t_{l}=2\unicode[STIX]{x03C0}l/(\unicode[STIX]{x1D6FA}N)$ ; $k=0,1/2,1,3/2$ , etc. such that $\hat{q}_{0}$ , $\hat{q}_{1/2}$ , $\hat{q}_{1}$ give the base flow, subharmonic mode and fundamental mode, respectively. The growth rate and phase velocity of the fundamental mode are given by

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}_{i}(x)=-{\displaystyle \frac{\mathit{Re}_{\unicode[STIX]{x1D6FF}}}{\mathit{Re }_{\infty }}}{\displaystyle \frac{1}{\hat{q}_{1}^{max}}}{\displaystyle \frac{\unicode[STIX]{x2202}\hat{q}_{1}^{max}}{\unicode[STIX]{x2202}x}}, & \displaystyle\end{eqnarray}$$

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle c(x)=\mathit{Re}_{\infty }F\left({\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}x}}\right)^{-1}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D719}_{k}=\arg \left(\hat{q}_{k}\right)$ is the phase angle of the perturbation and $F$ is the dimensionless frequency defined by (4.1).

The time integration is performed until the flow has reached a periodic solution. This is ensured by inspecting the subharmonic of the perturbation, whose amplitude is then at least one order of magnitude less than the fundamental perturbation.

4.1 The supercritical regime

Figure 3 shows the distribution of the local growth rate $(-\unicode[STIX]{x1D6FC}_{i})$ in the $F{-}\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ diagram. The contour levels start from the neutral curve ( $\unicode[STIX]{x1D6FC}_{i}=0$ ), up to the largest growth rate in the range of $\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ and $F$ considered. In each panel the canvas is divided into four quadrants with mirrored coordinates to help compare results of four Eckert numbers. The dotted lines represent the neutral curve for $\mathit{Ec}_{\infty }=0.05$ , which are also plotted in the quadrants with the higher Eckert numbers to highlight the compressibility effects. Among the results, figure 3(d) provides a reference for the ideal-gas regime with $T_{\infty }^{\ast }=800~\text{K}$ , where it can be seen that an increase in Eckert (Mach) number stabilizes the flow, as the growth rate reduces and the extent of the neutral curve decreases. This is expected and has been documented in the past by Mack (Reference Mack1984).

Figure 3.

Growth rates ( $-\unicode[STIX]{x1D6FC}_{i}$ ) of 2-D perturbations in the $F{-}\mathit{Re}_{\infty }$ stability diagram with supercritical free-stream temperatures: (a) $T_{\infty }^{\ast }=320~\text{K}$ , (b) $T_{\infty }^{\ast }=340~\text{K}$ , (c) $T_{\infty }^{\ast }=360~\text{K}$ and (d) $T_{\infty }^{\ast }=800~\text{K}$ . The coordinate in each panel has been mirrored in the centre in order to compare results for the four Eckert numbers.

Figure 4.

Phase velocities of perturbations (coloured contours) in the $F{-}\mathit{Re}$ stability diagram with supercritical free-stream temperatures: (a) $T_{\infty }^{\ast }=320~\text{K}$ , (b) $T_{\infty }^{\ast }=340~\text{K}$ , (c) $T_{\infty }^{\ast }=360~\text{K}$ and (d) $T_{\infty }^{\ast }=800~\text{K}$ . The solid lines indicate neutral curves of the corresponding cases.

However, if the free-stream temperature decreases towards values slightly above the pseudo-critical temperature, the non-ideal-gas effects considerably increase the stability of the base flow for high Eckert numbers. For example, the neutral curve for case T320E4 has become smaller, indicating a very narrow band of unstable frequencies. For low Eckert numbers (e.g. $\mathit{Ec}_{\infty }=0.05$ ), the size of the neutral curve is comparable, but the growth rate decreases with decreasing free-stream temperature.

Note that the results presented in figure 3 can be interpreted either by specifying $\mathit{Ec}_{\infty }$ and comparing results for different $T_{\infty }^{\ast }$ , or conversely setting $T_{\infty }^{\ast }$ and comparing different $\mathit{Ec}_{\infty }$ . Both show a stabilization of the base flow through non-ideal-gas (varying $T_{\infty }^{\ast }$ ) and compressibility effects (varying $\mathit{Ec}_{\infty }$ ). One may argue that $\mathit{Ma}_{\infty }$ can be used to show the compressibility effects. Here we clarify the advantage of using $\mathit{Ec}_{\infty }$ instead of $\mathit{Ma}_{\infty }$ . From table 1, it is clear that a higher $\mathit{Ec}_{\infty }$ corresponds to a higher $\mathit{Ma}_{\infty }$ , therefore the stabilization due to compressibility effects can be characterized by $\mathit{Ma}_{\infty }$ . However, $\mathit{Ma}_{\infty }$ is not a suitable parameter to measure the degree of compressibility when non-ideal-gas effects are present. In fact, the speed of sound for the non-ideal gas drops sharply near the critical point. For example, for $\text{CO}_{2}$ at 80 bar the speed of sound $a^{\ast }$ drops from $987.5$ to $178.9~\text{m}~\text{s}^{-1}$ and increases again to $433.3~\text{m}~\text{s}^{-1}$ when temperature increases from 220 K, to the pseudo-critical value of 307.7 K and further to 800 K, respectively. As a result, non-ideal-gas effects cannot be correctly quantified at fixed  $\mathit{Ma}_{\infty }$ .

The phase velocity ( $c=\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FC}_{r}$ ) shown in figure 4 is another important physical variable that characterizes the perturbation. The corresponding neutral curves are shown to indicate the region of instability. It is clear that the phase velocities within the neutral curves remain between 0.2 and 0.6. The distributions of $c$ in the $F{-}\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ diagram are similar for all the cases, i.e. an increase in $\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ or $F$ leads to a larger  $c$ . Non-ideal-gas effects, as well as compressibility effects, are seen to both increase the phase velocity. In accordance with an actual experimental/engineering set-up, the growth rate and phase velocity can be readily transferred into a physical $x^{\ast }-F^{\ast }$ diagram using (2.5) and (4.1).

We show profiles of perturbations in figure 5 for $T_{\infty }^{\ast }=320~\text{K}$ and $T_{\infty }^{\ast }=360$ . The results are also compared for two Eckert numbers of 0.05 and 0.20. The profiles correspond to the largest growth rate in the $F{-}\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ diagram shown in figure 3, and they are normalized by the amplitude of the streamwise velocity perturbation $\hat{u}$ . As shown in the figure, $\hat{u}$ and $\hat{v}$ are similar in all four cases. Perturbations are dominated by $\hat{u}$ except the T320E4 case, in which $\hat{\unicode[STIX]{x1D70C}}$ dominates the entire perturbation fields. Figure 5 demonstrates that compressibility (increase $\mathit{Ec}_{\infty }$ ) and non-ideal-gas effects (reduce $T_{\infty }^{\ast }$ ) both increase the amplitude of density perturbations.

Figure 5.

Profiles of the most amplified perturbations in the stability diagram of figure 4: (a) $T_{\infty }^{\ast }=320~\text{K}$ , $\mathit{Ec}_{\infty }=0.05$ , (b) $T_{\infty }^{\ast }=320~\text{K}$ , $\mathit{Ec}_{\infty }=0.20$ , (c) $T_{\infty }^{\ast }=360~\text{K}$ , $\mathit{Ec}_{\infty }=0.05$ and (d) $T_{\infty }^{\ast }=360~\text{K}$ , $\mathit{Ec}_{\infty }=0.20$ .

Figure 6.

Balance of the continuity equation with the same parameters as in figure 5. Legend shows the four terms in (4.2).

This can be explained by the balance of the linearized continuity equation,

(4.2) $$\begin{eqnarray}\underbrace{(\text{i}\unicode[STIX]{x1D6FC}u_{0}-\text{i}\unicode[STIX]{x1D714})\hat{\unicode[STIX]{x1D70C}}}_{1}+\underbrace{\text{i}\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{0}\hat{u} }_{2}+\underbrace{{\displaystyle \frac{\text{d}\unicode[STIX]{x1D70C}_{0}}{\text{d}y}}\hat{v}}_{3}+\underbrace{\unicode[STIX]{x1D70C}_{0}{\displaystyle \frac{\text{d}\hat{v}}{\text{d}y}}}_{4}=0.\end{eqnarray}$$

When non-ideal and compressibility effects are insignificant, see figure 5(c), the density gradient of the base flow is weak and, consequently, $\hat{\unicode[STIX]{x1D70C}}$ is small in amplitude. Therefore, equation (4.2) is mainly balanced by term 2 and term 4, as shown in figure 6(c). Either reducing $T_{\infty }^{\ast }$ or increasing $\mathit{Ec}_{\infty }$ results in an increase of $\text{d}\unicode[STIX]{x1D70C}_{0}/\text{d}y$ . This increase is fairly significant when $T_{\infty }^{\ast }$ is close to the pseudo-critical point (see figure 2). The increase of term 3 has to be balanced with the help of term 1, namely, a comparable density perturbation.

To account for possible non-parallel base flow effects and to validate the linear stability theory, DNS are performed for the T360 and T320 cases (the validation for the subcritical cases follows in the next section).

Figure 7.

DNS of the T360E1 case. Contour lines of the wall-normal velocity in five coloured regions show the laminar flow (1), receptivity stage (2), modal decay before branch-I of the neutral curve (3), followed by the modal growth (4) and modal decay (5) after branch-II of the neutral curve. The up/down arrows at $x=4$ show the introduced wall blowing/suction. A movie of the perturbation development is available as a supplementary file (movie 1) at https://doi.org/10.1017/jfm.2019.348.

Figure 7 shows an example of the DNS results for the T360E1 case. Results are displayed by the contour plot of the wall-normal velocity, which nicely shows the development of the perturbation. Five physical regions can be identified: before the sinusoidal wall blowing/suction is introduced at $x=4$ to excite the T–S wave, the flow remains laminar in region 1; the receptivity process takes places in region 2, where the forced external perturbation excites the T–S waves; the perturbation is well formed in the boundary layer and starts to follow the LST prediction in region 3. This region is yet ahead of branch-I of the neutral curve, therefore, perturbations goes through modal decay until region 4; in regions 4 and 5, perturbations follow the linear stability theory and goes through a modal growth and decay.

Figure 8.

DNS validation of the T320 (a1–a3) and T360 (b1–b3) cases. (a1) and (b1) show the neutral curve in the $F{-}x$ diagram. The blue solid line indicates the frequency of wall blowing/suction introduced to excite the T–S wave. (a2), (b2), (a3) and (b3) provide comparisons of the growth rate and phase velocity between DNS and LST. The arrows in (a2) and (b2) indicate the position where wall blowing/suction is introduced.

To quantitatively compare LST and DNS, figure 8(a1,b1) shows the neutral curve in the $F{-}x$ diagram. The arrow in figure 8(a2,b2) indicates the location where wall blowing/suction is introduced (upstream of branch-I of the neutral curves). The amplitude of the forced blowing/suction has been kept small ( $O(10^{-4})$ of the Mach number based on wall-normal momentum flux) so that the development of the perturbation stays in the linear regime. The frequencies $F=21\times 10^{-6}$ and $31\times 10^{-6}$ are chosen for the T320 and T360 cases respectively, which cut through the neutral curves (see figure 8 a1,b1). One would expect modal growth of the perturbation between the two branches of each neutral curve.

The DNS is thoroughly compared with LST in figure 8(a2,b2,a3,b3) in terms of the growth rate and phase velocity. A local average is applied to the DNS results to remove numerical oscillations. For the T360 cases (b1–b3), the DNS and LST match very well. Small differences can be seen for the high Eckert number cases, which can be attributed to the lower growth rate and the ensuing lower amplitude of the perturbation. To calculate the growth rate (2.14), a lower amplitude of perturbations can cause a relatively larger error.

On the other hand, the flow is considerably stabilized at $T_{\infty }^{\ast }=320~\text{K}$ (figure 8 a1–a3). As can be inferred from the neutral curve (figure 8 a1), the unstable region in $x$ is substantially reduced. Figure 8(a2,a3) shows that the DNS and LST match well at $\mathit{Ec}_{\infty }=0.05$ . The difference becomes larger with increasing Eckert number. In the two cases of $\mathit{Ec}_{\infty }=0.15$ and 0.20, the growth rate predicted by LST is very small, which is not completely captured by DNS. Note that we have shortened the domain of the DNS for the $\mathit{Ec}_{\infty }=0.15$ and 0.20 cases, the number of grid points both in $x$ and $y$ have been refined to make sure the results are grid independent. Differences in the phase velocity at $\mathit{Ec}_{\infty }=0.20$ also indicate that the modal instability predicted by LST (though the growth rate is small but positive) cannot be well captured in the DNS simulations. One of the possible reasons lies in the fact that in DNS the perturbed flow field is a result of multiple eigenmodes. Therefore, when the growth rate of the single unstable mode is exceedingly small, the flow field can be co-dominated by multiple decaying modes, and the post-processed growth rate and phase velocity from DNS cannot match a single mode from LST.

The above results have confirmed that the non-ideal-gas effects stabilize the flow in the supercritical regime. The stabilization is more significant when the free-stream temperature is closer to the pseudo-critical value and/or when coupled with strong compressibility effects.

4.2 The subcritical and transcritical regimes

Figure 9.

Growth rates of perturbations in the $F{-}\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ stability diagram with subcritical free-stream temperatures: (a) $T_{\infty }^{\ast }=220~\text{K}$ , (b) $T_{\infty }^{\ast }=240~\text{K}$ , (c) $T_{\infty }^{\ast }=260~\text{K}$ and (d) $T_{\infty }^{\ast }=280~\text{K}$ .

Figure 10.

Profiles of the most amplified perturbations in the stability diagram of figure 9: (a) $T_{\infty }^{\ast }=240~\text{K}$ , $\mathit{Ec}_{\infty }=0.05$ , (b) $T_{\infty }^{\ast }=240~\text{K}$ , $\mathit{Ec}_{\infty }=0.20$ , (c) $T_{\infty }^{\ast }=280~\text{K}$ , $\mathit{Ec}_{\infty }=0.05$ and (d) $T_{\infty }^{\ast }=280~\text{K}$ , $\mathit{Ec}_{\infty }=0.15$ .

In this section, we discuss cases with subcritical free-stream temperatures (T220, T240, T260 and T280). Recall the base flow in § 3, the laminar flow stays in the subcritical regime except for the transcritical case with $T_{\infty }^{\ast }=280~\text{K}$ and $\mathit{Ec}_{\infty }$ sufficiently large (T280E4).

The growth rates for these cases are shown in figure 9. Comparing the four quadrants in each panel, it can be seen that the flow is stabilized by increasing $\mathit{Ec}_{\infty }$ in the subcritical regime. Similar to the supercritical regime, non-ideal-gas effects (increase in $T_{\infty }^{\ast }$ ) further stabilize the subcritical flows. In the subcritical regime, the perturbation profiles for cases T240E1, T240E4, T280E1 and T280E3 are provided in figure 10. Recalling the discussion for the supercritical regime, the profiles in the subcritical regime follow a similar trend: the compressibility and non-ideal-gas effects increase the amplitude of the density perturbation due to the increase of the base flow’s density gradient $\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}/\unicode[STIX]{x2202}y$ and the balancing mechanism of the linearized continuity equation. To avoid repetition, the quantitative DNS validation of the stabilization of non-ideal-gas effects in the subcritical regime has been provided in appendix C.

More attention shall be paid to the T280 cases in figure 9(d). The stability of the base flow noticeably increases with growing Eckert number. For instance, at $\mathit{Ec}_{\infty }=0.15$ no modal instability can be seen for Reynolds numbers up to $\mathit{Re}_{\unicode[STIX]{x1D6FF}}=2000$ . After the base flow crosses the pseudo-critical point (at $\mathit{Ec}_{\infty }=0.20$ ), a co-existence of dual unstable modes is observed, which we term Mode I and Mode II. The neutral curve of Mode I has a similar shape to the subcritical regime, while Mode II has a much larger growth rate and unstable band. For this case, a typical eigenspectrum is provided in figure 11 with $\mathit{Re}_{\infty }=1500$ and $F=45\times 10^{-6}$ . As it can be clearly inferred, both modes are unstable with the eigenvalue of $\unicode[STIX]{x1D6FC}=0.1338{-}0.001293i$ (Mode I) and $\unicode[STIX]{x1D6FC}=0.2037{-}0.001887i$ (Mode II). As such, it can be concluded that with increasing Eckert number, the non-ideal-gas effects initially stabilize Mode I. However, after the base flow becomes transcritical, Mode II emerges and finally dominates the instability.

Figure 11.

Eigenspectrum of the transcritical case with $T_{\infty }^{\ast }=280~\text{K}$ , $\mathit{Ec}_{\infty }=0.20$ , $\mathit{Re}_{\unicode[STIX]{x1D6FF}}=1500$ and $F=45\times 10^{-6}$ . The dual modes are highlighted with arrows.

To further explore the transcritical regime for the T280 cases, we show in figures 12 and 13 the detailed evolution of the growth rate and phase velocity on gradually increasing the Eckert number from $\mathit{Ec}_{\infty }=0.11$ to $\mathit{Ec}_{\infty }=0.202$ . In order to show the compelling stabilizing effect as observed in figure 9(d) (at $\mathit{Ec}_{\infty }=0.15$ ), the range of $\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ has been extended to $\mathit{Re}_{\unicode[STIX]{x1D6FF}}=4000$ . Figure 12(a,b) shows the neutral curve of Mode I, while in figure 12(c) we show Mode II that becomes unstable at $\mathit{Ec}_{\infty }\geqslant 0.19$ . It appears that the maximum critical Reynolds number $\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ occurs at $\mathit{Ec}_{\infty }=0.16$ . The flow enters the transcritical regime (the temperature crosses the pseudo-critical point) at $\mathit{Ec}_{\infty }\geqslant 0.17$ , and the growth rate and the extent of the neutral curve of Mode I again increase. Figure 12(b) shows the evolution of Mode I in the transcritical regime. With an increase in $\mathit{Ec}_{\infty }$ , the range of unstable $\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ decreases, while the range for unstable $F$ increases. Most noteworthy is the growth rate of Mode II, shown in figure 12(c), which increases much faster with $\mathit{Ec}_{\infty }$ and becomes much larger than Mode I. This indicates that the flow in the transcritical regime is significantly destabilized by non-ideal-gas effects through Mode II. In addition, we also show the development of the phase velocity in figure 13. The phase velocity of Mode I remains between 0.2 and 0.54, while Mode II is between 0.3 and 0.35. Both increase with Eckert number.

Figure 12.

Growth rates of perturbations in the $F{-}\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ stability diagram with $T_{\infty }^{\ast }=280~\text{K}$ . (a) $\mathit{Ec}_{\infty }=0.11,0.12,\ldots ,0.19$ , (b) $\mathit{Ec}_{\infty }=0.190$ , $0.192,\ldots ,0.202$ (Mode I), (c)  $\mathit{Ec}_{\infty }=0.190$ , $0.192,\ldots ,0.202$ (Mode II).

Figure 13.

Phase velocities of perturbations in the $F{-}\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ stability diagram with $T_{\infty }^{\ast }=280~\text{K}$ . (a)  $\mathit{Ec}_{\infty }=0.11,0.12,\ldots ,0.19$ , (b)  $\mathit{Ec}_{\infty }=0.190$ , $0.192,\ldots ,0.202$ (Mode I), (c)  $\mathit{Ec}_{\infty }=0.190,0.192,\ldots ,0.202$ (Mode II).

Similar to the super- and subcritical cases discussed before, we also use DNS to validate the results from LST and to confirm the co-existence of the dual modes, as well as the dominance of Mode II. We show the neutral curve for $\mathit{Ec}_{\infty }=0.20$ (case T280E4) in the $F{-}x$ diagram in figure 14. Mode II has a much larger unstable region both in terms of $F$ and $x$ . In order to properly observe the two modes, the flow is perturbed with frequencies of $F_{1}=15\times 10^{-6}$ and $F_{2}=75\times 10^{-6}$ , respectively. From the LST prediction in figure 14, the forced frequency $F_{1}$ shall excite Mode II on entering the neutral curve, and the perturbation of frequency $F_{2}$ will be modulated by Mode I and II sequentially.

Figure 14.

Neutral curve of the T280E4 case in the $F{-}x$ diagram. The magenta and blue lines indicate the domain and frequency $F_{1}=15\times 10^{-6}$ and $F_{2}=75\times 10^{-6}$ simulated in the DNS.

Figure 15 shows the DNS validation, where the growth rate and phase velocity are compared with LST. In the $F_{1}$ case, the perturbation grows due to the positive growth rate of Mode II. In the $F_{2}$ case, the instability is initially dominated by Mode I and sequentially by Mode II, due to the cross-over of the growth rates of the two modes as shown in figure 15(c). In both cases, the DNS matches well with the LST predictions.

Figure 15.

Comparison of the growth rate and phase velocity between LST and DNS: $T_{\infty }^{\ast }=280~\text{K}$ ; $\mathit{Ec}_{\infty }=0.20$ , (a,b) $F=15\times 10^{-6}$ ; (c,d) $F=75\times 10^{-6}$ .

Figure 16.

Comparison of the profiles ( $\hat{u}$ , $\hat{v}$ , $\hat{\unicode[STIX]{x1D70C}}$ ) of the perturbations between LST and DNS. The amplitude has been normalized by $|\hat{u} |$ : $T_{\infty }^{\ast }=280~\text{K}$ , $\mathit{Ec}_{\infty }=0.20$ , $F=75\times 10^{-6}$ ; (a) $x=16.00$ , $\mathit{Re}_{\unicode[STIX]{x1D6FF}}=1264.9$ ; (b) $x=36.00$ , $\mathit{Re}_{\unicode[STIX]{x1D6FF}}=1894.7$ . The critical point ( $y=y_{c}$ , $u_{0}(y_{c})=c$ ) and the generalized inflection point ( $y=y_{i}$ ) are denoted with the blue dashed line and red dash-dotted line respectively.

The perturbation profiles with $F=75\times 10^{-6}$ obtained from LST and DNS are provided and compared in figure 16. The comparison is made at $x=16.00$ (Mode I dominates) and $x=36.00$ (Mode II dominates). Figure 16 shows that the perturbation profiles obtained from DNS match very well with the LST predictions of Mode I and Mode II at $x=16.00$ and $x=36.00$ , respectively. Due to the large gradient of $\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}/\unicode[STIX]{x2202}y$ , the density perturbation is the largest component for both modes. Note that $|\hat{\unicode[STIX]{x1D70C}}|$ is considerably larger for Mode II. On the other hand, the velocity components $\hat{u}$ and $\hat{v}$ remain similar for both modes, indicating that Mode II is caused by the dramatic variation of thermodynamic and transport properties. The above comparisons validated the co-existence of the dual modes and the dominance of Mode II in the transcritical regime.

4.3 Is Mode II comparable to Mack’s second mode?

Based on the findings in § 4.2, Mode II is generated when the fluid is transcritical although the flow is subsonic. In this section, using the inviscid theory, we clarify the relation to Mack’s second mode which appears in hypersonic flows.

The inviscid stability equations for the 2-D perturbations are given by

(4.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle (\unicode[STIX]{x1D6FC}u_{0}-\unicode[STIX]{x1D714})\hat{\unicode[STIX]{x1D70C}}+\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{0}\hat{u} =\text{i}\unicode[STIX]{x1D70C}_{0}D\hat{v}+\text{i}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}y}}\hat{v}\\ \displaystyle \unicode[STIX]{x1D6FC}\hat{p}+(\unicode[STIX]{x1D6FC}u_{0}-\unicode[STIX]{x1D714})\unicode[STIX]{x1D70C}_{0}\hat{u} =\text{i}\unicode[STIX]{x1D70C}_{0}{\displaystyle \frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}}\hat{v}\\ \displaystyle (\unicode[STIX]{x1D6FC}u_{0}-\unicode[STIX]{x1D714})\unicode[STIX]{x1D70C}_{0}\hat{v}-\text{i}D\hat{p}=0\\ \displaystyle (\unicode[STIX]{x1D6FC}u_{0}-\unicode[STIX]{x1D714})\unicode[STIX]{x1D70C}_{0}\hat{e}+\unicode[STIX]{x1D6FC}p_{0}\hat{u} =\text{i}p_{0}D\hat{v}+\text{i}\unicode[STIX]{x1D70C}_{0}{\displaystyle \frac{\unicode[STIX]{x2202}e_{0}}{\unicode[STIX]{x2202}y}}\hat{v}\end{array}\right\},\end{eqnarray}$$

where $e_{0}$ and $\hat{e}$ are the internal energy and its perturbation. Equation (4.3) can be reduced to a single equation,

(4.4) $$\begin{eqnarray}\displaystyle & & \displaystyle \left[\left(\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}e_{0}}}\right|_{p}{\displaystyle \frac{p_{0}}{\unicode[STIX]{x1D70C}_{0}^{2}}}-1\right)\unicode[STIX]{x1D6FC}^{2}+(\unicode[STIX]{x1D6FC}u_{0}-\unicode[STIX]{x1D714})^{2}\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}p_{0}}}\right|_{e}\right]\hat{p}-\left(\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}e_{0}}}\right|_{p}{\displaystyle \frac{p_{0}}{\unicode[STIX]{x1D70C}_{0}^{2}}}-1\right)D^{2}\hat{p}\nonumber\\ \displaystyle & & \displaystyle \quad +\left[\!\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}e_{0}}}\right|_{p}{\displaystyle \frac{p_{0}}{\unicode[STIX]{x1D70C}_{0}^{3}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}y}}+\!\left(\!\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}e_{0}}}\right|_{p}{\displaystyle \frac{p_{0}}{\unicode[STIX]{x1D70C}_{0}^{2}}}-1\!\right)\!{\displaystyle \frac{2\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6FC}u_{0}-\unicode[STIX]{x1D714}}}{\displaystyle \frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}}-{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}_{0}}}\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}e_{0}}}\right|_{p}\!{\displaystyle \frac{\unicode[STIX]{x2202}e_{0}}{\unicode[STIX]{x2202}y}}\!\right]D\hat{p}=0.\qquad \quad\end{eqnarray}$$

Following Mack (Reference Mack1984), neglecting the term related to $D\hat{p}$ results in

(4.5) $$\begin{eqnarray}\displaystyle D^{2}\hat{p}-\unicode[STIX]{x1D6FC}^{2}\left[1-{\displaystyle \frac{(u_{0}-c)^{2}}{1-\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}e_{0}}}\right|_{p}{\displaystyle \frac{p_{0}}{\unicode[STIX]{x1D70C}_{0}^{2}}}}}\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}p_{0}}}\right|_{e}\right]\hat{p}=0.\end{eqnarray}$$

We introduce the relative Mach number $M_{r}$ ,

(4.6) $$\begin{eqnarray}\displaystyle M_{r}={\displaystyle \frac{u_{0}-c}{\sqrt{\left(1-\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}e_{0}}}\right|_{p}{\displaystyle \frac{p_{0}}{\unicode[STIX]{x1D70C}_{0}^{2}}}\right)\left(\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}p_{0}}}\right|_{e}\right)^{-1}}}}.\end{eqnarray}$$

Using Maxwell relations for partial derivatives of thermodynamic properties, it can be mathematically proven (see appendix D) that the relative Mach number (4.6) for non-ideal gases equals to

(4.7) $$\begin{eqnarray}\displaystyle M_{r}={\displaystyle \frac{u_{0}-c}{\sqrt{\left.{\displaystyle \frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}}\right|_{s}}}}={\displaystyle \frac{u_{0}-c}{a_{0}}},\end{eqnarray}$$

which is equivalent to the result of Mack (Reference Mack1984). Equation (4.5) changes its behaviour from elliptic to hyperbolic when $(1-M_{r}^{2})$ changes sign (from positive to negative) and multiple solutions (modes) can be present. This has successfully explained Mack’s second mode in high-speed flows of ideal gas. Equation (4.7) indicates that the relative Mach number for non-ideal gases possesses the same physical nature, i.e. a local supersonic region relative to the phase velocity can give rise to multiple modes. We show in figure 17 the distribution of $M_{r}^{2}$ . The results are given for phase velocities $c=$ 0.30 and 0.54, which correspond to the lower and upper limits of the dual modes ( $0.19\leqslant \mathit{Ec}_{\infty }\leqslant 0.20$ ). The figure shows that $M_{r}$ stays always below 1 in the transcritical regime ( $T_{\infty }^{\ast }=280~\text{K}$ and $\mathit{Ec}_{\infty }=0.20$ ) where the dual modes co-exist. This implies that Mode II is different from Mack’s second mode.

Figure 17.

The relative Mach number with $T_{\infty }^{\ast }=280~\text{K}$ and $\mathit{Ec}_{\infty }=0.20$ .

Another evidence to highlight the difference between Mode II and Mack’s second mode is the co-existence of dual unstable modes with same parameters ( $F$ , $\unicode[STIX]{x1D6FD}$ , $\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ and $\mathit{Ec}_{\infty }$ ). In fact, according to the new terminology summarized in Fedorov & Tumin (Reference Fedorov and Tumin2011), Mack’s second mode is a result of synchronization between the fast mode and the slow mode, which stems from the continuous spectrum of the fast and slow acoustic waves. As such, only one unstable eigenmode is present for certain parameters ( $F$ , $\unicode[STIX]{x1D6FD}$ , $\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ and $\mathit{Ec}_{\infty }$ ) of the base flow and the perturbation. This leads to the conventional neutral curve for hypersonic boundary-layer flows (see for example figure 1 in Ren, Fu & Hanifi Reference Ren, Fu and Hanifi2016). However, the neutral curves of Mode I and Mode II presented here, indeed overlap with each other. Both modes can be unstable with the same frequency $F$ and wavenumber.

Figure 18.

Generalized derivatives of the base flow $\text{d}(\unicode[STIX]{x1D70C}_{0}\text{d}u_{0}/\text{d}y)/\text{d}y$ for cases (a) T260, (b) T280, (c) T320 and (d) T340.

Under the inviscid assumption, Lees & Lin (Reference Lees and Lin1946) derived the generalized inflection point criterion $D(Du_{0}/T_{0})$ , which gives the necessary condition for a compressible boundary-layer flow to be inviscidly unstable. They used the ideal-gas equation of state ( $\unicode[STIX]{x1D70C}_{0}T_{0}=1$ ) and below we show that a similar criterion applies for non-ideal gases. From (4.4) for $\hat{p}$ and using the relation $D\hat{p}=-\text{i}\unicode[STIX]{x1D70C}_{0}(\unicode[STIX]{x1D6FC}u_{0}-\unicode[STIX]{x1D714})\hat{v}$ , the equation for $\hat{v}$ then reads,

(4.8) $$\begin{eqnarray}(\unicode[STIX]{x1D6FC}^{2}+\unicode[STIX]{x1D703})\hat{v}+\unicode[STIX]{x1D709}D\hat{v}+\unicode[STIX]{x1D702}D^{2}\hat{v}=0,\end{eqnarray}$$

where

(4.9) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D703}={\displaystyle \frac{2M_{r}DM_{r}Du_{0}}{(1-M_{r}^{2})^{2}(u_{0}-c)}}+{\displaystyle \frac{D(\unicode[STIX]{x1D70C}_{0}Du_{0})}{\unicode[STIX]{x1D70C}_{0}(1-M_{r}^{2})(u_{0}-c)}},\\ \unicode[STIX]{x1D709}=-{\displaystyle \frac{2M_{r}DM_{r}}{(1-M_{r}^{2})^{2}}}-{\displaystyle \frac{D\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x1D70C}_{0}(1-M_{r}^{2})}},\\ \unicode[STIX]{x1D702}=-{\displaystyle \frac{1}{1-M_{r}^{2}}}.\end{array}\right\}\end{eqnarray}$$

It can be recognized that at the critical layer $y=y_{c}$ , where $u_{0}=c$ and $M_{r}=0$ , the term $D\left(\unicode[STIX]{x1D70C}_{0}Du_{0}\right)$ must vanish such that $y=y_{c}$ is a regular singular point of the $\hat{v}$ equation (4.8). Therefore, $D\left(\unicode[STIX]{x1D70C}_{0}Du_{0}\right)$ plays the same role as $D\left(Du_{0}/T_{0}\right)$ given by Lees & Lin (Reference Lees and Lin1946) and provides the generalized inflection point criterion for non-ideal gases.

In figure 18, we show the generalized derivatives $\text{d}(\unicode[STIX]{x1D70C}_{0}\text{d}u_{0}/\text{d}y)/\text{d}y$ for cases T260, T280, T320 and T340, each with four Eckert numbers. It indicates that like ideal gases, an inviscid mechanism exists in the supercritical regime, while the stability is viscous in nature in the subcritical regime. It is worth noting that a generalized inflection point suddenly appears once the base flow turns from the sub- to the transcritical regime. It suggests the inviscid nature of Mode II. Note in figure 16, if compared to Mode I, Mode II has the peak of fluctuations ( $|\hat{\unicode[STIX]{x1D70C}}|$ , $|\hat{u} |$ , $|\hat{p}|$ ) around the generalized inflection point, suggesting its inviscid nature (Lees & Lin Reference Lees and Lin1946). In figure 19 we compare the generalized derivative of the transcritical base flow for pressures $p_{0}^{\ast }=$ 78, 80 and 82 bar. It indicates that the closer to the critical point ( $p_{c}^{\ast }=73.9$ bar), the more inflectional the base flow becomes, and the more unstable the flow in the transcritical regime.

Figure 19.

Generalized derivative of the laminar base flow with $T_{\infty }^{\ast }=280~\text{K}$ , (a) $p_{0}^{\ast }=78$ bar, (b) $p_{0}^{\ast }=80$ bar and (c) $p_{0}^{\ast }=82$ bar.

4.4 Oblique perturbations

Figure 20.

Neutral surface along with slices of growth rate contours. $T_{\infty }^{\ast }=280~\text{K}$ , (a) $\mathit{Ec}_{\infty }=0.05$ ; (b) $\mathit{Ec}_{\infty }=0.20$ , Mode I; (c) $\mathit{Ec}_{\infty }=0.20$ , Mode II.

The LST and DNS studies in previous sections have both focused on the 2-D perturbations where $\unicode[STIX]{x1D6FD}=0$ . Here we comment on the oblique perturbations. Examples are provided in figure 20, which shows the neutral surface as well as slices of growth rate contours in the $\mathit{Re}_{\unicode[STIX]{x1D6FF}}{-}F{-}B$ diagram; $B=\unicode[STIX]{x1D6FD}/\mathit{Re}_{\unicode[STIX]{x1D6FF}}$ is defined as a global spanwise wavenumber. Figure 20(a) shows that with $T_{\infty }^{\ast }=280~\text{K}$ , the single mode at $\mathit{Ec}_{\infty }=0.05$ has a larger growth rate at $B=0$ . Increase in $B$ results in a diminished overall growth rate and a size-reduced neutral curve. This is however not true for Mode I at $\mathit{Ec}_{\infty }=0.20$ . It is shown in figure 20(b) that the largest growth happens at non-zero spanwise wavenumber $B$ . Another example we show is Mode II in figure 20(c). This mode reaches its maximum growth rate at $B=0$ .

A summary of the oblique effects is tabulated in table 2, where we show the results at $\mathit{Re}_{\unicode[STIX]{x1D6FF}}=2000$ . The optimal spanwise wavenumber at which the growth rate reaches its maximum (over all frequencies) is termed $\unicode[STIX]{x1D6FD}_{m}$ . The growth rate ratio shows the ratio between the maximum growth rate at $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{m}$ and the value with $\unicode[STIX]{x1D6FD}=0$ (over all frequencies). Table 2 clearly indicates that 2-D perturbations are the most unstable in the subcritical regime and in the transcritical regime for Mode II. In the supercritical regime for large $\mathit{Ec}_{\infty }$ and in the transcritical regime for Mode I, oblique perturbations become more unstable.

Table 2.

Summary of the 3-D perturbation effects at $\mathit{Re}_{\unicode[STIX]{x1D6FF}}=2000$ ; $\unicode[STIX]{x1D6FD}_{m}$ indicates the optimal spanwise wavenumber at which the growth rate reaches maximum. Growth rate ratio gives the ratio of maximum growth rate at $\unicode[STIX]{x1D6FD}_{m}$ to the maximum value of 2-D perturbations.

Figure 21.

Comparison of the 2-D ( $\unicode[STIX]{x1D6FD}=0$ ) and oblique ( $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{m}$ ) eigenfunctions at $\mathit{Re}_{\unicode[STIX]{x1D6FF}}=2000$ . (a) Mode I in the transcritical regime, $T_{\infty }^{\ast }=280~\text{K}$ , $\mathit{Ec}_{\infty }=0.20$ ; (b) the supercritical regime, $T_{\infty }^{\ast }=320~\text{K}$ , $\mathit{Ec}_{\infty }=0.15$ . $\unicode[STIX]{x1D714}$ is chosen such that each mode reaches its maximum growth rate.

In terms of eigenfunctions, figure 21 compares profiles of the oblique waves with their 2-D counterpart. Following table 2, we show case T280E4 (Mode I) and T320E3. Both reach maximum growth rate at a non-zero spanwise wavenumber. With $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{m}$ , the growth rate is considerably increased, and eigenfunctions are more mature and confined in the boundary layer. When the perturbation is oblique, the spanwise velocity perturbation $|{\hat{w}}|$ is present, while the density component $|\hat{\unicode[STIX]{x1D70C}}|$ becomes relatively smaller. Note that the density profile has been scaled by a factor of 4 in panel (a) for a better display of the results.