2.1 Flow conservation equations

The laws of conservation of mass, momentum and energy (the Navier–Stokes (N–S) equations), in dimensionless form, are given by

(2.1) $$\begin{eqnarray}\displaystyle & \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\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \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}}=F_{i}, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \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}}=u_{j}F_{j}, & \displaystyle\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 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, $F_{i}$ the body force and $p$ is the pressure. The viscous stress tensor, $\unicode[STIX]{x1D70F}_{ij}$ , and the heat flux, $q_{j}$ , are given by

(2.4a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{ij}=\frac{\unicode[STIX]{x1D707}}{Re}\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}}{Re}\unicode[STIX]{x1D6FF}_{ij}\frac{\unicode[STIX]{x2202}u_{k}}{\unicode[STIX]{x2202}x_{k}},\quad q_{j}=-\frac{\unicode[STIX]{x1D705}}{RePrEc}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}x_{j}}. & & \displaystyle\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. Results presented in the following sections are subject to $\unicode[STIX]{x1D707}_{b}=0$ . However, we will discuss the influence of the bulk viscosity on the linear stability in appendix C.

The above equations have been non-dimensionalized by reference values, as follows

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

which leads to the definition of the Reynolds number, $Re$ , Prandtl number, $Pr$ , Eckert number, $Ec$ , and Mach number, $Ma$ , which are given as

(2.6a-d ) $$\begin{eqnarray}\displaystyle Re=\frac{\unicode[STIX]{x1D70C}_{w}^{\ast }u_{r}^{\ast }h^{\ast }}{\unicode[STIX]{x1D707}_{w}^{\ast }},\quad Pr=\frac{\unicode[STIX]{x1D707}_{w}^{\ast }C_{pw}^{\ast }}{\unicode[STIX]{x1D705}_{w}^{\ast }},\quad Ec=\frac{u_{r}^{\ast 2}}{C_{pw}^{\ast }T_{w}^{\ast }},\quad Ma=\frac{u_{r}^{\ast }}{c_{w}^{\ast }}. & & \displaystyle\end{eqnarray}$$

The subscript $w$ denotes wall values, $h^{\ast }$ is the half-channel height, $c_{w}^{\ast }$ is the speed of sound at the wall, $u_{r}^{\ast }$ is the reference velocity. Note that for an ideal gas $Ec=(\unicode[STIX]{x1D6FE}-1)Ma^{2}$ , where $\unicode[STIX]{x1D6FE}$ is the heat capacity ratio. In this study, both walls are at the same temperature. Discussions on the choice of different reference scalings are provided in appendix D.

2.2 Fluid EoS

In order to find a closed form of the conservation equations, an EoS for the fluid has to be specified. As a representative example of non-ideal fluids, the study is performed with $\text{CO}_{2}$ at a pressure of $p^{\ast }=80$ bar, which is above the critical pressure, within the highly non-ideal thermodynamic region (see the isobar in figure 1). To account for the non-ideal gas effects, the NIST REFPROP library (Lemmon, Huber & McLinden Reference Lemmon, Huber and McLinden2013) has been used to obtain the most accurate thermodynamic and transport properties along with their gradients. The multi-parameter EoS (in functional forms) used in REFPROP are developed with an optimization algorithm. These EoS are suitable for a broad variety of fluids while high accuracy can be maintained. Readers may wish to refer to Span & Wagner (Reference Span and Wagner2003) for the derivation of the EoS. To build the linear stability equations (see appendix A), the temperature $T_{0}$ and density $\unicode[STIX]{x1D70C}_{0}$ are provided as input, while the required properties and their derivatives are obtained as output from REFPROP. Moreover, as a direct method to determine the thermodynamic properties, several cubic EoS (see appendix B), i.e. van der Waals (van der Waals Reference van der Waals1873), Redlich–Kwong (Redlich & Kwong Reference Redlich and Kwong1949) and Peng–Robinson (Peng & Robinson Reference Peng and Robinson1976), are used for the stability analysis as comparison. All results with the non-ideal EoS are also compared with an ideal gas (IG) model. A constant specific heat ratio $\unicode[STIX]{x1D6FE}=1.289$ is used for the IG model. All the fluid models are summarized in table 1.

Figure 2.

Thermodynamic and transport properties of $\text{CO}_{2}$ at $p^{\ast }=80$ bar for the fluid models in table 1. Distribution of (a) density $\unicode[STIX]{x1D70C}^{\ast }$ , (b) heat capacity at constant pressure $C_{p}^{\ast }$ , (c) viscosity $\unicode[STIX]{x1D707}^{\ast }$ and (d) thermal conductivity $\unicode[STIX]{x1D705}^{\ast }$ versus temperature $T^{\ast }$ . The pentagram shows the pseudo-critical temperature $T_{pc}^{\ast }$ (RP model). The shaded area indicates the pseudo-critical transition.

Table 1.

Fluid models studied in this paper. Gradients of the properties (with respect to temperature and density) are calculated analytically (using EoS, see appendix B) or numerically with a finite-difference algorithm (using REFPROP, see appendix A).

Figure 2 shows the thermodynamic and transport properties of $\text{CO}_{2}$ at a pressure of 80 bar. The pentagram in (a) shows the pseudo-critical temperature ( $T_{pc}^{\ast }=307.7~\text{K}$ , RP model), which is defined as the point on a supercritical isobar where $C_{p}^{\ast }$ reaches a maximum. Near $T_{pc}^{\ast }$ , all properties show large gradients, which do not exist in an ideal gas. As shown in figure 2(a,b), the Peng–Robinson (PR) EoS is closest to the highly accurate multiparameter EoS of $\text{CO}_{2}$ as implemented in REFPROP (RP). Principally, the cubic EoS do capture key features of the thermodynamic property variations. In figure 2(c,d), the power law (2.7) and Sutherland law (2.8), which fall on top of each other, are compared to the distributions from RP. The power and Sutherland laws for dynamic viscosity and thermal conductivity are given as

(2.7a,b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D707}^{\ast }}{\unicode[STIX]{x1D707}_{ref}^{\ast }}=\left(\frac{T^{\ast }}{T_{ref}^{\ast }}\right)^{n_{1}},\quad \frac{\unicode[STIX]{x1D705}^{\ast }}{\unicode[STIX]{x1D705}_{ref}^{\ast }}=\left(\frac{T^{\ast }}{T_{ref}^{\ast }}\right)^{n_{2}},\quad n_{1}=0.79,n_{2}=1.30,T_{ref}^{\ast }=273~\text{K},\qquad & & \displaystyle\end{eqnarray}$$

