Linear instability of Poiseuille flows with highly non-ideal fluids

Abstract

The objective of this work is to investigate linear modal and algebraic instability in Poiseuille flows with fluids close to their vapour–liquid critical point. Close to this critical point, the ideal gas assumption does not hold and large non-ideal fluid behaviours occur. As a representative non-ideal fluid, we consider supercritical carbon dioxide ($\text{CO}_{2}$) at a pressure of 80 bar, which is above its critical pressure of 73.9 bar. The Poiseuille flow is characterized by the Reynolds number ($Re=\unicode[STIX]{x1D70C}_{w}^{\ast }u_{r}^{\ast }h^{\ast }/\unicode[STIX]{x1D707}_{w}^{\ast }$), the product of the Prandtl ($Pr=\unicode[STIX]{x1D707}_{w}^{\ast }C_{pw}^{\ast }/\unicode[STIX]{x1D705}_{w}^{\ast }$) and Eckert numbers ($Ec=u_{r}^{\ast 2}/C_{pw}^{\ast }T_{w}^{\ast }$) and the wall temperature that in addition to pressure determine the thermodynamic reference condition. For low Eckert numbers, the flow is essentially isothermal and no difference with the well-known stability behaviour of incompressible flows is observed. However, if the Eckert number increases, the viscous heating causes gradients of thermodynamic and transport properties, and non-ideal gas effects become significant. Three regimes of the laminar base flow can be considered: the subcritical (temperature in the channel is entirely below its pseudo-critical value), transcritical and supercritical temperature regimes. If compared to the linear stability of an ideal gas Poiseuille flow, we show that the base flow is modally more unstable in the subcritical regime, inviscid unstable in the transcritical regime and significantly more stable in the supercritical regime. Following the principle of corresponding states, we expect that qualitatively similar results will be obtained for other fluids at equivalent thermodynamic states.

1 Introduction

Many processes in industrial applications constitute flows with fluids that do not follow the ideal gas law. For example, flows in vapour power systems, re-entry of spacecraft, supercritical dyeing, refrigeration and heat pump systems (examples in supercritical fluids can be found in Brunner (2010)). The non-ideality of fluids is especially significant in the thermodynamic region close to the vapour critical point. As such, it is of great importance to understand the fundamental physics that are related to flows with these fluids.

Recently, researchers have studied how non-ideal gas effects influence turbulence and heat transfer. For example, Kawai, Terashima & Negishi (2015), Kawai (2016) studied turbulent boundary layers with supercritical pressures and transcritical temperatures. They found that the mean velocity profiles (with density weighted van Driest transformation) coincide with the same log law as seen in an ideal gas. Sciacovelli et al. (2016), Sciacovelli, Cinnella & Gloerfelt (2017a ), Sciacovelli, Cinnella & Grasso (2017b ) comprehensively studied turbulence dynamics and near-wall turbulence of flows with molecularly complex fluids in the dense gas regime using direct numerical simulations. They found that dense gas flows with a heavy fluorocarbon exhibit almost negligible friction heating (in channel flows), weakening of compressive (and enhancement of expanding) structures (in homogeneous isotropic turbulence). Patel, Boersma & Pecnik (2016) studied the influence of variable properties on fully developed turbulent channel flows and derived a velocity transformation that allows us to collapse velocity profiles for heated or cooled non-ideal fluids. Moreover, Rinaldi et al. (2017) provided an explanation of near wall turbulence modulation, especially the intercomponent energy transfer that has been observed by, e.g. Morinishi, Tamano & Nakabayashi (2004), Pirozzoli, Bernardini & Grasso (2008), Duan, Beekman & Martin (2010). Nemati et al. (2016), Peeters et al. (2016) studied turbulent heat transfer to supercritical $\text{CO}_{2}$ , indicating that both the mean and instantaneous property variations have significant effects on turbulent structures and their self-regeneration processes in near-wall turbulence. Alferez & Touber (2017) have studied the refraction properties of compression shock waves in non-ideal gases. One of the new regimes found is that, due to the non-ideality of the fluid, it is possible that acoustic modes can be completely damped by a compression shock, leading to so-called ‘quiet shocks’.

