Skip to main content Accessibility help


  • Access
  • Open access



MathJax is a JavaScript display engine for mathematics. For more information see
      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Linear instability of Poiseuille flows with highly non-ideal fluids
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Linear instability of Poiseuille flows with highly non-ideal fluids
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Linear instability of Poiseuille flows with highly non-ideal fluids
        Available formats
Export citation


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}$$


(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}$$


(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 (ac), density (df) and velocity (gi) 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 gi). 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 (ac) as well as eigenfunctions (df) 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 df), 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 (ad, $PrEc=$ 0.01, 0.03, 0.05 and 0.07), transcritical case (eh, $PrEc=$ 0.01, 0.03, 0.05 and 0.06) and supercritical case (il, $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}$$


(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 1214, 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.


J.R. and R.P. thank the Netherlands Organization for Scientific Research (NWO), which funded this research through the grant with project no. 14711. J.R. acknowledges financial support from the National Natural Science Foundation of China through grants 11602127 and 11572176.

Appendix A. The stability equation

The non-zero elements in the stability equation (2.12) are given below. For simplicity, the derivative of a thermodynamic quantity with respect to $\unicode[STIX]{x1D70C}_{0}$ at constant $T_{0}$ (and vice versa) has been denoted as $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}$ , instead of $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}|_{T_{0}}$ . The elements are:

(A 1) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle {\mathcal{L}}_{t}(1,1)=1,\\ \displaystyle {\mathcal{L}}_{t}(2,1)=u_{0},\\ \displaystyle {\mathcal{L}}_{t}(2,2)={\mathcal{L}}_{t}(3,3)={\mathcal{L}}_{t}(4,4)=\unicode[STIX]{x1D70C}_{0},\\ \displaystyle {\mathcal{L}}_{t}(5,1)=e_{0}+\unicode[STIX]{x1D70C}_{0}\frac{\unicode[STIX]{x2202}e_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}},\\ \displaystyle {\mathcal{L}}_{t}(5,5)=\unicode[STIX]{x1D70C}_{0}\frac{\unicode[STIX]{x2202}e_{0}}{\unicode[STIX]{x2202}T_{0}},\end{array}\right\} & & \displaystyle\end{eqnarray}$$
(A 2) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle {\mathcal{L}}_{x}(1,1)=u_{0},\\ \displaystyle {\mathcal{L}}_{x}(1,2)=\unicode[STIX]{x1D70C}_{0},\\ \displaystyle {\mathcal{L}}_{x}(2,1)=u_{0}u_{0}+\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}},\\ \displaystyle {\mathcal{L}}_{x}(2,2)=2\unicode[STIX]{x1D70C}_{0}u_{0},\\ \displaystyle {\mathcal{L}}_{x}(2,3)=-\frac{1}{Re}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{x}(2,5)=\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}T_{0}},\\ \displaystyle {\mathcal{L}}_{x}(3,1)=-\frac{1}{Re}\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},\\ \displaystyle {\mathcal{L}}_{x}(3,2)=-\frac{1}{Re}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{x}(3,3)={\mathcal{L}}_{x}(4,4)=\unicode[STIX]{x1D70C}_{0}u_{0},\\ \displaystyle {\mathcal{L}}_{x}(3,5)=-\frac{1}{Re}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}T_{0}}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{x}(5,1)=e_{0}u_{0}+\unicode[STIX]{x1D70C}_{0}u_{0}\frac{\unicode[STIX]{x2202}e_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}},\\ \displaystyle {\mathcal{L}}_{x}(5,2)=\unicode[STIX]{x1D70C}_{0}e_{0}+p_{0},\\ \displaystyle {\mathcal{L}}_{x}(5,3)=-\frac{2}{Re}\unicode[STIX]{x1D707}_{0}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{x}(5,5)=\unicode[STIX]{x1D70C}_{0}u_{0}\frac{\unicode[STIX]{x2202}e_{0}}{\unicode[STIX]{x2202}T_{0}},\end{array}\right\} & & \displaystyle\end{eqnarray}$$
(A 3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}{\mathcal{L}}_{y}(1,3)=\unicode[STIX]{x1D70C}_{0},\\ \displaystyle {\mathcal{L}}_{y}(2,1)=-\frac{1}{Re}\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},\\ \displaystyle {\mathcal{L}}_{y}(2,2)=-\frac{1}{Re}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{y}(2,3)=\unicode[STIX]{x1D70C}_{0}u_{0},\\ \displaystyle {\mathcal{L}}_{y}(2,5)=-\frac{1}{Re}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}T_{0}}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{y}(3,1)=\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}},\\ \displaystyle {\mathcal{L}}_{y}(3,3)=-\frac{2}{Re}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}y}-\frac{1}{Re}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{y}(3,5)=\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}T_{0}},\\ \displaystyle {\mathcal{L}}_{y}(4,4)=-\frac{1}{Re}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{y}(5,1)=-\frac{1}{RePrEc}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D705}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}y}\right),\\ \displaystyle {\mathcal{L}}_{y}(5,2)=-\frac{2}{Re}\unicode[STIX]{x1D707}_{0}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{y}(5,3)=\unicode[STIX]{x1D70C}_{0}e_{0}+p_{0},\\ \displaystyle {\mathcal{L}}_{y}(5,5)=-\frac{1}{RePrEc}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D705}_{0}}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D705}_{0}}{\unicode[STIX]{x2202}T_{0}}\frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}y}\right),\end{array}\right\} & & \displaystyle\end{eqnarray}$$
(A 4) $$\begin{eqnarray}\displaystyle & \left.\begin{array}{@{}c@{}}\displaystyle {\mathcal{L}}_{z}(1,4)=\unicode[STIX]{x1D70C}_{0},\\ \displaystyle {\mathcal{L}}_{z}(2,4)=\unicode[STIX]{x1D70C}_{0}u_{0},\\ \displaystyle {\mathcal{L}}_{z}(3,4)=-\frac{1}{Re}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D706}_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{z}(4,1)=\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}},\\ \displaystyle {\mathcal{L}}_{z}(4,3)=-\frac{1}{Re}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{z}(4,5)=\frac{\unicode[STIX]{x2202}p_{0}}{\unicode[STIX]{x2202}T_{0}},\\ \displaystyle {\mathcal{L}}_{z}(5,4)=\unicode[STIX]{x1D70C}_{0}e_{0}+p_{0},\end{array}\right\} & \displaystyle\end{eqnarray}$$
(A 5) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle {\mathcal{L}}_{q}(1,3)=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{q}(2,1)=-\frac{1}{Re}\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}}-\frac{1}{Re}\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),\\ \displaystyle {\mathcal{L}}_{q}(2,3)=\unicode[STIX]{x1D70C}_{0}\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}+u_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{q}(2,5)=-\frac{1}{Re}\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}}-\frac{1}{Re}\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)\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{q}(3,1)=\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},\\ \displaystyle {\mathcal{L}}_{q}(3,5)=\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},\\ \displaystyle {\mathcal{L}}_{q}(5,1)=-\frac{1}{RePrEc}\left(\frac{\unicode[STIX]{x2202}^{2}T_{0}}{\unicode[STIX]{x2202}y^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D705}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}+\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D705}_{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]{x1D705}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x2202}T_{0}}\frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}y}\right)\frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}y}\right)\\ \displaystyle -\,\frac{1}{Re}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}\left(\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}\right)^{2},\\ \displaystyle {\mathcal{L}}_{q}(5,3)=e_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}y}+\unicode[STIX]{x1D70C}_{0}\frac{\unicode[STIX]{x2202}e_{0}}{\unicode[STIX]{x2202}y},\\ \displaystyle {\mathcal{L}}_{q}(5,5)=-\frac{1}{RePrEc}\left(\frac{\unicode[STIX]{x2202}^{2}T_{0}}{\unicode[STIX]{x2202}y^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D705}_{0}}{\unicode[STIX]{x2202}T_{0}}+\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D705}_{0}}{\unicode[STIX]{x2202}T_{0}^{2}}\frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D705}_{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)\frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}y}\right)\\ \displaystyle -\,\frac{1}{Re}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x2202}T_{0}}\left(\frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}\right)^{2},\end{array}\right\} & & \displaystyle\end{eqnarray}$$
(A 6) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle {\mathcal{V}}_{xx}(2,2)={\mathcal{V}}_{yy}(3,3)={\mathcal{V}}_{zz}(4,4)=-\frac{2\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D706}_{0}}{Re},\\ \displaystyle {\mathcal{V}}_{xx}(3,3)={\mathcal{V}}_{xx}(4,4)=-\frac{\unicode[STIX]{x1D707}_{0}}{Re},\\ \displaystyle {\mathcal{V}}_{yy}(2,2)={\mathcal{V}}_{yy}(4,4)=-\frac{\unicode[STIX]{x1D707}_{0}}{Re},\\ \displaystyle {\mathcal{V}}_{zz}(2,2)={\mathcal{V}}_{zz}(3,3)=-\frac{\unicode[STIX]{x1D707}_{0}}{Re},\\ \displaystyle {\mathcal{V}}_{xx}(5,5)={\mathcal{V}}_{yy}(5,5)={\mathcal{V}}_{zz}(5,5)=-\frac{\unicode[STIX]{x1D705}_{0}}{RePrEc},\\ \displaystyle {\mathcal{V}}_{xy}(2,3)={\mathcal{V}}_{xy}(3,2)=-\frac{\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D706}_{0}}{Re},\\ \displaystyle {\mathcal{V}}_{xz}(2,4)={\mathcal{V}}_{xz}(4,2)=-\frac{\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D706}_{0}}{Re},\\ \displaystyle {\mathcal{V}}_{yz}(3,4)={\mathcal{V}}_{yz}(4,3)=-\frac{\unicode[STIX]{x1D707}_{0}+\unicode[STIX]{x1D706}_{0}}{Re}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The second-order finite differences were used to determine the gradients of the properties. For instance, the gradients of viscosity

(A 7) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}(T_{0},\unicode[STIX]{x1D70C}_{0})}{\unicode[STIX]{x2202}T}=\frac{\unicode[STIX]{x1D707}(T_{0}+\unicode[STIX]{x0394}T,\unicode[STIX]{x1D70C}_{0})-\unicode[STIX]{x1D707}(T_{0}-\unicode[STIX]{x0394}T,\unicode[STIX]{x1D70C}_{0})}{2\unicode[STIX]{x0394}T}, & \displaystyle\end{eqnarray}$$
(A 8) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}(T_{0},\unicode[STIX]{x1D70C}_{0})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}=\frac{\unicode[STIX]{x1D707}(T_{0},\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x0394}\unicode[STIX]{x1D70C})-\unicode[STIX]{x1D707}(T_{0},\unicode[STIX]{x1D70C}_{0}-\unicode[STIX]{x0394}\unicode[STIX]{x1D70C})}{2\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}, & \displaystyle\end{eqnarray}$$
(A 9) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}(T_{0},\unicode[STIX]{x1D70C}_{0})}{\unicode[STIX]{x2202}T\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}=\frac{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}(T_{0}+\unicode[STIX]{x0394}T,\unicode[STIX]{x1D70C}_{0})-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}(T_{0}-\unicode[STIX]{x0394}T,\unicode[STIX]{x1D70C}_{0})}{2\unicode[STIX]{x0394}T}. & \displaystyle\end{eqnarray}$$