(2.8a,b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D707}^{\ast }}{\unicode[STIX]{x1D707}_{ref}^{\ast }}=\left(\frac{T^{\ast }}{T_{ref}^{\ast }}\right)^{3/2}\frac{T_{ref}^{\ast }+S_{1}^{\ast }}{T^{\ast }+S_{1}^{\ast }},\quad \frac{\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D705}_{ref}}=\left(\frac{T^{\ast }}{T_{ref}^{\ast }}\right)^{3/2}\frac{T_{ref}^{\ast }+S_{2}^{\ast }}{T^{\ast }+S_{2}^{\ast }}, & & \displaystyle\end{eqnarray}$$

where

(2.9) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle T_{ref}^{\ast }=273~\text{K},\quad \unicode[STIX]{x1D707}_{ref}^{\ast }=1.37\times 10^{-5}~\text{Pa}~\text{s},\quad \unicode[STIX]{x1D705}_{ref}^{\ast }=0.0146~\text{W}~(\text{m}~\text{K})^{-1},\\ \displaystyle n_{1}=0.79,\quad n_{2}=1.30,\quad S_{1}^{\ast }=222~\text{K},\quad S_{2}^{\ast }=1800~\text{K}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

As shown in figure 2(c,d), there is no discernible difference using the power or Sutherland law for the IG model, therefore results presented in this study for IG are based on the power law. In general, as temperature increases from subcritical to supercritical values, the fluid continuously transitions from compressed liquid to compressed vapour and finally reaches values that can be described by an ideal gas.

2.3 The linearized stability equations

Following the common procedure, the flow field is decomposed into the base flow and a perturbation, as

(2.10) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x1D70C}^{\prime },\\ \displaystyle u_{i}=u_{i0}+u_{i}^{\prime },\\ \displaystyle T=T_{0}+T^{\prime },\\ \displaystyle p=p_{0}+p^{\prime },\\ \displaystyle E=E_{0}+E^{\prime },\\ \displaystyle \unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D707}^{\prime },\\ \displaystyle \unicode[STIX]{x1D705}=\unicode[STIX]{x1D705}_{0}+\unicode[STIX]{x1D705}^{\prime }.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

It is known that for simple compressible systems (e.g. pure substances, uniform mixture of non-reacting gases), the thermodynamic state is defined by two independent thermodynamic properties. In this study, we keep $\unicode[STIX]{x1D70C}$ and $T$ as the two basic thermodynamic variables, while the other 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 pressure perturbation $p^{\prime }$ is expanded by a Taylor series with respect to $\unicode[STIX]{x1D70C}_{0}$ and $T_{0}$ in the following way

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

For the sake of brevity, the partial derivative of a quantity with respect to $T$ at constant $\unicode[STIX]{x1D70C}$ will be written as $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}T|_{\unicode[STIX]{x1D70C}_{0}}\equiv \unicode[STIX]{x2202}/\unicode[STIX]{x2202}T$ , and accordingly $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}|_{T_{0}}\equiv \unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}$ . The stability equation is derived by substituting (2.10) into the N–S equations (2.1)–(2.3), and then subtracting the governing equations of the base flow. With the nonlinear terms neglected, the linear stability equations are formulated as

(2.12) $$\begin{eqnarray}\displaystyle & & \displaystyle {\mathcal{L}}_{t}\frac{\unicode[STIX]{x2202}\boldsymbol{q}}{\unicode[STIX]{x2202}t}+{\mathcal{L}}_{x}\frac{\unicode[STIX]{x2202}\boldsymbol{q}}{\unicode[STIX]{x2202}x}+{\mathcal{L}}_{y}\frac{\unicode[STIX]{x2202}\boldsymbol{q}}{\unicode[STIX]{x2202}y}+{\mathcal{L}}_{z}\frac{\unicode[STIX]{x2202}\boldsymbol{q}}{\unicode[STIX]{x2202}z}+{\mathcal{L}}_{q}\boldsymbol{q}\nonumber\\ \displaystyle & & \displaystyle \quad +\,{\mathcal{V}}_{xx}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{q}}{\unicode[STIX]{x2202}x^{2}}+{\mathcal{V}}_{xy}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{q}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}y}+{\mathcal{V}}_{xz}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{q}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}z}+{\mathcal{V}}_{yy}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{q}}{\unicode[STIX]{x2202}y^{2}}+{\mathcal{V}}_{yz}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{q}}{\unicode[STIX]{x2202}y\unicode[STIX]{x2202}z}+{\mathcal{V}}_{zz}\frac{\unicode[STIX]{x2202}^{2}\boldsymbol{q}}{\unicode[STIX]{x2202}z^{2}}=0.\end{eqnarray}$$

Here $\boldsymbol{q}=(\unicode[STIX]{x1D70C}^{\prime },u^{\prime },v^{\prime },w^{\prime },T^{\prime })^{\text{T}}$ is the perturbation vector and ${\mathcal{L}}_{t}$ , ${\mathcal{L}}_{x}$ , ${\mathcal{L}}_{y}$ , ${\mathcal{L}}_{z}$ , ${\mathcal{L}}_{q}$ , ${\mathcal{V}}_{xx}$ , ${\mathcal{V}}_{yy}$ , ${\mathcal{V}}_{zz}$ , ${\mathcal{V}}_{xy}$ , ${\mathcal{V}}_{yz}$ and ${\mathcal{V}}_{xz}$ are matrices of size $5\times 5$ . The detailed expressions for these matrices are provided in appendix A. As can be seen, they are functions of the base flow, the thermodynamic and transport properties, $Re$ and $PrEc$ . The gradients of the properties are either calculated analytically using a cubic EoS (see appendix B) or numerically employing the finite-difference method within the REFPROP library.

2.4 Modal and algebraic stability

In modal stability, the perturbation is assumed to have the form

(2.13) $$\begin{eqnarray}\displaystyle \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.} & & \displaystyle\end{eqnarray}$$

where $\text{c.c.}$ stands for the complex conjugate. Substituting (2.13) into (2.12) results in

(2.14) $$\begin{eqnarray}\displaystyle & & \displaystyle (\!-\text{i}\unicode[STIX]{x1D714}{\mathcal{L}}_{t}+\text{i}\unicode[STIX]{x1D6FC}{\mathcal{L}}_{x}+{\mathcal{L}}_{y}D+\text{i}\unicode[STIX]{x1D6FD}{\mathcal{L}}_{z}+{\mathcal{L}}_{q}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\unicode[STIX]{x1D6FC}^{2}{\mathcal{V}}_{xx}+\text{i}\unicode[STIX]{x1D6FC}{\mathcal{V}}_{xy}D-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}{\mathcal{V}}_{xz}+{\mathcal{V}}_{yy}D^{2}+\text{i}\unicode[STIX]{x1D6FD}{\mathcal{V}}_{yz}D-\unicode[STIX]{x1D6FD}^{2}{\mathcal{V}}_{zz} )\hat{\boldsymbol{q}}=0,\end{eqnarray}$$