For ideal gases, the thermodynamic properties are associated with a simple equation of state (EoS). Additionally, the transport properties (namely, the viscosity and thermal conductivity) can be estimated as unary functions of the temperature (e.g. the widely used power law or Sutherland’s law). To assess to which degree the ideal gas law holds, it is possible to evaluate the compressibility factor, defined as $Z=p^{\ast }/\unicode[STIX]{x1D70C}^{\ast }R^{\ast }T^{\ast }$ . Figure 1 shows the $T{-}\unicode[STIX]{x1D717}$ diagram (temperature–specific-volume diagram, $\unicode[STIX]{x1D717}=1/\unicode[STIX]{x1D70C}$ ) of carbon dioxide ( $\text{CO}_{2}$ ). The white circle in each panel indicates the critical point, which for $\text{CO}_{2}$ is at a pressure and temperature of $p_{c}^{\ast }=73.9$ bar and $T_{c}^{\ast }=304.25~\text{K}$ . In this paper, we denote dimensional and critical quantities with the superscript ‘ $\ast$ ’ and subscript ‘ $c$ ’, respectively. Figure 1(a) shows the critical isobar (black thin dashed line), four isobars of 40 to 100 bar (black thin lines) and the compressibility factor $Z$ (coloured contour lines). Close to the critical point, the non-ideality is clearly indicated by low values of $Z$ , while the boundary between ideal and non-ideal gas behaviour is approximately indicated by the thick dashed line of $Z=0.99$ . The distribution of the thermodynamic and transport properties (specific heat capacity at constant volume $C_{v}^{\ast }$ , dynamic viscosity $\unicode[STIX]{x1D707}^{\ast }$ and thermal conductivity $\unicode[STIX]{x1D705}^{\ast }$ ) are shown in figure 1(b,c,d). In the ideal gas region, these contour lines become quasi-parallel to the $x$ -axis, indicating that they can be regarded as functions of temperature only. On the other hand, near the critical point, the gradients of these properties with respect to temperature and specific volume become significant.

Figure 1. $T{-}\unicode[STIX]{x1D717}$ diagram of $\text{CO}_{2}$ along with the critical point (white circle), the saturation curves (blue and red solid lines), the liquid–vapour region (grey area), the critical isobar (black thin dashed line), isobars of 40, 60, 80 and 100 bar (black thin lines), compressibility factor $Z=0.99$ (black thick dashed line) as well as coloured contour lines of (a) $Z=p^{\ast }/\unicode[STIX]{x1D70C}^{\ast }R^{\ast }T^{\ast }$ , (b) the specific heat capacity at constant volume $C_{v}^{\ast }$ , (c) the dynamic viscosity $\unicode[STIX]{x1D707}^{\ast }$ , and (d) the thermal conductivity $\unicode[STIX]{x1D705}^{\ast }$ .

In view of its great simplicity, most of the present knowledge on stability and laminar–turbulent transition is limited to the ideal gas (Fedorov 2011) or incompressible flows, where thermodynamic properties are constant. On the other hand, numerical simulations of real gas effects (high-enthalpy effects) in hypersonic flows have just gone through an initial stage (Zhong & Wang 2012; Marxen et al. 2013; Marxen, Iaccarino & Magin 2014). In fact, the well-known Orr–Sommerfeld equation (Orr 1907; Sommerfeld 1908, often termed the O-S equation) was derived by applying the linear stability theory (LST) to the incompressible parallel plane shear flow. Solved as an eigenvalue problem (in the time/space-asymptotic limit), the growth rate and profiles of the perturbation are obtained from the most unstable mode as its eigenvalue and eigenvector. This is known as the modal stability problem. The critical Reynolds number $Re_{c}$ , below which the flow is stable, regardless of the wavenumber and frequency of the perturbation, is often determined and emphasized in such modal stability analysis. For example, in plane Poiseuille flow, $Re_{c}$ is numerically determined to be 5772.22 (Thomas 1953; Orszag 1971). Here the Reynolds number is based on the half-channel height and the centreline flow velocity. Due to the non-normality of most practical linear systems, the modal stability analysis cannot cover the full behaviour of the linear instability (Schmid & Henningson 2001; Schmid & Brandt 2014). Instead of solving the eigenvalue problem, the stability equation can be formulated as an initial value problem under the framework of constrained optimization. Maximizing the energy growth in a finite domain of time or space, leads to the optimal perturbation, which grows transiently even below the critical Reynolds number $Re_{c}$ . This is termed the transient growth or algebraic growth. Accordingly, a ‘critical’ Reynolds number can be defined for the algebraic growth as well. For plane Poiseuille flow, this number is 49.6 (Busse 1969; Joseph & Carmi 1969).