An example of the sensitivity of $\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}^{\ast }/\unicode[STIX]{x2202}T^{\ast }\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\ast }$ to $\unicode[STIX]{x0394}T^{\ast }$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }$ is shown in figure 18. In fact, the gradients of the thermodynamic and transport properties became rather robust when $\unicode[STIX]{x0394}T^{\ast }\leqslant 1~\text{K}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }\leqslant 1~\text{kg}~\text{m}^{-3}$ . In this study, the results are obtained with $\unicode[STIX]{x0394}T^{\ast }=0.1~\text{K}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }=0.1~\text{kg}~\text{m}^{-3}$ .

Figure 18. Sensitivity of $\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D707}^{\ast }/\unicode[STIX]{x2202}T^{\ast }\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\ast }$ to $\unicode[STIX]{x0394}T^{\ast }$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }$ .

Appendix B. Cubic EoS

The material dependent parameters for CO $_{2}$ are provided in table 5. These parameters are necessary inputs for the cubic EoS detailed below and can be easily replaced for other fluids.

Table 5. The material dependent parameters for CO $_{2}$ .

B.1 The van der Waals EoS

The van der Waals (1873) EoS is the simplest cubic EoS that is capable of accounting phase separation and the critical point (see the introduction in Zappoli, Beysens & Garrabos 2015; Moran et al. 2012). The EoS can be written as

(B 1) $$\begin{eqnarray}\displaystyle p=\frac{R_{g}T}{\unicode[STIX]{x1D717}-b}-\frac{a}{\unicode[STIX]{x1D717}^{2}}, & & \displaystyle\end{eqnarray}$$

where $R_{g}$ is the specific gas constant, $a$ is a measure of the attraction forces between molecules and $b$ accounts for the finite volume occupied by the molecules. The constants $a$ and $b$ can be determined at the critical point where

(B 2) $$\begin{eqnarray}\displaystyle \left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D717}}\right|_{T_{c}}=\left.\frac{\unicode[STIX]{x2202}^{2}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D717}^{2}}\right|_{T_{c}}=0\quad \Rightarrow \quad a=\frac{27}{64}\frac{R_{g}^{2}T_{c}^{2}}{p_{c}},\quad b=\frac{R_{g}T_{c}}{8p_{c}}. & & \displaystyle\end{eqnarray}$$

Using the Maxwell relations and the departure function, it is possible to obtain the internal energy as

(B 3) $$\begin{eqnarray}\displaystyle e=C_{v}T-a\unicode[STIX]{x1D70C}. & & \displaystyle\end{eqnarray}$$

The required derivatives for stability equations are

(B 4a,b ) $$\begin{eqnarray}\displaystyle \left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}T}\right|_{\unicode[STIX]{x1D70C}}=\frac{\unicode[STIX]{x1D70C}R_{g}}{1-\unicode[STIX]{x1D70C}b},\quad \left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right|_{T}=\frac{R_{g}T}{(1-\unicode[STIX]{x1D70C}b)^{2}}-2a\unicode[STIX]{x1D70C}, & & \displaystyle\end{eqnarray}$$
(B 5a-c ) $$\begin{eqnarray}\displaystyle \left.\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right|_{T}\left(\left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}T}\right|_{\unicode[STIX]{x1D70C}}\right)=\frac{R_{g}}{(1-\unicode[STIX]{x1D70C}b)^{2}},\quad \left.\frac{\unicode[STIX]{x2202}^{2}p}{\unicode[STIX]{x2202}T^{2}}\right|_{\unicode[STIX]{x1D70C}}=0,\quad \left.\frac{\unicode[STIX]{x2202}^{2}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{2}}\right|_{T}=\frac{2R_{g}Tb}{(1-\unicode[STIX]{x1D70C}b)^{3}}-2a,\qquad & & \displaystyle\end{eqnarray}$$
(B 6a,b ) $$\begin{eqnarray}\displaystyle \left.\frac{\unicode[STIX]{x2202}e}{\unicode[STIX]{x2202}T}\right|_{\unicode[STIX]{x1D70C}}=C_{v}+\left.\frac{\unicode[STIX]{x2202}C_{v}}{\unicode[STIX]{x2202}T}\right|_{\unicode[STIX]{x1D70C}}T,\quad \left.\frac{\unicode[STIX]{x2202}e}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right|_{T}=-a. & & \displaystyle\end{eqnarray}$$

B.2 The Redlich–Kwong EoS

The Redlich–Kwong (Redlich & Kwong 1949) EoS is given as

(B 7) $$\begin{eqnarray}\displaystyle p=\frac{R_{g}T}{\unicode[STIX]{x1D717}-b}-\frac{a\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D717}(\unicode[STIX]{x1D717}+b)}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}=\sqrt{T_{c}/T}$ . Similarly, by satisfying the critical condition, the constants $a$ and $b$ are

