1. Introduction
High-pressure transcritical fluids operate within thermodynamic spaces in which supercritical gas- and liquid-like states can be differentiated across the pseudo-boiling region (Simeoni et al. Reference Simeoni, Bryk, Gorelli, Krisch, Ruocco, Santoro and Scopigno2010; Banuti Reference Banuti2015; Raju et al. Reference Raju, Banuti, Ma and Ihme2017; Jofre & Urzay 2020, 2021; Li Reference Li and Jin2024), across which a nonlinear phase transition is distinctively relevant (especially) within the range
$1 \lt P/P_c \lt 3$
(Banuti Reference Banuti2015). In particular, unlike subcritical two-phase fluids, single-component transcritical flows transition between phases in a continuous form as a consequence of the small Knudsen numbers encountered (
$K\!n \ll 1$
), resulting from reduced mean free paths at high pressures (Dahms & Oefelein Reference Dahms and Oefelein2013). The use of high-pressure supercritical fluids is a mature field in thermo-fluid engineering as they are utilised in a wide range of applications, such as liquid rocket engines, gas turbines and supercritical water-cooled reactors (Yoo Reference Yoo2013). Moreover, utilising direct numerical simulation (DNS) approaches, their inherent capacity to achieve microconfined turbulence (Bernades & Jofre Reference Bernades and Jofre2022; Bernades et al. Reference Bernades, Capuano and Jofre2023a
; Bernades et al. Reference Bernades, Duchanine, Capuano and Jofre2025) has been recently explored, which enables to significantly increase mixing and heat transfer rates in microfluidic applications. The resulting flow physics differs significantly from the typical behaviour of turbulent wall-bounded flows due to the presence of localised baroclinic torques responsible for remarkably increasing flow rotation. As a result, the flow becomes unstable and rotation is transformed into a wide range of scales (i.e. turbulent flow motions) through vortex stretching mechanisms. To that end, in this work, temperature differences are imposed between top and bottom walls (i.e. below and above the pseudo-boiling temperature) to force the fluid to undergo a transcritical trajectory by operating within a thermodynamic region across the pseudo-boiling line. However, the stability of such type of flow configurations is still not fully characterised. Thus, this work aims to conduct linear stability and transient growth analyses of relatively low-Mach/Reynolds-number wall-bounded flows at high-pressure transcritical fluid conditions to carefully identify and quantify the stability characteristics. For consistency, the same operating fluid as the DNS-based physics-characterisation analysis by Bernades et al. (Reference Bernades, Capuano and Jofre2023a
), nitrogen, is selected. In this regard, nitrogen is a substance typically used in studies of supercritical fluids. For example, among other DNS studies of high-pressure wall-bounded flows (Kawai, Terashima & Negishi Reference Kawai, Terashima and Negishi2015; Ma, Yang & Ihme Reference Ma, Yang and Ihme2018; Kawai Reference Kawai2019) and boundary layers (Boldini et al. Reference Boldini, Bugeat, Peeters, Kloker and Pecnik2024), Toki & Teramoto (Reference Toki and Teramoto2020) analysed the importance of thermodynamic variations on velocity and temperature distributions using nitrogen to accurately predict supercritical heat transfer in regenerative cooling channels.
Historically, the study of hydrodynamic stability in wall-bounded configurations was first established for incompressible parallel shear flows. In this context, linear stability theory (LST) gave rise to the well-known Orr–Sommerfeld equation (Orr Reference Orr1907; Sommerfeld Reference Sommerfeld1908) and related classical modal results, such as the critical Reynolds number
${\textit{Re}}_c = 5772.22$
for plane Poiseuille flow (Thomas Reference Thomas1953; Orzag Reference Orzag1971). Alternatively, energy stability-based methods yielded
${\textit{Re}}_c = 49.2$
(Busse Reference Busse1969; Joseph & Carmi Reference Joseph and Carmi1969), below which no energy growth was observed (Reddy & Henningson Reference Reddy and Henningson1993). Research on the stability of ideal gas compressible flows started later. For instance, Malik, Alam & Dey (Reference Malik, Alam and Dey2006) carried out LST selecting density and temperature as state variables along with the velocity vector and corresponding Jacobian matrices for compressible flows. They characterised the Y-shaped spectrum -- the so-called branches (Mack Reference Mack1976) -- and related even and odd modes. Over the past decades, variable-viscosity studies consisting of stratified, or Poiseuille flows with temperature dependency, have been performed based on a modified set of Orr–Sommerfeld equations (Potter & Graber Reference Potter and Graber1972; Govindarajan & Sahu Reference Govindarajan and Sahu2014). Related to viscosity effects, Wall (Reference Wall1996) investigated the stability limits at different conditions, whereas Malik, Dey & Alam (Reference Malik, Dey and Alam2008) demonstrated that viscosity stratification improves stability in compressible Couette flow, which was later confirmed by Saikia et al. (Reference Saikia, Ramachandran, Sinha and Govindarajan2017). Nonetheless, these studies were limited to either incompressible flow or ideal gas thermodynamics with temperature-dependent transport coefficients. In terms of wall-bounded compressible flow, Duck, Erlebacher & Hussaini (Reference Duck, Erlebacher and Hussaini1991) were some of the first authors to utilise LST in compressible plane Couette flow. In their work, destabilisation events were found even at the inviscid limit, which deviated noticeably from unbounded flows. Finally, recent studies (Gal et al. Reference Gal, Harlander, Borcia, Dizès, Chen and Favier2020; Cheng et al. Reference Cheng, Ma, Pullin and Luo2024) have leveraged LST to study the stability of generalised Couette-Poiseuille flows exposed to three-dimensional (3-D) oblique perturbations with buoyancy effects, as well as viscoelastic fluids in boundary layers (Cracco Reference Cracco2019). Moreover, by means of global stability analyses, Sarras et al. (Reference Sarras, Tayeh, Mons and Marquet2024) have investigated complex nonlinear events connected to flow hysteresis, low-frequency oscillations and 3-D stall cells.
In connection to complex thermophysical regimes, Ren et al. (Reference Ren, Marxen and Pecnik2019b ) have recently introduced a LST framework to study Poiseuille flows of non-ideal fluids. This study, in particular, studied the effects of varying (i) the dominant dimensionless numbers that characterise the flow, and (ii) the temperature of isothermal walls and bulk pressure. The isothermal limit (low Eckert and Prandtl numbers) results presented showed a good collapse with incompressible flow reference data. However, at larger Prandtl and Eckert numbers, the effects of large variations of the thermodynamic and transport properties driven by viscous heating of the flow became important. In particular, sub-/trans-/super-critical condition cases were compared against ideal gas scenarios. The results indicated that the base flow was modally more unstable in the subcritical regime, inviscid unstable in the transcritical regime and significantly more stable in the supercritical regime. Similar conclusions were extracted from algebraic growth analysis based on non-modal optimal energy growth with 3-D perturbations. It is important to highlight that the study, different to this work, focused on non-ideal gas destabilisation effects at relatively large Mach numbers for isothermal wall systems. Analogously, Zheng, Chen & Ding (Reference Zheng, Chen and Ding2024) carried out similar analyses for a Hagen--Poiseuille flow at supercritical regimes. More recently, Bugeat et al. (Reference Bugeat, Boldini, Hasan and Pecnik2024) have provided a quantitative characterisation of the stability diagram and dominant mechanisms for stratified plane Couette flow utilising supercritical fluids.
In particular, Ren et al. (Reference Ren, Fu and Pecnik2019a ) extended LST analyses to compressible boundary layers over adiabatic walls with fluids at supercritical pressure conditions. A second coexisting mode was found, causing flow destabilisation when crossing the pseudo-boiling line for two-dimensional (2-D) perturbations; prior, Ozgen & Kircali (Reference Ozgen and Kircali2008) validated the compressible LST framework utilising standard numerical methods against experimental data from the literature. In particular, the mode disappeared at supercritical conditions far from the pseudo-boiling region and at subcritical pressures. However, the effect of this mode on the value of the transition Reynolds number has not been characterised yet. Later, Ren & Kloker (Reference Ren and Kloker2022) expanded these studies to 3-D boundary layers. They found that at supercritical regimes, wall cooling stabilises the flow, whereas wall heating destabilises it. These effects were found to be reversed for the subcritical regime (same for ideal fluids). Finally, in the transcritical regime, wall heating stabilised the flow. In this regard, Gloerfelt et al. (Reference Gloerfelt, Robinet, Sciacovelli, Cinnella and Grasso2020) complemented the analysis for dense-gas effects, highlighting that complex thermophysical properties can lead to unconventional compressibility effects. In particular, the most important parameters were found to be the specific heat at constant pressure and the dynamic viscosity. Moreover, Wang et al. (Reference Wang, Liu, Ma and Yang2022) explored the difference between ideal and non-ideal gas LST frameworks for a mixing layer. They found that the most unstable frequency and the maximum spatial growth rate increase linearly with operating pressure, and its effect is diminished far from the critical point. In contrast, Guo et al. (Reference Guo, Gao, Jiang and Lee2020) found that existing correlations between time scales of the first mode and Reynolds numbers (based on the displacement thickness) are not suitable for supersonic flow regimes.
Therefore, the overall objective of this work is to perform linear stability and energy amplification analyses to characterise the stability conditions of high-pressure transcritical fluids wounded between walls at different temperatures. As observed in DNS studies at low-Mach-number flow regimes by Bernades et al. (Reference Bernades, Capuano and Jofre2023a
), notably higher levels of enstrophy were identified near the hot wall (
$50\times$
larger) in comparison to the cold wall. In addition, the enstrophy levels near the hot wall were
$100\times$
and
$10\times$
larger than equivalent laminar (iso-volumetric input power) and turbulent (iso-friction Reynolds number at cold/bottom wall) low-pressure systems, respectively. Likewise, Sahu & Matar (Reference Sahu and Matar2010) and Srivastava et al. (Reference Srivastava, Dalal, Sahu and Biswas2017) have explored the stability of low-pressure non-isothermal channel flow with viscous heating and found that increasing temperature across walls destabilised the flow. However, these investigations were limited to incompressible flow conditions using the Boussinesq approximation. Based on such approximation for stratified shear layers, Liu, Kaminski & Smyth (Reference Liu, Kaminski and Smyth2022) also performed a sensitivity analysis of this type of flow reporting maximum turbulence levels and mixing increase of approximately
$5\times$
and
$2\times$
, respectively, occurring when subharmonic and 3-D secondary instabilities were in strong competition. Furthermore, to identify optimum perturbation profiles, which yield turbulence enhancement and/or skin-friction drag reduction, modal and non-modal analyses are necessary. In this regard, focusing on temporal signals of incompressible Poiseuille flow, Massaro et al. (Reference Massaro, Martinelli, Schmid and Quadrio2023) explored the application of flow control techniques based on spanwise forcing. Typically, to achieve large destabilisation levels in isothermal wall flows, the base flow requires large Brinkman numbers as recently studied for Poiseuille flow using LST by Ren et al. (Reference Ren, Marxen and Pecnik2019b
). This, however, demands high-speed flows that are unfeasible in, for example, microconfined applications. Nonetheless, at non-isothermal conditions the system can operate at low velocities (and Brinkman numbers), while still exhibiting the destabilisation benefits of strong thermodynamic gradients occurring across the pseudo-boiling region. To this extent, this work is particularly focused on assessing the effects of subcritical, transcritical and supercritical thermodynamic regimes at isothermal and non-isothermal wall conditions operating at relatively low Reynolds and Mach numbers. In particular, with the aim of exploring the stability conditions at different wall temperatures and bulk pressures, the analyses consider results in terms of spectrum, unstable modes, range of instability, optimal perturbations and kinetic energy budget. Moreover, transient growth rates are also determined to quantify energy-based transient destabilisation regimes (Schmid Reference Schmid2007; Schmid & Brandt Reference Schmid and Brandt2014).
The paper, thus, is organised as follows. First, in § 2 the flow physics modelling is introduced. Next, the LST, linearised equations of supercritical fluids and discretisation methods utilised are described in § 3. Following, the set-up and parameters of the wall-bounded high-pressure transcritical DNS experiment performed are provided in § 4. The flow cases assessed by linear stability and modal and non-modal analyses are presented and discussed in § 5. Moreover, early instants of flow destabilisation are further examined in § 5.4 for non-isothermal channel flow configurations superimposing optimal perturbations from LST to a steady-state DNS solution. Finally, § 6 provides concluding remarks and proposes future research directions.
2. Flow physics modelling
The framework utilised in terms of (i) equations of fluid motion, (ii) real-gas thermodynamics and (iii) high-pressure transport coefficients is introduced below.
2.1. Equations of fluid motion
The flow motion of supercritical fluids is described by the following set of dimensionless equations of mass, momentum and total energy:
\begin{align} \frac {\partial \rho ^\star }{\partial t^\star } + \boldsymbol{\nabla} ^\star \boldsymbol{\cdot }\left ( \rho ^\star \boldsymbol{u}^\star \right ) & = 0, \end{align}
\begin{align} \frac {\partial \left ( \rho ^\star \boldsymbol{u}^\star \right )}{\partial t^\star } + \boldsymbol{\nabla} ^\star \boldsymbol{\cdot }\left ( \rho ^\star \boldsymbol{u}^\star \boldsymbol{u}^\star \right ) & = -\boldsymbol{\nabla} ^\star P^\star + \frac {1}{\textit{Re}} \boldsymbol{\nabla} ^\star \boldsymbol{\cdot }\boldsymbol{\tau ^\star } + \boldsymbol{F}^\star , \end{align}
\begin{align} \frac {\partial \left ( \rho ^\star E^\star \right ) }{\partial t^\star } + \boldsymbol{\nabla} ^\star \boldsymbol{\cdot }\left ( \rho ^\star \boldsymbol{u}^\star E^\star \right ) & = -\boldsymbol{\nabla} ^\star \boldsymbol{\cdot }(P^\star \boldsymbol{u}^\star ) - \frac {1}{{\textit{Re}} {\textit{Br}}} \boldsymbol{\nabla} ^\star \boldsymbol{\cdot }\boldsymbol{q}^\star \nonumber \\ & \quad + \frac {1}{\textit{Re}}\boldsymbol{\nabla} ^\star \boldsymbol{\cdot }(\boldsymbol{\tau ^\star } \boldsymbol{\cdot }\boldsymbol{u}^\star ) + \boldsymbol{F}^\star \boldsymbol{\cdot }\boldsymbol{u}^\star. \end{align}
Here superscript
$(\boldsymbol{\cdot })^\star$
denotes dimensionless quantities,
$\rho$
is the density,
$\boldsymbol{u}$
is the velocity vector,
$P$
is the pressure,
$\boldsymbol{\tau } = \mu ( \boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }\boldsymbol{u}^{T} ) + \lambda (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u})\boldsymbol{I}$
is the viscous stress tensor with
$\mu $
the dynamic viscosity,
$\lambda = - 2/3 \mu $
the bulk viscosity and
$\boldsymbol{I}$
the identity matrix,
$E = e + |\boldsymbol{u}|^{2}/2$
is the total energy with
$e$
the internal energy,
$\boldsymbol{q} = - {\kappa } \boldsymbol{\nabla }T$
is the Fourier heat flux with
$\kappa$
the thermal conductivity and
$T$
the temperature, and
$\boldsymbol{F}$
is a body force introduced to drive the flow in the streamwise direction.
The resulting set of scaled equations includes two dimensionless numbers: (i) the Reynolds number
${\textit{Re}} = \rho _{b}U_{r}\delta /\mu _{b}$
, where the subscript
$b$
refers to bulk quantities,
$\delta$
is the channel half-height and
$U_r$
is the reference streamwise velocity corresponding to its maximum value (i.e. centreline velocity for isothermal conditions), characterising the ratio between inertial and viscous forces; and (ii) the Brinkman number
${\textit{Br}} = \mu _b U_{r}^{2} / (\kappa T_b) = Pr Ec$
, quantifying the ratio of viscous heat generation to external heating through the walls (viz. larger
${\textit{Br}}$
values correspond to smaller heat conduction from viscous dissipation resulting in a temperature increase). The Brinkman number can also be expressed as the combination of Prandtl number
$Pr =\mu _{b} {c_{P}}_{b} / \kappa _{b}$
, where
$c_{P}$
is the isobaric heat capacity, expressing the ratio between momentum and thermal diffusivity, and Eckert number
$Ec = U_{r}^{2}/({c_{P}}_{b} T_{b})$
, which accounts for the ratio between advective mass transfer and heat dissipation potential. In this work, the Froude number, which represents the ratio between inertial and gravitational forces, is assumed to be large, and consequently buoyancy effects are not considered. The derivation of these dimensionless equations is based on the following set of inertial scalings (Jofre, del Rosario & Iaccarino Reference Jofre, del Rosario and Iaccarino2020, Reference Jofre, Bernades and Capuano2023b
):
\begin{align} \boldsymbol{x}^{\star } & = \frac {\boldsymbol{x}}{\delta }, \quad \boldsymbol{u}^{\star } = \frac {\boldsymbol{u}}{U_{r}}, \quad \rho ^{\star } = \frac {\rho }{\rho _{b}}, \quad T^{\star } = \frac {T}{T_{b}}, \quad P^{\star } = \frac {P}{\rho _{b} U_{r}^{2}}, \nonumber \\ & E^{\star } = \frac {E}{U_{r}^{2}}, \quad \mu ^{\star } = \frac {\mu }{\mu _{b}}, \quad \kappa ^{\star } = \frac {\kappa }{\kappa _{b}}, \quad \boldsymbol{F}^{\star } = \frac {\boldsymbol{F} U_{r} {\rho _b}^2}{\delta }. \end{align}
Here
$\boldsymbol{x}$
is position and
$H = 2 \delta$
is the total channel height.
2.2. Equation of state
The thermodynamic space of solutions for the state variables pressure
$P$
, temperature
$T$
and density
$\rho$
of a monocomponent substance is described by an equation of state. In the case of high-pressure systems, the Peng & Robinson (Reference Peng and Robinson1976) equation of state is a popular choice. In particular, high-pressure supercritical fluids are substances at temperatures and pressures above their critical values (
$T_c$
,
$P_c$
) that do not separate into different phases. Nonetheless, it is important to distinguish between supercritical gas-like and liquid-like fluids (Jofre & Urzay Reference Jofre and Urzay2021): (i) a supercritical liquid-like fluid is one whose density is large, whose isothermal compressibility is small and whose transport coefficients behave similar to a liquid; (ii) the density of supercritical gas-like fluids is smaller, their isothermal compressibility is larger and their transport coefficients vary similar to gases. At supercritical pressures, a smooth transition from liquid-like to gas-like supercritical behaviour occurs at the pseudo-boiling temperature
$T_{pb}$
, which corresponds approximately to the point where a second-order (continuous) phase transition occurs. In particular, such a smooth transition when crossing the pseudo-boiling line can be used to tune supercritical fluids to present liquid-like densities (
$\rho \sim 10^3$
kg m
$^{-3}$
) and gas-like viscosities (
$\mu \sim 10^{-5}$
Pa s), and therefore, favour inertial over viscous forces and resulting in turbulent flow (Bernades & Jofre Reference Bernades and Jofre2022). The rapid variations of thermophysical properties across the pseudo-boiling line are depicted in figure 28 (Appendix A). In particular, for a given pressure, the pseudo-boiling temperature corresponds to the point where isobaric specific heat capacity presents a peak-like maximum value. In the thermodynamic region around that temperature, large gradients in thermophysical properties are encountered.
Nonetheless, in this work, the CoolProp open-source library (Bell et al. Reference Bell, Wronski, Quoilin and Lemort2014) is used to describe the non-ideal thermophysical behaviour of high-pressure transcritical fluids. In particular, the thermodynamic quantities are derived from the Helmholtz energy equation of state. Their validation against NIST reference data (Linstrom & Mallard Reference Linstrom and Mallard2021) and sensitivity of the linear stability results to different thermodynamic frameworks (including ideal gas) is covered in Appendix A. In general form, the equation of state can be expressed in terms of the compressibility factor
$Z = P/(\rho R^{\prime } T)$
, where
$R^{\prime }$
is the specific gas constant. In dimensionless form, the generalised equation of state reads
\begin{equation} P^{\star } = \frac {Z \rho ^{\star } T^{\star }}{\hat {\gamma } M\!a_{b}^{2}}, \end{equation}
where
$\hat {\gamma } \approx Z (c_{P} / c_{V}) [(Z + T (\partial Z / \partial T)_\rho )/(Z+T(\partial Z / \partial T)_P)]$
is an approximated real-gas heat capacity ratio (Firoozabadi Reference Firoozabadi2016) with
$c_{V}$
the isochoric heat capacity. As it can be noted, the dimensionless bulk Mach number
$M\!a_{b} = u_b / c_b$
appears, where
$c_b$
is the bulk speed of sound, which represents the ratio of flow velocity to the local speed of sound. Real-gas equations of state need to be supplemented with the corresponding high-pressure thermodynamic variables (e.g. internal energy, heat capacities) based on departure functions (Reynolds & Colonna Reference Reynolds and Colonna2019) calculated as a difference between two states. These functions operate by transforming the thermodynamic variables from ideal gas conditions (low pressure – only temperature dependent) to supercritical conditions (high pressure). The ideal gas components are calculated by means of the NASA seven-coefficient polynomial (Burcat & Ruscic Reference Burcat and Ruscic2005), while the analytical departure expressions to high pressures are derived from the Peng-Robinson equation of state as detailed, for example, in Jofre & Urzay (Reference Jofre and Urzay2021).
2.3. High-pressure transport coefficients
The high pressures involved in the analyses conducted in this work prevent the use of simple relations for the calculation of dynamic viscosity
$\mu $
and thermal conductivity
$\kappa$
. Therefore, these quantities are also modelled by means of the CoolProp open-source library (Bell et al. Reference Bell, Wronski, Quoilin and Lemort2014). Alternatively, standard methods for computing these coefficients for Newtonian fluids are based on the correlation expressions proposed by Chung et al. (Reference Chung, Lee and Starling1984, Reference Chung, Ajlan, Lee and Starling1988). These correlations are mainly functions of the critical temperature
$T_c$
and density
$\rho _c$
, molecular weight
$W$
, acentric factor
$\omega$
, association factor
$\kappa _a$
and dipole moment
$\mathcal{M}$
, and the output from the NASA seven-coefficient polynomial (Burcat & Ruscic Reference Burcat and Ruscic2005); details can be found in dedicated works like, for example, Poling, Prausnitz & O’Connell (Reference Poling, Prausnitz and O’Connell2001) and Jofre & Urzay (Reference Jofre and Urzay2021). In this regard, similarly to the thermodynamic quantities, the differences in modelling the transport coefficients between CoolProp, NIST and Chung et al. correlations are presented in Appendix A.
3. Linear stability theory
The following subsections describe the linearised stability equations resulting from the flow model presented above and the corresponding discretisation method utilised to compute the results.
3.1. Linearised stability equations
The flow field can be decomposed into a base state and a perturbation denoted with subscript
$(\boldsymbol{\cdot })_0$
and superscript
$(\boldsymbol{\cdot })^\prime$
, respectively, yielding
\begin{equation} {} \boldsymbol{q} = \boldsymbol{q}_0 + \boldsymbol{q}^\prime , \end{equation}
where the vector
$\boldsymbol{q}$
is composed of five variables: three velocity components (
$u$
,
$v$
and
$w$
) and two independent thermodynamic variables (
$\rho$
and
$T$
). This decomposition yields a perturbation vector
$\boldsymbol{q}^\prime = (\rho ^\prime , u^\prime , v^\prime , w^\prime , T^\prime )^T$
superimposed to a base flow vector, which is assumed to be parallel to the walls, and consequently the only function of the wall-normal direction
$y$
in the form
$\boldsymbol{q}_0 = [\rho _0(y), u_0(y), 0, 0, T_0(y)]^T$
; the derivation of the reduced set of equations is detailed in Appendix B. From this point forward, the superscript
$(\boldsymbol{\cdot })^\star$
is omitted and it is assumed that all equations are in dimensionless form for the sake of exposition clarity.
The selection of the two main thermodynamic variables
$\rho$
and
$T$
allows the perturbations of the remaining thermodynamic variables (
$E, P, \mu , \kappa )$
to be expressed as a function of this pair of quantities by means of a Taylor expansion with respect to the base flow. For example, by neglecting higher-order terms, the pressure perturbation can be approximated by expanding the base flow pressure as
\begin{equation} P^\prime \approx \left .\frac {\partial P_0}{\partial \rho _0}\right \vert _{T_0} \rho ^\prime + \left .\frac {\partial P_0}{\partial T_0}\right \vert _{\rho _0} T^\prime , \end{equation}
with the resulting error of the approximation to be of the order of the amplitude (
$\epsilon$
) of the infinitesimal perturbations (Alves Reference Alves2016). Thus, by substituting (3.1) into the dimensionless equations of fluid motion ((2.1)–(2.3)), the linear stability equations are derived as a function of the perturbation vector
$\boldsymbol{q}^\prime$
and can be cast in compact form as
\begin{align} \boldsymbol{L_t} \frac {\partial \boldsymbol{q}^\prime }{{\partial t}} & + \boldsymbol{L_x} \frac {\partial \boldsymbol{q}^\prime }{{\partial x}} + \boldsymbol{L_y} \frac {\partial \boldsymbol{q}^\prime }{{\partial y}} + \boldsymbol{L_z} \frac {\partial \boldsymbol{q}^\prime }{{\partial z}} + \boldsymbol{L_q} \boldsymbol{q}^\prime \nonumber \\ & + \boldsymbol{V_{xx}} \frac {\partial ^2 \boldsymbol{q}^\prime }{{\partial x^2}} + \boldsymbol{V_{yy}} \frac {\partial ^2 \boldsymbol{q}^\prime }{{\partial y^2}} + \boldsymbol{V_{zz}} \frac {\partial ^2 \boldsymbol{q}^\prime }{{\partial z^2}} + \boldsymbol{V_{xy}} \frac {\partial ^2 \boldsymbol{q}^\prime }{{\partial x \partial y}} + \boldsymbol{V_{xz}} \frac {\partial ^2 \boldsymbol{q}^\prime }{{\partial x \partial z}} + \boldsymbol{V_{yz}} \frac {\partial ^2 \boldsymbol{q}^\prime }{{\partial y \partial z}} = 0, \end{align}
where all the nonlinear terms have been neglected and
$\boldsymbol{L_t}$
,
$\boldsymbol{L_x}$
,
$\boldsymbol{L_y}$
,
$\boldsymbol{L_z}$
,
$\boldsymbol{L_q}$
,
$\boldsymbol{V_{xx}}$
,
$\boldsymbol{V_{yy}}$
,
$\boldsymbol{V_{zz}}$
,
$\boldsymbol{V_{xy}}$
,
$\boldsymbol{V_{xz}}$
and
$\boldsymbol{V_{yz}}$
correspond to the Jacobian matrices of the base flow and thermophysical properties. In particular, these matrices are of size
$5\times 5$
, corresponding to each field of the perturbation vector
$\boldsymbol{q}^\prime$
. The components of these matrices can be obtained by inspection from the resulting linear stability equations presented in Appendix C. In particular, (C2) for continuity, (C4)–(C6) for momentum and (C9) for internal energy.
In linear modal stability analysis, the perturbation is assumed to have the ansatz
\begin{equation} \boldsymbol{q}^\prime (x,y,z,t) = \hat {\boldsymbol{q}}(y) e^{(i \alpha x + i \beta z - i \omega t)} + \textrm {c.c.}, \end{equation}
where
$\alpha$
and
$\beta$
are the streamwise and spanwise wavenumbers, respectively, while
$\omega$
is the temporal frequency whose real and imaginary parts correspond, respectively, to the wall-normal angular temporal frequency and its local growth rate, and c.c. stands for complex conjugate. Hence, by substituting (3.4) into (3.3), the following eigenvalue problem is obtained:
\begin{align} (- i \omega \boldsymbol{L_t} & + i \alpha \boldsymbol{L_x} + \boldsymbol{L_y} D + i \beta \boldsymbol{L_z} \boldsymbol{L_q} + \nonumber \\ & - \alpha ^2 \boldsymbol{V_{xx}} + i \alpha \boldsymbol{V_{xy}} D - \alpha \beta \boldsymbol{V_{xz}} + \boldsymbol{V_{yy}} D^2 + i \beta \boldsymbol{V_{yz}} D - \beta ^2 \boldsymbol{V_{zz}})\hat {\boldsymbol{q}} = 0. \end{align}
Here
$D$
is the derivative operator based on the Chebyshev discretisation presented in § 3.2. Equation (3.5) can be solved either on the temporal or spatial domain as
\begin{equation} \left . \begin{array}{ll} \displaystyle \boldsymbol{L_T} \hat {\boldsymbol{q}} = \omega \boldsymbol{L_t} \hat {\boldsymbol{q}},\\[8pt] \displaystyle \boldsymbol{L_S} \hat {\boldsymbol{q}} = \alpha (\beta \boldsymbol{V_{xz}} - i \boldsymbol{V_{xy}} D - i \boldsymbol{L_x})\hat {\boldsymbol{q}} + \alpha ^2 \boldsymbol{V_{xx}} \hat {\boldsymbol{q}}, \end{array}\right \} \end{equation}
where, by inspection from (3.5), the operator corresponding to the temporal (
$\boldsymbol{L_T}$
) and spatial (
$\boldsymbol{L_S}$
) matrices is written as
\begin{align} \boldsymbol{L_T} & = \alpha \boldsymbol{L_x} - i \boldsymbol{L_y} D + \beta \boldsymbol{L_z} - i \boldsymbol{L_q} \nonumber \\ &\quad + i \alpha ^2 \boldsymbol{V_{xx}} + \alpha \boldsymbol{V_{xy}} D + i \alpha \beta \boldsymbol{V_{xz}} - i \boldsymbol{V_{yy}} D^2 + \beta \boldsymbol{V_{yz}} D + i \beta ^2 \boldsymbol{V_{zz}}, \end{align}
\begin{align} \boldsymbol{L_S} & = - i \omega \boldsymbol{L_t} + \boldsymbol{L_y} D + i \beta \boldsymbol{L_z} + \boldsymbol{L_q} + \boldsymbol{V_{yy}} D^2 + i \beta \boldsymbol{V_{yz}} D - \beta ^2 \boldsymbol{V_{zz}}. \end{align}
For wall-bounded Poiseuille flow, the temporal problem is typically considered, and consequently the streamwise
$\alpha$
and spanwise
$\beta$
wavenumbers are prescribed, and the problem is solved for the eigenvalue
$\omega$
, with its real and imaginary parts corresponding, respectively, to the wall-normal wavenumber and its local growth rate. A comprehensive sensitivity analysis of the linear operator comparing ideal- and non-ideal gas frameworks is provided in Bernades (Reference Bernades2024). In brief, Appendix C covers the sensitivity of the linear operator with respect to maximum growth rate for modal analysis and transient growth for non-modal analysis.
3.2. Discretisation method
The discretisation of the linearised equations is based on Chebyshev collocation (Trefethen Reference Trefethen2000) with a domain spanning the interval
$0 \leqslant y/\delta \leqslant 2$
and discretised as
\begin{equation} y_j = \delta \left ( 1 - \textrm {cos} \frac {\pi j}{N} \right ), \quad j = 0, \ldots , N, \end{equation}
where
$N$
corresponds to the total number of collocation points. In this regard, Chebyshev differentiation matrices are utilised to obtain the discretised equations and define the LST eigenvalue problem operators. Specifically, the mesh size selected for this work is
$N = 200$
, which provides grid-independent results based on the convergence of the S-shaped Mack branches; the error scales with
$\mathcal{O}(h^2)$
, where
$h = H / N$
. Particularly, further increasing the grid size by
$50\,\%$
improves the accuracy by
$\mathcal{O}(10^{-4})$
; for brevity of exposition, the corresponding grid convergence results are not shown in this work. Moreover, the system of equations is subjected to
$u^\prime = v^\prime = w^\prime = T^\prime = 0$
boundary conditions for both walls.
3.3. Algebraic non-modal stability
The operator
$L_T$
in (3.6) is a non-normal matrix from which non-orthogonal eigenvectors and transient energy growth are also obtained. The individual eigenmodes of the modal system describe the behaviour of disturbances at a large asymptotic time, but they fail to capture the short-time dynamics. As a result, the linear superposition of the non-orthogonal eigenvectors exhibits substantial energy growth for a short period of time (Thummar, Bhoraniya & Narayanan Reference Thummar, Bhoraniya and Narayanan2024). Thus, non-modal stability theory, i.e. algebraic stability, based on the matrix exponential is used to capture the entire perturbation dynamics. Particularly, for algebraic stability analysis, an energy norm needs to be defined. In this regard, it is known that, for compressible flows, the Chu norm (Chu Reference Chu1965; Hanifi, Schmid & Henningson Reference Hanifi, Schmid and Henningson1996) is a well-suited candidate, which is defined as
\begin{equation} E(q) = \int \left [ ({u^\prime }^\dagger {u^\prime } + {v^\prime }^\dagger {v^\prime } + {w^\prime }^\dagger {w^\prime }) + m_\rho {\rho ^\prime }^\dagger {\rho ^\prime } + m_T {T^\prime }^\dagger {T^\prime } \right ] {\rm d}y, \end{equation}
where
$(\boldsymbol{\cdot })^\dagger$
denotes complex conjugate. For ideal gas thermodynamics, selecting the Mack’s energy norm with
$m_\rho = T_0 / [\rho _0^2 \gamma {M\!a}^2]$
and
$m_T = 1/(\gamma (\gamma - 1 )T_0 {M\!a}^2)$
is a common choice. However, it has been recently proposed that, for non-ideal fluids, the results are more robust when
$m_\rho = m_T = 1$
(Ren et al. Reference Ren, Marxen and Pecnik2019b
). Hence,
$E(\boldsymbol{q})$
represents the disturbance energy subjected to the eigenvector basis obtained from modal stability. Therefore, based on this norm, the algebraic growth rates and optimal energy amplification can be computed by performing a singular value decomposition (SVD) of the following quantity:
\begin{equation} G(t) = \text{max}_{\boldsymbol{q_0}} \frac {E[\boldsymbol{q}(t)]}{E(\boldsymbol{q_0})}. \end{equation}
It is important to highlight that, to increase the robustness and efficiency of the computations, when performing the SVD it is beneficial to exclude growth modes, i.e. eigenvalues, that could lead to non-transient growth (Hanifi et al. Reference Hanifi, Schmid and Henningson1996).
4. Computational DNS set-up
The DNS of high-pressure transcritical channel flow considered in § 5.4, conducted following the methodology outlined by Bernades et al. (Reference Bernades, Capuano and Jofre2023a
) using the RHEA flow solver (Jofre et al. Reference Jofre, Abdellatif and Oyarzun2023a
; Abdellatif et al. Reference Abdellatif, Ventosa-Molina, Grau, Torres and Jofre2023), is briefly introduced. The system operates with N
$_2$
at a supercritical bulk pressure of
$P_b/P_c = 2$
, confined between isothermal cold/bottom (
$cw$
) and hot/top (
$hw$
) walls; the system is periodic in the streamwise and spanwise directions. These walls are separated by a distance
$H = 2\delta$
, with
$\delta = 100\thinspace {\unicode{x03BC} \textrm {m}}$
representing the channel half-height. The wall temperatures are set to
$T_{cw}/T_c = 0.75$
and
$T_{hw}/T_c = 1.5$
, where the subscript
$(\boldsymbol{\cdot })_c$
denotes critical conditions. Hence, such a non-isothermal configuration forces the fluid to follow a transcritical thermodynamic trajectory across the pseudo-boiling region. The friction Reynolds number at the cold wall (
${Re}_{\tau ,cw}$
) is set to
$100$
(equivalent to
${Re}_b \approx 4000$
) to ensure fully developed turbulent conditions (Bernades & Jofre Reference Bernades and Jofre2022). The corresponding dimensional parameters and computational domain represent a DNS with spatial resolutions of
$\Delta x^{+} = 9.8$
and
$\Delta z^{+} = 3.3$
(normalised by the cold-wall friction velocity), and non-uniform stretching in the wall-normal direction. The vertical grid spacing ranges between
$0.2 \leqslant \Delta y^{+} \leqslant 2.1$
at the cold wall and
$0.5 \leqslant \Delta y^{+} \leqslant 5.6$
at the hot wall. The flow is initialised by steady-laminar flow to subsequently apply infinitesimal perturbations. The characteristic time or flow-through time is defined as
$t_{b} = L_x/u_{b} \sim \delta /u_{\tau ,bw}$
. The computational solution of the equations of fluid motion follows the numerical scheme developed by Bernades et al. (Reference Bernades, Jofre and Capuano2023b
). For time integration, a Courant--Friedrichs--Lewy number of approximately 0.9 is selected to ensure stability of the solution and proper extraction of data snapshots.
5. Results and discussion
The presentation and discussion of the results are covered in this section to characterise stability conditions. The flow cases studied are first introduced. Next, modal stability analyses for isothermal and non-isothermal conditions are investigated. Next, algebraic growth analyses are provided for 2-D and 3-D perturbations with the corresponding optimal inputs and responses. Finally, optimal linear perturbations are superimposed to a steady-state laminar DNS solution.
5.1. Flow cases
As mathematically introduced in § 3, wall-normal stability conditions are studied by means of a Poiseuille flow to isolate oblique 3-D effects from the analyses. In particular, as depicted in the sketch of figure 1, the Poiseuille flow is studied with two different base configurations: (i) isothermal wall conditions in which subcritical, transcritical or supercritical regimes are controlled by the bulk velocity and present a symmetric temperature profile with a maximum at the centreline of the channel; and (ii) non-isothermal (i.e. differentially heated) conditions by imposing a temperature difference between the walls that enforces the fluid to operate within transcritical trajectories. The complete list of cases considered and their flow parameters are provided in table 1. Particularly, figure 2 shows the converged base flows for the isothermal and non-isothermal cases studied. Details about the calculation of the base flows can be found in Appendix B, while a careful verification of the results at the isothermal limit is reported in Appendix D. It is important to highlight that (i) the velocity scaling is a function of
${\textit{Br}}$
and thermophysical quantities following the expression
$U_r = \sqrt {Br \thinspace \kappa _b \thinspace T_b/\mu _b}$
, and consequently it changes between cases as a result of the different levels of thermodynamic asymmetry between them; and (ii) case NI-5 has been designed to mimic the operation conditions of the microconfined high-pressure transcritical flow studied utilising DNS approaches by Bernades et al. (Reference Bernades, Capuano and Jofre2023a
), and consequently the Brinkman number is adjusted to
${\textit{Br}} = 5.6 \times 10^{-6}$
to obtain
$U_r = 1\thinspace \textrm {m}\, \rm s^{-1}$
. Subsequently, case NI-6 is the equivalent set-up at low pressure, which is known to be steady and laminar.
Figure 1. Sketch of the Poiseuille flow for isothermal (a) and non-isothermal (b) cases.
Figure 2. Base flow profiles of the isothermal cases listed in table 1 for dimensionless streamwise velocity (a) and reduced temperature (b) as a function of wall-normal direction.
Table 1. Base flow cases studied utilising LST. The first group of cases corresponds to symmetric Poiseuille flows with isothermal walls, whereas the second group considers non-isothermal (differentially heated) cases with different cold (
$cw$
) and hot (
$hw$
) wall temperatures. Note:
$\text{V-1}^{(\star )}$
is covered for isothermal limit analysis and verification in Appendix D.
5.2. Temporal modal stability analysis
This subsection aims at quantifying flow destabilisation when operating at high-pressure transcritical fluid conditions in comparison to subcritical and supercritical thermodynamic regimes under 2-D perturbations. In particular, the LST studies will consider cases with different base flows obtained from varying the Reynolds and Brinkman numbers as defined in table 1.
5.2.1. Isothermal wall cases
For the different isothermal wall cases studied (labelled I-1, I-2 and I-3 in table 1), the maximum normalised streamwise velocity and temperature of the base flows considered as a function of reduced wall temperature and Brinkman number are depicted in figure 3 (represented by solid, dashed and dashed–dotted black curves). The first observation to highlight is that, in the vicinity of the pseudo-boiling temperature at relatively large Brinkman numbers, the maximum normalised streamwise velocity is slightly reduced. However, it is important to note that, even if the normalised streamwise velocity decreases, the absolute streamwise velocity is characterised by large values. The second observation is that, for the isothermal wall set-ups, the only case operating across the pseudo-boiling region is the I-2 configuration.
Figure 3. Base flow in terms of maximum normalised streamwise velocity (a) and temperature (b) as a function of reduced wall temperature and Brinkman number for the cases listed in table 1.
Once the base flows have been characterised, the neutral curves for the different isothermal wall cases are shown in figure 4. First, the subcritical regime case I-1 (centreline temperature below
$T_{pb}$
) yields a reduced stability region with respect to the isothermal limit equivalent set-up when
${\textit{Br}}$
increases. In particular, the neutral curve expands to lower
${\textit{Re}}$
values for a wider range of wavenumbers and becomes unstable for
${\textit{Re}}_c \approx 1000$
at
$1.2 \lesssim \alpha \lesssim 1.4$
; see table 2. Next, the I-2 base flow undergoes a transcritical trajectory across the pseudo-boiling region and, as a result, both velocity and temperature become inflectional at
${\textit{Br}} \gtrsim 0.35$
. This results in a lower critical Reynolds number of value
${\textit{Re}}_c \approx 400$
(refer to table 2) and a significantly wider range of wavenumbers,
$1.2 \lesssim \alpha \lesssim 1.6$
, where early flow destabilisation may likely occur. Third, the I-3 case operates at supercritical conditions behaving similarly to the ideal gas solution (Ren et al. Reference Ren, Marxen and Pecnik2019b
), in which increasing
${\textit{Br}}$
enlarges the stability region. In particular, for
${\textit{Br}} \gt 0.1$
, the flow is stable for the entire
${\textit{Re}}-\alpha$
parameter space considered in this work.
Figure 4. Neutral curves for various Brinkman numbers at (a) subcritical, (b) transcritical and (c) supercritical regimes. The dashed–dotted line represents the isothermal limit (
${\textit{Br}} \xrightarrow {} 0$
) and the vertical dashed line indicates
${\textit{Re}}_c = 5772$
.
Table 2. Critical Reynolds numbers for the flow cases described in § 5.1 obtained from modal stability analysis. Here NI-5 and NI-6 are assessed at
${\textit{Br}} = 5.6 \times 10^{-6}$
, i.e. low
${\textit{Br}}$
similar to the isothermal limit (
${\textit{Br}} \xrightarrow {}0$
).
In connection to the neutral curves introduced above, figure 5 presents the perturbations of the most unstable mode for
${\textit{Re}} = 10\,000$
and
$\alpha = 1$
along the wall-normal direction for various Brinkman numbers at different thermodynamic conditions. As can be observed, unlike for the isothermal limit case, the subcritical thermodynamic regimes are dominated by density and temperature perturbations (thermophysical-driven mode). The effects of these perturbations, however, are diminished when the Brinkman number increases and, consequently, the streamwise velocity dominates. At transcritical conditions, instead, streamwise velocity perturbations dominate over the other fields the destabilisation of the flow in the vicinity of the walls, except for large Brinkman numbers where, again, the thermodynamic perturbations mandate. However, the wall-normal velocity governs flow destabilisation at the centreline independently of
${\textit{Br}}$
. In particular, figure 5(b) captures the perturbation dominance of this velocity component at the centre, where the other orthogonal velocity components and thermodynamic terms are not as important. Differently, for supercritical thermodynamic conditions, flow destabilisation is mainly dominated throughout the entire wall-normal direction by streamwise velocity perturbations (hydrodynamic-driven mode) with insignificant differences with Brinkman number increase regarding the hydrodynamic perturbation, but slowly rise in the vicinity of the walls for the thermodynamic perturbations.
Figure 5. Perturbation profiles of the most unstable mode at
${\textit{Re}} = 10\,000$
and
$\alpha = 1$
along the wall-normal direction for various Brinkman numbers at (a) subcritical, (b) transcritical and (c) supercritical regimes.
5.2.2. Non-isothermal cases
Figure 6 depicts the neutral curves for the non-isothermal cases. It is noted that in this case the fluid is forced to cross the pseudo-boiling region for all
${\textit{Br}}$
. In fact, the largest value of the Brinkman number was chosen to ensure that temperature remains below the hot-wall temperature, which corresponds to
${\textit{Br}} \leqslant 0.1$
. Particularly, while NI-1 and NI-3 cross the pseudo-boiling temperature, NI-2 is pressurised much beyond the critical pressure, and consequently it operates at supercritical conditions. The neutral curves of NI-1 are similar for all
${\textit{Br}}$
, where flow destabilisation is biased toward lower
${\textit{Re}}$
; especially, in comparison to the low-Brinkman-number cases at isothermal wall conditions. In particular, the wavenumber corresponding to the critical Reynolds number falls by roughly
$20\,\,\%$
. Nevertheless, the NI-3 case results in a larger
${\textit{Re}}_c$
value, which is above the isothermal limit transition when
${\textit{Br}} \lt 0.05$
. Therefore, as observed for the NI-1 case, by constraining the cold- and hot-wall temperatures closer to the pseudo-boiling temperature, the stability region is enhanced. Finally, when increasing the pressure of the system, the neutral curves display similar envelopes as for the subcritical isothermal case.
Figure 6. Neutral curves of cases (a) NI-1, (b) NI-2 and (c) NI-3 for different Brinkman numbers.
Figure 7. Perturbation profiles of the most unstable mode at
${\textit{Re}} = 10\,000$
along the wall-normal direction for various
${\textit{Br}}$
: (a) NI-1 (
$\alpha = 0.8$
), (b) NI-2 (
$\alpha = 1$
) and (c) NI-3 (
$\alpha = 0.8$
) cases.
Figure 8. Perturbation profiles of the most unstable mode for (a) NI-1 at
${\textit{Re}} = 4000$
and
$\alpha = 0.8$
, (b) NI-2 at
${\textit{Re}} = 4000$
and
$\alpha = 1$
, and (c) NI-3 at
${\textit{Re}} = 8000$
and
$\alpha = 0.8$
.
Analogously, figure 7 presents the perturbation profiles of the most unstable modes at the largest
${\textit{Re}}$
analysed. The perturbations are dominated by density and temperature variations near the pseudo-boiling region, but velocity is significantly high near the hot wall. On the other hand, at the cold wall, the thermodynamic and hydrodynamic modes have similar magnitudes, although they are
$50\,\,\%$
lower than at the hot wall. This behaviour is similar for both transcritical cases, NI-1 and NI-3, with small differences among
${\textit{Br}}$
numbers for NI-1. Although NI-3 is dominated by thermodynamic perturbations for
${\textit{Br}} = 0.01$
and
${\textit{Br}} = 0.05$
, which are roughly
$10\times$
and
$30\times$
greater in magnitude, compared with velocity perturbations. Nonetheless, at
${\textit{Br}} = 0.1$
, the system is governed again by streamwise perturbations in the vicinity of the hot wall. Instead, NI-2 is mainly dominated by the streamwise velocity near the walls and by the wall-normal velocity at the centre. Nevertheless, at large
${\textit{Br}}$
, the density and temperature perturbations near the cold wall become greater than those associated with the streamwise velocity. Additionally, the perturbations near the critical Reynolds value
${\textit{Re}}_c$
are depicted in figure 8. The NI-1 perturbations are similar to those obtained with larger
${\textit{Re}}$
values for low
${\textit{Br}}$
numbers. Instead, at
${\textit{Br}} = 0.05$
, thermodynamic modes dominate and at
${\textit{Br}} = 0.1$
they become negligible with respect to the streamwise velocity perturbation. However, for NI-2, the system is dominated by thermodynamics at
$y/\delta \approx 0.75$
at low
${\textit{Br}}$
. The larger the Brinkman number, the closer the effects of hydrodynamic and thermodynamic perturbations. Moreover, at
$y/\delta \approx 1.75$
, streamwise velocity perturbations dominate; in contrast, NI-3 is clearly dominated by thermodynamic modes at both peak locations. Finally, the NI-4 case operates at low-pressure conditions, and the linear stability results collapse to those corresponding to the ideal gas solution. Nonetheless, although it operates at different temperatures, the system is stable for all Reynolds numbers and wavenumbers considered.
Figure 9. Growth rate contours at
${\textit{Re}}-\alpha$
for (a) NI-5 and (b) NI-6 cases. The neutral curve for NI-5 is highlighted in red.
5.2.3. High-pressure transcritical versus superheated steam non-isothermal Poiseuille flow
This section compares the non-isothermal flow conditions at low- and high-pressure regimes operating at low
${\textit{Br}}$
, namely the NI-5 and NI-6 cases. As concluded in previous sections, non-isothermal flows are subject to larger destabilisation regions at low
${\textit{Br}}$
numbers. In particular, NI-5 is equivalent to the transcritical channel flow set-up studied via DNS by Bernades et al. (Reference Bernades, Capuano and Jofre2023a
), where the base flow is adjusted to obtain a reference velocity of
$U_r \approx 1\thinspace \textrm {m}\,\rm s^{-1}$
. The results show that, as expected, low Mach/Reynolds regimes are susceptible of an earlier flow destabilisation when operating at non-isothermal conditions, which confirms the hypothesis that flow destabilisation may occur earlier for non-isothermal wall-bounded systems. Figure 9(a) depicts the early modal-based destabilisation exhibited by NI-5 when reaching
${\textit{Re}}_c = 2000$
independently of the streamwise wavenumber. In particular, the perturbations near the pseudo-boiling point of solely the unstable mode (i.e.
${\textit{Re}}_b = 4000$
) are mainly dominated by the streamwise velocity, although the thermodynamic variables follow a similar trajectory (at lower magnitude) as shown in figure 10(a). At
${\textit{Re}} = 10\,000$
, the hydrodynamic components are still fully governing the destabilisation process (figure 10
c) specially near the hot wall. On the contrary, moving toward lower Reynolds numbers, all the analysed modes become stable. The resulting early stages of flow destabilisation are explored in § 5.4 under 2-D modal analysis perturbations and subsequent non-modal transient growth. In particular, based on DNS computations, the initial disturbance is given by the Fourier modes following
$u (\alpha ,y,0) = \epsilon \thinspace q (\alpha ,y)$
, where
$\epsilon$
is the infinitesimal amplitude (Yamamoto, Takahashi & Kambe Reference Yamamoto, Takahashi and Kambe2001) for the normal stability perturbations. This amplitude of the forcing is typically and properly chosen (between
$10^{-5}$
and
$10^{-3}$
of the free-stream Mach number) such that the excited perturbations stay in the linear regime. It is known that equivalent low-pressure set-ups become laminar, and therefore, are stable (Bernades, Capuano & Jofre Reference Bernades, Capuano and Jofre2024). As observed in the previous section, the non-isothermal superheated steam system becomes stable for
${\textit{Re}} \leqslant 10\,000$
. In particular, the NI-6 case is driven with the same
${\textit{Br}}$
number as NI-5 and its growth rate is depicted in figure 9(b), capturing the entire streamwise wavenumber range modal stability. The perturbations are dominated by the streamwise velocity for the least stable mode at
${\textit{Re}} = 10\,000$
and
${\textit{Re}} = 4000$
as depicted in figure 10(b,d), although thermodynamic components slowly grow their magnitudes at larger
${\textit{Re}}$
.
Figure 10. Perturbation profiles of the most unstable mode at
$\alpha = 0.8$
for (a) NI-5 at
${\textit{Re}} = 4000$
, (b) NI-6 at
${\textit{Re}} = 4000$
, (c) NI-5 at
${\textit{Re}} = 10\,000$
, and (d) NI-6 at
${\textit{Re}} = 10\,000$
.
5.2.4. Flow destabilisation summary
The critical Reynolds numbers for each case, extracted from the corresponding neutral curves, are summarised in table 2. These results are derived from modal stability analyses, and consequently are an indicator of the destabilisation trends. In addition, to characterise flow destabilisation thresholds, the algebraic stability needs to be computed.
5.2.5. Energy budgets
The decomposition in different terms of the kinetic energy growth further facilitates the characterisation of the destabilisation conditions. In particular, it identifies the main contributors to the temporal kinetic energy growth (
$K$
): production (
$P$
), thermodynamic (
$T$
) and viscous dissipation (
$V$
). It is important to highlight that the power of the external force introduced to drive the flow cancels out with the viscous dissipation and it has been consequently omitted from the budget, resulting in the growth balance expression
$K = P + T + V$
(Andreolli, Quadrio & Gatti Reference Andreolli, Quadrio and Gatti2021). The energy balance for the 2-D perturbations is derived by taking the inner product of the streamwise and spanwise momentum conservation equations with their respective velocity components for the most unstable mode (Boomkamp & Miesen Reference Boomkamp and Miesen1996). The temporal growth of density, embedded in the streamwise momentum, is obtained from the continuity equation. The resulting equation is typically averaged over the wavelength and integrated over the height of the channel (Moeleker Reference Moeleker1998; Sahu & Matar Reference Sahu and Matar2010; Ren et al. Reference Ren, Marxen and Pecnik2019b
), considering only the real part; the corresponding detailed expressions are presented in Appendix C.
The budgets are assessed at
${\textit{Re}} = 10\,000$
to represent the unstable modes of the spectra at their corresponding streamwise wavenumber ranges. It is well known that for single-phase Poiseuille flow, destabilisation is caused by a combination of the no-slip conditions at the boundaries and the viscous effects in the critical layer responsible for the generation of Reynolds stresses (Drazin & Reid Reference Drazin and Reid1981). In this regard, table 3 quantifies the different contribution terms for each case. Several observations can be extracted from this table: (i) production becomes the principal contributor to the kinetic energy growth scaling with
${\textit{Br}}$
except for the isothermal transcritical and supercritical base flows; (ii) in particular, the I-2 velocity perturbations are limited despite the large inflection points emerging in the vicinity of the walls at high-Brinkman-number base flows due to compressibility effects (Xie, Karimi & Girimaji Reference Xie, Karimi and Girimaji2017) and, as a result, the net production term decreases (the analysis of this result is extended in § 5.2.6); (iii) the thermodynamic term is negligible in comparison to the other terms and negative; (iv) the viscous term also decreases the temporal energy and barely changes with the Brinkman number (at low
${\textit{Br}}$
numbers, it weights approximately
$60\,\,\%$
of the production); (v) the viscous dissipation increases with
${\textit{Br}}$
at isothermal wall transcritical and supercritical conditions and for non-isothermal superheated steam; (vi) the stable flow cases result in a negative growth and an increase of the thermodynamic contribution; and (vii) the energy growth for the non-isothermal cases is barely sensitive to the Brinkman number in comparison to the isothermal wall set-ups. Similarly, Gloerfelt et al. (Reference Gloerfelt, Robinet, Sciacovelli, Cinnella and Grasso2020) noted that for boundary layers at subsonic regimes, the viscous mode was attenuated by compressibility effects favouring the presence of the first mode related to inviscid mechanisms. However, at supersonic speeds, the viscous mode was stable and independent of the thermodynamic conditions.
Table 3. Summary of kinetic energy budgets (
$\times 10^{-4}$
) for the 2-D perturbations at
${\textit{Re}} = 10\,000$
. Isothermal cases at
$\alpha = 1$
, non-isothermal NI-1–5 at
$\alpha = 0.6$
and NI-2–3-4–6 at
$\alpha = 0.8$
.
Figure 11 characterises the spatial growth along the wall-normal direction. First, the differences across the isothermal wall set-ups are depicted in the top row up to the half-plane symmetry. In particular, figure 11(a) shows the I-2 case at
${\textit{Br}} = 0.1$
(below the pseudo-boiling region) and
${\textit{Br}} = 0.5$
(crossing the pseudo-boiling region). As anticipated, it is found that the production term decreases (even in non-dimensionless form) due to the supersonic regimes achieved at large
${\textit{Br}}$
. Nonetheless, unlike at lower
${\textit{Br}}$
values, the production term rapidly grows away from the centreline achieving its maximum at the inflection point, where, in contrast, the thermodynamic term reaches its minimum. In general, the thermodynamic contribution dominates at the centreline and in the vicinity of the walls. Moreover, regardless of the
${\textit{Br}}$
value, the thermodynamic term introduces flow destabilisation near the walls, while viscous dissipation stabilises the flow at similar rates resulting in the cancellation of each other to provide a null net kinetic energy growth. In terms of subcritical versus supercritical behaviour, figure 11(b) compares I-1 against I-3 at large
${\textit{Br}}$
values. It can be clearly observed that cases at subcritical thermodynamic conditions tend to destabilise the flow, whereas the opposite is obtained at supercritical regimes. In particular, at subcritical conditions, production is important toward the walls and decreases toward the centre, whereas the supercritical cases present an opposite behaviour changing sign at the centreline and reaching its maximum (positive) in the vicinity of the walls. Moreover, the thermodynamic contribution is largest at the channel half-height for the subcritical set-ups, yielding negative values where production peaks and becoming once again the largest destabilisation contributor near the walls. Instead, for the I-3 case, the thermodynamic term is negative at the centreline and increases toward the boundaries becoming the main instability source. This behaviour, however, is modified in the wall-normal location where production peaks and the thermodynamic contribution reverses its effect following a double valley-peak trajectory. However, this change of behaviour is not important enough to obtain a positive kinetic energy growth. On the other hand, viscous dissipation stabilises the flow near the walls in both scenarios. Finally, the second row of the figure shows the behaviour of the non-isothermal cases. No important differences can be observed as the Brinkman number changes between cases. In this regard, figure 11(c)–(d) depicts the energy budget contributions at
${\textit{Br}} = 0.1$
for cases NI-1 and NI-3. Interestingly, production becomes the key mechanism in regions of relatively large density near the cold wall and achieves its negative peak near the pseudo-boiling region. This stable region is wider for NI-2, which presents larger temperature difference across walls. Moreover, case NI-3 (characterised by a smaller temperature difference between walls) yields lower production and thermodynamic peaks, but even so the production term is responsible for destabilising the flow near the hot wall. Thus, the thermodynamic term becomes the principal budget contributor in the vicinity of the walls, peaking at the pseudo-boiling region where it changes sign similar to the production term. Oppositely, viscous dissipation stabilises the flow in the vicinity of the walls and vanishes elsewhere. Hence, it is clear that for isothermal wall flow set-ups, the production term dominates the kinetic energy budget, which is similar to what happens in the cold region (near the wall) of the non-isothermal cases. Nonetheless, in the vicinity of the hot wall and pseudo-boiling region, flow destabilisation is dominated by the thermodynamic term.
Figure 11. Energy budget terms along the wall-normal direction at
${\textit{Re}} = 10\,000$
for (a) I-2 at
${\textit{Br}} = 0.1$
(left) and
${\textit{Br}} = 0.5$
(right), (b) I-1 (left) and I-3 (right) both at
${\textit{Br}} = 0.5$
, (c) NI-1 at
${\textit{Br}} = 0.1$
, and (d) NI-3 at
${\textit{Br}} = 0.1$
.
5.2.6. Production term
The growth in kinetic energy is mostly governed by the production term (Boomkamp & Miesen Reference Boomkamp and Miesen1996; Xie et al. Reference Xie, Karimi and Girimaji2017; Samanta Reference Samanta2020; Andreolli et al. Reference Andreolli, Quadrio and Gatti2021). It plays a crucial role in controlling disturbances by supplying energy to them by means of the Reynolds stresses. In this regard, it is helpful to examine the vorticity transport equation (Xie et al. Reference Xie, Karimi and Girimaji2017), and in particular, the components dominating the spanwise vorticity perturbations, which are relevant for shear flow instabilities. Thus, this section aims at carefully studying these effects.
The importance of the production term increases with larger energy inputs, i.e. higher Brinkman numbers. However, it is found that for the transcritical isothermal wall case, this term diminishes for
${\textit{Br}} \geqslant 0.1$
. In this regard, the singular behaviour of the production term in the case of isothermal wall transcritical high-speed flows is briefly characterised. Therefore, despite this work focusing mainly on the stability of high-pressure transcritical fluids at low
$M\!a$
and
${\textit{Re}}$
conditions, the large Brinkman numbers resulting in the isothermal wall cases exhibit
$M\!a \geqslant 1$
. In this regard, it is known that for incompressible flows, the instability mechanism related to streamwise perturbations is driven by the Tollmien–Schlichting (T-S) waves. However, for high-speed flows, the T-S waves are typically subjected to the effects of compressibility effects. This type of effect can be quantified by means of the gradient Mach number (Sarkar Reference Sarkar1995)
${M\!a}_g = [\partial u_0 / \partial y] / [c_0 \sqrt (\alpha ^2 + \beta ^2)]$
. In particular, given that obliqueness effects are null for the current 2-D modal assessment (
$\beta = 0$
),
${M\!a}_g$
is the dominant factor controlling the production rate. Generally, under such circumstances, the gradient Mach number exceeds unity near the walls, and consequently these regions are subjected to strong pressure waves and dilatational velocity fluctuations, exhibiting different flow behaviour compared with low-speed regimes. In this context, figure 12 depicts the behaviour of case I-2 varying from
${\textit{Br}} = 0.01$
to
${\textit{Br}} = 0.5$
. As can be observed, large
${M\!a}_g$
values are encountered near the boundaries, and as a result, the fluctuations in streamwise and wall-normal directions are significantly attenuated. Therefore, the production term of the energy budget is notably reduced yielding a decay of kinetic energy growth. Unlike for non-isothermal wall set-ups, isothermal wall cases are subjected to this attenuation on the production terms, which result in lower kinetic energy values to achieve flow destabilisation. This highlights the potential of transcritical regimes driven by temperature differences imposed across the boundaries. The inherently nonlinear nature of destabilisation in such flows emerges even at very low Mach and Reynolds numbers, leading to significantly different instability mechanisms and energy budgets compared with isothermal cases.
Figure 12. Energy budget compressibility effects for (a) the I-2 case at
${\textit{Br}} = 0.01$
(
${M\!a}_b = 0.1$
and
$U_r = 42.3\, \rm m\,\rm s^{-1}$
) and
${\textit{Br}} = 0.5$
(
${M\!a}_b = 1.4$
and
$U_r = 351\, \rm m\,\rm s^{-1}$
) depicting the evolution of
${M\!a}_g$
,
$u^\prime$
,
$v^\prime$
and
$P^\prime$
along the wall-normal direction. The perturbations are scaled by a factor of
$5\times$
for visualisation. For comparison, the subcritical I-1 case (b) is also reported.
The spanwise vorticity is decomposed into three main components (Xie et al. Reference Xie, Karimi and Girimaji2017): (i) the baroclinic effect
$(\partial \rho _0 / \partial y)(\partial \rho ^\prime / \partial x) / {\rho _0}^2$
, (ii) the compressible vortex production
$(\partial u_0 / \partial y)(\partial u^\prime / \partial x + \partial v^\prime / \partial y)$
, and (iii) the second derivative effect
$v^\prime (\partial ^2 u_0 / \partial ^2 y)$
. Among these contributions, shear flow destabilisation is typically critical for the second derivative of vorticity and production field from the velocity equation
$v^\prime \partial u / \partial y$
. These multiple mechanisms can be identified in figure 13. On the one hand, the isothermal wall transcritical case is depicted in figure 13(a) to clarify the importance of both the vorticity second derivative (especially at the channel centreline) and production terms. The baroclinic effect and compressible vortex production are one order of magnitude smaller. Interestingly, the larger the
${\textit{Br}}$
number the more important these two terms become in detriment of the dominant terms. On the other hand, figure 13(b) compares NI-1 at
${\textit{Br}} = 0.1$
and NI-5. For both cases, the second derivative of vorticity emerges as the principal destabilisation mechanism, especially near the pseudo-boiling region. In particular, in the case of low-speed conditions, the baroclinic effect achieves similar magnitude as for larger velocity flows at the pseudo-boiling region. It is important to note that, when the temperature difference between the walls is reduced (case NI-3), the magnitude of all terms become significantly smaller, especially the velocity production. The baroclinic effect, however, preserves similar magnitudes at the pseudo-boiling region.
Figure 13. Energy terms along the wall-normal direction at
${\textit{Re}} = 10\,000$
for (a) I-1 at
${\textit{Br}} = 0.1$
and
${\textit{Br}} = 0.5$
, and (b) NI-1 at
${\textit{Br}} = 0.1$
and NI-5. The baroclinic effect and compressible vortex production are scaled by a factor of
$10\times$
for visualisation.
5.3. Transient energy growth analysis
Linear modal stability is limited to the asymptotic behaviour of disturbances, omitting large transient amplifications at short times. This short-time energy growth can be significant compared with the initial perturbation level. Therefore, to infer more insights of flow stability characteristics, the non-normality of the linearised equations of the fluid motion operator needs to be taken into account (Schmid Reference Schmid2007; Ghosh et al. Reference Ghosh, Loiseau, Breugem and Brandt2019; Cabrit Reference Cabrit2020). In this regard, based on 3-D perturbations (
$\beta \gt 0$
), algebraic instabilities are quantified for the various flow cases studied. Particularly, this section characterises (i) the maximum growth rate and transient amplifications, and (ii) the optimal perturbation profiles that generate the maximum growth rate; for completeness, validation against incompressible (Reddy & Henningson Reference Reddy and Henningson1993) and isothermal wall (Ren et al. Reference Ren, Marxen and Pecnik2019b
) high-pressure Poiseuille flow cases are reported in Appendix E.
5.3.1. Maximum growth rate and transient amplification
The effects of the Brinkman number and thermodynamic regime at isothermal and non-isothermal flow conditions on the non-modal growth of the disturbance energy is investigated. The Reynolds number
${\textit{Re}} = 1000$
is chosen, far below the transition point based on the linear stability analyses summarised in table 1. First, the transient growth envelopes on the
$\alpha$
–
$\beta$
space for the isothermal wall set-ups are detailed in figure 14. They reveal that the kinetic energy of the perturbations grows by a factor of
$200$
when operating at subcritical and supercritical regimes. Furthermore, the kinetic energy can be enhanced by a factor of
$800\times$
when applying an invariant streamwise perturbation (
$\alpha = 0$
), and even become inviscidly unstable at
$\alpha = \beta = 1$
. In fact, a similar behaviour was observed by Schmid (Reference Schmid2007) for Poiseuille flow under 3-D perturbations at
${\textit{Re}} = 10\,000$
. Particularly, the largest (asymptotic exponential) transient growth occurs for streamwise independent perturbations at
$\alpha = 0$
and
$\beta = 2$
. Nonetheless, it is noted that case I-2 deviates from this trend at
${\textit{Br}} = 0.5$
, showcasing a region in which the maximum growth tends to infinity for
$\alpha = \beta = 1$
. In particular, the spectra yield an unstable mode responsible for the exponential growth rate over time, and consequently the critical Reynolds number is notably reduced. At
${\textit{Br}} = 0.25$
, this behaviour is no longer captured and smaller
${\textit{Br}}$
values result in maximum growth localised only at the invariant streamwise perturbation. Moreover, the colour maps of cases I-1 and I-3 at low
${\textit{Br}}$
values behave similarly to the results previously presented. In particular, they have an indistinguishable shape and maxima with a slight reduction of the maximum growth region, which is localised at
$\alpha = 0$
and
$\beta = 2$
. Hence, for the sake of brevity, only results for the largest Brinkman number
${\textit{Br}} = 0.5$
are depicted. Figure 15 depicts the evolution of the growth rate over time at the maximum growth rate on the
$\alpha -\beta$
space. As can be observed, the resulting growth rate presents a monotonically increasing trend reaching the maximum destabilisation energy after it starts attenuating. In particular, cases I-1 and I-3 develop a similar energy behaviour for the different base flows analysed. However, the transcritical case evolves to larger amplifications, which is exacerbated for
$\alpha = \beta = 1$
.
Figure 14. Transient growth envelopes at
${\textit{Re}} = 1000$
and
${\textit{Br}} = 0.5$
for (a) I-1, (b) I-2 and (c) I-3 cases. The resolution range for case I-2 has been limited to
$G_{\textit{max}} = 1000$
.
Figure 15. Amplification rates for
${\textit{Re}} = 1000$
at
$\alpha = 0$
and
$\beta = 2$
for (a) I-1, (b) I-2 and (c) I-3 cases. Inset of (b) corresponds to the exponential amplification rate at
$\alpha = \beta =1$
for case I-2.
Analogously to the discussion above focused on the isothermal wall cases, the transient growth of the non-isothermal NI-1, NI-2 and NI-3 cases are shown in figure 16. It can be observed that there is a clear range at
$\alpha = 0$
and
$1.0 \leqslant \beta \leqslant 3.0$
where amplification is large despite corresponding to fully stable spectra. Nonetheless, this range is notably reduced when operating at
$P_b/P_c = 5$
(case NI-2) for the optimum wavenumber
$\alpha = 0$
and
$\beta = 2$
in comparison to the isothermal wall case. In particular, the amplification is significantly reduced by a factor of roughly
$8\times$
. Moreover, narrowing the temperature operating window reduces the amplification by approximately
$5\times$
, but the energy boost from the non-modal instability is longer sustained within the system as quantified by the amplification rate evolution in figure 17. As previously introduced, Brinkman numbers are limited to
${\textit{Br}} \leqslant 0.1$
so that temperatures are constrained within the range imposed by the wall temperatures. Therefore, it is important to note that, when
${\textit{Br}}$
increases, the transient growth significantly increases along a large region of peak locations for
$\beta = 2$
and
$\alpha \geqslant 0$
. In this context, the results align with findings from the DNS study (Bernades et al. Reference Bernades, Capuano and Jofre2023a
), which reported flow destabilisation near the hot wall, where vorticity dissipation was approximately two orders of magnitude higher than in an equivalent laminar case. Notably, vortex stretching acted as a generator of turbulence by tilting and elongating vortices, shaping the topology of small-scale structures. The flow discriminant showed significantly larger positive values near the hot region (approximately 100 times higher), indicating that enstrophy dominates throughout the channel. More specifically, the analysis of flow invariants revealed a tendency toward a mild expansion across the channel, while becoming strongly negative near the hot wall. This behaviour corresponds to localised compression states in which the flow follows stable trajectories (Suman & Girimaji Reference Suman and Girimaji2010). In connection with the DNS results, the optimal output responses (figure 25
g,h,i) illustrate elongated structures and compression zones between the upper boundary and the pseudo-boiling region, with reversed flow directions toward the channel centre.
Figure 16. Transient growth envelopes at
${\textit{Re}} = 1000$
and
${\textit{Br}} = 0.1$
for (a) NI-1, (b) NI-2 and (c) NI-3 cases.
Figure 17. Amplification rates for
${\textit{Re}} = 1000$
at
$\alpha = 0$
and
$\beta = 2$
for (a) NI-1, (b) NI-2 and (c) NI-3 cases.
The different destabilisation properties of the non-isothermal cases with respect to the isothermal wall set-ups for low Reynolds and Mach number conditions has been discussed in § 5.2. Therefore, non-modal analyses are carried out only for cases NI-5 (
$P_b/P_c = 2$
) and NI-6 (
$P_b/P_c = 0.03$
). In this regard, figure 18 depicts the transient growth envelopes of these two cases. As can be observed, although the
${\textit{Br}}$
values are relatively small, the energy growth presents large rates when operating at transcritical conditions (NI-5) with only an approximate
$50\,\%$
reduction in maximum growth with respect to larger
${\textit{Br}}$
cases. Moreover, the optimum remains located at a similar wavelength region. Instead, if the system operates at low-pressure conditions, the growth rate is notably reduced.
Figure 18. Transient growth envelopes at
${\textit{Re}} = 1000$
and
${\textit{Br}} = 5.6 \times 10^{-6}$
for (a) NI-5 and (b) NI-6 cases.
The discussion below aims at detailing the stability maps obtained from the modal analyses. In particular, the focus is placed on the non-orthogonal effects of the operator. The neutral curves and the corresponding critical Reynolds numbers can be directly compared with the exponential energy growth. In this regard, figures 19 and 20 depict (i) the growth rate maps on the
$\alpha -Re$
space under 2-D perturbations (
$\beta = 0$
) at different contour levels, and (ii) the limit beyond which energy grows to infinity. Similar to the transient growth maps, equivalent behaviour is exhibited at low
${\textit{Br}}$
values as in the results from modal analysis. Consequently, colour maps for large
${\textit{Br}}$
cases are only depicted. As can be observed, since the plots present similar unstable regions, the results from energy growth validate the modal analysis results in the case of isothermal conditions. Furthermore, figure 4(c) confirms that transient growth also becomes stable at large
${\textit{Br}}$
values as a result of lower energy levels. The non-isothermal cases are also well captured, as well as the extended destabilisation region at low Reynolds numbers for case NI-1 at
$\alpha \approx 1$
. It is important to note that increasing the bulk pressure to
$P_b/P_c = 5$
results in a shift of the destabilisation region to larger streamwise wavenumbers, and consequently a narrower envelope is obtained with
$G_{max} \geqslant 200$
as depicted in figure 20(b). In addition, case NI-3 results in an even smaller region despite operating at transcritical conditions as shown in figure 20(c). It is worth noting that an additional large growth rate area is magnified at low streamwise wavenumbers. Although the energy remains stable over time, this narrower temperature difference set-up yields an important operating window of high energy growth at
$\alpha = 0.4$
. In particular, it confirms the trend observed in the energy growth rate map at
${\textit{Re}} = 1000$
(figure 16
c). In particular, this growth is within the same order of magnitude compared with the unstable wavenumber
$\alpha = 0.8$
at
${\textit{Re}} = 10\,000$
. Unlike
$\alpha = 0.8$
, the
$\alpha = 0.4$
region does not achieve an exponentially amplification growth, but the energy level remains similar. Moreover, figure 20(d) confirms that the low-velocity flow case NI-5 has a similar envelope to NI-1. Thus, this case is the most efficient candidate to enhance flow destabilisation. For brevity, cases NI-4 and NI-6 are not shown as they behave similarly to the supercritical flow case I-3 (figure 4
c).
Figure 19. Transient growth envelopes at
$\beta = 0$
and
${\textit{Br}} = 0.5$
for (a) I-1, (b) I-2 and (c) I-3 cases with
$10$
spaced contour levels. Grey-filled circles denote infinite energy growth. The resolution range for cases I-1 and I-2 has been limited to
$G_{\textit{max}} = 5000$
.
Figure 20. Transient growth envelopes at
$\beta = 0$
and
${\textit{Br}} = 0.1$
for (a) NI-1, (b) NI-2 and (c) NI-3 cases, and
${\textit{Br}} = 5.6 \times 10^{-6}$
for (d) NI-5. The resolution range has been limited to
$G_{\textit{max}} = 1000$
.
Finally, table 4 summarises the critical Reynolds number and maximum growth rate for the optimal
$\alpha -\beta$
parameter space found in the energy growth maps for the largest
${\textit{Br}}$
values studied. In particular,
$\alpha = 0$
and
$\beta =2$
result in a common region with a maximum growth rate for all flow cases, except for case I-2 (and I-1) in which a significantly large growth rate is also reported at
$\alpha = 1$
and
$\beta =1$
. Specifically, at higher
${\textit{Re}}$
this growth rises exponentially. These conditions result in early flow destabilisation for case NI-1–3 at
${\textit{Re}} \approx 2500$
, while case NI-5 at low Mach conditions achieves the earliest flow destabilisation at a similar regime. However, despite the energy increase, the behaviour is not exponential, and consequently the energy may eventually decay at larger times. Interestingly, the subcritical and transcritical cases I-1 and I-2 provide an unstable energy region only within the range
$900 \leqslant Re \leqslant 4000$
for I-2 and move toward large values
$1500 \leqslant Re \leqslant 5100$
for I-2. Nonetheless, I-1 presents rates that are roughly
$8\times$
smaller in comparison to the transcritical case. It is important to note that these amplification rates are achieved for non-isothermal cases at lower Reynolds numbers.
Table 4. Energy-based critical Reynolds numbers and energy growths for the flow cases listed in § 5.1 at two-parameter space
$(\alpha , \beta )$
levels. Cases I-1, I-2 and I-3 are assessed at
${\textit{Br}} = 0.5$
, NI-1, NI-2, NI-3 and NI-4 at
${\textit{Br}} = 0.1$
, and NI-5 and NI-6 at
${\textit{Br}} = 5.6 \times 10^{-6}$
. Transition criteria defined if energy grows beyond elapsed time (
$t \leqslant 400$
), highlighting with superscript
$(\boldsymbol{\cdot })^\star$
the exponentially algebraic growth.
5.3.2. Optimal input profile and output responses
Given that the growth rate over time describes the maximum energy amplification over all possible initial conditions, it is interesting to quantify which specific initial perturbation is responsible for the maximum amplification energy output at a given time. This initial condition is easily recovered by means of a SVD of the matrix exponential (Schmid Reference Schmid2007). In this regard, the optimal initial condition and the corresponding response are, respectively, shown in figures 21 and 22 at wavelengths
$\alpha = 0$
and
$\beta = 2$
providing the maximum amplification (infinitesimal perturbations scaled by the maximum perturbation value stated in the caption of the figures). On the one hand, for the isothermal wall cases, the streamwise vortices and velocity streaks are recovered for the optimal perturbation and its response, similarly to the case of incompressible flows. However, the thermal streaks are also important due to the crucial role of density and temperature at high-pressure transcritical conditions, whereas the hydrodynamic streaks remain of the same importance as the streamwise velocity. On the other hand, the non-isothermal cases yield a distinct behaviour. The optimal profile results in a parabolic wall-normal velocity perturbation. The response is prominent for thermal streaks located within the pseudo-boiling region (near the hot wall) and dominated by density over temperature variations. Moreover, for case NI-5, similar optimal input and output responses are exhibited. Instead, when operating at low-pressure conditions (case NI-6), the velocity streaks are again enhanced and more prominent in the vicinity of the cold wall. In particular, thermal streaks are significant within this region and more prominent than the isothermal wall subcritical and supercritical flow set-ups. It is important to note that the outstanding energy growth obtained for case I-2 at
$\alpha = \beta = 1$
is quantified in figure 23. Thus, the optimal input results in 3-D-like vortices with similar streamwise and spanwise components. Nonetheless, the wall-normal velocity is relatively lower for this case, and the response is enlarged for the density streaks whereas the streamwise velocity vortex becomes asymmetric and gradually decreases its magnitude toward the channel centreline.
Figure 21. Optimum input perturbation profiles at
${\textit{Re}} = 1000$
,
$\alpha = 0$
and
$\beta = 2$
for the isothermal flow cases (a) I-1, I-2 and I-3 at
${\textit{Br}} = 0.5$
, (b) NI-1, NI-2 and NI-3 at
${\textit{Br}} = 0.1$
, and (c) NI-4 and NI-5 at
${\textit{Br}} = 5.6 \times 10^{-6}$
. Results are normalised by
$w^\prime$
for all cases.
Figure 22. Optimum output response at
${\textit{Re}} = 1000$
,
$\alpha = 0$
and
$\beta = 2$
for the isothermal flow cases (a) I-1, I-2 and I-3 at
${\textit{Br}} = 0.5$
, (b) NI-1, NI-2 and NI-3 at
${\textit{Br}} = 0.1$
, and (c) NI-5 and NI-6 at
${\textit{Br}} = 5.6 \times 10^{-6}$
. Results are normalised by
$u^\prime$
for (a) and
$\rho ^\prime$
for (b,c).
Figure 23. Optimum (a) perturbation and (b) response at
${\textit{Re}} = 1000$
,
$\alpha = \beta = 1$
for the isothermal flow case I-2 at
${\textit{Br}} = 0.5$
. Results are normalised by
$w^\prime$
for (a) and
$\rho ^\prime$
for (b).
Figure 24 displays the optimal velocity patterns across the
$z-y$
plane at wavenumbers
$\alpha = 0$
and
$\beta = 2$
; in this case, invariant along the streamwise direction. The pattern resembles the streamwise vortices typical of shear flows (Malik et al. Reference Malik, Alam and Dey2006). These counter-rotating vortices tend to be compacted within the half-plane region at non-isothermal conditions for both Brinkman numbers. Nonetheless, at low velocity, the vortices are slightly biased to the hot wall. Likewise, the optimal responses recover the spanwise alternating streamwise velocities. The non-isothermal flow cases, however, report an interesting result. At large
${\textit{Br}}$
values for case NI-1, there is a single alternating streamwise velocity at wall-normal position
$0.5 \leqslant y/\delta \leqslant 1.5$
, presenting an elongated ellipsoid-like shape. The temperature field presents a similar behaviour, but its maximum is shifted to
$y/\delta \approx 1.6$
, which corresponds to the pseudo-boiling region. Interestingly, the density field presents a well-defined ellipsoid at
$y/\delta \approx 1.3$
. Furthermore, at lower Reynolds and Mach numbers (case NI-5), the counter-rotating vortices are clearly apparent and interact with the pseudo-boiling region that splits their rotation in two. Instead, the thermodynamic fields behave similarly to case NI-1, but their shapes are moderately stretched to the cold wall. It is important to highlight that DNS results of an equivalent 3-D turbulent flow set-up (Bernades, Capuano & Jofre Reference Bernades, Capuano and Jofre2022) indicated a similar behaviour of flow structures.
Figure 24. Optimum velocity perturbation at
${\textit{Re}} = 1000$
,
$\alpha = 0$
and
$\beta = 2$
for cases (a) I-1 at
${\textit{Br}} = 0.5$
, (b) NI-1 at
${\textit{Br}} = 0.1$
, and (c) NI-5 at
${\textit{Br}} = 5.6 \times 10^{-6}$
.
Figure 25. Optimum response at
${\textit{Re}} = 1000$
,
$\alpha = 0$
and
$\beta = 2$
for (a,b,c) I-1 at
${\textit{Br}} = 0.5$
, (d,e,f) NI-1 at
${\textit{Br}} = 0.1$
, and (g,h,i) NI-5 at
${\textit{Br}} = 5.6 \times 10^{-6}$
.
5.4. The DNS-based early evolution of modal and non-modal perturbations
Finally, with the aim of partially validating the linear stability predictions by solving the fully nonlinear equations, perturbations obtained from LST are superimposed onto a steady base flow field and evolved in time using DNS by means of the in-house flow solver RHEA (Jofre et al. Reference Jofre, Abdellatif and Oyarzun2023a ), described in § 4. In particular, the perturbations obtained from NI-5, viz. a similar set-up to the DNS study by Bernades et al. (Reference Bernades, Capuano and Jofre2023a ), corresponding to linearly unstable eigenmodes are superimposed to the steady solution to examine how the disturbances evolve nonlinearly in the early stages of flow destabilisation (Biau, Soueid & Bottaro Reference Biau, Soueid and Bottaro2008; Buffat et al. Reference Buffat, Penven, Cadiou and Montagnier2014; Kashyap, Duguet & Dauchot Reference Kashyap, Duguet and Dauchot2024). The initial perturbation field is constructed from the eigenfunctions of the least stable LST modes as
\begin{equation} \boldsymbol{q}^\prime (x,y,z,t) = \epsilon \thinspace {\textit{Re}} \thinspace \left[ \hat {\boldsymbol{q}}(y) e^{(i \alpha x + i \beta z)}\right]\!, \end{equation}
where a small amplitude
$\epsilon = 10^{-3}$
is used, consistent with the approach of Yamamoto et al. (Reference Yamamoto, Takahashi and Kambe2001). The perturbation is added to the base flow at
$t=0$
, and the resulting velocity field is evolved under the full equations of fluid motion.
In particular, two cases are analysed: (i) modal analysis with exponentially growing modes; and (ii) transient growth from algebraically amplifying perturbations, drawn from the most unstable modes obtained in §§ 5.2 and 5.3, respectively. Therefore, the eigenvectors corresponding to the least stable eigenvalues within their corresponding wavenumber parameter space are projected to the physical space, denormalised and imposed onto a steady flow field based on
$\boldsymbol{q}^\prime (x,y,z,0)$
. To monitor the evolution of disturbances, the perturbation energy disturbance has been tracked and computed for each time step as
\begin{equation} E(t) = \frac {1}{2} \thinspace \int _V (u^\prime + v^\prime + w^\prime ) {\rm d}V, \end{equation}
where the perturbations were obtained subtracting the base flow from the instantaneous fields. Hence, analogous to the decomposition of (3.1), they read
\begin{equation} \left . \begin{array}{ll} \displaystyle {u^\prime (x,y,z,t)} = u (x,y,z,t) - U(y),\\[8pt] \displaystyle {v^\prime (x,y,z,t)} = v (x,y,z,t),\\[8pt] \displaystyle {w^\prime (x,y,z,t)} = w (x,y,z,t). \end{array}\right \} \end{equation}
The energy growth
$E(t)/E(0)$
is compared against the growth predicted by LST, which for modal stability, corresponds to
$E(t) = e^{2 \omega _i t}$
, where
$\omega _i$
represents the imaginary part of the most unstable mode.
5.4.1. Modal stability perturbation
The modal linear stability growth rate is computed for the most unstable mode (
$\omega _i = 4.05 \times 10^{-3}$
) for a wavenumber pair
$(\alpha , \beta ) = (2, 0)$
based on the perturbation profile described in figure 10(a), which represents similar base flow conditions as the DNS by Bernades et al. (Reference Bernades, Capuano and Jofre2023a
) at
${\textit{Re}}_b \approx 4000$
. However, DNS results reveal a significantly faster disturbance amplification (not shown here) compared with LST. This behaviour suggests a dominant non-modal transient amplification mechanism. Indeed, DNS captures several effects beyond linear theory, including nonlinear interactions, finite-amplitude effects, mode coupling and secondary instabilities that are obviously absent in the LST framework. This reinforces the points raised in the modal analysis discussion, where caution is warranted due to the non-normality of the linear operator. Consequently, optimal input perturbations are evaluated using algebraic growth analysis in the following subsection to better characterise the early stage amplification evolution.
5.4.2. Non-modal stability optimal perturbation
A similar approach is employed to assess the predictions from transient growth theory. While DNS resolves all relevant nonlinear dynamics and multi-scale interactions, the linear transient growth framework captures energy amplification arising purely from the non-normality of the linearised equations of fluid motion operator. In this context, the transient growth mechanism results from the constructive interference of non-orthogonal eigenmodes, even in the absence of any unstable eigenvalues. This makes it a natural benchmark for early time flow evolution in scenarios dominated by non-modal mechanisms (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993; Jeong et al. Reference Jeong, Hussain, Schoppa and Kim1997; del Alamo & Jiménez Reference del Alamo and Jiménez2006; Schmid Reference Schmid2007).
Optimal perturbations (previously computed from the linear analysis and shown in figure 21(c) for
${\textit{Re}}_b = 1000$
) are introduced into the base flow as initial disturbances. These profiles remain qualitatively similar for
${\textit{Re}}_b = 4000$
, and to avoid redundancy, the reader is referred to the previous figures for representative shapes. The perturbations are imposed along the wall-normal and spanwise directions, corresponding to the streamwise-invariant structures (
$\alpha = 0$
), in accordance with (5.1). In this configuration, the disturbance energy growth observed in DNS agrees much more closely with the transient growth prediction than with the modal LST prediction, indicating that non-modal amplification mechanisms dominate the initial destabilisation of the flow. Figure 26 compares the normalised disturbance energy from DNS with the theoretical optimal growth
$G_{LST}(t)$
predicted by LST. At early times (
$t^\star \lt 40$
), the DNS evolution closely matches the linear prediction, validating the relevance of the optimal perturbation framework. However, as time progresses and perturbation amplitudes increase, nonlinear effects become significant, leading to divergence from the linear prediction. In the present case, figure 27 shows the evolution of vertical vorticity in the
$y{-}z$
plane (streamwise-invariant perturbation), where enhanced enstrophy levels appear just above the pseudo-boiling layer. This region coincides with strong baroclinic torque activity, suggesting a link between thermodynamic gradients and the amplification of vorticity structures. These observations are consistent with prior physics-based analyses of high-pressure transcritical channel flows (Bernades et al. Reference Bernades, Capuano and Jofre2023a
), and warrant further investigation of the transitional process using DNS.
Figure 26. (a) Disturbed energy over time based on initial optimal input perturbation from transient growth results at
${\textit{Re}}_b = 4000$
,
$\alpha = 0$
and
$\beta = 2$
comparing DNS with the non-modal energy growth, and (b) wall-normal ensemble average of normalised DNS disturbed energy.
Figure 27. Wall-normal vorticity on a
$y{-}z$
plane indicating the location of the pseudo-boiling line (purple dashed-doted line) and cross-velocity vector fields for (a)
$t^\star \approx 0$
, (b)
$t^\star \approx 1$
and (c)
$t^\star \approx 40$
.
In this regard, a comprehensive investigation of the full laminar-to-turbulent transition process, including detailed characterisation of coherent structures and transition pathways, is beyond the scope of the present study. Nonetheless, the present findings reinforce the dominant role of transient growth in the early dynamics of differentially heated high-pressure transcritical wall-bounded flows. Hence, the streak breakdown, the secondary instability effects and the emergence of vortical-like structures culminating in turbulence saturation when nonlinear effects become important are yet to be further explored.
6. Conclusions
A real-gas framework for linear stability analysis of plane Poiseuille flows has been developed and verified at the isothermal limit to study the stability of high-pressure transcritical flows. The subsequent modal analysis has focused on characterising the effects of streamwise perturbations for isothermal and non-isothermal flow set-ups at various Brinkman numbers. On the one hand, it has been found that for isothermal wall cases, the destabilisation occurs at lower Reynolds numbers when the Brinkman number increases. In particular, at transcritical conditions, the modal-based destabilisation is triggered at smaller Reynolds numbers than at subcritical thermodynamic regimes. Instead, similarly to ideal gas cases, the supercritical cases result in stability enhancement when the Brinkman number increases. On the other hand, non-isothermal set-ups have been studied by imposing temperature differences between walls such that the fluid undergoes a transcritical trajectory across the pseudo-boiling region. For these cases, it has been found that destabilisation effects arise at relatively lower Brinkman numbers. Moreover, componentwise energy budget analyses have revealed that their redistribution through pressure-velocity interactions plays a crucial role in driving energy growth/decay. In particular, kinetic energy balances have indicated that the production term is responsible for the destabilisation trends observed under isothermal subcritical conditions, whereas the transcritical cases are limited only to compressibility effects. Instead, the non-isothermal cases at low velocities are dominated by both the production and thermodynamic terms, especially the latter which mainly drives flow destabilisation. From a non-modal perspective, similar patterns are observed when treating unforced algebraic growth compared with modal analysis, where exponential energy growths are captured at similar
$\alpha -Re$
parameter values. This confirms that amplification rates are significantly large for the non-isothermal cases, even for low-Brinkman-number base flows. In addition, the algebraic growth of the 3-D perturbations provides a maximum amplification within a region similar to incompressible flows. Nonetheless, this region is significantly enlarged in the case of non-isothermal flows with amplification rates of roughly
$10\times$
in comparison to isothermal wall transcritical and equivalent sub-/supercritical non-isothermal set-ups. To this end, the non-isothermal (i.e. differentially heated) configuration demonstrates clear advantages in promoting earlier flow destabilisation, as evidenced from both modal and algebraic growth perspectives. Specifically, (i) the underlying destabilisation trends differ from those observed in isothermal cases, and (ii) the required input energy is significantly lower. Furthermore, both modal and non-modal optimal perturbations have been applied to a laminar channel flow DNS solution, revealing distinctly nonlinear energy dynamics even at early time instants. In particular, the disturbed energy rapidly diverges from linear predictions due to the onset of secondary instabilities. Nevertheless, there is evidence that flow destabilisation may inherently initiate just above the pseudo-boiling line, where preliminary streamwise structures form and enhance fluid rotation.
The modal stability analysis has also highlighted that, for the isothermal transcritical cases, the enhancement of destabilisation effects results in a critical Reynolds number that is
$5\times$
smaller than that of the isothermal limit case. Instead, at supercritical conditions at the same Brinkman number, the neutral curve is positively shifted by a factor of
$2\times$
with respect to the transcritical case. In particular, density and temperature perturbations are mostly responsible for such destabilisation for both regimes. However, supercritical flows enhance the stability when the Brinkman number increases. Alternatively, the non-isothermal set-ups exhibit flow destabilisation at even lower Brinkman numbers for wavenumbers that are roughly
$20\,\,\%$
smaller than for the isothermal cases. It is important to note that, by increasing the pressure of the system by
$5\times$
, similar neutral curves as for the isothermal subcritical case are recovered, although the critical Reynolds number remains slightly lower. Nonetheless, if the asymmetric wall temperature difference is reduced, destabilisation is significantly delayed and comparable to the isothermal subcritical cases. Moreover, the perturbations are dominated by thermodynamic modes, especially near the pseudo-boiling region. The additional stability provided by large Brinkman numbers has been shown to be due to a combination of decreased energy production and a more active velocity-pressure gradient term, which forces a component redistribution of energy in a way that significantly dampens wall-normal fluctuations. For isothermal wall flows, the velocity production and second derivative of the vorticity transport equation dominate, whereas the former is more important for the non-isothermal flows except near the pseudo-boiling region where the baroclinic effect becomes dominant and independent of the Brinkman number. Furthermore, the results of algebraic growth have been found to perform similar to those of the modal analysis. Additionally, under 3-D perturbations, it has been found that at
${\textit{Re}} = 1000$
a coexisting destabilisation region is obtained where amplifications tend to infinity. This region corresponds to
$\alpha = 0$
and
$\beta \approx 2$
, which for the non-isothermal transcritical flow cases, expands to a wider spanwise region at
$1 \leqslant \beta \leqslant 3$
. Interestingly, the isothermal wall transcritical case presents a coexisting region at
$\alpha = \beta \approx 1$
where amplifications grow exponentially. Moreover, it is noted that the maximum amplifications are exacerbated for the non-isothermal cases by a factor of approximately
$10\times$
in comparison to the isothermal wall transcritical cases. Nonetheless, this energy enhancement is lost when operating at low-pressure conditions or far beyond the critical point. Finally, fully nonlinear DNS computations of the amplification of both modal and non-modal perturbations were performed. Results suggest that the early evolution of disturbances in a high-pressure transcritical differentially heated channel flow is governed predominantly by non-modal mechanisms. When initialised with optimal transient growth perturbations, the DNS energy evolution closely follows the linear predictions, validating the relevance of non-modal analysis in capturing early time flow dynamics. In contrast, when modal LST eigenmodes are used as initial disturbances, the observed energy growth rate in DNS significantly exceeds the corresponding linear modal prediction. This discrepancy underscores the strong non-normality of the linearised operator, indicating that even linearly unstable modes can exhibit transient amplification beyond their asymptotic exponential growth.
Future work will extend these analyses to a broader range of non-isothermal configurations. In particular, energy budgets and growth rates will be further characterised by comparing isothermal and non-isothermal set-ups under 3-D perturbations. Additionally, the unstable modes arising from 2-D perturbations will be examined to identify the presence of additional instabilities beyond those covered in the current transcritical cases. While the present study provides valuable insight into the early destabilisation mechanisms of high-pressure transcritical channel flows using LST and DNS, a detailed investigation of the full laminar-to-turbulent transition process remains an important direction for future research. In particular, the emergence, evolution and breakdown of coherent structures (such as streaks, vortices and hairpin formations) warrant further analysis using high-fidelity DNS. To build upon the current findings, future work will focus on exploring optimal perturbations not only for energy amplification but also for practical objectives such as skin-friction reduction and transition advance/delay. The nonlinear evolution of base flows perturbed with both modal and non-modal disturbances will be systematically examined across a broader wavenumber spectrum to map the full landscape of transient growth behaviour. This will aid in identifying the critical modes and amplification pathways responsible for triggering transition under complex thermodynamic conditions. Ultimately, such efforts aim to uncover the underlying physical mechanisms driving transition in high-pressure, differentially heated flows and to bridge the gap between linear theory and fully developed turbulence. This progression from linear stability characterisation to nonlinear transition dynamics represents a natural and necessary extension of the present study. Finally, resolvent analyses will be conducted to assess the sensitivity of the operator to input forcing as well as to explore data-driven forcings and their corresponding flow responses.
Acknowledgements
The authors acknowledge support from the Formació de Professorat Universitari scholarship (FPU-UPC R.D 103/2019), and the Serra Húnter and SGR (2021-SGR-01045) programs of the Generalitat de Catalunya (Catalonia).
Funding
This work is funded by the European Union (ERC, SCRAMBLE, 101040379). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.
Declaration of interests
The authors report no conflict of interest.
Author contributions
Marc Bernades: conceptualisation, formal analysis, investigation, software, writing – original draft. Francesco Capuano: conceptualisation, investigation, supervision, writing – review & editing. Lluís Jofre: conceptualisation, funding acquisition, investigation, supervision, writing – review & editing.
Data availability statement
The data reported in this paper has been generated by means of the in-house MATLAB code named high-pressure linear stability analysis (HPLSA). The source code is openly accessible at https://github.com/marc-bernades/HPLSA. HPLSA has been designed to serve as a flexible tool to develop linear stability methods for ideal gas and real-gas thermodynamic frameworks, including wrappers for CoolProp and RefProp libraries. The code is equipped with several comments for readability. Additionally, the repository includes a description of the code, instructions and a guide for users. Consequently, the linear stability modal and non-modal analysis can be fully reproduced by the interested reader. Although the code is suitable for Poiseuille flows, it can be easily adapted to other wall-bounded cases, such as Couette flow, by adjusting the initial and boundary conditions. Furthermore, the operator is built for temporal eigenproblems (prescribing streamwise and spanwise wavenumbers and solving the eigenproblem for the angular frequency and growth rate of the perturbation), but it is prepared to be expanded to solve spatial problems also, which is typical in the case of external flows. This solver requires the previous installation of the high-pressure compressible flow solver (HPCFS) available at https://github.com/marc-bernades/HPCFS, which embeds some functionalities and thermodynamic models necessary for the complete usage of HPLSA.
Appendix A. Thermophysical properties validation
The thermodynamic quantities and transport coefficients obtained from the CoolProp library (Bell et al. Reference Bell, Wronski, Quoilin and Lemort2014) are validated against NIST data. They are also compared against the Peng-Robinson equation of state coupled with the high-pressure transport coefficients from Chung et al. (Reference Chung, Lee and Starling1984, Reference Chung, Ajlan, Lee and Starling1988) (labelled as `Model’). In this regard, figure 28 depicts the density, viscosity, isobaric heat capacity and thermal conductivity for nitrogen at
$P / P_c = 2$
and temperature range
$T / T_c = 0.75 {-} 1.5 $
. It can be observed that the CoolProp results match with the NIST reference data. However, the `Model’ solution results in differences in the vicinity of the pseudo-boiling region and at subcritical temperatures at a liquid-like state. The accuracy of such a model has been extensively analysed by Bernades & Jofre (Reference Bernades and Jofre2022) at different bulk pressures for different substances. Additionally, the results from the NIST-based library RefProp have also been validated obtaining a similar performance to the open-source CoolProp library (results not shown to facilitate visualisation).
Figure 28. Thermodynamic and transport properties of nitrogen at
$P/ P_c = 2$
with a temperature range from
$T / T_c = 0.75$
to
$T / T_c = 1.5$
for NIST, CoolProp and Model (Peng-Robinson equation of state with Chung et al. (Reference Chung, Lee and Starling1984, Reference Chung, Ajlan, Lee and Starling1988) high-pressure coefficients) for density
$\rho$
, dynamic viscosity
$\mu $
, isobaric heat capacity
$c_p$
and thermal conductivity
$\kappa$
.
Appendix B. Base flow equations of fluid motion
The base flow is driven by a body force in the streamwise direction and is obtained by solving the equations of fluid motion assuming that the flow is fully developed in a laminar state, spanwise and streamwise independent, steady and parallel, i.e.
$\partial (\boldsymbol{\cdot })/\partial x = \partial (\boldsymbol{\cdot })/\partial z = 0$
,
$\partial (\boldsymbol{\cdot })/\partial t = 0$
,
$ v = w = 0$
. The dimensionless equations of fluid motion ((2.1)–(2.3)) are thus simplified to
\begin{align} \overbrace {\dfrac {{\rm d} \mu }{{\rm d} y} \dfrac {{\rm d} u}{{\rm d} y} + \mu \dfrac {{\rm d}^2 u}{{\rm d} y^2} }^{\dfrac {\rm d}{{\rm d} y} \left ( \mu \dfrac {{\rm d} u}{{\rm d} y} \right )} & = - F, \end{align}
\begin{align} \frac {{\rm d} P}{{\rm d} y} & = 0, \end{align}
\begin{align} \underbrace {u \left ( \dfrac {{\rm d} \mu }{{\rm d} y} \dfrac {{\rm d} u}{{\rm d} y} + \mu \dfrac {{\rm d}^2 u}{{\rm d} y^2} \right ) + \mu \left ( \dfrac {{\rm d} u}{{\rm d} y} \right )^2 + \dfrac {\rm d \kappa }{{\rm d} y} \dfrac {{\rm d} T}{{\rm d} y} + \kappa \dfrac {{\rm d}^2 T}{{\rm d} y^2}}_{\dfrac {\rm d}{{\rm d} y} \left ( u \mu \dfrac {{\rm d} u}{{\rm d} y} \right ) + \dfrac {\rm d}{{\rm d} y} \left ( \kappa \dfrac {{\rm d} T}{{\rm d} y} \right )} & = - \textit{Fu}, \end{align}
where density is obtained from the equation of state. In this regard, the streamwise velocity is calculated from (B1) and the temperature from (B3). At this stage, density, viscosity and thermal conductivity can be determined with the updated temperature. Based on the initial conditions, the body force
$F$
can be calculated from (B1) and controls the velocity profile to achieve the desired Reynolds number.
B.1. Base flow verification
The base flow of the linear stability solver is verified against the solution of a laminar steady-state channel flow. This case operates with nitrogen at low-pressure conditions
$P/P_c = 0.03$
and the set-up is such that it matches the same input pressure as the high-pressure transcritical equivalent case at
$P/P_c = 2$
(Bernades et al. Reference Bernades, Capuano and Jofre2023a
). The resulting dimensionless numbers and set-up for the iterative LST base flow are imposed similar to the DNS solution for
${\textit{Br}} = 7.7 \times 10^{-4}$
, where the reference velocity scaling is based on the bulk velocity. For consistency, instead of imposing a constant normalised body force, the LST incorporates a feedback control that adjusts the value of
$\hat {F}$
. Moreover, for this verification, the LST uses the same thermodynamic model used for the DNS computation, i.e. the Peng-Robinson equation of state and high-pressure transport coefficients from Chung et al. (Reference Chung, Lee and Starling1984, Reference Chung, Ajlan, Lee and Starling1988). The base flow iterative process for the linear stability analysis is independent of the Reynolds number. In this regard, figure 29 depicts the velocity, temperature and transport properties compared between the LST and DNS frameworks. As can be observed, good agreement between both approaches is found, which consequently verifies the LST algorithm.
Figure 29. Low-pressure
$P/P_c = 0.03$
non-isothermal base flow with a temperature range from
$T / T_c = 0.75$
to
$T / T_c = 1.5$
utilising nitrogen from (i) ensemble-averaged DNS high-pressure transcritical channel flow, and (ii) linear stability solver for velocity, temperature, dynamic viscosity and thermal conductivity. Both frameworks utilise the thermodynamic model of Peng-Robinson equation of state and Chung et al. (Reference Chung, Lee and Starling1984, Reference Chung, Ajlan, Lee and Starling1988) high-pressure transport coefficients.
B.2. Grid convergence
A grid convergence study is briefly summarised below with results for modal maximum growth rate and non-modal transient growth. The results shown in figure 30 are calculated based on the
$L_1$
-norm with respect to the finest grid assessed (
$N = 450$
) for a channel height of
$L/\delta = 2$
. As can be observed, (i) the spatial error is second-order (note that
$L/N$
is varying across the wall-normal direction based on the Chebyshev collocation), and (ii) the value
$N = 200$
provides an optimal balance between accuracy and computational cost for the eigenproblems and SVDs conducted within this work.
Figure 30. Grid convergence results for case I-1 of maximum growth rate modal analysis at
$\alpha = 1.0$
,
${\textit{Br}} = 0.1$
and
${\textit{Re}} = 10\,000$
(a), and maximum transient growth non-modal analysis at
$\alpha = 0.0$
,
$\beta = 2.0$
,
${\textit{Br}} = 0.1$
and
${\textit{Re}} = 1000$
(b).
Appendix C. Linear stability equations
This section presents the linearised equations of fluid motion. They are derived by substituting (3.1) into the dimensionless conservation ((2.1)–(2.3)) with appropriate algebra developments. The derived discrete equations can then be rewritten in matrix form as formulated in (3.3). Therefore, by inspection, the corresponding matrices can be obtained for each of the
$5\times 5$
terms of
$\boldsymbol{L_t}$
,
$\boldsymbol{L_x}$
,
$\boldsymbol{L_y}$
,
$\boldsymbol{L_z}$
,
$\boldsymbol{L_q}$
,
$\boldsymbol{V_{xx}}$
,
$\boldsymbol{V_{yy}}$
,
$\boldsymbol{V_{zz}}$
,
$\boldsymbol{V_{xy}}$
,
$\boldsymbol{V_{xz}}$
and
$\boldsymbol{V_{yz}}$
.
C.1. Continuity equation
Substituting the decomposition
$\rho = \rho _0 + \rho ^\prime$
into (2.1) reads
\begin{align} & \frac {\partial (\rho _0 + \rho ^\prime )}{\partial t} = - \boldsymbol{\nabla }\boldsymbol{\cdot }[ (\rho _0 + \rho ^\prime ) (\boldsymbol{u_0} + \boldsymbol{u}^\prime ) ] = \nonumber \\ & -\rho _0 \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^\prime - \underbrace {\rho _0 \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u_0}}_{=0} - \rho ^\prime \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^\prime - \underbrace {\rho ^\prime \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u_0}}_{=0} - \underbrace {\boldsymbol{u_0} \boldsymbol{\nabla }\rho _0}_{=0} - \underbrace {\boldsymbol{u_0} \boldsymbol{\nabla }\rho ^\prime }_{u_0 \dfrac {\rm d \rho ^\prime }{{\rm d}x}} - \underbrace {\boldsymbol{u}^\prime \boldsymbol{\nabla }\rho _0}_{v^\prime \dfrac {\rm d \rho _0}{\rm {\rm d}y}} - \underbrace {\boldsymbol{u}^\prime \boldsymbol{\nabla }\rho ^\prime }_{= 0}, \end{align}
which, by further manipulation, results in
\begin{align} \frac {\partial \rho ^\prime }{\partial t} = - \rho _0 \left ( \frac {{\rm d}u^\prime }{{\rm d}x} + \frac {{\rm d}v^\prime }{{\rm d}y} + \frac {{\rm d}w^\prime }{{\rm d}z}\right ) - u_0 \frac {{\rm d} \rho ^\prime }{{\rm d}x} - v^\prime \frac {{\rm d} \rho _0}{{\rm d}y}. \end{align}
C.2. Momentum equations
Analogously to the continuity equation, substituting the decomposition of variables into (2.2) yields
\begin{align} & \frac {\partial [ (\rho _0 + \rho ^\prime ) (\boldsymbol{u_0} + \boldsymbol{u}^\prime ) ]}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }[ (\rho _0 + \rho ^\prime ) (\boldsymbol{u_0} + \boldsymbol{u}^\prime ) (\boldsymbol{u_0} + \boldsymbol{u}^\prime )] + \boldsymbol{\nabla }((P_0 + P^\prime )) = \nonumber \\ & \frac {1}{\textit{Re}} \boldsymbol{\nabla }\boldsymbol{\cdot }\{ (\mu _0 + \mu ^\prime ) [\boldsymbol{\nabla }(\boldsymbol{u_0} + \boldsymbol{u}^\prime ) + \boldsymbol{\nabla }(\boldsymbol{u_0} + \boldsymbol{u}^\prime )^{T} ] + (\lambda _0 + \lambda ^\prime )[\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u_0} + \boldsymbol{u}^\prime )]\boldsymbol{I} \} + \boldsymbol{F}. \end{align}
Based on this expression, the linearised momentum in each direction is detailed as follows.
Momentum streamwise direction:
\begin{align} & \rho _0 \frac {\partial u^\prime }{\partial t} + u_0 \frac {\partial \rho ^\prime }{\partial t} + u_0 u_0 \frac {\partial \rho ^\prime }{\partial x} + \frac {\partial P_0}{\partial \rho _0} \frac {\partial \rho ^\prime }{\partial x} + 2 \rho _0 u_0 \frac {\partial u^\prime }{\partial x} + \rho _0 u_0 \frac {\partial v^\prime }{\partial y} \nonumber \\ &\quad + \rho _0 u_0 \frac {\partial w^\prime }{\partial z} + \rho _0 \frac {\partial u_0}{\partial y} v^\prime + u_0 \frac {\partial \rho _0}{\partial y} v^\prime = \frac {1}{\textit{Re}} \frac {\partial \mu _0}{\partial y}\frac {\partial v^\prime }{\partial x} + \frac {1}{\textit{Re}} \frac {\partial \mu _0}{\partial y}\frac {\partial u^\prime }{\partial y} \nonumber \\ &\quad + \frac {2 \mu _0 - \lambda _0}{\textit{Re}} \frac {\partial ^2 u^\prime }{\partial ^2 x} + \frac {\mu _0}{\textit{Re}} \frac {\partial ^2 u^\prime }{\partial ^2 y} + \frac {\mu _0}{\textit{Re}} \frac {\partial ^2 u^\prime }{\partial ^2 z} + \frac {\mu _0 + \lambda _0}{\textit{Re}} \frac {\partial ^2 v^\prime }{\partial x \partial y} + \frac {\mu _0 + \lambda _0}{\textit{Re}} \frac {\partial ^2 w^\prime }{\partial x \partial z} \nonumber \\ &\quad + \frac {\partial P_0}{\partial T_0}\frac {\partial T^\prime }{\partial x} + \frac {1}{\textit{Re}} \frac {\partial \mu _0}{\partial \rho _0}\frac {\partial u_0}{\partial y}\frac {\partial \rho ^\prime }{\partial y} + \frac {1}{\textit{Re}} \frac {\partial \mu _0}{\partial T_0}\frac {\partial u_0}{\partial y}\frac {\partial T^\prime }{\partial y} + \frac {1}{\textit{Re}} \frac {\partial \mu _0}{\partial \rho _0}\frac {\partial ^2 u_0}{\partial ^2 y} \rho ^\prime \nonumber \\ &\quad + \frac {1}{\textit{Re}} \frac {\partial \mu _0}{\partial T_0}\frac {\partial ^2 u_0}{\partial ^2 y} T^\prime + \frac {1}{\textit{Re}} \left ( \frac {\partial ^2 \mu _0}{\partial ^2 \rho _0} \frac {\partial \rho _0}{\partial y} + \frac {\partial ^2 \mu _0}{\partial \rho _0 \partial T_0} \frac {\partial T_0}{\partial y} \right ) \frac {\partial u_0}{\partial y} \rho ^\prime \nonumber \\ &\quad + \frac {1}{\textit{Re}} \left ( \frac {\partial ^2 \mu _0}{\partial ^2 T_0} \frac {\partial T_0}{\partial y} + \frac {\partial ^2 \mu _0}{\partial T_0 \partial \rho _0} \frac {\partial \rho _0}{\partial y} \right ) \frac {\partial u_0}{\partial y} T^\prime + \frac {\hat {F}}{\textit{Re}}. \end{align}
Momentum wall-normal direction:
\begin{align} & \rho _0 \frac {\partial v^\prime }{\partial t} + \rho _0 u_0 \frac {\partial v^\prime }{\partial x} = \frac {1}{\textit{Re}} \frac {\partial \lambda _0}{\partial y}\frac {\partial u^\prime }{\partial x} + \frac {1}{\textit{Re}} \frac {\partial \lambda _0}{\partial y}\frac {\partial v^\prime }{\partial y} + \frac {2}{\textit{Re}} \frac {\partial \mu _0}{\partial y}\frac {\partial v^\prime }{\partial y} + \frac {1}{\textit{Re}} \frac {\partial \lambda _0}{\partial y}\frac {\partial w^\prime }{\partial z} \nonumber \\ &\quad + \frac {\mu _0}{\textit{Re}} \frac {\partial ^2 v^\prime }{\partial ^2 x} + \frac {2 \mu _0 + \lambda _0}{\textit{Re}} \frac {\partial ^2 v^\prime }{\partial ^2 y} + \frac {\mu _0}{\textit{Re}} \frac {\partial ^2 v^\prime }{\partial ^2 z} + \frac {\mu _0 + \lambda _0}{\textit{Re}} \frac {\partial ^2 u^\prime }{\partial x \partial y} + \frac {\mu _0 + \lambda _0}{\textit{Re}} \frac {\partial ^2 w^\prime }{\partial y \partial z} \nonumber \\ &\quad + \frac {1}{\textit{Re}} \frac {\partial \mu _0}{\partial \rho _0}\frac {\partial u_0}{\partial y}\frac {\partial \rho ^\prime }{\partial x} + \frac {1}{\textit{Re}} \frac {\partial \mu _0}{\partial T_0}\frac {\partial u_0}{\partial y}\frac {\partial T^\prime }{\partial x} + \frac {\partial P_0}{\partial \rho _0}\frac {\partial \rho ^\prime }{\partial y} + \frac {\partial P_0}{\partial T_0}\frac {\partial T^\prime }{\partial y} + \frac {\partial ^2 P_0}{\partial ^2 \rho _0} \frac {\partial \rho _0}{\partial y} \rho ^\prime \nonumber \\ &\quad + \frac {\partial ^2 P_0}{\partial T_0 \partial \rho _0} \frac {\partial T_0}{\partial y} \rho ^\prime + \frac {\partial ^2 P_0}{\partial ^2 T_0} \frac {\partial T_0}{\partial y} T^\prime + \frac {\partial ^2 P_0}{\partial \rho _0 \partial T_0} \frac {\partial \rho _0}{\partial y} T^\prime . \end{align}
Momentum spanwise direction:
\begin{align} &\rho _0 \frac {\partial w^\prime }{\partial t} + \rho _0 u_0 \frac {\partial w^\prime }{\partial x} = \frac {1}{\textit{Re}} \frac {\partial \mu _0}{\partial y}\frac {\partial w^\prime }{\partial y} + \frac {1}{\textit{Re}} \frac {\partial \mu _0}{\partial y}\frac {\partial v^\prime }{\partial z} + \frac {\mu _0}{\textit{Re}} \frac {\partial ^2 w^\prime }{\partial ^2 x} + \frac {\mu _0}{\textit{Re}} \frac {\partial ^2 w^\prime }{\partial ^2 y} \nonumber \\ &\quad + \frac {2 \mu _0 + \lambda _0}{\textit{Re}} \frac {\partial ^2 w^\prime }{\partial ^2 z} + \frac {\mu _0 + \lambda _0}{\textit{Re}} \frac {\partial ^2 u^\prime }{\partial x \partial z} + \frac {\mu _0 + \lambda _0}{\textit{Re}} \frac {\partial ^2 v^\prime }{\partial y \partial z} + \frac {\partial P_0}{\partial \rho _0}\frac {\partial \rho ^\prime }{\partial z} + \frac {\partial P_0}{\partial T_0}\frac {\partial T^\prime }{\partial z}. \end{align}
C.3. Internal energy equation
Expressing (2.3) in terms of internal energy by removing the kinetic energy part, the corresponding transport equation reads as
\begin{equation} \frac {\partial \left ( \rho e \right ) }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \rho \boldsymbol{u} e \right ) = -P \thinspace \boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{u}) - \frac {1}{{\textit{Re}} {\textit{Br}}} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q} + \frac {1}{{\textit{Re}}}\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{\tau } \boldsymbol{\cdot }\boldsymbol{u}) + \boldsymbol{F} \boldsymbol{\cdot }\boldsymbol{u}. \end{equation}
Similar to the previous subsections, by substituting the expressions
$\rho = \rho _0 + \rho ^\prime$
,
$\boldsymbol{u} = \boldsymbol{u_0} + \boldsymbol{u}^\prime $
and
$e = e_0 + e^\prime$
into (C7), the equation reads as
\begin{align} \frac {\partial \left ( (\rho _0 + \rho ^\prime ) (e_0 + e^\prime )\right ) }{\partial t^\star } & + \boldsymbol{\nabla} ^\star \boldsymbol{\cdot }\left ( (\rho _0 + \rho ^\prime ) (\boldsymbol{u_0} + \boldsymbol{u}^\prime ) (e_0 + e^\prime ) \right ) \nonumber \\ & + (P_0 + P^\prime ) \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u_0} + \boldsymbol{u}^\prime ) = +\frac {1}{{\textit{Re}} {\textit{Br}}} \boldsymbol{\nabla }\boldsymbol{\cdot }\left ((\kappa _0 + \kappa ^\prime ) \boldsymbol{\nabla }(T_0 + T^\prime ) \right ) \nonumber \\ & + \frac {1}{{\textit{Re}}}\boldsymbol{\nabla }\boldsymbol{\cdot }\left [\left ( (\mu _0 + \mu ^\prime ) \left ( \boldsymbol{\nabla }(\boldsymbol{u_0} + \boldsymbol{u}^\prime ) + \boldsymbol{\nabla }(\boldsymbol{u_0} + \boldsymbol{u}^\prime )^{T} \right ) \right . \right . \nonumber \\ & + \left . \left . \vphantom{\left ( \boldsymbol{\nabla }(\boldsymbol{u_0} + \boldsymbol{u}^\prime ) + \boldsymbol{\nabla }(\boldsymbol{u_0} + \boldsymbol{u}^\prime )^{T} \right )} (\lambda _0 + \lambda ^\prime )(\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u_0} + \boldsymbol{u}^\prime ))\boldsymbol{I} \right ) \boldsymbol{\cdot }(\boldsymbol{u_0} + \boldsymbol{u}^\prime )\right ] + \boldsymbol{F}\boldsymbol{\cdot }(\boldsymbol{u_0} + \boldsymbol{u}^\prime ), \end{align}
which, by (i) substituting the terms
$e^\prime$
,
$P^\prime$
,
$\kappa ^\prime$
and
$\mu ^\prime$
with their corresponding Taylor series expansion approximation as a function of the independent thermodynamic pair of variables selected (
$\rho ^\prime$
,
$T^\prime$
), and (ii) linearising the equations, the following expression is obtained:
\begin{eqnarray} && \left ( e_0 + \rho _0 \frac {\partial e_0}{\partial \rho _0} \right ) \frac {\partial \rho ^\prime }{\partial t} + \rho _0 \frac {\partial e_0}{\partial T_0} \frac {\partial T^\prime }{\partial t} + \left ( e_0 u_0 + \rho _0 u_0 \frac {\partial e_0}{\partial \rho _0} \right ) \frac {\partial \rho ^\prime }{\partial x} \nonumber \\ && + \left ( \rho _0 e_0 + P_0 \right ) \left (\frac {\partial u^\prime }{\partial x} + \frac {\partial v^\prime }{\partial y} + \frac {\partial w^\prime }{\partial z} \right ) + \rho _0 u_0 \frac {\partial e_0}{\partial T_0} \frac {\partial T^\prime }{\partial x} = \frac {2}{\textit{Re}} \mu _0 \frac {\partial u_0}{\partial y} \frac {\partial v^\prime }{\partial x} + \frac {2}{\textit{Re}} \mu _0 \frac {\partial u_0}{\partial y} \frac {\partial u^\prime }{\partial y} \nonumber \\ && +\frac {1}{\textit{Re}\,\textit{Br}} \frac {\partial \kappa _0}{\partial \rho _0} \frac {\partial T_0}{\partial y} \frac {\partial \rho ^\prime }{\partial y} + \frac {1}{\textit{Re}\,\textit{Br}} \left (\frac {\partial \kappa _0}{\partial y} + \frac {\partial \kappa _0}{\partial T_0} \frac {\partial T_0}{\partial y} \right ) \frac {\partial T^\prime }{\partial y} + \frac {\kappa _0}{\textit{Re}\,\textit{Br}} \left ( \frac {\partial ^2 T^\prime }{\partial ^2 x} + \frac {\partial ^2 T^\prime }{\partial ^2 y} + \frac {\partial ^2 T^\prime }{\partial ^2 z}\right ) \nonumber \\ && + \frac {1}{\textit{Re}\,\textit{Br}} \left ( \frac {\partial ^2 T_0}{\partial y^2} \frac {\partial \kappa _0}{\partial \rho _0} + \left ( \frac {\partial ^2 \kappa _0}{\partial {\rho _0}^2} \frac {\partial \rho _0}{\partial y} + \frac {\partial ^2 \kappa _0}{\partial \rho _0 \partial T_0} \frac {\partial T_0}{\partial y} \right ) \frac {\partial T_0}{\partial y} \right ) \rho ^\prime + \frac {1}{\textit{Re}} \frac {\partial \mu _0}{\partial \rho _0} \left ( \frac {\partial u_0}{\partial y} \right )^2 \rho ^\prime \nonumber \\ && + \frac {1}{\textit{Re}\,\textit{Br}} \left ( \frac {\partial ^2 T_0}{\partial y^2} \frac {\partial \kappa _0}{\partial T_0} + \left ( \frac {\partial ^2 \kappa _0}{\partial {T_0}^2} \frac {\partial T_0}{\partial y} + \frac {\partial ^2 \kappa _0}{\partial T_0 \partial \rho _0} \frac {\partial \rho _0}{\partial y} \right ) \frac {\partial T_0}{\partial y} \right ) T^\prime + \frac {1}{\textit{Re}} \frac {\partial \mu _0}{\partial T_0} \left ( \frac {\partial u_0}{\partial y} \right )^2 T^\prime \nonumber \\ && + \frac {\hat {F}}{\textit{Re}} u^\prime . \end{eqnarray}
It is important to note that the Jacobian terms of the equation of state fields are obtained from the Coolprop library. However, the derivatives of the high-pressure transport coefficients
$\mu $
and
$\kappa$
are computed numerically by means of second-order finite differences with a derivation step of
$\Delta T = 10^{-2}\, \rm K$
and
$\Delta \rho = 10^{-2}\, \rm kg\, \rm m^{-3}$
, respectively. This relatively small step is robust enough to accurately calculate the gradients of the thermodynamic and transport coefficient fields at each grid point.
C.1. Energy balance equations
The turbulent kinetic energy can be expressed as the difference between the averaged Reynolds stress terms, which determines the production of energy rate due to transfer of energy from the base flow to the disturbances, and the dissipation energy. Therefore, the equation isolates the mechanisms by which energy is transferred from the base flow to the disturbances. In this regard, the contributors of the kinetic energy balance equation described in § 5.2.5 are detailed below. For ease of exposition, the notation `
$\hat {\boldsymbol{\cdot }}$
’ denoting perturbation vectors is dropped from the equations, yielding
\begin{align}&\qquad\qquad\qquad K = - i \omega \int \rho _0 (u u^\dagger + v v^\dagger ) {\rm d}y, \end{align}
\begin{align}&\qquad\qquad\qquad \varTheta = - i \alpha \int \rho _0 u_0 (u u^\dagger + v v^\dagger ) {\rm d}y, \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\quad P = - \int \rho _0 \frac {\partial u_0}{\partial y} v u^\dagger {\rm d}y, \end{align}
\begin{align} T = & - \int \left [ i \alpha \left ( \frac {\partial P_0}{\partial \rho _0} \rho + \frac {\partial P_0}{\partial T_0} T \right ) u^\dagger + \left (\frac {\partial P_0}{\partial \rho _0} \frac {\partial \rho }{\partial y} + \frac {\partial P_0}{\partial T_0} \frac {\partial T}{\partial y} \right ) v^\dagger \right . \nonumber \\ & + \left . \left ( \frac {\partial ^2 P_0}{\partial {\rho _0}^2} \frac {\partial \rho _0}{\partial y} + \frac {\partial ^2 P_0}{\partial {\rho _0} \partial T_0} \frac {\partial T_0}{\partial y} \right ) \rho v^\dagger + \left ( \frac {\partial ^2 P_0}{\partial {T_0}^2} \frac {\partial T_0}{\partial y} + \frac {\partial ^2 P_0}{\partial {T_0} \partial T_0} \frac {\partial \rho _0}{\partial y} \right ) T v^\dagger \right ] {\rm d}y, \end{align}
\begin{align} V = & \frac {1}{\textit{Re}} \int \left [ - \alpha ^2 (2 \mu _0 + \lambda _0) u u^\dagger + \mu _0 \frac {\partial ^2}{\partial y^2} u^\dagger + i \alpha (\mu _0 + \lambda _0) \frac {\partial v}{\partial y} u^\dagger \right . \nonumber \\ & + i \alpha \frac {\partial \mu _0}{\partial y} v u^\dagger + \frac {\partial \mu _0}{\partial \rho _0} \frac {\partial u_0}{\partial y} \frac {\partial \rho }{\partial y} u^\dagger + \frac {\partial mu_0}{\partial y} \frac {\partial u}{\partial y} u^\dagger + \frac {\partial \mu _0}{\partial T_0} \frac {\partial u_0}{\partial y} \frac {\partial T}{\partial y} u^\dagger \nonumber \\ & + \frac {\partial \mu _0}{\partial \rho _0} \frac {\partial ^2 u_0}{\partial y^2} \rho u^\dagger + \frac {\partial u_0}{\partial y} \left ( \frac {\partial ^2 \mu _0}{\partial {\rho _0}^2} \frac {\partial \rho _0}{\partial y} + \frac {\partial ^2 \mu _0}{\partial \rho _0 \partial T_0} \frac {\partial T_0}{\partial y} \right ) \rho u^\dagger \nonumber \\ & + \frac {\partial \mu _0}{\partial T_0} \frac {\partial ^2 u_0}{\partial y^2} T u^\dagger + \frac {\partial u_0}{\partial y} \left ( \frac {\partial ^2 \mu _0}{\partial {T_0}^2} \frac {\partial T_0}{\partial y} + \frac {\partial ^2 \mu _0}{\partial T_0 \partial \rho _0} \frac {\partial \rho _0}{\partial y} \right ) T u^\dagger \nonumber \\ & - \alpha ^2 \mu _0 v v^\dagger + (2 \mu _0 + \lambda _0) \frac {\partial ^2 v}{\partial y^2} v^\dagger + i \alpha (\mu _0 + \lambda _0) \frac {\partial u}{\partial y} v^\dagger \nonumber \\ & \left . + i \alpha \frac {\partial \mu _0}{\partial \rho _0} \frac {\partial u_0}{\partial y} \rho v^\dagger + i \alpha \frac {\partial \lambda _0}{\partial y} u v^\dagger + i \alpha \frac {\partial \mu _0}{\partial T_0} \frac {\partial u_0}{\partial y} T v^\dagger + \left ( 2 \frac {\partial \mu _0}{\partial y} + \frac {\partial \lambda _0}{\partial y} \right ) \frac {\partial v}{\partial y} v^\dagger \right ] {\rm d}y. \end{align}
C.2. Operator sensitivity
The sensitivity of the linear operator is briefly assessed for the isothermal I-1 case. In particular, the effect of the boundary conditions, wall temperature and bulk pressure on the maximum growth rate is quantified. Specifically, figure 31(a) highlights that the wall temperature linearly decreases the maximum growth at a constant rate of
$-0.0055$
units of maximum growth for a rise of a unit of
$T/T_c$
. Nevertheless, when approaching the pseudo-boiling region, the sensitivity reverses sign with an increase of maximum growth with temperature. Moreover, figure 31(b) shows that the growth rate monotonically decreases with the system pressure by a rate of
$-10^{-4}$
units of growth for a unit of normalised pressure. To this extent, Bernades (Reference Bernades2024) has extensively characterised the sensitivity of the linear operator comparing ideal gas and real-gas frameworks under non-isothermal conditions.
Figure 31. Operator sensitivity of the wall temperature (a) and bulk pressure (b) for the maximum growth rate for case I-1 at
$\alpha = 1.0$
,
${\textit{Br}} = 0.1$
and
${\textit{Re}} = 10\,000$
.
Figure 32. (a) Eigenspectrum at
${\textit{Re}} = 10\,000$
and wavenumber
$\alpha = 1$
and (b) neutral curve. Real-gas framework with CoolProp thermodynamic and transport properties model (RG), ideal gas with power law (IG) and incompressible framework (IC). In (a) the red highlighted eigenvalue corresponds to an unstable mode (
$\omega = 0.2375 + 0.0037i$
), whereas dark yellow corresponds to a stable mode (
$\omega = 0.4164 - 0.1382i$
) whose perturbations are depicted in figure 33.
Appendix D. Linear stability verification at the isothermal limit
The linear stability solver is verified with respect to the incompressible reference (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993) and extended replications (Trefethen Reference Trefethen2000). In this regard, the base flow is solved from the stability equations at the isothermal limit with (i) the CoolProp real-gas and transport coefficients framework, and (ii) the ideal gas equation of state and power law for the transport coefficients. For this comparison, wall scalings have been used with carbon dioxide (CO
$_2$
) as the operating fluid to match the Poiseuille flow conditions from Ren et al. (Reference Ren, Marxen and Pecnik2019b
) at
$T_{cw} = T_{hw} = 290\thinspace \textrm {K}$
and
$P = 8\thinspace \textrm {MPa}$
. The spectrum is computed at
${\textit{Re}} = 10\,000$
and wavenumber
$\alpha = 1$
solved for 2-D perturbations, i.e.
$\beta = 0$
. Figure 32(a,b) depicts the corresponding eigenspectrum and neutral curve, respectively. It can be observed that the
$A, P$
and
$S$
branches of Mack (Reference Mack1976) are reproduced by the incompressible case. At the isothermal limit, the ideal gas result collapses to these branches. However, the real-gas solution adds additional modes due to thermodynamic effects, i.e. modes driven by thermodynamics. The same resulting unstable mode is captured by the studied frameworks as it is a hydrodynamic mode independent of thermodynamic effects. Hence, the thermodynamic and transport property models do not influence the hydrodynamic modes in the isothermal limit. Moreover, the neutral curve shows good agreement for the different models considered and, in fact, they all provide the same critical Reynolds number
${\textit{Re}}_c = 5772$
determined for Poiseuille flow by Thomas (Reference Thomas1953) and Orzag (Reference Orzag1971). Next, figure 33(a,b) depicts the perturbations of the unstable mode and stable mode (dark yellow highlighted eigenvalue in figure 32
a), respectively. Thus, the statement claimed above can be verified since the unstable modes of the ideal gas and real-gas frameworks are the same. Instead, some small differences emerge for the stable mode since the compressibility effects on the density perturbation are represented by the real-gas framework as shown in the inset of figure 33(b).
Figure 33. Perturbation profiles of (a) the unstable mode (
$\omega = 0.2375 + 0.0037i$
) and (b) the stable mode (
$\omega = 0.4164 - 0.1382i$
) normalised by
$| u^{\prime }|$
. Real-gas framework with CoolProp thermodynamic and transport properties model (RG) depicted by solid lines, and ideal gas with power law (IG) with markers.
Figure 34. Transient growth map at incompressible conditions (
${\textit{Br}} \approx 0$
) for
${\textit{Re}} = 2000$
.
Appendix E. Transient growth verification
The transient growth framework is firstly verified against the incompressible reference (Reddy & Henningson Reference Reddy and Henningson1993); refer to figure 15. Figure 34 highlights that the largest transient growth is achieved at
$\alpha = 0$
and
$\beta \approx 2$
. The classic optimal perturbation with streamwise vortices, and the corresponding outputs with streaks, are also recovered as depicted in figure 35. Hence, no significant differences are observed when operating at the isothermal limit using the real-gas framework with respect to the incompressible flow reference solution.
Figure 35. Optimum eigenvector profiles at incompressible conditions (
${\textit{Br}} \approx 0$
) for
${\textit{Re}} = 2000$
at maximum growth (
$\alpha = 0$
and
$\beta = 2$
) for (a) input and (b) output. Results are normalised by (a)
$w^\prime$
and (b)
$u^\prime$
.
Figure 36. Transient growth maps at
${\textit{Re}} = 1000$
and
${\textit{Br}} = 0.07$
for CO
$_2$
at isothermal conditions with (a)
$T = 290\thinspace \textrm {K}$
, (b)
$T = 300\thinspace \textrm {K}$
and (c)
$T = 310\thinspace \textrm {K}$
.
Finally, the transient growth framework is also verified against the high-pressure transcritical results presented by Ren et al. (Reference Ren, Marxen and Pecnik2019b
) at isothermal wall conditions with CO
$_2$
as the working fluid (refer to figures 12, 13 and 14). Despite the utilisation of different thermodynamic frameworks, the largest amplifications and their corresponding wavenumbers are properly recovered. It is important to note that the optimum perturbation profiles are also accurately obtained. The small differences observed can be attributed to the large sensitivity of the operator to the thermodynamic framework utilised.
Open access