where $D=\text{d}/\text{d}y$ . Equation (2.14) is solved as an eigenvalue problem, which describes the development of the perturbations in temporal or spatial domain, i.e.

(2.15) $$\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}ll@{}}L_{T}\hat{\boldsymbol{q}}=\unicode[STIX]{x1D714}{\mathcal{L}}_{t}\hat{\boldsymbol{q}}\quad & \text{(temporal)},\\ L_{S}\hat{\boldsymbol{q}}=\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FD}{\mathcal{V}}_{xz}-\text{i}{\mathcal{V}}_{xy}D-\text{i}{\mathcal{L}}_{x})\hat{\boldsymbol{q}}+\unicode[STIX]{x1D6FC}^{2}{\mathcal{V}}_{xx}\hat{\boldsymbol{q}}\quad & (\text{spatial}),\end{array}\right. & & \displaystyle\end{eqnarray}$$

where

(2.16) $$\begin{eqnarray}\displaystyle L_{T} & = & \displaystyle \unicode[STIX]{x1D6FC}{\mathcal{L}}_{x}-\text{i}{\mathcal{L}}_{y}D+\unicode[STIX]{x1D6FD}{\mathcal{L}}_{z}-\text{i}{\mathcal{L}}_{q}\nonumber\\ \displaystyle & & \displaystyle +\,\text{i}\unicode[STIX]{x1D6FC}^{2}{\mathcal{V}}_{xx}+\unicode[STIX]{x1D6FC}{\mathcal{V}}_{xy}D+\text{i}\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}{\mathcal{V}}_{xz}-iD^{2}{\mathcal{V}}_{yy}+\unicode[STIX]{x1D6FD}D{\mathcal{V}}_{yz}+\text{i}\unicode[STIX]{x1D6FD}^{2}{\mathcal{V}}_{zz},\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle L_{S}=-\text{i}\unicode[STIX]{x1D714}{\mathcal{L}}_{t}+{\mathcal{L}}_{y}D+\text{i}\unicode[STIX]{x1D6FD}{\mathcal{L}}_{z}+{\mathcal{L}}_{q}+D^{2}{\mathcal{V}}_{yy}+\text{i}\unicode[STIX]{x1D6FD}D{\mathcal{V}}_{yz}-\unicode[STIX]{x1D6FD}^{2}{\mathcal{V}}_{zz}. & & \displaystyle\end{eqnarray}$$

Here we consider the temporal problem only, therefore $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ are the prescribed streamwise and spanwise wavenumbers. $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{r}+\text{i}\unicode[STIX]{x1D714}_{i}$ is solved as the eigenvalue, where $\unicode[STIX]{x1D714}_{r}$ and $\unicode[STIX]{x1D714}_{i}$ give the angular frequency and growth rate of the perturbation. The domain $0\leqslant y\leqslant 2$ is discretized with Chebyshev collocation points, defined by:

(2.18) $$\begin{eqnarray}\displaystyle y_{j}=1-\cos \frac{\unicode[STIX]{x03C0}j}{N},\quad j=0,1,\ldots ,N-1,N. & & \displaystyle\end{eqnarray}$$

The derivatives with respect to $y$ (the $D$ operator) in (2.15) are evaluated with the matrix form of the Chebyshev collocation derivatives. Numerical tests indicate that typically $N=200$ (used here) is sufficient to give a grid-independent solution of the physical modes.

With regard to the algebraic stability, following Schmid & Henningson (Reference Schmid and Henningson2001), the optimal energy amplification is defined as:

(2.19a,b ) $$\begin{eqnarray}\displaystyle G(t)=\underset{\boldsymbol{q}_{0}}{\max }\frac{E(\boldsymbol{q}(t))}{E(\boldsymbol{q}_{0})},\quad G(x)=\underset{\boldsymbol{q}_{0}}{\max }\frac{E(\boldsymbol{q}(x))}{E(\boldsymbol{q}_{0})}. & & \displaystyle\end{eqnarray}$$

Here $E(\boldsymbol{q})$ is the disturbance energy with the definition as given in (5.1). The perturbation $\boldsymbol{q}$ is expanded using the eigenvector basis obtained from the modal stability. The calculation of the optimal energy amplification $G$ , the corresponding optimal perturbation (the input), as well as the resulting perturbation (the output), are performed using singular value decomposition method described in Schmid & Henningson (Reference Schmid and Henningson2001), Schmid & Brandt (Reference Schmid and Brandt2014).

The (modal and algebraic) perturbations are solved subjected to the boundary condition: $u^{\prime }=v^{\prime }=w^{\prime }=T^{\prime }=0$ at the lower ( $y=0$ ) and upper wall ( $y=2$ ).

3.1 The isothermal limit

When $PrEc\rightarrow 0$ , the viscous heating is negligible if compared to the thermal conduction. Therefore, the temperature, as well as the other thermodynamic properties, remain constant, namely $T_{0}=1,\unicode[STIX]{x1D70C}_{0}=1,\unicode[STIX]{x1D707}_{0}=1,\unicode[STIX]{x1D705}_{0}=1$ . The flow is thus simply governed by $\unicode[STIX]{x2202}^{2}u/\unicode[STIX]{x2202}y^{2}=-\hat{F}$ . Choosing $u_{r}^{\ast }$ as the centreline velocity leads to setting $\hat{F}=2$ . As a result, the dimensionless base flow is independent of any parameters (e.g.  $T_{w}$ , $\hat{F}$ and $PrEc$ ) and is given by $u_{0}=y(2-y)$ . A sketch of this base flow, which is free from any non-ideal gas effects, is shown in figure 3 (dashed lines).

Figure 3.

Sketch of the laminar base flow. Dashed lines show the isothermal limit with $PrEc\rightarrow 0$ , such that $u_{0}=y(2-y)$ , $T_{0}=1$ , $\unicode[STIX]{x1D70C}_{0}=1$ . Solid lines represent a transcritical case, in which $T_{0}$ crosses $T_{pc}$ and $u_{0}$ is inflectional.

3.2 The compressible base flow