(B 8a,b ) $$\begin{eqnarray}\displaystyle a=\frac{0.42748R_{g}T_{c}^{2}}{p_{c}},\quad b=\frac{0.08664R_{g}T_{c}}{p_{c}}. & & \displaystyle\end{eqnarray}$$

The internal energy is

(B 9) $$\begin{eqnarray}\displaystyle e=C_{v}T+\frac{3}{2}\frac{a}{b}\unicode[STIX]{x1D6FC}\ln \frac{1}{1+b\unicode[STIX]{x1D70C}}. & & \displaystyle\end{eqnarray}$$

The derivatives in the stability equations are

(B 10) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}T}\right|_{\unicode[STIX]{x1D70C}}=\frac{\unicode[STIX]{x1D70C}R_{g}}{1-b\unicode[STIX]{x1D70C}}+\frac{1}{2}T^{-(3/2)}T_{c}^{1/2}\frac{a\unicode[STIX]{x1D70C}^{2}}{1+b\unicode[STIX]{x1D70C}},\\ \displaystyle \left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right|_{T}=\frac{R_{g}T}{(1-\unicode[STIX]{x1D70C}b)^{2}}-T^{-(1/2)}T_{c}^{1/2}\frac{2a\unicode[STIX]{x1D70C}+ab\unicode[STIX]{x1D70C}^{2}}{(1+\unicode[STIX]{x1D70C}b)^{2}},\end{array}\right\} & & \displaystyle\end{eqnarray}$$
(B 11) $$\begin{eqnarray}\displaystyle \left.\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right|_{T}\left(\left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}T}\right|_{\unicode[STIX]{x1D70C}}\right)=\frac{R_{g}}{(1-b\unicode[STIX]{x1D70C})^{2}}+\frac{1}{2}T^{-(3/2)}T_{c}^{1/2}\frac{2a\unicode[STIX]{x1D70C}+ab\unicode[STIX]{x1D70C}^{2}}{(1+b\unicode[STIX]{x1D70C})^{2}}, & & \displaystyle\end{eqnarray}$$
(B 12) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \left.\frac{\unicode[STIX]{x2202}^{2}p}{\unicode[STIX]{x2202}T^{2}}\right|_{\unicode[STIX]{x1D70C}}=-\frac{3}{4}T^{-(5/2)}T_{c}^{1/2}\frac{a\unicode[STIX]{x1D70C}^{2}}{1+b\unicode[STIX]{x1D70C}},\\ \displaystyle \left.\frac{\unicode[STIX]{x2202}^{2}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{2}}\right|_{T}=\frac{2bR_{g}T}{(1-\unicode[STIX]{x1D70C}b)^{3}}-T^{-(1/2)}T_{c}^{1/2}\frac{2a}{(1+\unicode[STIX]{x1D70C}b)^{3}},\end{array}\right\} & & \displaystyle\end{eqnarray}$$
(B 13) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \left.\frac{\unicode[STIX]{x2202}e}{\unicode[STIX]{x2202}T}\right|_{\unicode[STIX]{x1D70C}}=C_{v}+\left.\frac{\unicode[STIX]{x2202}C_{v}}{\unicode[STIX]{x2202}T}\right|_{\unicode[STIX]{x1D70C}}T-\frac{3}{4}\frac{a}{b}T^{-(3/2)}T_{c}^{1/2}\ln \frac{1}{1+b\unicode[STIX]{x1D70C}},\\ \displaystyle \left.\frac{\unicode[STIX]{x2202}e}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right|_{T}=-\frac{3}{2}T^{-(1/2)}T_{c}^{1/2}\frac{a}{1+b\unicode[STIX]{x1D70C}}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

B.3 The Peng–Robinson EoS

The Peng–Robinson (Peng & Robinson 1976) EoS modifies the original RK and SRK (RK modified by Soave (1972)) EoS, giving better results regarding the liquid density, vapour pressure and equilibrium ratios. It is one of the most used EoS. It is given as

(B 14) $$\begin{eqnarray}\displaystyle p=\frac{R_{g}T}{\unicode[STIX]{x1D717}-b}-\frac{a\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D717}^{2}+2b\unicode[STIX]{x1D717}-b^{2}}. & & \displaystyle\end{eqnarray}$$

The constants $a$ , $b$ and parameter $\unicode[STIX]{x1D6FC}$ are given by

(B 15a-c ) $$\begin{eqnarray}\displaystyle a=\frac{0.457235R_{g}^{2}T_{c}^{2}}{p_{c}},\quad b=\frac{0.077796R_{g}T_{c}}{p_{c}},\quad \unicode[STIX]{x1D6FC}=(1+K(1-\sqrt{T/T_{c}}))^{2}.\qquad & & \displaystyle\end{eqnarray}$$

Here $K=0.37464+1.54226\unicode[STIX]{x1D714}-0.26992\unicode[STIX]{x1D714}^{2}$ , $\unicode[STIX]{x1D714}$ is the acentric factor of the species. The internal energy

(B 16) $$\begin{eqnarray}\displaystyle e=C_{v}T+\frac{a}{2\sqrt{2}b}[(1+K)^{2}-K(1+K)\sqrt{T/T_{c}}]\ln \frac{1+b(1-\sqrt{2})\unicode[STIX]{x1D70C}}{1+b(1+\sqrt{2})\unicode[STIX]{x1D70C}}. & & \displaystyle\end{eqnarray}$$

