2.1. Fluid models

Four fluids will be considered throughout this paper, each of them being associated with a different equation of state and viscosity law. However, they all share the common property of assuming an extremum of kinematic viscosity $\nu$ at a given temperature. Recalling that $\nu = \mu / \rho$, where $\mu$ and $\rho$ are the dynamic viscosity and the density, respectively, different ways to generate a minimum of $\nu$ can be imagined. Three theoretical fluid models will be used to control and study a restricted number of parameters. A fourth more realistic model for supercritical fluids, based on the van der Waals equation of state, will also be considered. A summary of the different fluids is provided in figure 1 while a detailed description is given in the next subsections.

Figure 1.

The four fluid models considered in this paper. The definition of $\nu _\infty$, used to normalise $\nu$, is given in Appendix B in the particular case of fluid VdW.

2.1.1. Fluid VB: bump of dynamic viscosity with constant density

In this model, density is assumed constant, while the viscosity is chosen to locally exhibit a bump at a temperature $T_m^*$, using a Gaussian function (in this paper, all dimensional quantities are noted using the superscript ‘*’). Note that, in this case, the non-dimensional kinematic viscosity is equal to the non-dimensional dynamic viscosity, both reading

(2.1)\begin{equation} \mu = \nu = 1 + A_{\nu} \exp\left({ - \left( \frac{T - 1 }{ \delta T_\nu } \right)^2 }\right), \end{equation}

using $T_m^*$ as the reference temperature. The reference viscosities are $\mu _\infty ^*$ and $\nu _\infty ^*$, which are the asymptotic values away from the bump. The parameter $A_\nu$ controls the amplitude of the bump, and its sign determines whether the kinematic viscosity admits a minimum ($A_\nu < 0$) or a maximum ($A_\nu > 0$). The characteristic width of the bump is set through $\delta T_{\nu }$, which is again made non-dimensional using $T_m^*$. Finally, the thermal conductivity $\lambda$ is assumed to be constant.

2.1.2. Fluid DB: bump of density with constant dynamic viscosity

Inversely to fluid VB, $\mu$ is kept constant in fluid DB while a bump is introduced in the density profile. This bump is chosen such that the resulting kinematic viscosity has the same expression as in (2.1). Hence, the density law simply reads

(2.2)\begin{equation} \rho = \frac{1}{ 1 + A_{\nu} \exp\left({ - \left( \dfrac{T - 1 }{ \delta T_\nu } \right)^2 }\right) }, \end{equation}

and since, in this case, $\nu = 1/ \rho$, $\nu$ is the same as in fluid VB. The conductivity $\lambda$ is again chosen to be constant.

2.1.3. Fluid HT: hyperbolic tangent laws

In fluid HT, thermal conductivity is also kept constant, while dynamic viscosity and density are now both allowed to vary according to hyperbolic tangent laws. In order to generate an extremum of kinematic viscosity, a small shift is introduced between the two hyperbolic tangents, controlled by the non-dimensional parameter $\varepsilon _T$. This choice is inspired by supercritical fluids and represents an attempt to mimic some of their features in the vicinity of the pseudo-boiling region. This will be discussed in more detail in § 2.1.4 after the supercritical fluid laws are introduced. The non-dimensional governing laws for fluid HT are formally written as

(2.3)\begin{gather} \rho = 1 - \gamma \tanh \left[ \frac{ T - 1 }{ \delta T_{\mu,\rho} } \right], \end{gather}

(2.4)\begin{gather}\mu = 1 - \gamma \tanh \left[ \frac{ T - (1+\varepsilon_T) }{ \delta T_{\mu,\rho} } \right]. \end{gather}

The reference temperature $T_m^*$ is here defined as the point of anti-symmetry of the density profile. The density at $T=1$ and the viscosity at $T=1+\varepsilon _T$ are used as the reference scales. The parameter $\gamma$ controls the jump of density and dynamic viscosity while $\delta T_{\mu,\rho }$ sets the temperature range over which this jump takes place. The fluid properties are shown in figure 1, where it is verified that the kinematic viscosity admits a minimum around $T=1$. By analogy with fluids VB and DB, it is possible to estimate the amplitude $A_\nu$ of this minimum, as well as the characteristic width $\delta T_\nu$ of the temperature range onto which it occurs. The following relations will be used in this paper:

(2.5)\begin{gather} A_\nu =\frac{\gamma \varepsilon_T}{\delta T_{\mu,\rho}} , \end{gather}

(2.6)\begin{gather}\delta T_\nu = 1.2 \delta T_{\mu,\rho} . \end{gather}

The derivation and verification of these expressions are detailed in Appendix A.

2.1.4. Fluid VdW: a model for supercritical fluids

The reduced van der Waals equation of state is used and reads

(2.7)\begin{equation} \check{p} = \frac{8 \rho \check{T}}{3 - \check{\rho}} - 3 \check{\rho}^2, \end{equation}

where the reduced variables $\check {p} = p^* / p_c^*$, $\check {\rho }= \rho ^* / \rho _c^*$ and $\check {T}= T^* / T_c^*$ have been introduced, with $p_c^*$, $\rho _c^*$ and $T_c^*$ being the critical pressure, density and temperature, respectively. A choice of diffusion laws is required in addition to this equation of state. The models proposed by Jossi, Stiel & Thodos (Reference Jossi, Stiel and Thodos1962) and Stiel & Thodos (Reference Stiel and Thodos1964) are used for the dynamic viscosity and the thermal conductivity, respectively. They provide analytical expressions for non-polar supercritical fluids based on theoretical scalings and experimental fittings. In supercritical fluids, these diffusion laws depend both on $\check {T}$ and $\check {\rho }$. The density, dynamic and kinematic viscosity profiles are plotted in figure 1. Note that the reference temperature, here again noted $T_m^*$ to maintain consistency with the previous fluids, is usually termed pseudo-boiling or pseudo-critical temperature in supercritical fluids (Banuti Reference Banuti2015). Density and dynamic viscosity are strongly correlated, and both exhibit strong gradients in the pseudo-boiling region. This motivated the choice of fluid HT, where density and viscosity are both defined using a hyperbolic tangent function, aiming at capturing these gradients while neglecting other variations away from them. The kinematic viscosity admits a minimum around the pseudo-critical temperature but is not localised, as opposed to the other fluids. The relatively simple model of fluid HT is found to decently reproduce this minimum as a result of the shift $\varepsilon _T$ introduced between the hyperbolic tangent laws, but differs away from the point where $\nu$ remains strictly constant in fluid HT.

In analogy to the previous fluids, we would like to extract the characteristic scales $\delta T_\nu$ and $A_\nu$ from the kinematic viscosity law. However, while $\nu$ does have a minimum in fluid VdW, it is not clear that this minimum is localised over a finite, identified range $\delta T_\nu$. Still, it can be observed, after calculation, that $\nu (T)$ admits two inflection points in the vicinity of the pseudo-boiling temperature – one below and one above this temperature. This can be used to define the scale $\delta T_\nu$ as the width between these two inflection points. From this, an amplitude $A_\nu$ can be naturally defined. The procedure is thoroughly described in Appendix B. Finally, note that the reduced pressure is the control parameter of the kinematic viscosity seen as a function of the temperature. In other words, $\nu (T)$ is different for each $\check {p}$ and, consequently, so are $\delta T_\nu$ and $A_\nu$.

2.2. Base flow