Equations (3.1)–(3.3) show that the compressible base flow is determined by $PrEc$ , $\hat{F}$ and $T_{w}^{\ast }$ . Either by increasing $PrEc$ or $\hat{F}$ , the compressibility effects become more significant. Without loss of the generality, a constant body force $\hat{F}=2$ is specified in this work, while $PrEc$ is varied from the isothermal limit (we assume $PrEc=10^{-5}$ ) to a typical compressible state with $PrEc=0.1$ . For example, setting $PrEc=0.1$ and $T_{w}^{\ast }=290,300$ or $310~\text{K}$ , the Mach number is $Ma=0.40,0.58$ or $1.35$ , respectively. In this work, the wall temperature $T_{w}^{\ast }$ is considered in a range from 265 to 320 K. Note, given our non-dimensionalization, the base flow is free from the choice of the Reynolds number.

Figure 4 shows the contours of the centreline temperature $T_{centre}^{\ast }$ and velocity $u_{centre}=u_{centre}^{\ast }/u_{r}^{\ast }$ as a function of wall temperature $T_{w}^{\ast }$ and $PrEc$ . Regardless of the wall temperature, an increase of $PrEc$ is accompanied by an increase of $T_{centre}^{\ast }$ and $u_{centre}$ as compressible effects become more prominent. Interestingly, a distinct right-angled triangular area emerges in each panel of figure 4. At the hypotenuse of this triangle, the centreline temperature, $T_{centre}^{\ast }$ , and velocity, $u_{centre}$ , suddenly increase, forming a discontinuity in the $PrEc$ – $T_{w}^{\ast }$ plane.

It is also interesting to note that the hypotenuse of the triangle almost coincides with the line where $T_{centre}^{\ast }$ reaches the pseudo-critical temperature $T_{pc}^{\ast }=307.7~\text{K}$ (shown with the dot-dashed line). Likewise, the upper boundary of the triangle coincides with the dotted line where $T_{w}^{\ast }=T_{pc}^{\ast }$ .

Table 2.

Cases investigated in this study.

Figure 4.

Centreline (a) temperature $T_{centre}^{\ast }$ and (b) velocity $u_{centre}$ as functions of wall temperature $T_{w}^{\ast }$ and $PrEc$ . Model RP, $p^{\ast }=80$ bar, $\hat{F}=2$ .

Figure 5.

Temperature (a–c), density (d–f) and velocity (g–i) profiles of the base flow. The wall temperature is (a,d,g) $T_{w}^{\ast }=290~\text{K}$ , (b,e,h) $T_{w}^{\ast }=300~\text{K}$ , (c,f,i) $T_{w}^{\ast }=310~\text{K}$ respectively. $PrEc$ increases uniformly from 0.01 to 0.1. The black and blue lines on the left and right half denote the non-ideal (RP) and ideal (IG) gases respectively. The dashed lines in each panel show the isothermal limit. The REFPROP library is used for the transport and thermodynamic properties of the non-ideal gas. The orange and red lines in panels (b,e,h) show the profiles at $PrEc=0.05115$ and 0.05116 respectively. The dash-dotted lines (the triangle) in (g,h,i) show the lines of constant gradient $|\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y|=\hat{F}$ .

For a more detailed discussion, we will now define three cases with different wall temperatures that are summarized in table 2 and highlighted by dashed lines in figure 4. These cases will also be used in the subsequent sections regarding the linear modal and algebraic instability analysis. The wall temperature for these cases has been set to 290, 300 and 310 K, such that their temperature profile in the considered range of $PrEc$ is either subcritical, transcritical or supercritical, respectively. Their base flow profiles are plotted in figure 5, together with the incompressible limit, indicated by the dashed line in each panel. The profiles on the left half (black lines) and right half (blue lines) represent the base flow of the non-ideal (RP) and ideal (IG) gases, respectively. As $PrEc$ uniformly increases from 0.01 to 0.1 it can be seen that the temperature and velocity increase, while the density decreases. For the transcritical case, however, a sudden jump of the base flow profiles can be observed. This jump occurs between $PrEc=$ 0.05115 and 0.05116, as highlighted by the orange and red lines in figure 5(b,e,h). Note, the jump is caused by an inflectional velocity profile as highlighted by the red line in figure 5(h).

The discontinuous behaviour with respect to $PrEc$ can be explained as follows. Integrating (3.1) gives $\unicode[STIX]{x1D707}\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y=-\hat{F}y+C.$ Applying the symmetry condition at the channel centre ( $y=1$ ), it follows that $C=\hat{F}$ . Therefore, (3.1) can be written as

(3.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}^{2}u}{\unicode[STIX]{x2202}y^{2}}=-\hat{F}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}y}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}=-\hat{F}\left(1+\frac{1}{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}y}(1-y)\right). & & \displaystyle\end{eqnarray}$$

Based on (3.4), it can be seen that an inflectional velocity profile occurs if the viscosity gradient is large enough to change the sign within the parenthesis in (3.4), namely if

(3.5) $$\begin{eqnarray}\displaystyle \frac{1}{\unicode[STIX]{x1D707}}\left|\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}y}\right|>1. & & \displaystyle\end{eqnarray}$$

In the cases considered herein, it appears that the viscosity gradient at the wall is large enough to cause an inflectional profile to occur when the temperature in the channel centre reaches $T_{pc}^{\ast }$ . Recall figure 2(c), a sharp gradient of the viscosity ( $\unicode[STIX]{x2202}\unicode[STIX]{x1D707}/\unicode[STIX]{x2202}T\ll 0$ ) is seen close to the pseudo-critical point. As $PrEc$ increases, $\unicode[STIX]{x2202}T/\unicode[STIX]{x2202}y$ increases at the wall, such that $\unicode[STIX]{x2202}\unicode[STIX]{x1D707}/\unicode[STIX]{x2202}y\cong (\unicode[STIX]{x2202}\unicode[STIX]{x1D707}/\unicode[STIX]{x2202}T)(\unicode[STIX]{x2202}T/\unicode[STIX]{x2202}y)$ can drop below $-1$ at the wall, leading to inflectional velocity profiles. The jump of the base flow solution can thus be explained by referring to figure 5(h). Since, $\unicode[STIX]{x1D707}|\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y|$ at the wall is equal to the constant forcing $\hat{F}$ , regardless of $PrEc$ and wall temperature, the velocity profiles with/without inflectional points are isolated by the line of constant gradient $\hat{F}$ (the dash-dotted lines that form a triangle in figure 5 g–i). Therefore, a velocity profile without an inflection point cannot reach the apex of the triangle ( $|\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y|$ decreases towards the channel centre) and the sudden increase of the centreline velocity appears once an inflection point is formed.