The derivatives used in the linear stability equations are given by

(B 17) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}T}\right|_{\unicode[STIX]{x1D70C}}=\frac{\unicode[STIX]{x1D70C}R_{g}}{1-\unicode[STIX]{x1D70C}b}+K\sqrt{\frac{\unicode[STIX]{x1D6FC}}{TT_{c}}}\frac{a\unicode[STIX]{x1D70C}^{2}}{1+2b\unicode[STIX]{x1D70C}-b^{2}\unicode[STIX]{x1D70C}^{2}}, & \displaystyle\end{eqnarray}$$
(B 18) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right|_{T}=\frac{R_{g}T}{(1-\unicode[STIX]{x1D70C}b)^{2}}-\unicode[STIX]{x1D6FC}\frac{2a\unicode[STIX]{x1D70C}+2ab\unicode[STIX]{x1D70C}^{2}}{(1+2b\unicode[STIX]{x1D70C}-b^{2}\unicode[STIX]{x1D70C}^{2})^{2}}, & \displaystyle\end{eqnarray}$$
(B 19) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right|_{T}\left(\left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}T}\right|_{\unicode[STIX]{x1D70C}}\right)=\frac{R_{g}}{(1-\unicode[STIX]{x1D70C}b)^{2}}+K\sqrt{\frac{\unicode[STIX]{x1D6FC}}{TT_{c}}}\frac{2a\unicode[STIX]{x1D70C}+2ab\unicode[STIX]{x1D70C}^{2}}{(1+2b\unicode[STIX]{x1D70C}-b^{2}\unicode[STIX]{x1D70C}^{2})^{2}}, & \displaystyle\end{eqnarray}$$
(B 20) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}^{2}p}{\unicode[STIX]{x2202}T^{2}}\right|_{\unicode[STIX]{x1D70C}}=-\frac{K(1+K)}{2\sqrt{T^{3}T_{c}}}\frac{a\unicode[STIX]{x1D70C}^{2}}{1+2b\unicode[STIX]{x1D70C}-b^{2}\unicode[STIX]{x1D70C}^{2}}, & \displaystyle\end{eqnarray}$$
(B 21) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}^{2}p}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{2}}\right|_{T}=\frac{2R_{g}bT}{(1-\unicode[STIX]{x1D70C}b)^{3}}-\unicode[STIX]{x1D6FC}\frac{2a(2b^{3}\unicode[STIX]{x1D70C}^{3}+3b^{2}\unicode[STIX]{x1D70C}^{2}+1)}{(1+2b\unicode[STIX]{x1D70C}-b^{2}\unicode[STIX]{x1D70C}^{2})^{3}}, & \displaystyle\end{eqnarray}$$
(B 22) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}e}{\unicode[STIX]{x2202}T}\right|_{\unicode[STIX]{x1D70C}}=C_{v}+\left.\frac{\unicode[STIX]{x2202}C_{v}}{\unicode[STIX]{x2202}T}\right|_{\unicode[STIX]{x1D70C}}T+\frac{a}{4\sqrt{2}b}[-K(1+K)\sqrt{1/TT_{c}}]\ln \frac{1+b(1-\sqrt{2})\unicode[STIX]{x1D70C}}{1+b(1+\sqrt{2})\unicode[STIX]{x1D70C}}, & \displaystyle\end{eqnarray}$$
(B 23) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\frac{\unicode[STIX]{x2202}e}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}\right|_{T}=-\frac{a}{1+2b\unicode[STIX]{x1D70C}-b^{2}\unicode[STIX]{x1D70C}^{2}}[\left(1+K\right)^{2}-K(1+K)\sqrt{T/T_{c}}]. & \displaystyle\end{eqnarray}$$

Appendix C. Influence of the bulk viscosity

The dynamics of a fluid is described by the Navier–Stokes equation, which in its simplest form contains a linear relation between deformation of a fluid element and the resulting stress, with the shear viscosity $\unicode[STIX]{x1D707}$ the coefficient of proportionality. Phenomenologically, another coefficient is possible, the second viscosity $\unicode[STIX]{x1D706}$ , which was introduced by Stokes (1845). Stokes anticipated that the second viscosity might play a role in compressible fluids. However, for the cases he considered, the fluids can be assumed incompressible with negligible dilatational effects, such that the bulk viscosity within the second viscosity can be ignored. This is known as the Stokes approximation. Consequently, setting the bulk viscosity $\unicode[STIX]{x1D707}_{b}$ to zero, has been broadly adopted in numerical simulations of compressible flows (see a succinct review by Graves & Argrow 1999).

Cramer’s (2012) numerical estimates indicate that $\unicode[STIX]{x1D707}_{b}/\unicode[STIX]{x1D707}$ of some common gases can reach $O(10^{3})$ . To investigate the influence of $\unicode[STIX]{x1D707}_{b}$ on the results of the linear stability, we performed simulations with $\unicode[STIX]{x1D707}_{b}=1000\unicode[STIX]{x1D707}$ . The results are shown in figures 19 and 20, which show the comparison of the linear stability results for $\unicode[STIX]{x1D707}_{b}=1000\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D707}_{b}=0$ , using the RP model (the other parameters are kept the same). Figure 19 shows that the neutral curves are barely affected. A discernible difference only exists in the transcritical case ( $T_{w}^{\ast }=300~\text{K}$ , $PrEc=0.06$ ), where the neutral curve with $\unicode[STIX]{x1D707}_{b}=1000\unicode[STIX]{x1D707}$ becomes slightly more expanded. On the other hand, the algebraic instability does not vary with bulk viscosity. Only the modal growth region ( $G_{max}=\infty$ ) in figure 20(b) becomes larger with $\unicode[STIX]{x1D707}_{b}=1000\unicode[STIX]{x1D707}$ and is consistent with the results shown in figure 19(b).