Linear stability analysis requires the knowledge of a base flow, defined as a steady solution of the nonlinear Navier–Stokes equations. After recasting the nonlinear Navier–Stokes equations given the physical assumptions associated with this flow, the equations are numerically solved. Plane Couette flow occurs between two plates and is driven by the upper plate moving at speed $u_1^*$, which is used as the reference velocity scale. The streamwise and wall-normal directions are noted with $x$ and $y$, respectively. The flow is assumed to be parallel: the streamwise velocity $u$ does not depend on $x$, and the wall-normal and spanwise components of the velocity $v$ and $w$, respectively, are zero. The lower plate is fixed and, given the no-slip conditions, the non-dimensional streamwise velocity at the boundaries verifies $u(0)=0$ and $u(1)=1$. The distance $h^*$ between the two plates is used as the reference length scale. The lower plate is kept at temperature ${T_0}^*$, chosen as the reference temperature. We choose to consider the non-dimensional temperature gradient $\tau$ between the two plates as an input parameter, which in turn sets the temperature of the upper plate. The boundary conditions for the temperature are then $T(0)=1$ and $T(1)=1 + \tau$. Under the assumption of a steady flow without pressure gradient – the flow is driven by the top wall – the non-dimensional Navier–Stokes equations reduce to a system of ordinary differential equations:

(2.8)\begin{gather} ( \bar{\mu} \bar{u}' )' = 0, \end{gather}

(2.9)\begin{gather}( \bar{\lambda} \bar{T}' )' = 0, \end{gather}

where the superscript $'$ denotes the wall-normal derivative and the overbars identify base-flow variables. The inertial terms are zero given the parallel flow assumption, and the problem does not depend on the Reynolds and Prandtl numbers. Besides, the temperature is decoupled from the velocity field. When $\lambda$ is constant, as it is supposed to be in fluids VB, DB and HT, the temperature profile is readily obtained as $T(y) = 1 + \tau y$. As for fluid VdW, (2.9) is solved using Newton's method by setting the initial guess as the aforementioned linear profile. Once $T$ is obtained, the density profile is also known via the equation of state. The velocity profile is finally obtained by integration of (2.8) with the knowledge of the dynamic viscosity profile as a function of $T$ and $\rho$. Finally, note that we make the arbitrary choice to locate the extremum of kinematic viscosity at the centre line of the flow, $y=1/2$. This is achieved by accordingly setting $T_{m}^*/{T_0}^* = 1 + \tau /2$ under the assumption that the temperature profile is linear – which is indeed the case for fluids VB, DB and HT.

2.3. Inviscid linear stability theory

2.3.1. Rayleigh equation with density gradients

Assuming infinitely small, inviscid, two-dimensional perturbations, the linearised Navier–Stokes equations in the low-Mach approximation (Rehm & Baum Reference Rehm and Baum1978; Paolucci Reference Paolucci1982) can be written

(2.10)\begin{gather} \frac { \partial u } { \partial x } + \frac { \partial v } { \partial y } = 0, \end{gather}

(2.11)\begin{gather}\bar{\rho} \left( \frac { \partial u } { \partial t } + \bar{u} \frac { \partial u } { \partial x } \right) ={-} \frac { \partial p } { \partial x } - \bar{\rho} \frac { \partial \bar{u} } { \partial y } v, \end{gather}

(2.12)\begin{gather}\bar{\rho} \left( \frac { \partial v } { \partial t } + \bar{u} \frac { \partial v } { \partial x } \right) ={-} \frac { \partial p } { \partial y } , \end{gather}

(2.13)\begin{gather}\frac { \partial T } { \partial t } + \bar{u} \frac { \partial T } { \partial x } ={-} \frac { \partial \bar{T} } { \partial y } v. \end{gather}

Perturbations of the form $q(x,y,t) = \Re \{ \hat {q}(y) \exp ({{\rm i}(\alpha x - \omega t)})\}$ are now considered, with ${q = [u, v, T]}$ being the state vector of the perturbations and $\Re \{ \}$ the real part. These linearised equations can then be recast into the Rayleigh equation governing the linear dynamics of incompressible flows with density gradients (see also Fontane & Joly Reference Fontane and Joly2008):

(2.14)\begin{equation} \bar{\rho} ( \bar{u} - c) \left[ \hat{v}'' + \frac{\bar{\rho}'}{\bar{\rho}} \hat{v}' - \alpha^2 \hat{v} \right] - (\bar{\rho} \bar{u}')' \hat{v} = 0, \end{equation}

where $c = \omega / \alpha$ is the complex phase velocity. Note that the disturbance temperature does not appear in (2.14) since the linearised mass and momentum equations are decoupled from the energy equation (2.13). Temperature disturbances are deduced from the hydrodynamic disturbances, which can be calculated independently. Thermal effects are, however, at play in the velocity and density profiles of base flow, which the momentum equations (2.11) and (2.12), and ultimately the Rayleigh equation (2.14), depend on. A temporal framework is adopted: the wavenumber $\alpha$ is a real parameter while the frequency $\omega$ is a complex number that is to be determined. The temporal growth rate is given by its imaginary part, $\omega _i$. A positive value corresponds to an inviscid instability. The (real) phase velocity $c_\varphi$ of the perturbation is simply $c_r$, the real part of $c$.

Equation (2.14) can be classically solved numerically as an eigenvalue problem. The boundary condition $\hat {v} = 0$ is used at the wall. A pseudo-spectral method is employed to discretise the system and to obtain the derivative matrices (Orszag Reference Orszag1971). In order to avoid the singularity at the critical layer for neutral modes, a parabolic complex mapping is used, following Boyd (Reference Boyd1985). This allows the growth rate to be computed even when it reaches small values, while a real mapping would produce spurious numerical oscillations.

2.3.2. Vorticity

In order to interpret some results, it can be useful to consider an alternative formulation of the problem in terms of the disturbance vorticity $\xi = \partial v / \partial x - \partial u / \partial y$. For a parallel base flow without pressure gradients, $\xi$ is governed, in the physical space, by the linear equation

(2.15)\begin{equation} \frac { \partial \xi } { \partial t } + \bar{u} \frac { \partial \xi } { \partial x } = \underbrace{- \varOmega' v }_{S_\xi} \underbrace{ - \frac{\bar{\rho}'}{\bar{\rho}^2} \frac { \partial p } { \partial x } }_{B_\xi}, \end{equation}

where $\varOmega$ is the vorticity of the base flow. The left-hand side represents the material derivative of $\xi$ by the base flow. The right-hand side corresponds to vorticity sources, which may induce an instability. The term $S_\xi$ is the production of vorticity responsible for shear flow instabilities. The second term, $B_\xi$, is the inertial baroclinic torque, which may generate vorticity when the density and pressure gradients are not aligned. In the absence of density gradients, this term is evidently zero.

5.1. Piecewise linear base flows

Piecewise linear base flows have been extensively used to study a variety of stability problems with constant density as well as variable density – usually in the framework of the Taylor–Goldstein equation, under the Boussinesq approximation (Drazin & Howard Reference Drazin and Howard1966). While being simple approximations, useful analytical predictions can be derived from these models, especially predicting the linear stability of long waves (Gallaire Reference Gallaire2015). We will here consider arbitrary large variations of density.

The base flows of fluids VB, DB and HT are divided into three layers. The central layer, centred around $y=1/2$, has a width $\delta$ (figure 3). This approach follows that proposed by Rayleigh (Reference Rayleigh1887) for bounded, constant-density flows. Velocity profiles are continuous at the interfaces between layers, but their gradient may not be; a vorticity jump may occur at the interface. This is the case of fluid VB, where the viscosity bump is modelled by a discontinuous jump in the central layer (figure 3a). This generates a stronger shear rate (i.e. vorticity) in the central layer. Note that this is precisely the configuration studied by Rayleigh (Reference Rayleigh1887). In fluid DB, the shear rate is constant throughout the flow, but density exhibits a jump in the central layer (figure 3b). In fluid HT, density linearly varies in the central layer but remains constant in the two other ones (figure 3c). The same profile is used for the dynamic viscosity. All fluids feature an excess $\Delta \varPhi >0$ of $|\varPhi |$ in the central layer in order to model the smooth profile of $|\varPhi |$ that was observed in the previous section (figure 3d).

Figure 3.

Piecewise linear base-flow profiles for fluids VB, DB and HT. The central layer (layer 2) is centred around $y=1/2$.

5.2. Expressions of $\delta$ and $\Delta \varPhi$

The relations between the parameters $\delta$ and $\Delta \varPhi$ of the sought stability model and the physical input parameters of our system $A_\nu$, $\delta T_\nu$ and $\tau$ are now examined. The region of the smooth base-flow profiles across which $\varPhi$ varies is the same as that across which kinematic viscosity varies, as expressed in (3.4). The characteristic length of this region in the base-flow profile is proportional to $\delta T_\nu$, which is related to the fluid property, and inversely proportional to the temperature gradient $\tau$ of the flow. The thickness of the layer $\delta$ thus follows the proportionality relation:

(5.1)\begin{equation} \delta \propto (1 + \tau/2) \frac{\delta T_\nu}{\tau}, \end{equation}

where the factor $(1 + \tau /2)$ results from the factor $T_m^* / T_0^*$ (see § 2.2) that appears when $\delta T_\nu$ is made dimensionless with $T_0^*$. A choice of a prefactor is ultimately required in order to assign a definitive value to $\delta$ in (5.1), and will be specified for each fluid.

We define the quantity $\Delta \varPhi$ as the jump of $|\varPhi |$ at the interface: $\Delta \varPhi = \varDelta ( |- \bar {\rho } \bar {u}' | )$. The following convention is used: $\Delta \varPhi > 0$ corresponds to configurations in which the magnitude of $|\varPhi |$ is larger in the central layer than in the other layers. Given that $\bar {\rho } \bar {u}' = \bar {\mu } \bar {u}' / \bar {\nu }$ and that integrating the momentum equation (2.8) yields $\bar {\mu } \bar {u}' = K$, with $K$ a positive constant, we can express $\Delta \varPhi$ as

(5.2)\begin{equation} \Delta \varPhi = K \varDelta \left( \frac{1}{\bar{\nu}} \right) . \end{equation}

The characteristic value of $\bar {\nu }$ being $1+A_\nu$ in the central region and $1$ elsewhere, the jump of $1/\bar {\nu }$ at the interface reads

(5.3)\begin{equation} \varDelta \left( \frac{1}{\bar{\nu}} \right) ={-} \frac{A_\nu}{1 + A_\nu} . \end{equation}

An excess of $|\varPhi |$ in the central layer ($\Delta \varPhi > 0$) is associated with a minimum of $\nu$ ($A_\nu < 0$), consistent with § 3. The derivation of the different expressions of the constant $K$ and the final expression of $\Delta \varPhi$ associated with each fluid model is detailed in Appendix E.

5.3. Stability calculations in the long-wave approximation

5.3.1. Derivation

The Rayleigh equation (2.14) is solved for the three piecewise linear base flows introduced in § 5.1. We will restrict our analysis to long waves ($\alpha \ll 1$). Following Rayleigh (Reference Rayleigh1887) (see also Drazin & Reid Reference Drazin and Reid2004; Charru Reference Charru2011), (2.14) is first solved separately in each of the three layers of the piecewise linear base flows. As $\varPhi$ is constant in each layer, the last term of the Rayleigh equation vanishes. Furthermore, assuming $\alpha \ll 1$ and writing $\hat {v}$ and $c$ as a power series of $\alpha$, (2.14) reduces, at the order $\alpha ^0$, to

(5.4)\begin{equation} (\bar{\rho} \hat{v}') ' = 0 . \end{equation}

This equation can be solved in each layer. When density is constant across a layer, the solution is simply

(5.5)\begin{equation} \hat{v}_i = A_{i} y + B_{i}, \end{equation}

where the index ‘$i$’ refers to the layer 1, 2 or 3 (see figure 3). In the central layer of fluid HT, where density varies, only the first derivative of $\hat {v}$ will be needed. This is because interface conditions, described hereafter, set the value of $\hat {v}_2$ using $\hat {v}_1$ and $\hat {v}_3$ (readily obtained from (5.5)). The following expression is immediately found:

(5.6)\begin{equation} \hat{v}_2' = \frac{ A_{2} y }{\bar{\rho}}. \end{equation}

At each interface between the layers, the kinematic and dynamic conditions (Charru Reference Charru2011) read

(5.7)\begin{gather} \varDelta \left[ \frac{\hat{v}}{\bar{u} - c} \right] = 0 , \end{gather}

(5.8)\begin{gather}\varDelta [ (\bar{u} - c ) \hat{v}' + \varPhi \hat{v}] = 0 . \end{gather}

Note that (5.7) reduces to $\Delta \hat {v} = 0$ as $\bar {u}$ is continuous. Using (5.7) and (5.8) as well as the boundary conditions $\hat {v}=0$ at the walls lead to a linear system on the coefficients $A_i$ and $B_i$. Equating the determinant to zero provides an expression of $c^2$. The system is unstable for $c^2<0$ and stable for $c^2>0$. In this section, derivations are carried out in the frame of reference moving with $\bar {u}(y=1/2) = \bar {u}_{1/2}$.

5.3.2. Results for fluid VB

As previously mentioned, the dispersion relation for fluid VB corresponds to that derived by Rayleigh (Reference Rayleigh1887). We will use his result in the long-wave regime. Using geometrical reasoning on the piecewise linear base flow of fluid VB (figure 3), $\Delta u$ can be written as $\Delta u = \delta [ 1 + \Delta \varPhi ( 1 - \delta ) ]$, recalling that $\Delta \varPhi = \Delta \varOmega$ for this fluid. Injecting this into Rayleigh's result yields

(5.9)\begin{equation} c^2 = \frac{\delta}{4} (\delta - \Delta \varPhi (1 - \delta)^2). \end{equation}

Note that $\Delta \varPhi$ was defined as $\Delta \varPhi = \varDelta (|\varPhi |)$ in § 5.2, and that $\Delta \varPhi >0$ corresponds to an excess of $|\varPhi |$ in the central layer. The instability criterion $\Delta \varPhi > \Delta \varPhi _c$ is deduced, with the instability threshold $\Delta \varPhi _c$:

(5.10)\begin{equation} \Delta \varPhi_c = \frac{\delta}{(1 - \delta)^2}. \end{equation}

5.3.3. Results for fluid DB

The expression for $c^2$ in fluid DB is solved by first remarking that $\Delta u = \delta$, given that $\bar {u}(y) = y$ in this case. It is also noticed that $\Delta \rho = \Delta \varPhi$. Analytical calculations lead to the relation

(5.11)\begin{equation} c^2 = \frac{\delta [\delta - \Delta \varPhi (1 - \delta)]}{4 [1 + \Delta \varPhi( 1 - \delta)]}. \end{equation}

Because $\delta < 1$ and $\Delta \varPhi >-1$ (since $\bar {\rho }>0$), the denominator is always strictly positive. The instability criterion is thus deduced from the sign of the numerator, leading to $\Delta \varPhi > \Delta \varPhi _c$, with

(5.12)\begin{equation} \Delta \varPhi_c = \frac{\delta}{(1 - \delta)}. \end{equation}

5.3.4. Results for fluid HT

In fluid HT, the transformation of $\Delta \rho$ into the parameters $\delta$ and $\Delta \varPhi$ cannot be found. We will then use $\Delta \rho = 2 \gamma / (1+\gamma )$, which comes from the definition ${\Delta \rho = \bar {\rho }(0) - \bar {\rho }(1)}$ introduced in figure 3 while using $\bar {\rho }^*(0)$ as the reference scale. The stability model will then depend on $\gamma$ in addition to $\delta$ and $\Delta \varPhi$, but will become independent of $\gamma$ in some regimes of interest. The approximation $\Delta u = \delta$ will be used, which is valid when the central layer or the viscosity jump are small. This will be shown to be of practical interest for the more realistic fluid VdW. Under these considerations, the following expression can be derived:

(5.13)\begin{equation} c^2 = \tfrac{1}{4} \{ \sqrt{ \delta [ 1 - \Delta \varPhi_0 (1 - \delta)] [ \delta - \Delta \varPhi_0(1-\delta) + \gamma^2 \Delta \varPhi_0 (1 - \delta)^2] } + \delta \gamma \Delta \varPhi_0 (1 - \delta)\}^2, \end{equation}

where the quantity

(5.14)\begin{equation} \Delta \varPhi_0 = \frac{ \Delta \varPhi }{1 - \gamma} , \end{equation}

is introduced to simplify (5.13). It will also be shown to be of practical interest for fluid VdW as $\Delta \varPhi _0$ only depends on the jump of kinematic viscosity $\varDelta (1/\bar {\nu })$ (Appendix E). A criterion of instability is obtained by examining the sign of the expression under the square root in (5.13). Given that the factor ${1 - \Delta \varPhi _0 (1 - \delta )}$ is always positive, the criterion of instability is given by the third factor, reading $\Delta \varPhi > \Delta \varPhi _c$ with

(5.15)\begin{equation} \Delta \varPhi_c = \frac{\delta ( 1 - \gamma)}{(1 - \delta)[1 + \gamma^2 (1 - \delta)]}. \end{equation}

5.3.5. Comments and limiting case

The three different criteria all state that a certain excess of $|\varPhi |$ is required in the central layer ($\Delta \varPhi > \Delta \varPhi _c$) to generate an instability. This can be seen as an improvement of the Fjørtoft criterion, which states that an instability may occur only if $\Delta \varPhi >0$, but does not specify the magnitude of the excess of $|\varPhi |$ that is sufficient to make the system unstable. It can be noticed that the derived instability thresholds always increase with $\delta$. In the limit of small $\gamma$, fluid DB and fluid HT possess the same criterion of instability, (5.15) reducing to (5.12). Furthermore, in this regime, these equations differ from that of fluid VB (5.10) only at order $\mathcal {O}(\delta ^2)$. At order $\mathcal {O}(\delta )$, all three fluids share the common criterion of instability:

(5.16)\begin{equation} \Delta \varPhi > \delta . \end{equation}

Moreover, the growth rate near the instability threshold can be calculated from the different expressions of $c^2$, using $\Delta \varPhi = \mathcal {O}(\delta ) \ll 1$. A general expression is obtained for the unstable modes of all fluids:

(5.17)\begin{equation} \omega_i = \alpha \sqrt{ \delta ( \Delta \varPhi - \delta ) } . \end{equation}

This shows the fundamental role the quantities $\Delta \varPhi$ and $\delta$ play in modelling these instabilities, regardless of the types of stratifications.

5.3.6. Phase velocity

The equations on $c^2$ obtained for each fluid also provide interesting results regarding the phase velocity of the unstable modes. For fluids VB and DB, if $c^2 < 0$, then $c$ is purely imaginary. Therefore, an unstable mode will have a phase velocity $\bar {u}_{1/2}$, which, by symmetry, is equal to $1/2$ for these fluids. This does not hold for fluid HT, for which the phase velocity is shifted from $\bar {u}_{1/2}$ by ${\delta \gamma \Delta \varPhi _0 (1 - \delta )}$. By integrating the conservation of shear stress in the three layers, the following expression can be obtained for fluid HT:

(5.18)\begin{equation} \bar{u}_{1/2} = \frac{1}{2} - \frac{\gamma (1 - \delta)}{2} . \end{equation}

Note that the departure from $1/2$ in (5.18) is of order $\gamma$, while the aforementioned additional shift is at most $\gamma ^2$ for small $\Delta \varPhi _0$ given that ${\Delta \varPhi _0 \sim \gamma | \varepsilon _T / \delta T_{\mu,\rho } | < \gamma }$ (Appendix E). Equation (5.18) is therefore expected to be a good approximation of the phase velocity.

6.1. Growth rate and phase velocity

Numerical stability calculations are carried out for fluids VB, DB and HT using three different thicknesses of the central layer $\delta$. Note that $\delta$ is initially not an input parameter of the problem: it is calculated from (5.1), which requires a prefactor. This prefactor, which is a priori different for each fluid, is set so as to yield the best agreement between the calculated and predicted stability diagram, which will be presented in the next subsection. The value 1.12 is used for fluids VB and DB. While not fully predictive – this value is not obtained by the model and requires one calculation point in order to be calibrated, it remains close to one: the model can provide order-of-magnitude predictions even without further knowledge. A prefactor equal to 1 is used for HT, requiring no external data.

First, a constant value of $A_\nu =-0.04$, which sets the magnitude of $\Delta \varPhi$ for small $\delta$ (E9), is chosen for all fluids. Table 1 indicates the corresponding values of $\Delta \varPhi /\Delta \varPhi _c$, providing a useful reference when comparing with figure 5, which will be presented later. As shown in figure 4(a–c), similar behaviours as well as close quantitative values of the growth rate are found for all fluids, despite the fundamentally different stratifications in each fluid. All fluids exhibit a long-wave instability: low wavenumbers are always unstable, while a cutoff wavenumber $\alpha _c$ exists beyond which the system is stable. This observation justifies the restricted analysis to long waves to predict the stability of the system (§ 5.3). As $\delta$ increases, $\alpha _c$ decreases, as it would for a constant-density, unbounded shear layer for which ${\alpha _c \sim 1/\delta }$ (Charru Reference Charru2011). At constant $A_\nu$, the maximum growth rate increases as $\delta$ is reduced. However, this does not hold for the growth rate of long waves, for which confinement can have a destabilising effect – increasing $\delta$ being equivalent to approaching the walls closer to the central layer. This behaviour is not unexpected given that $c^2$ depends on polynomials of $\delta$ of degree larger than 1 (§ 5.3). Differentiating (5.17) with respect to $\delta$, a general estimate of the value $\delta _m$ that yields the maximum growth rate is found to be ${\delta _m = \Delta \varPhi /2}$. Using (E9) for $\Delta \varPhi$, we find that $\delta _m \simeq 0.02$, which is consistent with the observations. Note that a different behaviour is observed in a bounded, constant-density mixing layer, in which Healey (Reference Healey2009) found that confinement reduces the temporal growth rate of the instability.

Table 1.

Values of $\Delta \varPhi /\Delta \varPhi _c$ for the different fluids and values of $\delta$ presented in figure 4.

Figure 4.

Growth rate, as a function of $\alpha$, of the fluid models VB, DB and HT obtained for $A_\nu = -0.04$ and different thicknesses of the central layer $\delta$.

The growth rate of long waves is presented in figure 5(a–c). The quantity ${\omega _i / \alpha }$ is plotted, which corresponds to the slope of the growth rate at $\alpha =0$. For all fluids, the instability threshold is well predicted. The behaviour of the growth rate past the threshold is also reasonably well captured, but piecewise linear models do not yield exact quantitative matches. The validity of the general approximation of the growth rate (5.17) is verified for $\delta =0.003$. At $\delta =0.01$, this approximation still captures the threshold well. Yet, significant departures from the full theoretical prediction are observed for fluids DB and HT; for these fluids, more terms are indeed neglected in the derivation leading to (5.17). At larger $\delta$, noticeable differences appear for all fluids. Finally, the phase velocity of the unstable mode in each fluid is very well predicted by theoretical models (figure 5d–f). It is equal to $1/2$ for fluids VB and DB, and does not depend on $\delta$ and $\Delta \varPhi$. The phase velocity markedly differs from $1/2$ in fluid HT as the velocity of the base flow $\bar {u}_{1/2}$ depends on both $\delta$ and the viscosity ratio $\gamma$ (on which $\Delta \varPhi$ depends), as discussed § 5.3.6.

Figure 5.

(a–c) Growth rate, normalised by $\alpha$, of long waves, obtained at $\alpha = 10^{-2}$. (d–f) Associated phase velocities. Solid lines are the numerical results and dashed lines are the theoretical predictions given, for each fluid, in (5.9), (5.11) and (5.13). Black dashed lines correspond to the theoretical general growth rate of (5.17), obtained for $\delta \ll 1$ near the instability threshold. Instability thresholds $\Delta \varPhi _c$ are given in (5.10), (5.12) and (5.15). Theoretical phase velocities are given in the frame of reference associated with the lower wall using the values of $\bar {u}_{1/2}$ discussed in § 5.3.6. For all figures, each colour is associated with a size $\delta$ of the central layer and each column is associated with a fluid (VB, DB or HT).

6.2. Stability diagrams

The stability diagram of each fluid is represented in the space $(\delta,\Delta \varPhi )$ in figure 6. Neutral curves are calculated by detecting, for each $\Delta \varPhi$, the value of $\delta$ for which $\omega _i / \alpha$ becomes close to zero – the threshold being set to $10^{-4}$. Calculations are carried out for long waves at $\alpha =10^{-2}$. Theoretical predictions are generally in good agreement with the numerical results for all fluids. For fluid VB, good predictions are observed at low $\delta$, but important discrepancies appear for $\delta >0.2$. Such a mismatch is not observed in fluid DB, for which the neutral curve is still accurately predicted for $\delta =0.5$, corresponding to a configuration in which the central layer occupies half of the flow. As for fluid HT, the prediction is also excellent. Note that, for this fluid, the range of $\Delta \varPhi$ that can be studied is limited by the range of $A_\nu$ (Appendix E), which is itself limited by $\gamma <1$ and the set value of $\varepsilon _T / \delta T_{\mu,\rho }$ (Appendix A). As predicted in (5.16), the neutral curve of all fluids collapses in the limit of $\delta \ll 1$, following $\Delta \varPhi = \delta$. As $\delta$ increases, higher-order terms in $\delta$ destabilise fluid HT: the magnitude of $\Delta \varPhi$ required to produce an instability becomes smaller than $\delta$. Conversely, higher-order terms stabilise fluids VB and DB.

Figure 6.

Neutral curves in the space $(\delta, \Delta \varPhi )$, obtained for fluids VB, DB and HT. Solid lines are numerical calculations. Dashed lines with circles are theoretical predictions, given in (5.10), (5.12) and (5.15).

6.3. The different sources of vorticity production

The generation of disturbance vorticity can be examined from the structure of the unstable modes in the physical space. The wall-normal velocity perturbations are made of a plane progressive wave along $x$ (figure 7a). A peak is observed in the central region of the flow, consistent with the linear increase from zero at the wall predicted in the outer regions in (5.5). The associated pressure field is a plane progressive wave with a phase shift of ${\rm \pi} /2$ with respect to $v$. Note that results are here presented for fluid HT, but nearly identical fields are obtained for fluids VB and DB.

Figure 7.

Wall-normal velocity (a) and pressure (b) of the unstable mode in fluid HT at $\alpha =10^{-2}$, $\delta = 10^{-2}$ and ${A_\nu = -0.04}$, plotted in the physical space $(x,y)$ at an arbitrary time $t$.

In fluid VB, vorticity production only results from the shear term $S_\xi$ in (2.15), given that $\bar {\rho }'=0$. Given its mathematical expression, $S_\xi$ follows the same wave structure as $v$, multiplied by the factor $-\varOmega '$. The vorticity profile $|\varOmega |$ of the base flow, which is equal to $|\varPhi |$ for this fluid, is only non-zero around the central layer. It features strong positive and negative gradients around the lower and upper interfaces, respectively (figure 2d). As a result, the structure of $S_\xi$ contains two out-of-phase waves that are localised around each interface of the central layer, as seen in figure 8(a). The field of total vorticity production, ${S_\xi + B_\xi }$, has the same structure as $S_\xi$ since $B_\xi =0$.

Figure 8.

Terms of the vorticity equation (2.15), plotted in the physical space $(x,y)$ at an arbitrary time $t$, for the unstable modes presented in figure 7. The shear term, the baroclinic term, their sum and the resulting vorticity are shown from top to bottom, in that order. Each column corresponds to one fluid. In each panel, fields are normalised by their maximum absolute value, and colour bars are the same as in figure 7. The central layer of thickness $\delta =10^{-2}$ is shown in (a) by the red dashed line; its location is identical in each figure, albeit not reproduced in order to ease visualisation. The location $y=1/2$ is indicated by the black solid line in (c,f), revealing an offset of the wave below and above this line, respectively.

A similar reasoning can be applied to fluid DB, in which only inertial baroclinic effects are at play given that the base-flow vorticity is constant. The structure of the inertial baroclinic term $B_\xi$ (2.15) is deduced from that of $p$ and the profile of $- \bar {\rho }'/\bar {\rho }^2$. Given that $\bar {\rho } = \varPhi$ for this fluid, the $B_\xi$ field presented in figure 8(e) is found to be similar to that of $S_\xi$ observed for fluid VB (figure 8a). As $S_\xi =0$, it follows that the total vorticity production ${S_\xi + B_\xi }$ also resembles that of fluid VB (figure 8g,h).

As for fluid HT, both $S_\xi$ and $B_\xi$ are active in the generation of disturbance vorticity, and their structure (figure 8c,f) is markedly different from that of two previous fluids. Both feature a peak around the central region. This is again a result of the profiles of $- \varOmega$ and $- \bar {\rho }'/\bar {\rho }^2$. Moreover, $S_\xi$ and $B_\xi$ exhibit a phase difference of ${\rm \pi}$. This is readily understood from the phase difference of ${\rm \pi} /2$ between $v$ and $p$, another phase shift of ${\rm \pi} /2$ being added to $B_\xi$ as it contains $\partial p / \partial x$ (2.15). Despite being out of phase and having a similar structure, the sum of $S_\xi$ and $B_\xi$ is not zero. Instead, $S_\xi +B_\xi$ is composed of two out-of-phase waves around each interface (figure 8i), similar to what was observed for fluids DB and VB. This behaviour can be linked to the profile of $\varPhi$ with the following arguments. Around the central region, the $x$-momentum equation (2.11) can be approximated to ${{\varPhi v \simeq \partial p / \partial x }}$. This results from $\partial u / \partial t + \bar {u} \partial u / \partial x$ being much smaller than $\varPhi v$ in this region, given that the phase velocity of $u$ is here close to $\bar {u}$, and that the growth rate $\omega _i$ is small (this is more evident in the spectral space, where this term is simply ${{\rm i} \alpha (\bar {u} - c) u}$). Under this approximation, the linearised vorticity equation (2.15) can be recast as

(6.1)\begin{equation} \frac { \partial \xi } { \partial t } + \bar{u} \frac { \partial \xi } { \partial x } \simeq{-} \frac{\varPhi'}{\bar{\rho}} v, \end{equation}

where the right-hand side corresponds to $S_\xi +B_\xi$; this shows how the quantity $\varPhi$ encapsulates both shear and inertial baroclinic effects. The profile of $\varPhi '$ and the structure of $v$ ultimately explain the structure of $S_\xi +B_\xi$ in fluid HT (figure 8i). Note that the denominator ${\bar {\rho }}$, attached to $\varPhi '$ in (6.1), modulates the amplitude of the minimum and maximum of $\varPhi '$; it can be seen, in figure 8(i), that larger peaks observed around the upper interface than those around the lower interface, given that $\bar {\rho }$ is smaller in the upper region.

Overall, the structure of the total vorticity production $S_\xi +B_\xi$ is similar for each fluid, regardless of the fluid stratification. The associated vorticity fields are finally displayed in figure 8( j–l). Their structure follows that of $S_\xi +B_\xi$, with an alteration resulting from advection effects (left-hand side of the vorticity equation (2.15)). The final picture is two vorticity waves travelling along each interface, with a phase difference of ${\rm \pi}$ minus a phase shift induced by advection. It can also be noted that these waves are generated around the critical layer $y_c$, defined as $\bar {u}(y_c)= c_\varphi$. In fluids VB and DB, $c_\varphi =1/2$ (§ 5.3.6) which, by symmetry of the base flow, leads to $y_c=1/2$. In fluid HT, the phase speed of the mode presented in figure 8 is $c_\varphi \simeq 0.45$; we have verified that $\bar {u}(1/2)\simeq 0.45$.

The relative amplitudes of the terms plotted in figure 8 are hidden by the normalisation of each field. The maximum absolute value that each field reaches within the physical space $(x,y)$ is used in the normalisation procedure. The relative values are indicated in table 2, noting by $M_S$, $M_B$ and $M_{B+S}$ the maxima reached by $S_\xi$, $B_\xi$ and $(S_\xi +B_\xi )$, respectively. Results are straightforward for fluids VB and DB for which $B_\xi =0$ and $S_\xi =0$, respectively. As for fluid HT, it is shown that the shear and baroclinic effects act with similar strengths. Furthermore, the total source of vorticity $(S_\xi +B_\xi )$ is approximately six times weaker than these effects as it proceeds from the interference of the two waves, cancelling out a large part of their amplitude.

Table 2.

Relative maximum amplitudes of the fields shown in figure 8.

6.4. Interpretation

These results can be interpreted within the wave-interaction theory as reviewed by Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011). In this framework, instabilities are seen as a result of vorticity waves developing along two interfaces that are located close enough so that each wave modifies the velocity field of the other. The modified velocity field further deforms each interface, which yields additional vorticity production. This forms a positive feedback loop in which the two vorticity waves are amplified, constituting an instability. In order to achieve amplification, these waves have to be phase locked, which was indeed observed for the three fluids in figure 8( j–l). The appearance of vorticity waves along each interface has been shown to be driven by the structure of $\varPhi$ (6.1). The physical mechanisms giving rise to these vorticity waves are now further examined.

In fluid VB, the instability arises from an excess of vorticity in a localised layer. This is the essential ingredient of the Kelvin–Helmholtz instability, which is well known and will not be further discussed. The mechanism of this instability has indeed been given following either kinematic (Batchelor Reference Batchelor2000) or dynamic (Charru Reference Charru2011) arguments. Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2011) also examined this instability from a wave-interaction perspective.

In fluid DB, only inertial baroclinic effects produce vorticity disturbance. This originates from misalignments between gradients of density and pressure (Soteriou & Ghoniem Reference Soteriou and Ghoniem1995). More precisely, given the flow assumptions, this misalignment can only be generated via the density stratification of the base flow and the gradient of pressure perturbations in the streamwise direction. Physically, two parcels of fluids at different wall-normal locations, having two different densities (that of the base flow), do not have the same streamwise acceleration when submitted to a streamwise perturbation of pressure (Reinaud et al. Reference Reinaud, Joly and Chassaing1999). This induces a wall-normal gradient of streamwise velocity, i.e. vorticity. In fluid DB, the central region has an excess of density. The above mechanism of vorticity production, therefore, occurs between the lower and central layer, as well as between the upper and central layer. This results in the two vorticity waves that have been observed in figure 8(h).

As for fluid HT, the generation of the two vorticity waves is not as straightforward. A sketch of the mechanism is shown in figure 9(a–c). On the one hand, disturbance vorticity is generated following the inertial baroclinic mechanism described in the previous paragraph. However, contrary to fluid DB, this occurs only between two regions of the flow, the lower and the upper ones, which have different densities (figure 2c) – the central layer playing the role of an interface between them. This idealised representation is illustrated in figure 9(b). Therefore, only one interface is felt by the baroclinic perturbations, instead of two as in fluid DB. This results in only one vorticity wave, which is generated along the central region, as previously observed in figure 8(f). On the other hand, disturbance vorticity is also produced by a shear mechanism, which consists in the wall-normal transport of base-flow vorticity by the perturbations. The base flow can also be divided into two regions of vorticity – smaller and larger magnitudes of the shear rate are indeed observed in the lower and upper part of the flow, respectively (figure 2c). The central layer again plays the role of an interface between these two regions, and only one interface is felt by the shear perturbations, as sketched in figure 9(a). As a result, only one wave is generated by the shear mechanism, as previously observed in figure 8(c). This is a consequence of the plane wave structure of $v$, which, as it takes positive and negative values along $x$, alternatively transports parcels of fluid that contain smaller and larger magnitudes of $\varOmega$ towards the central region. At this point, each mechanism generates one wave that is localised in the central layer. These two waves are out of phase, as discussed in the previous section. Moreover, because of the existing shift between the viscosity and density profiles in fluid HT (figure 2c), a small shift $\Delta y \sim \varepsilon _T / \tau$ also exists between the interface at play in each mechanism (figure 9). As a result, the two central waves are in fact slightly shifted from each other, such that their superposition gives rise to two waves that are localised along each side of the central layer (figure 9c). Ultimately, it is these two phase-locked, interacting vorticity waves that produce the instability in fluid HT, following the interpretation of the wave-interaction theory.

Figure 9.

Sketch illustrating the mechanism generating the two vorticity waves in fluid HT. Red and blue vortices indicate positive and negative values of the vorticity disturbance $\xi$, respectively. These disturbances are generated by shear (a) and inertial baroclinic (b) mechanisms. Solid black lines represent the vorticity (a) and density (b) profiles of the base flow, $\varOmega$ and $\bar {\rho }$, respectively. Note that $\varOmega <0$ and that a shift $\Delta y$ exists between the interfaces of $\varOmega$ and $\bar {\rho }$. Vertical green arrows in (a) represent the transport of $\varOmega$ by the disturbance $v$ in and out of a control volume centred around the interface. Horizontal orange arrows in (b) show the disturbance streamwise velocity of a fluid parcel that is induced by a streamwise gradient of disturbance pressure. The magnitude of this velocity depends on $\bar {\rho }$, which modifies the inertia of the fluid parcel whether it is located in the lower or upper region. The linear superposition of the waves generated in (a) and (b) is shown in (c).

7.1. Base flow

A typical base flow of fluid VdW is shown in figure 10(a). The general behaviour is similar to that observed for fluid HT in figure 2(c), the latter being, indeed, an attempt to model some important features of the former. Two regions of markedly different shear rate can be identified. The density and dynamic viscosity profiles are not as simple as in fluid HT: strong gradients are present in the central region, but properties also exhibit weaker variations away from it. The resulting profile of $|\varPhi |$ figure 10(b) follows that of the kinematic viscosity profile, presented in § 2.1.4. As noted, the extremum of $\nu$ (and therefore $|\varPhi |$) is seemingly not as localised as in the other fluid models. Modelling it with a piecewise view of the $\varPhi$ profile, as proposed in figure 3(d), does not appear as an obvious choice. However, it was also noted that scales associated with the width and amplitude of the minimum of $\nu$ could be introduced for fluid VdW (Appendix B). By using them, we will show that the piecewise linear model of fluid HT can indeed provide a good representation of fluid VdW, allowing some stability features to be predicted, and thereby providing useful elements to interpret the instability. Note that, contrary to the previous fluids, the minimum of $\nu$ is not exactly reached at $y=1/2$ but is shifted upwards. This is linked to the procedure used to set the location of the minimum through the reference temperature at the wall and based on the assumption of a constant temperature gradient (see § 2.2). This assumption is not exact for fluid VdW as thermal conductivity varies, causing the observed shift.

Figure 10.

Base-flow profile of fluid VdW at $\check {p}=1.06$ and $\tau =0.5$. (a) Velocity, density and dynamic viscosity profiles. (b) Normalised density-weighted vorticity profile.

7.2. Numerical stability calculations: comparison with the model of fluid HT

The stability diagram of fluid VdW can be calculated in the space $(\delta, \Delta \varPhi )$ by varying both the reduced pressure $\check {p}$ and the temperature gradient $\tau$. Results are compared with the following theoretical predictions obtained for fluid HT in § 5.3.4. For small $\gamma$, (5.15) leads to the theoretical criterion of instability

(7.1)\begin{equation} \Delta \varPhi_0 > \frac{ \delta }{1 - \delta}. \end{equation}

This equation is of interest for fluid VdW as it does not involve $\gamma$ (since $\Delta \varPhi _0$ only depends on $A_\nu$, as shown in Appendix E), which would have had to be defined for such a fluid. Note that (7.1) has an error of order $\mathcal {O}(\gamma ^2)$ (and not $\mathcal {O}(\gamma )$, as might have been expected), widening the validity of this approximation. A generally good agreement is obtained between the numerical calculations and the theoretical prediction (figure 11). This shows the robustness of the model based on fluid HT with respect to the property variations outside the central layer and to the upward shift of this layer. Furthermore, the definition of $\delta T_\nu$ and $A_\nu$ by the inflection point of $\nu (T)$ (Appendix B) is proved to be relevant, leading to a quantitative prediction of the neutral curve. Note that the prefactor used in the definition of $\delta$ (5.1) is here equal to 1.

Figure 11.

Stability diagram of fluid VdW for $\check {p} \in [1.03;1.095]$ and $\tau \in [0.005;0.045]$. Red circles and blue squares are numerical results and indicate instability and stability, respectively. The black line is the theoretical prediction given in (7.1).

The different terms of the vorticity equation (2.15) for the unstable mode in fluid VdW are now examined. The shear term $S_\xi$ (figure 12a) and the inertial baroclinic term $B_\xi$ (figure 12b) are each composed of a unique wave, which reaches a maximum in the central region (following the discussion of § 7.1, the central region is here defined as being centred around the minimum of kinematic viscosity). This is similar to the observations made for fluid HT in figures 8(c) and 8(f). But because non-zero gradients of $\bar {\mu }$ and $\bar {\rho }$ persist away from the central region in fluid VdW (figure 10a), these terms are not as localised in the centre as in fluid HT. This is particularly noticeable for $B_\xi$, which extends further in the upper region, as $\bar {\rho }'$ is non-zero and $\bar {\rho }$ is much smaller than in the lower region. The sum $S_\xi +B_\xi$ contains two out-of-phase waves (figure 12c) given that as $S_\xi$ and $B_\xi$ are themselves out of phase and are slightly shifted from each other in the wall-normal direction; see the discussion for fluid HT in § 6.4. However, because of the aforementioned asymmetrical structure of $B_\xi$, the upper wave is not localised around the upper interface of the central layer. This constitutes a significant difference with fluid HT (figure 8i). Nevertheless, the final picture is essentially identical: after advection is accounted for, the vorticity field contains two waves localised around each interface, with an additional phase shift leading to a phase difference smaller than ${\rm \pi}$ between them (figure 12d). The role played by advection into the localisation of $S_\xi + B_\xi$ in the central region can be understood from the vorticity equation (2.15), which can be recast, in the spectral space, as

(7.2)\begin{equation} |\xi| = \frac{|S_\xi + B_\xi|}{\alpha | \bar{u} - c |}. \end{equation}

Therefore, $|\xi |$ increases as the phase velocity of the disturbance approaches that of the base flow – which occurs in the central region, as will be shown in the next subsection.

Figure 12.

Terms of the vorticity equation (2.15) for the unstable mode in fluid VdW at $\alpha =10^{-2}$, $\check {p}=1.06$ and $\tau =0.5$. Normalisation, colour bars and annotations are identical to those detailed in figure 8. Note that, in fluid VdW, the central region is defined as being centred around the minimum of kinematic viscosity. (a) Shear production term. (b) Baroclinic production term. (c) Sum of the shear and baroclinic production terms. (d) Vorticity.

Overall, the theoretical stability model developed for fluid HT predicts important features of the stability fluid VdW. This indicates that the main ingredients of the inviscid instability developing in supercritical fluids are indeed contained in fluid HT. These ingredients are the presence of strong, localised gradients of density and viscosity and the associated existence of a localised minimum of viscosity – whose characteristic scales control the instability through the more general parameters $\delta$ and $\Delta \varPhi$. The inviscid instability in supercritical fluids can ultimately be interpreted, like in fluid HT, as a result of the combination of shear and inertial baroclinic mechanisms which generate two interacting waves around the central layer (§ 6.4).

7.3. Growth rate and phase velocity

Having shown the correspondence between the stability of fluids HT and VdW, further numerical results of fluid VdW can now be examined through the lens of the stability model. The growth rate of the unstable mode is shown in figure 13(a) for different supercritical pressures ($\check {p}>1$) and temperature gradients $\tau$. The long-wave nature of the instability is recovered. At constant $\check {p}$, the maximum of the growth rate and the cutoff wavenumber increase with $\tau$. This is consistent with the observations in figure 4(c) for fluid HT, since increasing $\tau$ alone only affects $\delta$, decreasing it. Note that the potential destabilisation of long waves by confinement effects (§ 6.1) is not observed in the present range of parameters. When $\tau$ is fixed, the growth rate is reduced as the supercritical pressure increases. The interpretation is not straightforward since $\check {p}$ affects both $A_\nu$ and $\delta T_\nu$ (Appendix B), and therefore both $\Delta \varPhi$ and $\delta$ in the stability model. Equation (5.17), which can reasonably be invoked given that the magnitudes of $\Delta \varPhi$ and $\delta$ are here of the order of $10^{-2}$, can shed light on this behaviour. Both $\Delta \varPhi$, which can be approximated by $|A_\nu |$ in this regime, and $\delta$, which is proportional to $\delta T_\nu$, increase with $\check {p}$. Therefore, the reduction of the growth rate with $\check {p}$ proceeds from a faster decrease of ${(\Delta \varPhi - \delta )}$ than the increase of $\delta$.

Figure 13.

Growth rate (a) and phase velocity (b) of the unstable mode in fluid VdW. Results in solid lines are obtained for different values of $\tau$ at $\check {p}=1.03$. For $\tau =1$, results at $\check {p}=1.06$ and $\check {p}=1.09$ are also presented by dashed lines. The circles, the square and the triangle in (b) indicate the value of $\bar {u}_m$ (see text) at $\check {p}=1.03$, $\check {p}=1.06$ and $\check {p}=1.09$, respectively (these values are not a function of $\alpha$ and are added to this plot for comparison purposes).

Note that negative growth rates are obtained in figure 13(a) despite carrying out an inviscid stability analysis, in which only neutral modes are usually expected in the absence of an instability. This is a consequence of the use of a complex mapping for $y$ (see § 2.3.1), which is used to remove the singularity of the critical layer for neutral modes. Whilst these negative growth rates do not carry any physical significance, this approach allows us to compute accurate cutoff wavenumbers, presented later in this section.

The phase velocity of the unstable mode is displayed in figure 13(b). It is found to be reasonably constant for all wavenumbers. The velocity of the base flow at the minimum of $\bar {\nu }$, $\bar {u}_{m}$, provides a good prediction of the phase velocity. This is consistent with the results of § 5.3.6, substituting $\bar {u}_{1/2}$ by $\bar {u}_{m}$ because of the aforementioned upward shift of the central layer. This is also consistent with defining the location of the central layer around the minimum of $\bar {\nu }$, as used in the previous subsection. A noticeable departure from $\bar {u}_{m}$ can, however, be observed as $\tau$ increases, i.e. $\delta$ decreases. This behaviour is qualitatively unexpected as the additional shift predicted by the theoretical model should decrease with $\delta$.

The cutoff wavenumber $\alpha _c$ can be plotted from the growth rate calculations by detecting the value of $\alpha \ne 0$ where $\omega _i = 0$. Results at two pressures are presented in figure 14. For each pressure, the cutoff wavenumber is calculated for a range of temperature gradients $\tau$ in order to vary $\delta$. The scaling $\alpha _c \sim \delta ^{-1}$ is found. While this scaling was not derived from the stability models – the long-wave approximation does not allow $\alpha _c$ to be examined, this result is analogous to that of constant-density mixing layers (Charru Reference Charru2011). Note that the scaling is not valid as $\delta$ increases towards the stability threshold. Indeed, as $\delta$ is finite, it would predict finite values of $\alpha _c$ near the threshold, when $\alpha _c$ tends to zero.

Figure 14.

Cutoff wavenumber obtained for fluid VdW at two reduced pressures.

8.1. Summary

The stability of strongly density- and viscosity-stratified plane Couette flow was examined. It was shown that a minimum of kinematic viscosity in the base-flow profile produces a GIP, which furthermore satisfies the generalised Fjørtoft criterion. We first consider three fluid models (termed VB, DB and HT) which were designed to all exhibit a minimum of $\nu$, while having different stratifications of density and dynamic viscosity. Modelling their base flow via piecewise linear profiles, the Rayleigh equation that accounts for strong density gradients was solved for each of them in the long-wave approximation. Expressions of the growth rate and phase velocity were derived, as well as a criterion of instability for each fluid. All these criteria express that a sufficient excess of $|\varPhi |$, the density-weighted vorticity of the base flow, shall be localised in a layer of thickness $\delta$. This improves, for the specific flow studied in this paper, the generalised Fjørtoft criterion – which does not provide such an instability threshold. A general criterion was obtained for all fluids in the limit of small $\delta$, regardless of the type of stratification.

Theoretical predictions were compared with numerical stability calculations. A qualitatively good agreement was found for the growth rate, while an excellent match was observed for the phase velocity. Stability diagrams are also generally well predicted by the models, and the collapse of the neutral curves for small $\delta$ was indeed observed in the calculations. The physical mechanisms responsible for the instability in each fluid were then examined. For all fluids, two phase-locked vorticity waves travelling along each interface of the central layer were found. The growth of the instability is ultimately interpreted as a result of the interaction between these two waves. These waves are generated by either shear effects in fluid VB or inertial baroclinic effects in fluid DB. In fluid HT, which features strong stratifications of both density and viscosity, the two waves were shown to result from a combination of shear and inertial baroclinic mechanisms.

The stability of a fluid governed by the van der Waals equation of state at supercritical pressure was finally examined. The stability model was found to quantitatively predict the neutral curve of this more realistic fluid. The concurrent action of shear and inertial baroclinic effects in the vorticity production was shown to produce two vorticity waves around the central layer. This suggests that the physical interpretation proposed for fluid HT also provides the essential mechanism of the instability in supercritical fluids. Ultimately, our study provides evidence that the minimum of kinematic viscosity, reached at the Widom line, is the leading cause of this instability.

8.2. Supercritical fluids and beyond

We finally address the link between the present paper and the unstable mode recently found in supercritical fluids by Ren et al. (Reference Ren, Marxen and Pecnik2019b), before discussing how this study may be relevant to other types of fluid and flow. It has been shown that, in plane Couette flow, a heated supercritical fluid features a similar instability to that developing in simpler fluids, whose only common property was to assume a minimum of kinematic viscosity. This indicates that other non-ideal thermodynamic effects, such as non-unity compressibility factor $Z$, large values of $c_p$ exhibited near the Widom line or non-monotonic large variations of the speed of sound, are unlikely to play a decisive role in the instability mechanism. In the original study of Ren et al. (Reference Ren, Marxen and Pecnik2019b), the existence of the additional unstable mode was reported for a flat-plate boundary layer flow, considering supercritical $\mathrm {CO}_2$ at non-negligible Mach numbers. We aimed at simplifying their configuration by considering a simpler flow (plane Couette), a simpler supercritical fluid model (VdW equation of state and simple analytical diffusion laws) and by neglecting acoustic phenomena. In this more canonical framework, the instability could be recovered and further analysed. The role played by the minimum of kinematic viscosity remains to be clearly demonstrated in flat-plate boundary layer flows as studied by Ren et al. (Reference Ren, Marxen and Pecnik2019b), since the theoretical development presented in § 3 cannot directly apply to this flow. However, the authors observed the additional unstable mode only when the temperature profile crosses the Widom line, about which the presence of a minimum of $\nu$ is a common property (see figure 15a). This indicates strong links with our results. Bugeat et al. (Reference Bugeat, Boldini and Pecnik2022) furthermore pointed out how, in some cases, inertial effects can be neglected around the pseudo-boiling region, resulting in the correct prediction of a GIP in the boundary layer. In this case, the momentum equation reduces to that of Couette flow, making the aforementioned theoretical result valid for a flat-plate boundary layer flow. However, a rigorous understanding of the conditions of the existence of a GIP is so far missing in this flow, and further efforts could be directed towards the search for a sufficient condition which would factor in inertial effects. The link with an excess layer of $\varPhi$ should also be analysed in this case, as the criterion of instability in (7.1), involving $\Delta \varPhi$ and $\delta$, may be altered for a different flow configuration.

Figure 15.

(a) Kinematic viscosity of different real fluids at supercritical pressure $\check {p}= 1.1$ as a function of temperature. The subscript ‘pb’ refers to the pseudo-boiling point. (b) Kinematic viscosity of a mixture of methanol and toluene as a function of the mole fraction for different temperatures, from McAllister (Reference McAllister1960).

Other types of flow may exhibit a localised minimum of viscosity, hence, potentially, a similar instability. In the present study, the scalar quantity that controls the spatial distribution of $\nu$ is the temperature. Other quantities could play an analogous role. The case of a shear flow made of two miscible fluids, mentioned in the introduction, is an interesting example. The mole fraction is constant in each fluid but varies through the interface, which can be assumed to have a small, finite thickness. Ern et al. (Reference Ern, Charru and Luchini2003) studied the instability developing in this system at low Reynolds numbers. It can be noted that certain fluid mixtures may exhibit a minimum of kinematic viscosity for intermediate mole fraction (between 0 and 1), as shown in figure 15(b). For these fluids, a minimum of $\nu$ would then exist within the diffused interface, as it does for supercritical fluids in the region of the Widom line. An inviscid instability can be expected in this case, which may compete with other instability at low Reynolds numbers if viscous damping does not stabilise it too quickly. As pointed out by Ern et al. (Reference Ern, Charru and Luchini2003), a stability analysis of such a system is valid as long as the time scale associated with the diffusion of this interface is large compared with that of the instability. Similarly, in a single (not necessarily supercritical) fluid flow, heating or cooling a small layer of liquid or gas, respectively, would produce a localised minimum of $\nu$. An instability could then appear, provided that its associated time scale is small compared with that of thermal diffusion of the cooled or heated layer. Note that these time scale considerations were avoided in the present study, as the temperature field strictly was a steady solution of the Navier–Stokes equations. A comparison between the viscous time scale of the perturbations and that of the instability could, however, evaluate the potential of prediction of the inviscid analysis at finite Reynolds numbers.

Finally, the case of non-Newtonian fluids can be mentioned as dynamic viscosity, and hence kinematic viscosity, may vary as a function of the stress profile. In a canonical plane Couette, no variation of stress would be observed, and the flow would be linearly stable. Introducing temperature gradients through viscous heating can alter the stability of the flow (Eldabe, El-Sabbagh & El-Sayed Reference Eldabe, El-Sabbagh and El-Sayed2007). In an isothermal configuration, Nouar & Frigaard (Reference Nouar and Frigaard2009) added a streamwise pressure gradient, leading to a plane Couette–Poiseuille flow. The presence of a maximum (minimum) of stress for a shear-thinning (shear-thickening) fluids may then generate a minimum of kinematic viscosity. However, the criterion of instability derived in the present paper would not apply because of the presence of a pressure gradient, and an isothermal parallel flow of non-Newtonian fluid may therefore not experience this instability.