In general, the base flow solutions can be summarized as follows:

1. (i) In the subcritical case, the wall temperature is much lower than $T_{pc}^{\ast }$ , and in the range of $PrEc$ considered, $T_{centre}^{\ast }$ is always less than $T_{pc}^{\ast }$ . Hence, the velocity profile is not inflectional.

2. (ii) In the transcritical case, the wall temperature is close to $T_{pc}^{\ast }$ , such that for large enough $PrEc$ , $T_{centre}^{\ast }$ reaches $T_{pc}^{\ast }$ . Consequently, a jump of the solution with respect to $PrEc$ occurs and the velocity profile becomes inflectional. From figure 4(b), it can be inferred that the lower the wall temperature $T_{w}^{\ast }$ , the larger the discontinuity will be.

3. (iii) In the supercritical case, the wall temperature is higher than $T_{pc}^{\ast }$ . The properties of the fluid are gas-like (compressed vapour) and the velocity is not inflectional.

4.1 The isothermal limit

Figure 6.

Eigenspectrum (a) and neutral curve (b) for the isothermal limit. The eigenspectrum is subject to $\unicode[STIX]{x1D6FC}=1$ , $\unicode[STIX]{x1D6FD}=0$ and $Re=10\,000$ . The neutral curve is solved for two-dimensional perturbations ( $\unicode[STIX]{x1D6FD}=0$ ). Symbols show results using different fluid models (RP, PR, RK, VW, IG) and incompressible equations (IC). IC in (b) shows the results given by Schmid & Henningson (Reference Schmid and Henningson2001, p. 71). In (c) and (d), profiles of the unstable mode ( $\unicode[STIX]{x1D714}=0.2375+0.0037\text{i}$ ) and one of the stable modes ( $\unicode[STIX]{x1D714}=0.4164-0.1382\text{i}$ , highlighted in orange in the spectrum) are shown. The perturbations are normalized by $|u^{\prime }|$ . An offset of $-$ 0.1 and $-$ 0.2 is applied to $|\unicode[STIX]{x1D70C}^{\prime }|$ and $|T^{\prime }|$ . The solid lines are results with fluid model IG.

With the base flow obtained in § 3.1, we solve the stability equations (2.12) for the isothermal limit with different fluid models (RP, PR, RK, VW, IG), as well as for the incompressible equations (IC). As shown in figure 6(a), at $Re=10\,000$ , $\unicode[STIX]{x1D6FC}=1$ and $\unicode[STIX]{x1D6FD}=0$ , the A-, P- and S-branches (originally named by Mack (Reference Mack1976)) are reproduced by incompressible equations. Comparing the results using different EoS, the eigenvalues fall on top of the incompressible counterparts, verifying the correct behaviour of the compressible models at low Eckert (Mach) numbers. One of the modes (highlighted in red) is exclusively unstable. Despite being solved with different thermodynamic models, this mode is shown to be a dynamic mode, which leads to identical neutral curve and eigenfunctions as shown in figure 6(b,c). The contour lines in figure 6(b) show the growth rate $\unicode[STIX]{x1D714}_{i}$ (RP model). In fact, inspecting the stability equations (2.12) (see appendix A), it can be shown that the thermodynamic and transport properties do not influence the dynamic modes in the isothermal limit. For instance, gradients of properties, which vary among different models, are multiplied with thermodynamic components of the perturbations.

On the other hand, more stable modes emerge when the compressible equations are solved. By looking into the corresponding eigenfunctions (not shown), thermodynamic components become important in these modes, and as such dependent on the non-ideal gas properties. We plot one of the stable modes in figure 6(d), where density perturbations are captured by compressible equations (indicated by blue ellipses).

4.2 Compressible flows

To achieve a first impression of the non-ideal gas effects, the problem is first studied with the RP model, where thermodynamic and transport properties are taken from the REFPROP library. Figure 7 shows the neutral curves (a–c) as well as eigenfunctions (d–f) at representative parameters. As discussed in § 3.2, the temperature is subcritical, transcritical and supercritical with $T_{w}^{\ast }=290~\text{K}$ , 300 K and 310 K respectively. The results are compared with those of an ideal gas (IG).

By increasing $PrEc$ , the base flow of the ideal gas becomes more stable as the critical Reynolds number increases, regardless of the $T_{w}$ specified. In fact, despite the difference in wall temperature, the dimensionless thermodynamic and transport properties (scaled with wall values) remain much the same. On the other hand, the behaviour of the non-ideal cases is different for the three cases investigated. In the subcritical case, the flow becomes more unstable when $PrEc$ is increased. This is manifested by the enlargement of the neutral curve. Similarly, the transcritical case becomes more unstable as $PrEc$ increases. However, once $PrEc$ reaches the critical value (in this case $PrEc=0.05115$ ), the base flow becomes inflectional. The flow is thus inviscid unstable and the critical Reynolds number is substantially reduced. For instance, the flow is unstable for $Re<1000$ and $PrEc=0.06$ . In the supercritical case, the increase of $PrEc$ stabilizes the base flow and the non-ideal gas is even more stable than the ideal gas. In this case, when $PrEc$ reaches 0.03, the modal instability is found after $Re>8000$ . Interestingly, a weak influence of $PrEc$ on the velocity perturbations is observed (see figure 7 d–f), while the amplitudes of density and temperature perturbations are considerably larger if $PrEc$ increases. For the transcritical case, when the flow enters the triangular zone, the density perturbations are dominant.

Figure 7.

Neutral curves and profiles of perturbations for the non-ideal gas (RP model) and ideal gas (IG model). (a,d) $T_{w}^{\ast }=290~\text{K}$ , (b,e) $T_{w}^{\ast }=300~\text{K}$ , (c,f) $T_{w}^{\ast }=310~\text{K}$ . The neutral curves are obtained for two-dimensional perturbations ( $\unicode[STIX]{x1D6FD}=0$ ). The profiles shown are subject to $\unicode[STIX]{x1D6FC}=1$ , $\unicode[STIX]{x1D6FD}=0$ and $Re=10\,000$ (close to the most unstable area of the neutral curves), and they are normalized with $|u^{\prime }|_{max}$ . The left and right halves show the non-ideal and ideal gas respectively.

Below, we compare the fidelity of the cubic EoS models with the EoS model from REFPROP. The solutions for the ideal EoS are also shown to highlight the difference with respect to the results obtained with the non-ideal EoS models. Figure 8 shows the growth rate of the unstable modes for all EoS models. Recall the discussion in § 4.1, all these curves collapse under the isothermal limit. As can be inferred from each row of figure 8, the differences between these models magnify when $PrEc$ is increased. In all three cases, the cubic EoS models predict the correct trend that the flow becomes more unstable in sub-/transcritical cases, and more stable in supercritical cases as $PrEc$ increases.