Figure 19. Influence of bulk viscosity on neutral curves in (a) subcritical (b) transcritical and (c) supercritical case.

Figure 20. Influence of bulk viscosity on the algebraic growth. $PrEc=0.07$ , $Re=1000$ . (a) Subcritical (b) transcritical and (c) supercritical case.

The above comparisons support the Stokes hypothesis used in this study. In fact, $\unicode[STIX]{x1D707}_{b}$ is frequency dependent, this means that, depending on the perturbation one prescribes in the stability analysis, the bulk viscosity will have different values. Therefore for a more rigorous investigation we would need reliable frequency resolved data for the bulk viscosity, either from theories, experiments (Karim & Rosenhead 1952) or molecular dynamics simulations (Hoover et al. 1980).

Appendix D. Influence of the reference scaling

Previous studies have shown that the scaling of the governing equations may have a large influence on the results. For example, if the viscosity at the cold wall is used as a reference value, Sahu & Matar (2010) concluded that increasing the temperature difference between both walls destabilizes the flow, while Sameen & Govindarajan (2007) concluded the opposite behaviour if the viscosity at the hot wall is used. On the other hand, Rinaldi et al. (2018) investigated the edge state solutions of viscosity stratified flows where they showed that a different reference value for viscosity does not qualitatively change the results. In this appendix, we show how the definition of the non-dimensional quantities will influence the results presented in the paper.

Figure 21. Influence of the reference scaling on the neutral curve for the (a) subcritical, (b) transcritical and (c) supercritical cases. Solid lines show the results with wall scaling (in the $\unicode[STIX]{x1D6FC}{-}Re$ diagram), while dashed lines indicate the average scaling (in the $\unicode[STIX]{x1D6FC}{-}\overline{Re}$ diagram).

We introduce the averaged values of the thermodynamic and transport properties:

(D 1a,b ) $$\begin{eqnarray}\displaystyle T_{av}^{\ast }=\frac{1}{h^{\ast }}\int T^{\ast }\,\text{d}y,\quad \unicode[STIX]{x1D70C}_{av}^{\ast }=\frac{1}{h^{\ast }}\int \unicode[STIX]{x1D70C}^{\ast }\,\text{d}y, & & \displaystyle\end{eqnarray}$$
(D 2a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}_{av}^{\ast }=\frac{1}{h^{\ast }}\int \unicode[STIX]{x1D707}^{\ast }\,\text{d}y,\quad \unicode[STIX]{x1D705}_{av}^{\ast }=\frac{1}{h^{\ast }}\int \unicode[STIX]{x1D705}^{\ast }\,\text{d}y. & & \displaystyle\end{eqnarray}$$

When the governing equations are scaled by the above averaged values, one obtains the averaged Reynolds number, $\overline{Re}$ , and the product of the averaged Prandtl and Eckert numbers, $\overline{PrEc}$ :

(D 3a,b ) $$\begin{eqnarray}\displaystyle \overline{Re}=\frac{\unicode[STIX]{x1D70C}_{av}^{\ast }u_{r}^{\ast }h^{\ast }}{\unicode[STIX]{x1D707}_{av}^{\ast }},\quad \overline{PrEc}=\frac{\unicode[STIX]{x1D707}_{av}^{\ast }u_{r}^{\ast 2}}{\unicode[STIX]{x1D705}_{av}^{\ast }T_{av}^{\ast }}. & & \displaystyle\end{eqnarray}$$

We name it the average scaling, to distinguish from the wall scaling presented in § 2.1 of the paper. Note that the reference velocity is not independent, and is given by

(D 4) $$\begin{eqnarray}\displaystyle u_{r}^{\ast }=\sqrt{PrEc\frac{\unicode[STIX]{x1D705}_{w}^{\ast }T_{w}^{\ast }}{\unicode[STIX]{x1D707}_{w}^{\ast }}}=\sqrt{\overline{PrEc}\frac{\unicode[STIX]{x1D705}_{av}^{\ast }T_{av}^{\ast }}{\unicode[STIX]{x1D707}_{av}^{\ast }}}. & & \displaystyle\end{eqnarray}$$

Using either scalings resulted in a qualitatively similar conclusion as shown in figure 21 for the modal instability. That is, the flow becomes more unstable in the subcritical regime, inviscid unstable in the transcritical regime and more stable in the supercritical regime.

Figure 22. Influence of the reference scaling on the maximum algebraic growth $G_{max}$ for the (a) subcritical, (b) transcritical and (c) supercritical cases. The left and right halves show the results with wall scaling and average scaling respectively. The non-ideal gas model is RP, $PrEc=0.07$ , $Re=1000$ .

Table 6. Maximum algebraic growth $G_{max}$ over $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ . $PrEc=0.07$ , $Re=1000$ (wall scaling) or $\overline{Re}=1000$ (average scaling). The values for the $T_{w}^{\ast }=300$ case are for $\unicode[STIX]{x1D6FC}=0$ where modal instability is absent.