Studies of viscosity stratified flows, where the viscosity depends on temperature, has recently received attention, readers may refer to Govindarajan & Sahu (2014) for a review. Based on the modified O-S equations, early studies show that including a linear temperature profile destabilizes the Poiseuille flow (Potter & Graber 1972) and stabilizes/destabilizes the water boundary layer flow (depend on wall heating/cooling) (Wazzan et al. 1972). However, viscosity and temperature perturbations were ignored in both studies and were later examined by Pinarbasi & Liakopoulos (1995). Wall & Wilson (1996, 1997) investigated the effects of different viscosity models, indicating that the flow can either be more stable or unstable. The study on wall heating and viscosity stratification has also been extended to transient growth, secondary instability as well as instabilities in other types of flows (Chikkadi, Sameen & Govindarajan 2005; Sameen & Govindarajan 2007; Sahu 2011; Sameen, Bale & Govindarajan 2011; Sahu & Govindarajan 2014). For compressible plane Couette flow, Malik, Dey & Alam (2008) showed that the flow is more stable with viscosity stratification, while recently, a further study on this flow is given by Saikia et al. (2017), in which the effects of individual/combined viscosity–thermal conductivity stratification are elucidated. The influence of viscosity gradients on the edge state has recently been studied by Rinaldi, Schlatter & Bagheri (2018), showing that in minimal channel flows, the kinetic energy level and the driving force of self-sustained cycle of the edge state depend on viscosity distribution. The above studies are based on the incompressible flow assumption or the ideal gas EoS; at the same time, transport properties are estimated as functions of temperature only.

Since there is very limited knowledge on flow stability with highly non-ideal fluids, we investigate Poiseuille flows with fluids close to the thermodynamic vapour–liquid critical point. In § 2, the gas model, the formulation of the stability analysis and the related numerical methods are outlined in detail. The results and discussions on the base flow are provided in § 3, followed by the modal growth and algebraic instability in §§ 4 and 5, respectively. The paper is concluded in § 6.

2 Governing equations

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 2013) 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 (2003) 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 1873), Redlich–Kwong (Redlich & Kwong 1949) and Peng–Robinson (Peng & Robinson 1976), 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 (2001), 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 (2001), Schmid & Brandt (2014).

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 The laminar base flow

The base flow is driven by a constant body force in the streamwise direction and is obtained by solving (2.1)–(2.3) with the assumption that the flow is fully developed, spanwise and streamwise independent, steady and parallel, i.e.  $\unicode[STIX]{x2202}(\,)/\unicode[STIX]{x2202}x=0,\unicode[STIX]{x2202}(\,)/\unicode[STIX]{x2202}z=0,\unicode[STIX]{x2202}(\,)/\unicode[STIX]{x2202}t=0,v=w=0$ . The N–S equations are thus simplified as

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left(\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}\right)=-ReF=-\hat{F}, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}y}=0, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left(\frac{\unicode[STIX]{x1D705}}{PrEc}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}y}+\unicode[STIX]{x1D707}u\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}y}\right)=-ReF\cdot u=-\hat{F}\cdot u. & \displaystyle\end{eqnarray}$$

