Linear instability of Poiseuille flows with highly non-ideal fluids

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 (CO$_2$) at pressure of 80 bar, which is above its critical pressure of 73.9 bar. The Poiseuille flow is characterized by the Reynolds number ($Re=\rho_{w}^{*}u_{r}^{*}h^{*}/\mu_{w}^{*}$), the product of Prandtl ($Pr=\mu_{w}^{*}C_{pw}^{*}/\kappa_{w}^{*}$) and Eckert number ($Ec=u_{r}^{*2}/C_{pw}^{*}T_{w}^{*}$), and the wall temperature that in addition to pressure determines 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, subcritical (temperature in the channel is entirely below its pseudo-critical value), transcritical, and supercritical temperature regime. If compared to the linear stability of an ideal gas Poiseuille flow, we show that the base flow is more unstable in the subcritical regime, inviscid unstable in the transcritical regime, while significantly more stable in the supercritical regime. Following the corresponding states principle, we expect that qualitatively similar results will be obtained for other fluids at equivalent thermodynamic states.


Introduction
Many processes in industrial applications constitute of flows with fluids that do not follow the ideal gas law. For example, flows in vapour power systems, re-entry of spacecrafts, 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 et al. (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. (2017aSciacovelli et al. ( ,b, 2016 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 et al. (2016) studied the influence of variable properties on fully developed turbulent channel flows and derived a velocity transformation that allows 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 et al. (2004), Pirozzoli et al. (2008), Duan et al. (2010). Nemati et al. (2016); Peeters et al. (2016) studied turbulent heat transfer to supercritical 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 * /ρ * R * T * . Figure 1 shows the T − ϑ diagram (temperature -specific volume diagram, ϑ = 1/ρ) of carbon dioxide CO 2 . The white circle in each subplot indicates the critical point, which for CO 2 is at a pressure and temperature of p * c = 73.9 bar and T * c = 304.25 K. In this paper, we denote dimensional and critical quantities with superscript ' * ' 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 (colored 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 behavior is roughly 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 , dynamic viscosity µ * and thermal conductivity κ * ) 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.
In view of its great simplicity, most of the present knowledge on stability and laminarturbulent 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 has just gone through an initial stage (Zhong & Wang 2012;Marxen et al. 2013Marxen et al. , 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 Figure 1. T − ϑ diagram of CO2 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 colored contour lines of (a) Z = p * /ρ * R * T * , (b) the specific heat capacity at constant volume C * v , (c) the dynamic viscosity µ * , and (d) the thermal conductivity κ * . 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 centerline flow velocity. Due to the non-normality of most practical linear systems, the modal stability analysis cannot cover the full behavior 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. Also for plane Poiseuille flow, this number is 49.6 (Joseph & Carmi 1969;Busse 1969).
Studies of viscosity stratified flows, where the viscosity depends on temperature, has recently received attention, readers may refer to  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 instabilites in other types of flows (Chikkadi et al. 2005;Sameen & Govindarajan 2007;Sameen et al. 2011;Sahu 2011;. For compressible plane Couette flow, Malik et al. (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 is recently studied by Rinaldi et al. (2018), showing that in minimal channel flows, the kinetic energy level and the driving force of self-sustained cycle of the edge state depends on viscosity distribution. The above studies are based on the incompressible flow assumption or the ideal gas equation-of-state (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.

Flow conservation equations
The laws of conservation of mass, momentum and energy (the Navier-Stokes (N-S) equations), in dimensionless form, are given by 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, ρ 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, τ ij , and the heat flux, q j , are given by Here µ is the dynamic viscosity, λ = µ b − 2/3µ the second viscosity, µ b the bulk viscosity, and κ is the thermal conductivity. Results presented in the following sections are subject to µ 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.6) The subscript w denotes wall values, h * is the half channel height, c * w is the speed of sound at the wall, u * r is the reference velocity. Note that for an ideal gas Ec = (γ −1)Ma 2 , where γ 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.

Fluid equations of state
In order to find a closed form of the conservation equations, an equation of state for the fluid has to be specified. As a representative example of non-ideal fluids, the study is performed with CO 2 at a pressure of p * = 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 et al. 2002) 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. The EoS are suitable for a broad variety of fluids while high accuracy can be maintained. Readers shall 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 ρ 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 model (IG). A constant specific heat ratio γ = 1.289 is used for the IG model. All the fluid models are summarized in table 1. Figure 2 shows the thermodynamic and transport properties of CO 2 at a pressure of 80 bar. The pentagram in subplot (a) shows the pseudo-critical temperature (T * pc = 307.7 K, RP model), which is defined as the point on a supercritical isobar where C * p reaches a maximum. Near T * pc , 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 CO 2 as implemented in REFPROP (RP). In general, the cubic EoSs 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.9) 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.

The linearized stability equations
Following the common procedure, the flow field is decomposed into the base flow and a perturbation, as It is known that for simple compressible systems (e.g. pure substances, uniform mixture of nonreacting gases), the thermodynamic state is defined by two independent thermodynamic properties. In this study, we keep ρ and T as the two basic thermodynamic variables, while the other thermodynamic and transport properties (e.g. E, p, µ, κ) are determined as functions of ρ and T . For example, the pressure perturbation p is expanded by a Taylor-series with respect to ρ 0 and T 0 in the following way For the sake of brevity, the partial derivative of a quantity with respect to T at constant ρ, will be written as ∂/∂T | ρ0 ≡ ∂/∂T , and accordingly ∂/∂ρ| T0 ≡ ∂/∂ρ. The stability equation is derived by substituting (2.10) into the N-S equations (2.1), (2.2) and (2.3), and then subtracting the governing equations of the base flow. With the nonlinear terms neglected, the linear stability equations are formulated as and V xz are matrices of size 5 × 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 cubic EoS (see Appendix B) or numerically employing finite-difference method within the REFPROP library.

Modal and algebraic stability
In modal stability, the perturbation is assumed to have the form ( 2.13) where c.c. stands for the complex conjugate. Substituting (2.13) into (2.12) results in where D = d/dy. The equation (2.14) is solved as an eigenvalue problem, which describes the development of the perturbations in temporal or spatial domain, i.e.
Here we consider the temporal problem only, therefore α and β are the prescribed streamwise and spanwise wave numbers. ω = ω r + iω i is solved as the eigenvalue, where ω r and ω i give the angular frequency and growth rate of the perturbation. The domain 0 y 2 is discretized with Chebyshev collocation points, defined by The differentiation of (2.15) is accomplished with the matrix form of Chebyshev collocation derivatives. Numerical tests indicate that, typically N = 200 (used here) are 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: Here E (q) is the disturbance energy with the definition as given in (5.1). The perturbation q is expanded by the eigenvector obtained from the modal stability. The calculation of the optimal energy amplification G and the corresponding optimal perturbation (the input), as well as the resulting perturbation (the output), lead to a singular value problem, which is solved with the same Chebyshev differentiation method as in the modal growth (Schmid & Henningson 2001;Schmid & Brandt 2014). The (modal and algebraic) perturbations are solved subjected to the boundary condition: u = v = w = T = 0 at the lower (y = 0) and upper wall (y = 2).

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.2) and (2.3) with the assumption that the flow is fully developed, spanwise and streamwise independent, steady and parallel, i.
pseudo-critical temperature inflection-point It is worth noting that the above equations are independent of density, therefore, ρ 0 can be separately determined by the EoS. We assume the body force,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 = constant, µ and κ 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). µ and κ are then updated using the obtained temperature. This procedure is repeated until the solution is converged.

The isothermal limit
When PrEc → 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, ρ 0 = 1, µ 0 = 1, κ 0 = 1. The flow is thus simply governed by ∂ 2 u/∂y 2 = −F . Choosing u * r as the centerline velocity, leads to settingF = 2. As a result, the dimensionless base flow is independent of any parameters (e.g. T w ,F and P rEc) 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).

The compressible base flow
Equations (3.1), (3.2) and (3.3) show that the compressible base flow is determined by PrEc,F , and T * w . Either by increasing PrEc orF , the compressibility effects become more significant. Without loss of the generality, a constant body forceF = 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 = 290, 300, or 310 K, the Mach number is Ma = 0.40, 0.58, or 1.35, respectively. In this work, the wall temperature T * w 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 centerline temperature T * center and velocity u center = u * center /u * r as a function of wall temperature T * w and PrEc. Regardless of the wall temperature, an increase of PrEc is accompanied with an increase of T * center and u center as compressible effects become more prominent. Interestingly, a distinct right-angled triangular area emerges in each subplot of figure 4. At the hypotenuse of this triangle, the centerline temperature, T * center , and velocity, u center , suddenly increase, forming a discontinuity in the P rEc-T * w plane. It is also interesting to note that the hypotenuse of the triangle almost coincides with  the line where T * center reaches the pseudo critical temperature T * pc = 307.7 K (shown with the dot-dashed line). Likewise, the upper boundary of the triangle coincides with the dotted line where T * w = T * pc . 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 P rEc 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 subplot. 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 P rEc can be explained as follows. Integrating (3.1) gives µ∂u/∂y = −F y + C. Applying the symmetry condition at the 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 center reaches T * pc . Recall figure 2(c), a sharp gradient of the viscosity (∂µ/∂T 0) is seen close to the pseudo-critical point. As PrEc increases, ∂T /∂y increases at the wall, such that ∂µ/∂y ∼ = (∂µ/∂T )(∂T /∂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, µ|∂u/∂y| at the wall is equal to the constant forcingF , regardless of PrEc and wall temperature, the velocity profiles with/without inflectional points are isolated by the line of constant gradientF (the dash-dotted lines that form a triangle in figure 5(g-i)). Therefore, a velocity profile without an inflection point cannot reach the apex of the triangle (|∂u/∂y| decreases towards channel center) and the sudden increase of the centerline velocity appears once an inflection point is formed.
In general, the base flow solutions can be summarized as follows: • In the subcritical case, the wall temperature is much lower than T * pc , and in the range of PrEc considered, T * center is always less than T * pc . Hence, the velocity profile is not inflectional.
• In the transcritical case, the wall temperature is close to T * pc , such that for large enough PrEc, T * center reaches T * pc . 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 , the larger the discontinuity will be. • In the supercritical case, the wall temperature is higher than T * pc . The properties of the fluid are gas-like (compressed vapour) and the velocity is not inflectional.

Linear modal instability
Depending on the cases discussed below, we will use the definition of dynamic and thermodynamic modes, as ρ = 0, and T = 0 (dynamic modes) ρ = 0, or T = 0 (thermodynamic modes). (4.1)

The isothermal limit
With the base flow obtained in section 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 = 10000, α = 1 and β = 0, the A-, P-and S-branches (originally named by Mack 1976) are reproduced by incompressible equations. Comparing the results using different equations of state, 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 ω 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). Figure 6. Eigen spectrum (a) and neutral curve (b) for the isothermal limit. The eigen spectrum is subject to α = 1, β = 0, and Re = 10000. The neutral curve is solved for 2-D perturbations (β = 0). Symbols show results using different fluid models (RP, PR, RK, VW, IG) and incompressible equations (IC). IC in subplot (b) shows the results given by Schmid & Henningson (2001, pp. 71). In subplot (c) and (d), profiles of the unstable mode (ω = 0.2375 + 0.0037i) and one of the stable modes (ω = 0.4164−0.1382i, highlighted in orange in the spectrum) are shown. The perturbations are normalized by |u |. An offset of -0.1 and -0.2 is applied to |ρ | and |T |. The solid lines are results with fluid model IG.

Compressible flows
To achieve a first impression of the non-ideal gas effects, the problem is first studied with the RP model, where thermodynamic and transport properties are taken from the REFPROP library. Figure 7 shows the neutral curves (a-c) as well as eigenfunctions (d-f) at representative parameters. As discussed in § 3.2, the temperature is subcritical, transcritical and supercritical with T * w = 290K, 300K and 310K respectively. The results are compared with ideal gas (IG).
By increasing PrEc, the base flow of the ideal gas becomes more stable as the critical Reynolds number increases, regardless of T w specified. In fact, despite the difference in wall temperature, the dimensionless thermodynamic and transport properties (scaled with wall values) remain much the same. On the other hand, the behaviour of the non-ideal cases is different for the three cases investigated. In the subcritical case, the flow becomes more unstable when PrEc is increased. This is manifested by the enlargement of the neutral curve. Similarly, the transcritical case becomes more unstable as PrEc increases. However, once PrEc reaches the critical value (in this case PrEc = 0.05115), the base flow becomes inflectional. The flow is thus inviscid unstable and the critical Reynolds number is substantially reduced. For instance, the flow is unstable for Re < 1000 and PrEc = 0.06. In the supercritical case, the increase of PrEc stabilizes the base flow and the non-ideal gas is even more stable than the ideal gas. In this case, when PrEc reaches 0.03, the modal instability is found after Re > 8000. Interestingly, a weak influence of PrEc on the velocity perturbations is observed (see figure 7 (d-f)), while the amplitudes of density and temperature perturbations are considerably larger if PrEc increases. For the transcritical case, when the flow enters the triangular zone, the density perturbation are the most dominant.
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 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 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 state the fidelity of the cubic EoS models in terms of the growth rate.

The kinetic energy budget
To further understand the instability mechanism of the non-ideal fluids, we perform a kinetic energy budget analysis for the 2D perturbation. The energy balance equation is the sum of the x-momentum perturbation equation, multiplied withû † , and the y-equation, multiplied withv † . Here, dagger stands for the complex conjugate. The continuity equation is used to substitute the temporal growth of density, which appears in the x-momentum equation. This gives the following kinetic energy balance equation: (4.7) The real part of the equation (4.2) describes the balance of the kinetic energy growth. In particular, K r is the temporal growth of the kinetic energy, Θ 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 α = 1 and Re = 10000. 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 = 300 K, PrEc = 0.06), which has a considerably larger growth rate (as also shown in figure 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 = 300 K. It can be inferred that the inflectional velocity profile (PrEc = 0.06) has caused a larger ρ 0 ∂u 0 /∂y near both walls, the amplitude of the velocity perturbation vû † is larger as well. Therefore, a large production term − ρ 0 ∂u0 ∂yvû † dy and accordingly the large growth rate can be explained.

Choice of the energy norm
The Mack's energy norm (Mack 1969;Hanifi et al. 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 components. In fact, Mack's norm is equivalent to Chu's norm (Chu 1965;George & Sujith 2011). In the current non-ideal gas flows, the equation of states can be different (PR, RK, VW, IG), or even implicit (look-up table) as in the case of the RP EoS model. Therefore, we choose a general form of the norm: where † 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 ρ , for PrEc = 0.05 (thermodynamic components become important) and a wall temperature of T * w = 290 K. When m ρ 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 ρ = 1. Therefore, the results presented in this section are mainly obtained for m ρ = m T = 1. A comparison with Mack's energy norm (m ρ = T 0 /(ρ 2 0 γMa 2 ), m T = 1/(γ(γ − 1)T 0 Ma 2 )) is proivded at the end of this section.

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 spectrum 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 α − β 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 α = 0 and 2 β 2.1, which is well-known for ideal gas. The optimal perturbation and the corresponding output are shown in figure  11(b,c) for α = 0, β = 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.

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 IG model in figure 12, 13 and 14, respectively. The three cases actually start from the same results at the isothermal limit (figure 11a). Regardless of the wall temperature and PrEc, the largest transient growth occurs at α = 0 and 2 β 2.1 for both ideal and non-ideal gases. In the subcritical and transcritical cases (figure 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 → ∞ 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 β = 0) in the α − β diagram is observed (where G max → ∞). For better display of the results, we have limited the color 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 case, and more stable for the supercritical case. On the whole, the nonideal gas effects increase the algebraic instability in all regimes, most prominently in the transcritical regime.
The typical optimal perturbation and the resulting output are shown in figure 15 at PrEc = 0.07, α = 0 and β = 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 (ρ and T ) 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 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, α = 0 and β = 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. The main results presented in this section are based on the energy norm: m ρ = 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, α = 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.

Conclusion
Linear stability of highly non-ideal plane Poiseuille flows is studied. We have chosen carbon dioxide (CO 2 ) at supercritical pressures (p * =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 P r and Eckart 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 → 0) to typical compressible flows with PrEc = 0.1. Three cases with wall temperatures in the vicinity of the pseudo-critical point (T * pc = 307.7 K) have been investigated in more detail. In particular, the wall temperatures are set such that the temperature profile is subcritical (T * w = 290 K), transcritical (T * w = 300 K) and supercritical (T * w = 310 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 laminarturbulent transition is more likely to be dominated by modal growth 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. Additionally, we show that the linear stability analysis with simple cubic equations of state give qualitatively similar results than using the more accurate multi-parameter equation of state 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.

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 ρ 0 at constant T 0 (and viceversa) has been denoted as ∂ ∂ρ0 , instead of ∂ ∂ρ0 T0 . The elements are, The second-order finite differences were used to determine the gradients of the properties. For instance, the gradients of viscosity An example of the sensitivity of ∂ 2 µ * ∂T * ∂ρ * to ∆T * and ∆ρ * is shown in figure 18. In fact, the gradients of the thermodynamic & transport properties became rather robust when ∆T * 1 K and ∆ρ * 1 Kg/m 3 . In this study, the results are obtained with ∆T * = 0.1 K and ∆ρ * = 0.1 Kg/m 3 . Figure 18. Sensitivity of ∂ 2 µ * ∂T * ∂ρ * to ∆T * and ∆ρ * . gas constant heat capacity ratio acentric factor critical pressure critical temperature R * g = 188.9 J/(Kg K) γ = 1.289 ω = 0.228 p * c = 73.9 bar T * c = 304.1 K Table 5. The material dependent parameters for CO2.

Appendix B. Cubic equation of state
The material dependent parameters for CO2 are provided in table 5. These parameters are necessary inputs for the cubic equation of states detailed below and can be easily replaced for other fluids.

B.1. The van der Waals equation of state
The van der Waals (1873) equation of state (EoS) is the simplest cubic equation of state that is capable of accounting phase separation and the critical point (see the introduction in Zappoli et al. 2015;Moran et al. 2012). The EoS can be written as 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 ∂p ∂ϑ Tc = ∂ 2 p ∂ϑ 2 Tc = 0 ⇒ a = 27 64 Using the Maxwell relations and the departure function, it is possible to obtain the internal energy as Here K = 0.37464 + 1.54226ω − 0.26992ω 2 , ω is the acentric factor of the species. The internal energy The derivatives used in the linear stability equations are give by

Appendix C. Influence of the bulk viscosity
The dynamics of a fluid are described by the Navier-Stokes equation, which in its simplest form contain a linear relation between deformation of a fluid element and the resulting stress, with the shear viscosity µ the coefficient of proportionality. Phenomenologically, another coefficient is possible, the second viscosity λ, 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 µ b to zero, has been broadly adopted in numerical simulations of compressible flows (see a succinct review by Graves & Argrow 1999).
Cramer (2012)'s numerical estimates indicate that µ b /µ of some common gases can reach O(10 3 ). To investigate the influence of µ b on the results of the linear stability, we performed simulations with µ b = 1000µ. The results are shown in figure 19 and 20, which show the comparison of the linear stability results for µ b = 1000µ and µ 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 = 300 K, PrEc=0.06), where the neutral curve with µ b = 1000µ becomes slightly more expanded. On the other hand, the algebraic instability does not vary with bulk viscosity. Only the modal growth region (G max = ∞) in figure 20(b) becomes larger with µ b = 1000µ and is consistent with the results shown in figure 19(b).
The above comparisons support the Stokes' hypothesis used in this study. In fact, µ 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) .  Table 6. Maximum algebraic growth Gmax over α and β. PrEc = 0.07, Re = 1000 (wall scaling) or Re = 1000 (average scaling). The values for the T * w = 300 case are for α = 0 where modal instability is absent.