Regarding the algebraic instability using the average scaling, as can be seen from figure 22, the maximum growth shows a minor reduction in the subcritical regime. Increases in $G_{max}$ are noticed for the trans- and supercritical regimes. Comparisons with the ideal gas have been summarized in table 6. The ideal gas is not sensitive to the wall temperature under either scaling. With average scaling, the conclusion for the algebraic instability will not change.


Alferez, N. & Touber, E. 2017 One-dimensional refraction properties of compression shocks in non-ideal gases. J. Fluid Mech. 814, 185221.
Brunner, G. 2010 Applications of supercritical fluids. Annu. Rev. Chem. Biom. Engng 1, 321342.
Busse, F. H. 1969 Bounds on the transport of mass and momentum by turbulent flow between parallel plates. Z. Angew. Math. Phys. 20 (1), 114.
Chikkadi, V., Sameen, A. & Govindarajan, R. 2005 Preventing transition to turbulence: a viscosity stratification does not always help. Phys. Rev. Lett. 95 (26), 264504.
Chu, B.-T. 1965 On the energy transfer to small disturbances in fluid flow (part i). Acta Mechanica 1 (3), 215234.
Cramer, M. S. 2012 Numerical estimates for the bulk viscosity of ideal gases. Phys. Fluids 24 (6), 066102.
Duan, L., Beekman, I. & Martin, M. P. 2010 Direct numerical simulation of hypersonic turbulent boundary layers. Part 2. Effect of wall temperature. J. Fluid Mech. 655, 419445.
Fedorov, A. 2011 Transition and stability of high-speed boundary layers. Annu. Rev. Fluid Mech. 43 (1), 7995.
George, K. J. & Sujith, R. I. 2011 On Chu’s disturbance energy. J. Sound Vib. 330 (22), 52805291.
Govindarajan, R. & Sahu, K. C. 2014 Instabilities in viscosity-stratified flow. Annu. Rev. Fluid Mech. 46, 331353.
Graves, R. E. & Argrow, B. M. 1999 Bulk viscosity: past to present. J. Thermophys. Heat Transfer 13 (3), 337342.
Hanifi, A., Schmid, P. J. & Henningson, D. S. 1996 Transient growth in compressible boundary layer flow. Phys. Fluids 8 (3), 826837.
Hoover, W. G., Ladd, A. J. C., Hickman, R. B. & Holian, B. L. 1980 Bulk viscosity via nonequilibrium and equilibrium molecular dynamics. Phys. Rev. A 21 (5), 17561760.
Joseph, D. D. & Carmi, S. 1969 Stability of Poiseuille flow in pipes, annuli, and channels. Q. Appl. Maths 26 (4), 575599.
Karim, S. M. & Rosenhead, L. 1952 The second coefficient of viscosity of liquids and gases. Rev. Mod. Phys. 24, 108116.
Kawai, S. 2016 Direct numerical simulation of transcritical turbulent boundary layers at supercritical pressures with strong real fluid effects. In 54th AIAA Aerospace Sciences Meeting, p. 1934.
Kawai, S., Terashima, H. & Negishi, H. 2015 A robust and accurate numerical method for transcritical turbulent flows at supercritical pressure with an arbitrary equation of state. J. Comput. Phys. 300, 116135.
Lemmon, E. W., Huber, M. L. & McLinden, M. O.2013 NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties – REFPROP, Version 9.1. National Institute of Standards and Technology. Available at:
Mack, L. M.1969 Boundary layer stability theory. JPL Rep. 900-277. Jet Propulsion Laboratory.
Mack, L. M. 1976 A numerical study of the temporal eigenvalue spectrum of the Blasius boundary layer. J. Fluid Mech. 73 (03), 497520.
Malik, M., Dey, J. & Alam, M. 2008 Linear stability, transient energy growth, and the role of viscosity stratification in compressible plane Couette flow. Phys. Rev. E 77 (3), 036322.
Marxen, O., Iaccarino, G. & Magin, T. E. 2014 Direct numerical simulations of hypersonic boundary-layer transition with finite-rate chemistry. J. Fluid Mech. 755, 3549.
Marxen, O., Magin, T. E., Shaqfeh, E. S. & Iaccarino, G. 2013 A method for the direct numerical simulation of hypersonic boundary-layer instability with finite-rate chemistry. J. Comput. Phys. 255, 572589.
Moran, M. J., Shapiro, H. N., Boettner, D. D. & Bailey, M. B. 2012 Principles of Engineering Thermodynamics, 7th edn. Wiley.
Morinishi, Y., Tamano, S. & Nakabayashi, K. 2004 Direct numerical simulation of compressible turbulent channel flow between adiabatic and isothermal walls. J. Fluid Mech. 502, 273308.
Nemati, H., Patel, A., Boersma, B. J. & Pecnik, R. 2016 The effect of thermal boundary conditions on forced convection heat transfer to fluids at supercritical pressure. J. Fluid Mech. 800, 531556.
Orr, W. M’F. 1907 The stability or instability of the steady motions of a perfect liquid and of a viscous liquid. Part II: A viscous liquid. Proc. R. Irish Acad. 27, 69138.
Orszag, S. A. 1971 Accurate solution of the Orr–Sommerfeld stability equation. J. Fluid Mech. 50 (4), 689703.
Patel, A., Boersma, B. J. & Pecnik, R. 2016 The influence of near-wall density and viscosity gradients on turbulence in channel flows. J. Fluid Mech. 809, 793820.
Peeters, J. W. R., Pecnik, R., Rohde, M., van der Hagen, T. H. J. J. & Boersma, B. J. 2016 Turbulence attenuation in simultaneously heated and cooled annular flows at supercritical pressure. J. Fluid Mech. 799, 505540.
Peng, D.-Y. & Robinson, D. B. 1976 A new two-constant equation of state. Ind. Engng Chem. Fundam. 15 (1), 5964.
Pinarbasi, A. & Liakopoulos, A. 1995 The role of variable viscosity in the stability of channel flow. Intl Commun. Heat Mass Transfer 22 (6), 837847.
Pirozzoli, S., Bernardini, M. & Grasso, F. 2008 Characterization of coherent vortical structures in a supersonic turbulent boundary layer. J. Fluid Mech. 613, 205231.
Potter, M. C. & Graber, E. 1972 Stability of plane Poiseuille flow with heat transfer. Phys. Fluids 15 (3), 387391.
Redlich, O. & Kwong, J. N. S. 1949 On the thermodynamics of solutions. V. An equation of state. Fugacities of gaseous solutions. Chem. Rev. 44 (1), 233244.
Rinaldi, E., Patel, A., Schlatter, P. & Pecnik, R. 2017 Linear stability of buffer layer streaks in turbulent channels with variable density and viscosity. Phys. Rev. Fluids 2 (11), 113903.
Rinaldi, E., Schlatter, P. & Bagheri, S. 2018 Edge state modulation by mean viscosity gradients. J. Fluid Mech. 838, 379403.
Sahu, K. C. 2011 The instability of flow through a slowly diverging pipe with viscous heating. Trans. ASME J. Fluids Engng 133 (7), 071201.
Sahu, K. C. & Govindarajan, R. 2014 Instability of a free-shear layer in the vicinity of a viscosity-stratified layer. J. Fluid Mech. 752, 626648.
Sahu, K. C. & Matar, O. K. 2010 Stability of plane channel flow with viscous heating. Trans. ASME J. Fluids Engng 132 (1), 011202.
Saikia, B., Ramachandran, A., Sinha, K. & Govindarajan, R. 2017 Effects of viscosity and conductivity stratification on the linear stability and transient growth within compressible Couette flow. Phys. Fluids 29 (2), 024105.
Sameen, A., Bale, R. & Govindarajan, R. 2011 The effect of wall heating on instability of channel flow – CORRIGENDUM. J. Fluid Mech. 673, 603605.
Sameen, A. & Govindarajan, R. 2007 The effect of wall heating on instability of channel flow. J. Fluid Mech. 577, 417442.
Schmid, P. J. & Brandt, L. 2014 Analysis of fluid systems: stability, receptivity, sensitivity. Lecture notes from the FLOW-NORDITA summer school on advanced instability methods for complex flows, Stockholm, Sweden, 2013. Appl. Mech. Rev. 66 (2), 024803.
Schmid, P. J. & Henningson, D. S. 2001 Stability and Transition in Shear Flows. Springer.
Sciacovelli, L., Cinnella, P., Content, C. & Grasso, F. 2016 Dense gas effects in inviscid homogeneous isotropic turbulence. J. Fluid Mech. 800, 140179.
Sciacovelli, L., Cinnella, P. & Gloerfelt, X. 2017a Direct numerical simulations of supersonic turbulent channel flows of dense gases. J. Fluid Mech. 821, 153199.
Sciacovelli, L., Cinnella, P. & Grasso, F. 2017b Small-scale dynamics of dense gas compressible homogeneous isotropic turbulence. J. Fluid Mech. 825, 515549.
Soave, G. 1972 Equilibrium constants from a modified Redlich–Kwong equation of state. Chem. Engng Sci. 27 (6), 11971203.
Sommerfeld, A.1908 Ein Beitrag zur hydrodynamischen Erklaerung der turbulenten Fluessigkeitsbewegungen. Atti del IV Congresso Internationale dei Matematici, Rome, pp. 116–124. Accademia dei Lincei.
Span, R. & Wagner, W. 2003 Equations of state for technical applications. I. Simultaneously optimized functional forms for nonpolar and polar fluids. Intl J. Thermophys. 24, 139.
Stokes, G. G. 1845 On the theories of the internal friction of fluids in motion, and of the equilibrium and motion of elastic solids. Trans. Camb. Phil. Soc. 8 (22), 287319.
Thomas, L. H. 1953 The stability of plane Poiseuille flow. Phys. Rev. 91, 780783.
van der Waals, J. D.1873 Over de continuiteit van den gas- en vloeistoftoestand. PhD thesis, Leiden, Netherlands.
Wall, D. P. & Wilson, S. K. 1996 The linear stability of channel flow of fluid with temperature-dependent viscosity. J. Fluid Mech. 323, 107132.
Wall, D. P. & Wilson, S. K. 1997 The linear stability of flat-plate boundary-layer flow of fluid with temperature-dependent viscosity. Phys. Fluids 9 (10), 28852898.
Wazzan, A. R., Keltner, G., Okamura, T. T. & Smith, A. M. O. 1972 Spatial stability of stagnation water boundary layer with heat transfer. Phys. Fluids 15 (12), 21142118.
Zappoli, B., Beysens, D. & Garrabos, Y. 2015 Heat Transfers and Related Effects in Supercritical Fluids. Springer.
Zhong, X. & Wang, X. 2012 Direct numerical simulation on the receptivity, instability, and transition of hypersonic boundary layers. Annu. Rev. Fluid Mech. 44 (1), 527561.