It is worth noting that the above equations are independent of density, therefore, $\unicode[STIX]{x1D70C}_{0}$ can be separately determined by the EoS. We assume the body force, $\hat{F}$ , which drives the flow, to be uniform. To obtain a solution of the base flow, an initial temperature field is assumed, e.g.  $T=T_{w}=$ const., $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D705}$ are determined from REFPROP according to the temperature and pressure. First, the velocity is solved using equation (3.1), followed by an update of temperature by solving (3.3). Values of $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D705}$ are then updated using the obtained temperature. This procedure is repeated until the solution is converged.

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 Linear modal instability

Depending on the cases discussed below, we will use the definition of dynamic and thermodynamic modes, as

(4.1) $$\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70C}^{\prime }=0,\text{ and}\quad T^{\prime }=0\quad & \text{(dynamic modes)}\\ \unicode[STIX]{x1D70C}^{\prime }\neq 0,\text{ or}\quad T^{\prime }\neq 0\quad & \text{(thermodynamic modes).}\end{array}\right. & & \displaystyle\end{eqnarray}$$

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 (2001, 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 (1976)) 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 Algebraic growth

5.1 Choice of the energy norm

Mack’s energy norm (Mack 1969; Hanifi, Schmid & Henningson 1996) 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 1965; George & Sujith 2011). 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$ .

6 Conclusion

Linear stability of highly non-ideal plane Poiseuille flows is studied. We have chosen carbon dioxide ( $\text{CO}_{2}$ ) at supercritical pressures ( $p^{\ast }=80$ bar) as an example of a fluid in a highly non-ideal thermodynamic region. The investigation is based on the fully compressible Navier–Stokes equations in which the product of two dimensionless parameters, namely the Prandtl $Pr$ and Eckert $Ec$ numbers, determines the viscous heating and consequently the temperature increase between the two isothermal walls. The investigated range of $PrEc$ is from the isothermal limit ( $PrEc\rightarrow 0$ ) to typical compressible flows with $PrEc=0.1$ . Three cases with wall temperatures in the vicinity of the pseudo-critical point ( $T_{pc}^{\ast }=307.7~\text{K}$ ) have been investigated in more detail. In particular, the wall temperatures are set such that the temperature profile is subcritical ( $T_{w}^{\ast }=290~\text{K}$ ), transcritical ( $T_{w}^{\ast }=300~\text{K}$ ) and supercritical ( $T_{w}^{\ast }=310~\text{K}$ ). In all cases, the thermodynamic and transport properties are strongly dependent on the thermodynamic state of the fluid (e.g. temperature, density) and they influence the stability in a coupled way through the base flow and the linear stability operator.

In the isothermal limit, the three cases with different wall temperatures have the same base flow as the ideal gas. When $PrEc$ increases, the base flows of the three cases deviate from the ideal gas solution. In the subcritical regime, as $PrEc$ increases, the flow becomes more unstable with regard to both the modal and algebraic growth, while for ideal gas the trend is opposite. When $PrEc$ is large enough, or $T_{w}$ is closer to (but lower than) $T_{pc}$ , the flow falls in the transcritical regime. Due to the large gradient of the viscosity near $T_{pc}$ , the base flow becomes inflectional and inviscid unstable. As a consequence, the stability of the non-ideal gas flow is very different from the ideal gas. The neutral curve is expanded, which results in a very low critical Reynolds number. Moreover, the algebraic growth is also enhanced. It should be expected that the laminar–turbulent transition is considerably enhanced in this regime. When $T_{w}>T_{pc}$ , the fluid is in the thermodynamic supercritical regime. In this case, the results of the modal growth show that the non-ideal gas is substantially more stable than the ideal gas. However, the transient growth shows only a weak dependence on the non-ideal gas effects and is slightly more unstable than ideal gas flows. Additionally, we show that the linear stability analysis with simple cubic EoS gives qualitatively similar results to using the more accurate multi-parameter EoS implemented in the REFPROP library. Discussions on the reference scaling indicate that the conclusion is not influenced by the choice of the reference variables. This investigation constitutes the first systematic study of linear stability with highly non-ideal fluids close to the thermodynamic critical point. Future studies will focus on the validation of the results using direct numerical simulations.