Specifically, the van der Waals EoS shows a good agreement with the RP EoS model in the subcritical case, while both the Peng–Robinson and Redlich–Kwong EoS predict a lower growth rate (shown in figure 8). In the transcritical case, the van der Waals and Redlich–Kwong EoS give acceptable growth rates if compared to RP. When the base flow becomes inflectional ( $PrEc=0.06$ ), the Peng–Robinson EoS shows the best approximation. In the supercritical case, Redlich–Kwong produces the best results, while the van der Waals EoS gives a much lower growth rate. Given these observations, it can be concluded that all non-ideal EoS models give the same trends. However, it is not possible to conclude on the fidelity of the cubic EoS models in terms of the growth rate.

Figure 8.

Growth rates of the perturbation for different gas models. Results shown are at $Re=10\,000$ , $\unicode[STIX]{x1D6FD}=0$ for the subcritical case (a–d, $PrEc=$ 0.01, 0.03, 0.05 and 0.07), transcritical case (e–h, $PrEc=$ 0.01, 0.03, 0.05 and 0.06) and supercritical case (i–l, $PrEc=$ 0.01, 0.015, 0.02 and 0.03). Note that the $y$ -coordinate of (h) is different from the others.

4.3 The kinetic energy budget

To further understand the instability mechanism of the non-ideal fluids, we perform a kinetic energy budget analysis for the two-dimensional perturbation. The energy balance equation is the sum of the $x$ -momentum perturbation equation, multiplied with $\hat{u} ^{\dagger }$ , and the $y$ - equation, multiplied with $\hat{v}^{\dagger }$ . Here, dagger stands for the complex conjugate. The continuity equation is used to substitute the temporal growth of the density, which appears in the $x$ -momentum equation. This gives the following kinetic energy balance equation:

(4.2) $$\begin{eqnarray}\displaystyle K=\unicode[STIX]{x1D6E9}+P+T+V, & & \displaystyle\end{eqnarray}$$

where

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle K=-\text{i}\unicode[STIX]{x1D714}\int \unicode[STIX]{x1D70C}_{0}(\hat{u} \hat{u} ^{\dagger }+\hat{v}\hat{v}^{\dagger })\,\text{d}y, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E9}=-\text{i}\unicode[STIX]{x1D6FC}\int \unicode[STIX]{x1D70C}_{0}u_{0}(\hat{u} \hat{u} ^{\dagger }+\hat{v}\hat{v}^{\dagger })\,\text{d}y, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle P=-\int \unicode[STIX]{x1D70C}_{0}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}\hat{v}\hat{u} ^{\dagger }\,\text{d}y, & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle T & = & \displaystyle -\int \left[\text{i}\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}\hat{\unicode[STIX]{x1D70C}}\hat{u} ^{\dagger }+\text{i}\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}T_{0}}\hat{T}\hat{u} ^{\dagger }+\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}y}\hat{v}^{\dagger }+\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}T_{0}}\frac{\unicode[STIX]{x2202}\hat{T}}{\unicode[STIX]{x2202}y}\hat{v}^{\dagger }\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\left(\frac{\unicode[STIX]{x2202}^{2}p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}^{2}p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}T_{0}}\frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}y}\right)\hat{\unicode[STIX]{x1D70C}}\hat{v}^{\dagger }+\left(\frac{\unicode[STIX]{x2202}^{2}p_{0}}{\unicode[STIX]{x2202}T_{0}^{2}}\frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}^{2}p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}T_{0}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}y}\right)\hat{T}\hat{v}^{\dagger }\right]\text{d}y,\qquad\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle V & = & \displaystyle \frac{1}{Re}\int \left[-\unicode[STIX]{x1D6FC}^{2}(2\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D706}_{0})\hat{u} \hat{u} ^{\dagger }+\unicode[STIX]{x1D707}_{0}\frac{\unicode[STIX]{x2202}^{2}\hat{u} }{\unicode[STIX]{x2202}y^{2}}\hat{u} ^{\dagger }+\text{i}\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D706}_{0})\frac{\unicode[STIX]{x2202}\hat{v}}{\unicode[STIX]{x2202}y}\hat{u} ^{\dagger }\right.\nonumber\\ \displaystyle & & \displaystyle +\,\text{i}\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}y}\hat{v}\hat{u} ^{\dagger }+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}\frac{\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}y}\hat{u} ^{\dagger }+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}y}\frac{\unicode[STIX]{x2202}\hat{u} }{\unicode[STIX]{x2202}y}\hat{u} ^{\dagger }+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}T_{0}}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}\frac{\unicode[STIX]{x2202}\hat{T}}{\unicode[STIX]{x2202}y}\hat{u} ^{\dagger }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}^{2}u_{0}}{\unicode[STIX]{x2202}y^{2}}\hat{\unicode[STIX]{x1D70C}}\hat{u} ^{\dagger }+\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}T_{0}}\frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}y}\right)\hat{\unicode[STIX]{x1D70C}}\hat{u} ^{\dagger }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}T_{0}}\frac{\unicode[STIX]{x2202}^{2}u_{0}}{\unicode[STIX]{x2202}y^{2}}\hat{T}\hat{u} ^{\dagger }+\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}T_{0}^{2}}\frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}T_{0}\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}y}\right)\hat{T}\hat{u} ^{\dagger }\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D707}_{0}\hat{v}\hat{v}^{\dagger }+(2\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D706}_{0})\frac{\unicode[STIX]{x2202}^{2}\hat{v}}{\unicode[STIX]{x2202}y^{2}}\hat{v}^{\dagger }+\text{i}\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D706}_{0})\frac{\unicode[STIX]{x2202}\hat{u} }{\unicode[STIX]{x2202}y}\hat{v}^{\dagger }\nonumber\\ \displaystyle & & \displaystyle +\left.\text{i}\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}\hat{\unicode[STIX]{x1D70C}}\hat{v}^{\dagger }+\text{i}\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}_{0}}{\unicode[STIX]{x2202}y}\hat{u} \hat{v}^{\dagger }+\text{i}\unicode[STIX]{x1D6FC}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}T_{0}}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}\hat{T}\hat{v}^{\dagger }+\left(\!2\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}_{0}}{\unicode[STIX]{x2202}y}\right)\frac{\unicode[STIX]{x2202}\hat{v}}{\unicode[STIX]{x2202}y}\hat{v}^{\dagger }\right]\text{d}y.\qquad\end{eqnarray}$$

The real part of (4.2) describes the balance of the kinetic energy growth. In particular, $K_{r}$ is the temporal growth of the kinetic energy, $\unicode[STIX]{x1D6E9}$ is purely imaginary and does therefore not contribute to the temporal growth, $P_{r}$ is the production term, $T_{r}$ is the thermodynamic term and $V_{r}$ is the viscous dissipation.

The results of the kinetic energy budget analysis are summarized in table 3. The analysis is performed for all three cases at $\unicode[STIX]{x1D6FC}=1$ and $Re=10\,000$ . It clearly shows that, for all the cases, the energy growth $K_{r}$ originates from the production term $P_{r}$ . The thermodynamic term $T_{r}$ slightly reduces the growth. The viscous dissipation $V_{r}$ is not sensitive to the parameters and remains almost constant, except in the transcritical case ( $T_{w}^{\ast }=300~\text{K}$ , $PrEc=0.06$ ), which has a considerably larger growth rate (as also shown in figures 7 and 8). The reason for this lies in a much larger production and a smaller viscous dissipation. Figure 9 compares the production of the two cases with $PrEc=0.05$ and 0.06 at $T_{w}^{\ast }=300~\text{K}$ . It can be inferred that the inflectional velocity profile ( $PrEc=0.06$ ) has caused a larger $\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}u_{0}/\unicode[STIX]{x2202}y$ near both walls, the amplitude of the velocity perturbation $\hat{v}\hat{u} ^{\dagger }$ is larger as well. Therefore, a large production term $-\int \unicode[STIX]{x1D70C}_{0}(\unicode[STIX]{x2202}u_{0}/\unicode[STIX]{x2202}y)\hat{v}\hat{u} ^{\dagger }\,\text{d}y$ and accordingly the large growth rate can be explained.

Table 3.

Kinetic energy budget analysis for two-dimensional perturbations. $\unicode[STIX]{x1D6FC}=1$ , $Re=10\,000$ .

Figure 9.

Production of the kinetic perturbation energy with $T_{w}^{\ast }=300~\text{K}$ , $\unicode[STIX]{x1D6FC}=1$ , $Re=10\,000$ . (a) $PrEc=0.05$ , (b) $PrEc=0.06$ .

5.1 Choice of the energy norm

Mack’s energy norm (Mack Reference Mack1969; Hanifi, Schmid & Henningson Reference Hanifi, Schmid and Henningson1996) has been extensively used in compressible flows. The norm is designed under the ideal gas assumption, therefore the pressure-related energy transfer terms can be eliminated by choosing suitable coefficients for each component. In fact, Mack’s norm is equivalent to Chu’s norm (Chu Reference Chu1965; George & Sujith Reference George and Sujith2011). In the current non-ideal gas flows, the EoS can be different (PR, RK, VW, IG), or even implicit (the look-up table method) as in the case of the RP EoS model. Therefore, we choose a general form of the norm:

(5.1) $$\begin{eqnarray}\displaystyle E(\boldsymbol{q})=\int (u^{\prime \dagger }u^{\prime }+v^{\prime \dagger }v^{\prime }+w^{\prime \dagger }w^{\prime })+m_{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D70C}^{\prime \dagger }\unicode[STIX]{x1D70C}^{\prime }+m_{T}T^{\prime \dagger }T^{\prime }\,\text{d}V, & & \displaystyle\end{eqnarray}$$

where $\dagger$ denotes the complex conjugate. This norm has been tested for the compressible ideal/non-ideal gas flows at various conditions. Figure 10 shows the optimal energy growth $G_{max}$ (the maximum of $G(t)$ over time $t$ ) as a function of $m_{T}$ and $m_{\unicode[STIX]{x1D70C}}$ , for $PrEc=0.05$ (thermodynamic components become important) and a wall temperature of $T_{w}^{\ast }=290~\text{K}$ . When $m_{\unicode[STIX]{x1D70C}}$ is set to 0, $G_{max}$ converges to a constant value when $m_{T}$ is large enough. On the other hand, the energy norm is shown to be rather robust when the density component is properly accounted for, e.g.  $m_{\unicode[STIX]{x1D70C}}=1$ . Therefore, the results presented in this section are mainly obtained for $m_{\unicode[STIX]{x1D70C}}=m_{T}=1$ . A comparison with Mack’s energy norm ( $m_{\unicode[STIX]{x1D70C}}=T_{0}/(\unicode[STIX]{x1D70C}_{0}^{2}\unicode[STIX]{x1D6FE}Ma^{2})$ , $m_{T}=1/(\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D6FE}-1)T_{0}Ma^{2})$ ) is provided at the end of this section.

Figure 10.

Maximum energy growth $G_{max}$ versus $m_{T}$ using the energy norm (5.1).  $m_{\unicode[STIX]{x1D70C}}=0$ , 1, 2 and 10. The fluid (with RP model) is at $PrEc=0.05$ , $Re=2000$ , $\unicode[STIX]{x1D6FC}=1.0$ , $\unicode[STIX]{x1D6FD}=0.25$ and $T_{w}^{\ast }=290~\text{K}$ .

5.2 The isothermal limit

Although all EoS considered in this work give the same most unstable mode in the isothermal limit (discussed in § 4.1), their eigenvalue spectra can be rather different (see figure 6 a). Their corresponding eigenfunctions form the basis of the optimal perturbation and the algebraic growth. We show the contour plot of $G_{max}$ in $\unicode[STIX]{x1D6FC}$ – $\unicode[STIX]{x1D6FD}$ diagram in figure 11(a). Lines and circles show results of RP and IG models, respectively. It is evident that they fall on top of each other. In fact, all five models (RP, PR, RK, VW, IG) show the same results, and correspond to the results using incompressible equations. The largest transient growth occurs at $\unicode[STIX]{x1D6FC}=0$ and $2\leqslant \unicode[STIX]{x1D6FD}\leqslant 2.1$ , which is well known for an ideal gas. The optimal perturbation and the corresponding output are shown in figure 11(b,c) for $\unicode[STIX]{x1D6FC}=0$ , $\unicode[STIX]{x1D6FD}=2$ . The classic streamwise vortices (the optimal perturbation) and streaks (the corresponding output) are recovered. There is no discernible difference between the non-ideal and ideal gases under the isothermal limit.

Figure 11.

Transient growth in the isothermal limit. $PrEc\rightarrow 0$ , $Re=2000$ , $T_{w}^{\ast }=290~\text{K}$ . (a) Contour plot of $G_{max}$ . (b) The optimal perturbation (input) for $\unicode[STIX]{x1D6FC}=0$ , $\unicode[STIX]{x1D6FD}=2$ . (c) The corresponding output. Lines and circle symbols show results of non-ideal (RP) and ideal gas (IG) respectively.

5.3 Compressible flows

The algebraic growth has been studied for the subcritical, transcritical and supercritical cases at $Re=1000$ and $PrEc=0.01$ , 0.03, 0.05, 0.07. The optimal energy growth $G_{max}$ for RP model is compared with the IG model in figures 12–14, respectively. The three cases actually start from the same results at the isothermal limit (figure 11 a). Regardless of the wall temperature and $PrEc$ , the largest transient growth occurs at $\unicode[STIX]{x1D6FC}=0$ and $2\leqslant \unicode[STIX]{x1D6FD}\leqslant 2.1$ for both ideal and non-ideal gases. In the subcritical and transcritical cases (figures 12 and 13), when $PrEc$ is increased, the ideal gas tends to be slightly more stable, while the non-ideal gas becomes more unstable. In fact, due to the power/Sutherland law (for the transport properties), the results for the ideal gas are weakly dependent on the wall temperature. Notably in figure 13(d), where $PrEc=0.07$ , an area of $G_{max}\rightarrow \infty$ stands out. Recall the discussion in § 4, the base flow has entered the triangular zone (see figure 5) and becomes inflectional. Hence, the flow is inviscid unstable and the critical $Re$ is reduced considerably (see figure 7 c). As a result, a sub-zone of modal growth (near $\unicode[STIX]{x1D6FD}=0$ ) in the $\unicode[STIX]{x1D6FC}{-}\unicode[STIX]{x1D6FD}$ diagram is observed (where $G_{max}\rightarrow \infty$ ). For better display of the results, we have limited the colour band to $G_{max}=450$ in figure 13. In the supercritical case (figure 14), the plots are almost symmetrical, indicating the non-ideal gas effects are rather insignificant. The non-ideal gas is only slightly more unstable than the ideal gas. Table 4 summarizes the above maximum transient growth $G_{max}$ . With the increase in $PrEc$ , a similar trend as for the modal growth can be observed. Namely, the ideal gas becomes more stable, while the non-ideal gas tends to be more unstable for the subcritical and transcritical cases, and more stable for the supercritical case. On the whole, the non-ideal gas effects increase the algebraic instability in all regimes, most prominently in the transcritical regime.

Figure 12.

Contour plot of $G_{max}$ at $T_{w}^{\ast }=290~\text{K}$ . $Re=1000$ . On the left- and right-hand side of each panel, we show the results for non-ideal (RP) and ideal (IG) gases respectively. (a) $PrEc=0.01$ , (b) $PrEc=0.03$ , (c) $PrEc=0.05$ , (d) $PrEc=0.07$ .

Figure 13.

Same as figure 12 but for $T_{w}^{\ast }=300~\text{K}$ .

Figure 14.

Same as figure 12 but for $T_{w}^{\ast }=310~\text{K}$ .

Table 4.

Maximum transient growth $G_{max}$ of the perturbations at $Re=1000$ .

Figure 15.

Optimal perturbations (a) and the resulting output (b). $PrEc=0.07$ , $\unicode[STIX]{x1D6FC}=0$ and $\unicode[STIX]{x1D6FD}=2$ . Only significant components are plotted, namely in (a) $|v^{\prime }|$ , $|w^{\prime }|$ , (b) $|u^{\prime }|$ , $|\unicode[STIX]{x1D70C}^{\prime }|$ and $|T^{\prime }|$ .

The typical optimal perturbation and the resulting output are shown in figure 15 at $PrEc=0.07$ , $\unicode[STIX]{x1D6FC}=0$ and $\unicode[STIX]{x1D6FD}=2$ . Similar to an incompressible flow, the streamwise vortices and velocity streaks are recovered as the optimal perturbation and the output, respectively. For compressible flows, thermal streaks ( $\unicode[STIX]{x1D70C}^{\prime }$ and $T^{\prime }$ ) also become significant. Considering the non-ideal gas effects, the subcritical and supercritical cases share similar optimal perturbations as the ideal gas. In the transcritical case, the profiles are strongly influenced by the inflectional base flow and the strong property variations. On the other hand, the output perturbations are almost the same with regard to the $u^{\prime }$ component, indicating similar dynamic streaks are being generated. The amplitude of the thermal streak is much larger in the transcritical case close to the wall.

We have shown in § 4.2 that cubic EoS cannot guarantee accurate results for the growth rate if compared to results obtained with the accurate REFPROP EoS. This is also true for the algebraic instability as shown in figure 16, depicting $G$ – $t$ curves of the three cases with different EoS at $PrEc=0.07$ , $\unicode[STIX]{x1D6FC}=0$ and $\unicode[STIX]{x1D6FD}=2$ . For example, the van der Waals EoS over-predicts $G_{max}$ by 270 % for the transcritical case. In the supercritical case, the non-ideal gas effects are less significant, and the results of all considered EoS are close to each other.

Figure 16.

The transient amplification curve $G(t)$ at $PrEc=0.07$ , $\unicode[STIX]{x1D6FC}=0$ and $\unicode[STIX]{x1D6FD}=2$ . (a) $T_{w}^{\ast }=290~\text{K}$ , (b) $T_{w}^{\ast }=300~\text{K}$ , (c) $T_{w}^{\ast }=310~\text{K}$ .

The main results presented in this section are based on the energy norm: $m_{\unicode[STIX]{x1D70C}}=m_{T}=1$ . When Mack’s energy norm is used, figure 17 provides a comparison for all three regimes with highly non-ideal gas effects ( $PrEc=0.07$ , $\unicode[STIX]{x1D6FC}=0$ ). Indeed, the non-ideal gas has a larger algebraic growth in all three cases with Mack’s energy norm, while on the other hand, the ideal gas are rather insensitive to different norms. As a result, the conclusion on algebraic growth will not change, the flow is more unstable compared with the ideal gas in all three regimes (sub-, trans- and supercritical).

Figure 17.

Comparison of maximum algebraic growth using Mack’s energy norm. Cases with (a) $T_{w}^{\ast }=290~\text{K}$ , (b) $T_{w}^{\ast }=300~\text{K}$ and (c) $T_{w}^{\ast }=310~\text{K}$ . The other parameters are kept constant: $\unicode[STIX]{x1D6FC}=0$ , $Re=1000$ , $PrEc=0.07$ .