Non-Oberbeck–Boussinesq effects on the linear stability of a vertical natural convection boundary layer

Abstract

The non-Oberbeck–Boussinesq effects on the stability of a vertical natural convection boundary layer are investigated using the linearised disturbance equations for air flows up to a temperature difference of $\Delta T=100\,{\rm K}$. Based on the linear stability results, the neutral curve is shown to be sensitive to the choice of reference temperature. When evaluated using the film temperature $T_f$, a lower film Grashof number is required to trigger the linear instability for larger $\Delta T$. The relative contributions of shear and buoyant production to the perturbation kinetic energy budget reveals that the marginally unstable modes are amplified based on different mechanisms: for lower wavenumbers at relatively small Grashof number, the instability is driven by buoyancy; whereas for higher wavenumbers and larger Grashof number, the flow becomes unstable due to a shear instability. The use of reference temperature is found to scale the shear- and buoyant-driven instabilities differently so that no single reference temperature definition would collapse the neutral curves. The linear stability result further demonstrates that at a given Grashof number a higher temperature difference would give a larger amplification rate of the perturbation, which then leads to an earlier onset of the nonlinearities when evaluated at $T_f$. Finally, by comparing the amplification rates obtained from direct numerical simulation and the linear stability results, the extent of the linear regime is determined for $\Delta T = 100\,{\rm K}$.

1. Introduction

Natural convection boundary layer (NCBL) flow is ubiquitous in a vast variety of industrial and geophysical applications where the temperature difference between a vertical surface and the ambient fluid induces a spontaneous fluid motion. The majority of past efforts have been dedicated to investigating such convective flows with small temperature differences where the fluid property variations are considered negligible and the buoyancy force is modelled by a linear variation with local temperature difference (e.g. Versteegh & Nieuwstadt 1999; Abedin, Tsuji & Hattori 2009; Nakao, Hattori & Suto 2017). However, many natural convection flows operate at large temperature differences – for instance, the temperature difference between the surface of a solar thermal receiver and the ambient air usually exceeds $650\,{\rm K}$ (Ho & Iverson 2014), in which cases, the local fluid density and properties would vary significantly. In an early study, Siebers, Moffatt & Schwind (1985) experimentally investigated the heat transfer characteristics of the NCBL along a vertical isothermal plate up to a temperature difference of $500\,{\rm K}$ for air. Based on their measurements, they found that the critical Grashof number which marks the onset of the laminar–turbulent transition, when evaluated at the film temperature $T_f = ( T_w+ T_\infty )/2$, decreases with increasing temperature ratio $T_w/T_\infty$ (where $T_w$ is the isothermal wall temperature and $T_\infty$ is the ambient temperature). Such an observation implies that the stability results obtained from the Oberbeck–Boussinesq (OB) limit may not be directly applied to the NCBL with large temperature differences using the film temperature.

While the linear stability for NCBLs under OB conditions has been studied extensively over decades (see, e.g. Xin & Le Quéré 2001; Tumin 2003; Xin & Le Quéré 2012), the non-Oberbeck–Boussinesq (NOB) effects on the stability for NCBL flows have received less attention. Carey & Mollendorf (1978, 1980) analytically demonstrated the applicability of similarity analysis to NCBLs in liquids wherein the viscosity is assumed to vary linearly with temperature while treating all other properties as constants. Although their first-order viscosity variation approximation is only valid when the temperature difference is small (the viscosity variation itself being small enough), their results suggest that the temperature-dependent viscosity would have significant effects on the stability and transition of the NCBL flows by adjusting the mean velocity and temperature profiles. Carey & Mollendorf (1978) further pointed out that the flow similarity exists only when the temperature distribution of the flow is fixed (i.e. a self-similar temperature profile). Sabhapathy & Cheng (1986) numerically investigated the effects of temperature-dependent viscosity and coefficient of thermal expansion on the linear stability properties of NCBLs along a vertical isothermal wall for liquids with Prandtl numbers $7\leq {{Pr}}\leq 10$. In their study, the variation of viscosity is approximated by a linearised Taylor series expansion about the film temperature. Using Orr–Sommerfeld eigenvalue formulation, they found that the variable viscosity stabilises the liquid flow near heated walls and destabilises the flow adjacent to cooled walls, which appears consistent with the findings of Piau (1974); while the variable coefficient of thermal expansion would stabilise the flow near a heated wall locally but destabilise the flow farther downstream. A similar approach was also adopted by Jiin-Yuh & Mollendorf (1988) who numerically studied the linear stability of a vertical plate NCBL in ethylene glycol (${{Pr}}=100$ at film temperature) with a temperature difference of $23\,{\rm K}$. In their analysis, the linear variation of viscosity with temperature is incorporated into the Orr–Sommerfeld equation based on self-similar velocity and temperature base flow profiles. Their numerical results confirmed the findings of Sabhapathy & Cheng (1986), and showed that the glycol flow is more stable when heated than the cooled case with same temperature difference and film temperature. Chenoweth & Paolucci (1985) investigated the vertical NCBL in differentially heated enclosed slots for gases. Using Sutherland's law, they analytically obtained the steady-state exact solution to the laminar flow up to $\epsilon =1$ with both viscosity and thermal conductivity considered temperature-dependent, where $\epsilon = 2(T_h-T_c)/(T_h+T_c)$ is a dimensionless temperature difference given by the differentially heated walls $T_h$ and $T_c$. In their study, they showed that although the variable property only has minimal effect on the wall heat transfer, the temperature and velocity distributions are sensitive to the property variations. Their study was extended by Chenoweth & Paolucci (1986) who numerically investigated the laminar vertical NCBL in cavities. Using a weakly compressible formulation, Chenoweth & Paolucci (1986) showed that flow instability may follow different mechanisms with increasing temperature difference. Suslov & Paolucci (1995b) further investigated NOB effects on the linear stability. In their analysis, the authors demonstrated that the cavity NCBL becomes unstable to a shear instability at relatively low temperature difference and the critical Rayleigh number which marks the onset of flow instability increases with increasing temperature difference (i.e. stabilising the flow) until a critical temperature difference $\epsilon ^*$. For greater temperature differences $\epsilon >\epsilon ^*$, the flow instability is switched to a buoyancy-driven one, resulting in a sudden reduction in the critical Rayleigh number (i.e. destabilising the flow) and the most unstable wavenumber. A similar behaviour is also seen by Suslov & Paolucci (1995a) for mixed convection flows in tall channel. More recently, Rajamanickam, Coenen & Sánchez (2017) studied the NOB effects on the linear stability of an inclined NCBL using a simplified low-Mach-number formulation and the fluid properties are approximated using a power-law. Their numerical results show that the NOB air flow can have travelling waves over a significantly wider inclination angle than those predicted by the Boussinesq limit. Liu et al. (2018) investigated the linear stability in Rayleigh–Bénard convection, where the NOB effects are found to stabilise (destabilise) the flow for large (small) aspect ratios. For natural convection in cavities (including Rayleigh–Bénard systems) and differentially heated channels, the strong influences of the NOB effects on the flow instability are found mainly due to the breakdown of symmetry when heated and cooled (e.g. Ahlers 1980; Suslov & Paolucci 1995a,b; Liu et al. 2018); whereas such a symmetry does not exist for NCBL flows along a vertical plate. Despite the existing efforts, a fundamental understanding of the NOB effects, induced by the bulk temperature difference, on the onset of such flow instability remains unclear.

In this paper, the linear stability of a two-dimensional laminar, periodic vertical NCBL in air is investigated by employing the linearised disturbance equations as an initial value problem. The periodic parallel flow is practically important to those start-up transient flows with large downstream distances where the one-dimensional conductive regime could exist for a significant amount of time before the arrival of the leading edge signal. The existence of such a one-dimensional conductive flow regime has been confirmed both experimentally and numerically (e.g. Siegel 1958; Miyamoto 1977; Gebhart 1988; Sammakia, Gebhart & Qureshi 1982) with analytical solutions made available for the OB case (Illingworth 1950; Schetz & Eichhorn 1962). In particular, Goldstein & Briggs (1964) and Joshi & Gebhart (1987) showed that the flow could become unstable and therefore transition to turbulence before the arrival of the leading edge signal. The linear stability of such parallel transient flow has been investigated under OB conditions for NCBL in cavities (Brooker et al. 2000) and along a vertical flat plate (Joshi & Gebhart 1987; Krane & Gebhart 1993; Ke et al. 2019). However, with increasing characteristic temperature difference, the stability results obtained from the OB studies become largely unreliable in predicting the flow instability and transition in real-world applications. Using a weakly compressible formulation (Batchelor 1953; Paolucci 1982), for the first time, we are able to obtain the marginal stability curve of such a periodic NCBL flow with large temperature differences (up to $100\,{\rm K}$); and provide a detailed insight of the NOB effects on its instability mechanisms.

The rest of paper is organised as follows. An overview and the mathematical formulation for this initial-value problem is given in § 2. Based on the linear stability results, the marginal stability curves for various temperature differences are obtained in § 3.1. The instability dynamics are further explored in § 3.2 by inspecting the energy budgets for two representative marginally unstable modes. The linear stability results are re-normalised using the film temperature in § 3.3 to allow comparisons with observations in the literature. In § 3.4, a direct numerical simulation is carried out to incorporate the unsteady and nonlinear effects – which is then compared with the linear results in § 3.1 to determine the linear stability regime. Section 4 briefly summarises the findings in this paper.

2. Numerical formulation

The problem under consideration is a vertical natural convection flow adjacent to an isothermally heated plate in air. With the use of Squire's theorem (Squire 1933), the transient stability problem is reduced to a two-dimensional one. In the analysis that follow, equations and results are made dimensionless using the intrinsic length scale $L_{s,\infty }={\xi ^*}^{2/3}/(g^*\beta ^*\Delta T)^{1/3}$ and velocity scale $U_{s,\infty }=(\xi ^*g^*\beta ^*\Delta T)^{1/3}$ since the temporally developing parallel flow has no nature length scale. Here, $\xi$ represents the thermal diffusivity, $g$ is the gravitational acceleration, $\Delta T$ is the temperature difference given by the heated wall $T_w$ and the ambient $T_\infty$, and $^*$ indicates the dimensional quantities that are evaluated at the ambient temperature $T_\infty$. Fluid density, coefficient of thermal expansion and specific heats are normalised using the ambient temperature quantities so that $\rho^* = \rho \rho_\infty^*$, $\beta^* = \beta \beta_\infty^*$ and $C_p^* = C_p C_{p\infty}^*$.

The NCBL flow is governed by the conservation equations with low-Mach-number approximation in which the sound waves are filtered while allowing arbitrary density variations (Batchelor 1953; Paolucci 1982). The dimenionless form of these governing equations read

(2.1a)$$\begin{gather} \frac{\partial \rho}{\partial t}+ \frac{\partial \rho u_i}{\partial x_i} = 0, \end{gather}$$

(2.1b)$$\begin{gather}\frac{\partial \rho u_i}{\partial t} + \frac{\partial \rho u_i u_j}{\partial x_j}={-}\frac{\partial p}{\partial x_i} +\frac{\partial \tau_{ij}}{\partial x_j} + \frac{(\rho-1)}{\epsilon}\boldsymbol{n}_i, \end{gather}$$

(2.1c)$$\begin{gather}\rho C_p\left(\frac{\partial \theta}{\partial t} + \frac{\partial u_i\theta }{\partial x_i}\right) = \frac{\partial}{\partial x_i}\left(\kappa\frac{\partial \theta}{\partial x_i}\right)+ \varGamma \frac{\partial \varpi}{\partial t}, \end{gather}$$

(2.1d)$$\begin{gather}\varpi = \rho \left(1+\theta \epsilon\right), \end{gather}$$

where the subscript $i,j=1,2$ denotes the streamwise ($x$) and wall-normal ($y$) directions, and $u_1, u_2 = u, v$ are the streamwise and wall-normal velocity, respectively; $\kappa$ is the thermal conductivity, $\mu$ is the dynamic viscosity, $p$ represents the hydrodynamic pressure; $n_i$ is the unit gravitational vector pointing in the $-x$ direction. The viscous stress tensor $\tau _{ij}$ is given by

(2.2)\begin{equation} \tau_{ij} = \mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right) - \frac{2}{3}\mu\frac{\partial u_k}{\partial x_k}\delta_{ij}, \end{equation}

and $\delta _{ij}$ is the Kronecker delta. The thermal field is made dimensionless using $\theta = (T-T_\infty ) /\Delta T$, whereas the bulk temperature difference is given by $\epsilon =\Delta T/T_\infty$; $\varpi$ represents the spatially uniform thermostatic pressure, and $\varGamma = (\gamma _\infty -1)/\gamma _\infty$ where $\gamma _\infty = C_p/C_v$ is the ratio of the specific heats under constant pressure and volume at $T_\infty$. In the present study, both dynamic viscosity and thermal conductivity of the air flow are assumed to vary with the local temperature $\theta$ following Sutherland's law (Sutherland 1893), which in dimensionless form reads

(2.3a)$$\begin{gather} {\mu} = {{Pr}}_\infty\left(1+\epsilon\theta\right)^{3/2}\frac{1+S_\mu/T_\infty}{1+\epsilon\theta+S_\mu/T_\infty}, \end{gather}$$

(2.3b)$$\begin{gather}{\kappa} = \left(1+\epsilon\theta\right)^{3/2}\frac{1+S_\kappa/T_\infty}{1+\epsilon\theta+S_\kappa/T_\infty}, \end{gather}$$

where ${{Pr}}_\infty \equiv C_p\mu _\infty /\kappa _\infty$ is the Prandtl number evaluated at $T_\infty$, and $S_\mu =111\,{\rm K}$ and $S_\kappa =194\,{\rm K}$ are Sutherland constants for air (Hilsenrath 1955; White 1988). The use of Sutherland's law was shown to be able to accurately predict both viscosity and conductivity for air up to $\epsilon =0.6$, beyond which its accuracy degrades rapidly (Chenoweth & Paolucci 1985). For $T_\infty = 288\,{\rm K}$, this limit corresponds to a maximum $\Delta T = 172\,{\rm K}$. All other properties are considered constant and evaluated at $T_\infty$ (Bergman et al. 2011, e.g. the specific heat $C_p$ varies by ${\sim }1.5\,\%$ for a temperature difference of $\Delta T =172\,{\rm K}$ from $T_\infty =288\,{\rm K}$) – similar assumptions have also been invoked in early studies (e.g. Chenoweth & Paolucci 1985, 1986; Emery & Lee 1999, for NCBL in vertical slots and cavities).

Since the laminar NCBL flow under consideration is perfectly parallel to the isothermal wall and homogeneous in the streamwise direction (see figure 1), the wall-normal velocity and streamwise gradients must vanish $v=\partial /\partial x=0$, with which system (2.1) is simplified to (cf. Suslov & Paolucci 1995a,b)

(2.4a)$$\begin{gather} \frac{\partial U_b}{\partial t} = \frac{1}{\rho_b}\frac{\partial}{\partial y}\left(\mu\frac{\partial U_b}{\partial y}\right) + \frac{\left(1-\rho_b\right)}{\epsilon}, \end{gather}$$

(2.4b)$$\begin{gather}\frac{\partial \theta_b}{\partial t} = \frac{1}{\rho_bC_p}\frac{\partial}{\partial y}\left(\kappa\frac{\partial \theta_b}{\partial y}\right) + \frac{\varGamma }{\rho_bC_p}\frac{\partial \varpi}{\partial t}, \end{gather}$$

(2.4c)$$\begin{gather}\rho_b = \frac{\varpi}{1+\theta_b\epsilon}, \end{gather}$$

(2.4d)$$\begin{gather}\varpi = \left(\int_\mathcal{V} \frac{1}{1+\theta_b\epsilon} \,\mathrm{d} \mathcal{V}\right)^{{-}1}, \end{gather}$$

where $\mathcal {V}$ is the computational volume; $U_b$, $\theta _b$ and $\rho _b$ are the streamwise velocity, temperature and density profiles of the laminar flow. For the purpose of linear stability analysis, an artificial temperature perturbation $\tilde {\theta }$ is superimposed on the laminar base flow ($U_b, \theta _b$) at a designated Grashof number $Gr_\infty$ – at which the base flow is ‘frozen’ (treated as quasi-steady) since the temporal evolution of the base flow is considered much slower than that of the perturbations (Otto 1993; Brooker et al. 2000; Ke et al. 2019). The validity of quasi-steady treatment for the transient one-dimensional periodic NCBL flows has been extensively discussed in the limit of both short and long waves under OB conditions (Daniels & Patterson 1997, 2001), where they showed that the quasi-steady approximation is most accurate for short waves, and the accuracy for the long-wave disturbances will increase with the amplification rate (or, with increase Reynolds number, Tromans 1978; Dwoyer & Hussaini 1987). The temporal evolution of the perturbations is then examined in a similar sense to Ke et al. (2019) as an initial value problem.

Figure 1. Systematic sketch of the periodic NCBL problem, where $u$ represents the streamwise velocity and $\theta$ denotes the temperature distribution of the flow; and $g$ is the gravitational acceleration (not to scale).

Here, the ambient Grashof number $Gr_\infty$ and the integral velocity boundary layer thickness is defined by

(2.5a,b)\begin{equation} Gr_\infty = \frac{{\rho_\infty^*}^2g^* \Delta T {\delta^*}^3}{T_\infty {\mu^*_\infty}^2}= \frac{\rho_\infty^2\delta^3}{\mu_\infty^2}, \quad \delta = \int_0^{\infty} \frac{U_b}{U_{m}} \,\mathrm{d} y, \end{equation}

where $U_m$ is the maximum streamwise velocity of the laminar flow. For a parallel periodic NCBL, both $\delta$ and $Gr_\infty$ are functions of time $t$ only as the flow is invariant in the streamwise ($x$) direction but unsteady in time. Notably, the exact solutions to the boundary layer in (2.4) are analytically difficult to obtain since the flow fields lose their self-similarity due to the temperature-dependent thermal conductivity (Carey & Mollendorf 1978). Although the laminar streamwise velocity profile of the NCBL under different temperature differences do not appear self-similar, the use of the ambient Grashof number (2.5) allows a direct comparison of these flows with the same dimensionless boundary layer thickness $\delta$. In this study, the laminar base flow profiles are numerically obtained using a precursor simulation which solves (2.4) simultaneously with the boundary conditions $U_b=0$, $\theta _b=1$ at $y=0$ and $\partial U_b/\partial y=\theta _b=0$ at $y=\infty$ (see Ke et al. 2019, 2020, for numerical details of the precursor DNS), allowing the flow to develop from quiescence until a designated $Gr_\infty$ is reached. We also note that for larger temperature difference, a shorter development time $t$ is required to reach a given $Gr_\infty$ than that of the smaller $\epsilon$ case, resulting in a lower maximum streamwise velocity $U_m$. For the cases considered in the present study, the maximum streamwise velocities at $Gr_\infty =600$ are $U_m=3.67, 3.62,3.38, 3.14$ for $\epsilon =0$ (Oberbeck–Boussinesq case), $0.035, 0.174, 0.347$. The streamwise velocity magnitude difference is not as obvious as for the NCBL flows in cavities and vertical channels: for example, in an open cavity, Juárez et al. (2011) showed the maximum velocity is decreased by ${\sim }30\,\%$ and the boundary layer thickness $\delta$ is increased by $\sim 70\,\%$ for a $0.3$ increase in the temperature difference $\epsilon$ at a given cavity height Rayleigh number ($Ra_H$) in the steady state; in a fully enclosed cavity, Zhong, Yang & Lloyd (1985) showed that there is a ${\sim }20\,\%$ decrease in the $U_m$ for $\epsilon =0.1$ (although the fluid density is considered constant in their study) when compared with the OB flows. Notably, the steady-state boundary layer thickness decreases for the NCBL in cavities at larger $Ra_H$ since the geometric length scale $H$ and cavity aspect ratio limit the development time to reach steady state for different $\epsilon$ values (Juárez et al. 2011); whereas in our temporally developing flow, the boundary layer thickness increases with the development time. Typical laminar base flow profiles for our periodic NCBL are shown in figure 2 at $Gr_\infty =600$, where the OB results are shown to agree well with the low-temperature-difference case ($\Delta T = 10\,{\rm K}$, $\epsilon =0.035$) indicating the NOB effects on the base flow profiles are negligible at this temperature difference. With increasing $\epsilon$ (equivalently, $\Delta T$ in dimensional arguments), the laminar base flow profiles deviate from the OB approximation: as depicted by figure 2(a), with increasing $\epsilon$, the velocity profile shifts away from the wall (see inset). This trend is consistent with the observations of Zhong et al. (1985) for the natural convection in a square cavity that the variable-fluid-property effect would reduce the velocity in the near-wall region but increase it in the outer region. Temperature profiles, shown in figure 2(b), exhibit a slight difference from the OB case, where the local temperature is increased with increasing $\epsilon$ in the boundary layer (see inset).

Figure 2. Laminar base flow (a) streamwise velocity and (b) temperature profiles for the NCBL obtained by the precursor direct numerical simulation at $Gr_\infty =600$. Insets show the magnified view of grey boxes; black arrows in insets indicate increasing temperature difference $\epsilon$.

The temporal response of the two-dimensional (as per Squire's theorem) artificial perturbations, however, is characterised by the linearised disturbance equations, which can be obtained by substituting $\phi = \phi _b+\tilde {\phi }$ (where $\phi = u,v,\theta,p, \rho, \mu, \kappa$) into (2.1), dropping out the higher-order nonlinear terms and neglecting the base flow variation in time

(2.6a)$$\begin{gather} \frac{\partial \tilde{\rho}}{\partial t}+ \rho_b\frac{\partial \tilde{u}}{\partial x} +U_b\frac{\partial \tilde{\rho}}{\partial x}+\frac{\partial \rho_b \tilde{v}}{\partial y}= 0, \end{gather}$$

(2.6b)$$\begin{gather} \frac{\partial \tilde{u}}{\partial t} + U_b \frac{\partial \tilde{u}}{\partial x} + \tilde{v}\frac{\partial U_b}{\partial y} = \frac{1}{\rho_b}\left[-\frac{\partial \tilde{p}}{\partial x} + {\mu}_b \triangledown^2\tilde{u} + \tilde{\mu} \triangledown^2 U_b + \frac{1}{3}\mu_b \frac{\partial }{\partial x}\left(\frac{\partial \tilde{u}}{\partial x}+\frac{\partial \tilde{v}}{\partial y}\right)\right.\nonumber\\ \left.+\frac{\partial {\mu}_b}{\partial y}\frac{\partial \tilde{v}}{\partial x}+ \frac{\partial {\mu}_b}{\partial y}\frac{\partial \tilde{u}}{\partial y} + \frac{\partial \tilde{\mu}}{\partial y}\frac{\partial U_b}{\partial y}+\frac{\tilde{\rho}}{\epsilon} \right], \end{gather}$$

(2.6c)$$\begin{gather} \frac{\partial \tilde{v}}{\partial t} + U_b\frac{\partial \tilde{v}}{\partial x} = \frac{1}{\rho_b}\left[-\frac{\partial \tilde{p}}{\partial y}+{\mu}_b\triangledown^2 \tilde{v} +\frac{1}{3}\mu_b\frac{\partial}{\partial y}\left(\frac{\partial \tilde{u}}{\partial x}+\frac{\partial \tilde{v}}{\partial y}\right) \right.\nonumber\\ \left.+ \frac{\partial \tilde{\mu} }{\partial x }\frac{\partial U_b}{\partial y}+\frac{4}{3}\frac{\partial {\mu}_b }{\partial y}\frac{\partial \tilde{v}}{\partial y} - \frac{2}{3}\frac{\partial \mu_b}{\partial y} \frac{\partial \tilde{u}}{\partial x}\right], \end{gather}$$

(2.6d)$$\begin{gather}\frac{\partial \tilde{\theta}}{\partial t} + U_b\frac{\partial \tilde{\theta}}{\partial x} + \tilde{v}\frac{\partial \theta_b}{\partial y}= \frac{1}{\rho_b C_p}\left({\kappa}_b \triangledown^2\tilde{\theta} + \tilde{\kappa} \triangledown^2\theta_b + \frac{\partial \tilde{\kappa} }{\partial y}\frac{\partial \theta_b}{\partial y} + \frac{\partial {\kappa}_b}{\partial y} \frac{\partial \tilde{\theta} }{\partial y}\right), \end{gather}$$

where the Laplacian operator $\triangledown ^2=\partial ^2/\partial x^2+\partial ^2/\partial y^2$. The temperature perturbation, which takes the form

(2.7)\begin{equation} \tilde{\theta} = \theta_b A_n\sin\left(2{\rm \pi} n \frac{x}{L_x}\right), \end{equation}

is fed into the base flow at time $t=0$ and interacts linearly with the base flow following (2.6). Here, $n$ is the mode number – a real positive integer indicating the wavelength $L_x/n$ (the dimensionless wavenumber evaluated at ambient temperature $T_\infty$ is therefore $k_\infty =2{\rm \pi} n /L_x$); $A_n=10^{-3}$ is the initial magnitude of the perturbation and $L_x=L_x^*/L_{s,\infty }$ the dimensionless streamwise domain length. Perturbations must vanish at the non-slip wall and decay towards zero in the far-field, giving the perturbation boundary conditions $\tilde {u}=\tilde {v}=\tilde {\theta } = 0$ at both wall ($y=0$) and far-field ($y=\infty$). No velocity perturbation is added at $t=0$ in the current study as the temperature perturbation will directly disturb the velocity field via the buoyancy term in (2.6b). The linearised system (2.6) is solved numerically for a given wavenumber $k_\infty$ (or mode number $n$) with a base flow at $Gr_\infty$ and temperature difference $\epsilon$ as an initial-value problem. In the present study, the NOB effects on the stability properties are explored for $400\leq Gr _\infty \leq 2000$, $4\leq n\leq 32$ ($0.02 \leq k_\infty \leq 0.19$) and $10K\leq \Delta T\leq 100\,{\rm K}$ (equivalently, $0.035\leq \epsilon \leq 0.347$). The ambient temperature is set to $T_\infty =288\,{\rm K}$ with an ambient Prandtl number of ${{Pr}}_\infty = 0.71$ (based on Bergman et al. 2011) for all cases considered. The calculations are carried out in a computational domain of size $L_x^* \times L_y^* = 1035L_{s,\infty } \times 210L_{s,\infty }$, using a fractional step method (Xia et al. 2016; Demou, Frantzis & Grigoriadis 2018) with a collocated finite volume grid ($N_x\times N_y =1024 \times 256$). The domain is made sufficiently large so that both gradients and magnitudes of the flow variables are reduced to their ambient values at the far-field $y=L_y$. A uniform mesh is employed in the streamwise direction while the grids in the wall-normal direction are stretched using a Gamma function with a maximum stretching rate of $1.24\,\%$.

3. Results and discussion

3.1. Linearised disturbance equations

With the laminar base flow at a given ambient Grashof number $Gr_\infty$, the flow is unstable (or stable) to the perturbation $k$ if the magnitude of its temporal response grows (or decreases) with time. Notably, due to the parallel nature of the unsteady periodic flow, perturbations have no streamwise dependence as would have been seen in the spatially developing flows (Huerre & Monkewitz 1990). As a result, the perturbations would become unstable (or stable) simultaneously at the same amplification (or damping) rate in the flow anywhere, in which sense the convective instability (or stability) of the flow is global rather than local (see, e.g. later in figures 6 and 8). In this study, the temporal response to the initial perturbation $k$ is tracked by the time trace of the velocity and temperature signals ($\tilde {u}$ and $\tilde {\theta }$) at an arbitrary location $(x_p,y_p)$ within the boundary layer (see also figure 6 of Ke et al. 2019). Figure 3 shows the typical temporal responses of the perturbation signals for $k_\infty =0.0486$ ($n=8$) and $k_\infty =0.0971$ ($n=16$) at $Gr_\infty =600$ and $\epsilon = 0.174$. As seen in figure 3, at $Gr_\infty =600$, the oscillatory temporal signals quickly adjust their magnitude to a fully periodic state within $1\sim 2$ temporal periods (i.e. receptivity as the flow receives the external perturbation at $t=0$). In this fully periodic state, the envelope of the temporal signals follow an exponential growth (for $k_\infty =0.0486$) or decay (for $k_\infty =0.0971$) with time. The growth rate of the perturbation signals, termed the amplification rate, is then given by

(3.1)\begin{equation} \sigma_k = \frac{1}{A_t}\frac{\partial A_t}{\partial t}, \end{equation}

where $A_t$ is the instantaneous magnitude of the perturbation signal at time $t$ (i.e. the envelope magnitude in figure 3), which is obtained by taking the fast Fourier transform in the streamwise ($x$) direction at $y=y_p$. For a given combination of $k_\infty$ and $Gr_\infty$, a positive $\sigma _k$ indicates that the temporal response is amplified by the laminar base flow while a negative $\sigma _k$ represents a damped signal. The neutral curve, which characterises the minimal Grashof numbers at which the flow is linearly unstable to a given perturbation $k_\infty$ (i.e. $\sigma _k=0$), can be traced by bisectioning the amplification rates in the $k$–$Gr$ plane.

Figure 3. Temporal responses of the temperature signals at $(x_p,y_p)=(781.6,0.25)$ for base flow at ${Gr_\infty =600}$ and $\epsilon =0.174$ ($\Delta T=50\,{\rm K}$) and (a) $k_\infty =0.0486$ ($n=8$); (b) $k_\infty =0.0971$ ($n=16$). Black dotted lines are the envelopes $\pm A_t$ of the signal that follow the exponential growth or decay; red markers in panel (b) highlight two arbitrarily chosen neighbouring peaks which describes the temporal period $\Delta \tau$ of the signal.

Figure 4 compares the neutral curves ($\sigma _k=0$) for $\epsilon =0.035,0.174$ and $0.347$ with the OB case (Ke et al. 2019). While the neutral curves of all the NOB cases show apparent qualitative similarity to the OB case, these curves demonstrate a clear $\epsilon$-dependence: with increasing $\epsilon$, the neutral curve shifts towards higher Grashof numbers at both lower branch ($k_\infty < O(0.1)$) and upper branch ($k_\infty >O(0.1)$). As seen in figure 4, although the critical mode $k_\infty =0.049$ in the lower branch has minimal $\epsilon$-dependency, the critical Grashof number is increased to $Gr_\infty = 473, 536$ and $622$ for $\epsilon = 0.035,0.174$ and $\epsilon =0.347$, respectively, when compared with $Gr_{\infty }=454$ for the OB case. At larger wavenumbers in the upper branch ($k_\infty >O(0.1)$), such a dependency on the temperature difference $\epsilon$ is considerably less than that of the lower branch, where the marginally unstable Grashof number of $k_\infty =0.12$ is increased to $Gr_\infty =948, 959$ and $980$ for $\epsilon = 0.035,0.174$ and $\epsilon =0.347$, respectively. A similar $\epsilon$ dependency of the critical Grashof number has also been reported by Suslov & Paolucci (1995b) for an NCBL in a differentially heated cavity when using a reference temperature obtained by averaging the hot and cold walls. Based on their eigenvalue analysis, Suslov & Paolucci (1995b) found that the critical point (i.e. critical $Gr$ and $k$) for the NCBL in a cavity initially locates in the upper branch for relatively low $\epsilon$ and the critical Grashof number would increase with increasing $\epsilon$ – until $\epsilon$ is sufficiently large ($\epsilon >0.54$) such that the critical wavenumber would then shift to the lower branch. However, in the present study, such a drastic change in the critical wavenumber is not observed – we will show later in § 3.2 that this is due to a different instability mechanism in the cavity.

Figure 4. Neutral curve of the NCBL in air with non-Oberbeck–Boussinesq (NOB) effects under various temperature differences, compared with OB results obtained by Ke et al. (2019); dimensionless wavenumber $k_\infty$ and Grashof number $Gr_\infty$ are based on the fluid properties at ambient temperature $T_\infty$.

Since the perturbation wavelength $\lambda$ is fully determined by the mode number $n$ and the streamwise domain length $L_x$, the temporal response signal in figure 3 allows further examination of the wave speed of the perturbation travelling wave $v_p$,

(3.2)\begin{equation} v_p \equiv \frac{\omega}{k}= \frac{\lambda}{\Delta \tau} = \frac{L_x}{n\Delta \tau}, \end{equation}

where $\omega$ is the oscillation frequency and $\Delta \tau$ is the dimensionless time interval between two neighbouring peaks of the response signal, as shown in figure 3(b). Based on linear theory, the interaction between the perturbation and its base flow is linear so that $\Delta \tau$ is constant for a given choice of $Gr_\infty$, $k_\infty$ and $\epsilon$. Figure 5(a,b) compare the maximum base flow velocity $U_m$ with the wave velocity obtained by (3.2) at $Gr_\infty =600$ and $Gr_{\infty }=1000$, respectively. For all temperature differences and base flows considered, the wave velocity decreases with increasing wavenumber and temperature difference $\epsilon$; and the NOB effects on the phase velocity appear more apparent for higher wavenumbers. Notably, Daniels & Patterson (1997, 2001) showed that for vertical NCBL flows, the ‘less accurate’ long-wave perturbations will always have a phase speed greater than the maximum base flow velocity ($v_p/U_m\approx 2$); whereas the phase velocity for the short-waves are always less than the maximum base flow velocity ($v_p< U_m$). In this sense, perturbations considered in our study are mostly short-wave disturbances, with a small portion of the ‘less accurate’ intermediate–long waves (since $v_p/U_m\approx 1.2 < 2$ in figure 5a,b) whose accuracy is improved with increasing Grashof (Reynolds) number. At $Gr_\infty =600$, the critical mode $k_\infty =0.049$ (cf. figure 4) is found to have a wave speed close to the maximum base flow velocity $U_m$. This observation is consistent with the findings of Armfield & Patterson (1992) in which the authors showed that the travelling waves amplified by the OB base flow have a wave speed close to the maximum base flow velocity; and with the eigenvalue analysis for the OB case (Ke et al. 2019) from which the phase velocity of the critical point is found to be equal to $U_m$. Such a result is seen more clearly in figure 5(c), where the amplification rate $\sigma _k$ is plotted against the wave speed $v_p$. From figure 5(c), only a small band of wavenumbers are amplified by the base flow ($\sigma _k>0$) at $Gr_\infty =600$ – all of which are shown to have a wave speed close to the velocity maximum in the base flow $U_m$, regardless of the temperature differences. Figure 5(c) also shows apparent NOB effects on the amplification rates: the base flow is unstable to $k_\infty =0.049$ ($v_p/U_m=0.98$) and $k_\infty =0.061$ ($v_p/U_m=0.89$) for $\varepsilon \leq 0.174$, but appears stable for $\varepsilon =0.347$ (see inset). At $Gr_{\infty }=1000$, the base flow becomes unstable to a broader band of wavenumbers at the lower phase velocity range at $v_p/U_m\approx 0.55$ (i.e. shorter waves with higher wavenumbers). As shown in figure 5(d), while the critical mode $k_\infty =0.049$ remains the most amplified mode at $Gr_{\infty }=1000$, its phase velocity is decreased to $v_p/U_m=0.88$.

Figure 5. Wave speed of the perturbations $v_p$ against (a,b) wavenumber $k_\infty$ at $Gr_\infty = 600$ and $Gr_{\infty }=1000$; and (c,d) amplification rate $\sigma _k$ at $Gr_{\infty }=600$ and $Gr_{\infty }=1000$. Inset in panel (b) shows the magnified view for the amplified modes; black dotted lines indicates $\sigma _k=0$ and $v_p/U_m=1$.

3.2. Energy budgets for the buoyancy-driven and shear-driven modes

In figure 4, the lower wavenumbers demonstrate a greater sensitivity to the temperature difference $\epsilon$ compared with the upper branch (i.e. larger wavenumbers $k_\infty >O(0.1)$), indicating potential differences in their instability dynamics. In this subsection, we examine the instantaneous perturbation field as well as the perturbation kinetic energy budget of the representative modes in the lower and upper branches to identify and discern the underlying driving mechanisms for these two branches.

Figure 6 depicts the instantaneous contours of the temperature and velocity perturbations of the critical mode $k_\infty =0.049$ at $Gr_\infty =600$, which is representative of the lower branch instabilities. For brevity, only $\epsilon =0.174$ is shown here since all other temperature differences show qualitatively similar structures. From figure 6, it is seen that the velocity fluctuations extend farther out from the wall than the temperature fluctuations, and both instantaneous fluctuation fields have their peaks located beyond the maximum base flow velocity location – indicating that the instability is strongest in the outer part of the boundary layer. There also exists a secondary emerging peak in the velocity perturbation field at $y\approx 16$, which can be seen more clearly in figure 6(c), where the Reynolds normal stress $\overline {\tilde {u}\tilde {u}}$ profile is normalised by its instantaneous maximum (denoted by subscript $_m$). Here, $\bar {\phi }$ denotes the spatial average of flow variable $\phi$ over the homogeneous direction. Although it is difficult to make a direct quantitative comparison due to flow differences, a similar velocity perturbation structure has also been observed in earlier OB stability results for spatially developing natural convection in air (Versteegh & Nieuwstadt 1998; Aberra et al. 2012) and water (Brooker, Patterson & Armfield 1997; Brooker et al. 2000).

Figure 6. (a) Temperature and (b) velocity perturbation contours of the critical mode $k_\infty =0.049$ at ${Gr_\infty =600}$ and $\varepsilon =0.174$. Contour lines are equally spaced ($10\,\%$, $40\,\%$ and $70\,\%$ of the maximum), and the red solid, blue dashed and black dotted contour lines indicate positive, negative and zero levels. (c) Streamwise Reynolds normal stress $\overline {\tilde {u}\tilde {u}}$ profile normalised by its instantaneous maximum $\overline {\tilde {u}\tilde {u}}_m$. Vertical dash-dotted lines indicate the maximum base flow velocity location $\delta _m$ (green).

To further understand the energy contributions to the perturbation pattern shown in figure 6, we inspect the two-dimensional perturbation energy $E=(\tilde {u}^2+\tilde {v}^2)/2$ budget for this critical mode at $Gr_{\infty }=600$. With streamwise homogeneity, the energy budget equation can be obtained by summing (2.6b) multiplied by $\tilde {u}$ and (2.6c) multiplied by $\tilde {v}$. The resulting budget equation for our temporally developing NCBL reads

(3.3)\begin{equation} \frac{\partial {E}}{\partial t} = \mathcal{P}+\mathcal{B} + \mathcal{D}_b+ \varepsilon_b + \varPi + \varSigma, \end{equation}

where $\mathcal {P}$ is the shear production, $\mathcal {B}$ is the buoyancy production; $\mathcal {D}_b$ is the base flow diffusion; $\varepsilon _b$ is the base flow viscous pseudo dissipation; $\varPi$ is the pressure transport; $\varSigma$ represents the effects due to spatial variation in fluid density and properties. These terms read

(3.4a–c)\begin{gather} \mathcal{P} ={-} \overline{\tilde{u}\tilde{v}}\frac{\partial {U}_b}{\partial y}, \quad \mathcal{B} = \frac{\overline{\tilde{u}\tilde{\rho}}}{\rho_b\epsilon}, \quad \mathcal{D}_b = \frac{1}{2}\frac{\mu_b}{\rho_b}\left(\frac{\partial^2 \overline{\tilde{u}^2} }{\partial y^2}+\frac{4}{3}\frac{\partial^2 \overline{\tilde{v}^2} }{\partial y^2}\right), \end{gather}

(3.4d)\begin{gather} \varepsilon_b ={-}\frac{\mu_b}{\rho_b} \left\{ \overline{ \left(\frac{\partial \tilde{u}_i}{\partial x_j}\right)^2 } +\frac{1}{3} \left[ \delta_{ij}\overline{\left(\frac{\partial \tilde{u}_i}{\partial x_i}\right)^2} +\overline{\frac{\partial \tilde{u}}{\partial x}\frac{\partial \tilde{v}}{\partial y}} +\overline{\frac{\partial \tilde{v}}{\partial x}\frac{\partial \tilde{u}}{\partial y}} \right] \right\}, \end{gather}

(3.4e,f)\begin{gather} \varPi = \frac{1}{\rho_b} \left(\overline{\tilde{p}\frac{\partial \tilde{u}_i}{\partial x_i}}-\frac{\partial \overline{\tilde{v}\tilde{p}}}{\partial y} \right), \quad \varSigma = \mathcal{G}_b +\mathcal{G}_\mu +\mathcal{D}_\mu. \end{gather}

Here, $\mathcal {G}_b$ and $\mathcal {G}_\mu$ are the viscosity gradient effects due to the base flow and local instantaneous viscosity fluctuation, respectively; and $\mathcal {D}_\mu$ is diffusion due to viscosity fluctuation:

(3.5a)\begin{gather} \mathcal{G}_b = \frac{1}{2\rho_b}\frac{\partial \mu_b}{\partial y}\left(\frac{\partial \overline{\tilde{u}^2}}{\partial y}+\frac{4}{3}\frac{\partial \overline{\tilde{v}^2}}{\partial y}+2\overline{\tilde{u}\frac{\partial \tilde{v} }{\partial x}}-\frac{4}{3}\overline{\tilde{v}\frac{\partial \tilde{u}}{\partial x}}\right), \end{gather}

(3.5b,c)\begin{gather} \mathcal{G}_\mu = \frac{1}{\rho_b}\left(\overline{\tilde{u}\frac{\partial \tilde{\mu}}{\partial y}}+\overline{\tilde{v}\frac{\partial \tilde{\mu}}{\partial x}}\right)\frac{\partial U_b}{\partial y},\quad \mathcal{D}_\mu = \frac{1}{\rho_b}\overline{\tilde{\mu}\tilde{u}}\frac{\partial^2 U_b}{\partial y^2}. \end{gather}

Figure 7(a) demonstrates the individual contributors to the energy budget (3.3) for the critical mode $k_\infty =0.049$ at $Gr_\infty =600$. In linear analysis, the perturbations are allowed to grow/decay indefinitely, resulting in a time-dependent scale for the energy budget. However, we note that it is the relative contributions of the terms in (3.4) to the energy budget, rather than the scales, that is of crucial importance to the understanding of the instability mechanisms. In the present study, the energy budgets are normalised in a way such that the base flow viscous diffusion $\mathcal {D}_b$ at the wall ($y=0$) is one. As shown in figure 7(a), the buoyancy production $\mathcal {B}$, being the dominant term in driving the instability, correlates well with the peak of the $\overline {\tilde {u}\tilde {u}}$ seen in figure 6(c); whereas the shear production, $\mathcal {P}$, makes only a small contribution at the peak location and acts as a sink for $y>7$. The production terms shown here behave differently to what is observed in fully turbulent flow where the buoyancy production is usually negligible and negative, and the turbulence is predominantly maintained by shear production (for ${{Pr}}<1$, see Janssen & Armfield 1996; Versteegh & Nieuwstadt 1998; Ke et al. 2023). The strong buoyancy production is often seen in the onset of the buoyant instabilities of the unstable atmospheric surface layers (e.g. Leclerc et al. 1990, although the gravitational vector aligns differently to the mean flow direction). The high level buoyancy production in figure 7(a) suggests that the instability of the critical mode $k_\infty =0.049$ at $Gr_\infty =600$ is buoyancy driven. However, the effects of local viscosity gradient and fluctuations $\varSigma$ only have minimal effect on the balance. The contributors to $\varSigma$, given by (3.5), are shown in figure 7(b). It is seen that the viscosity fluctuation gradient effect $\mathcal {G}_\mu$ is negligible; while other terms are two orders of magnitude smaller than those in figure 7(a). The buoyancy production is then balanced by base flow viscous diffusion $\mathcal {D}_{b}$ and pressure transport $\varPi$, redistributing the production to the near-wall region and the outer free-shear layer. The sum of the contributors, giving the temporal growth of the perturbation kinetic energy $E$, is shown in figure 7(c). The net energy gain peaks beyond the velocity maximum location $\delta _m$ with a secondary peak at $y\approx 16$. These energy peaks correlate well with those of the instantaneous $\overline {\tilde {u}\tilde {u}}$ profile in figure 6(c).

Figure 7. (a) Perturbation kinetic energy budget for the critical mode in the lower branch $k_\infty =0.049$ at $Gr_\infty =600$ and $\varepsilon =0.174$; (b) contributors involving spatial variation in fluid properties (note change in scale); (c) rate of change $\partial E/\partial t$. All quantities are normalised in such a way that $\mathcal {D}_b=1$ at $y=0$. Vertical dash-dotted lines indicate the maximum base flow velocity location $\delta _m$ (green).

At $Gr_\infty =1000$ and $k_\infty =0.12$, which is representative of the upper branch instabilities, a fundamentally different perturbation field is seen in figure 8. An additional peak in the velocity contour is apparently observed in the near-wall shear layer ($y\approx 2$), as shown in figure 8(b,c). The strong oscillations in the velocity perturbation correlate well with the critical layer locations (marked in violet in figure 8), where the local base flow velocity $U_b\approx 0.55U_m$ (cf. figure 5d) matches the phase velocity $v_p$ of $k_\infty =0.12$ (see also in, e.g. Mack 1984; Wall & Wilson 1996). The existence of such a near-wall peak gives rise to a stronger shear stress and marks the onset of a shear instability for the wall-bounded NCBL (see also figure 8 of Janssen & Armfield 1996), rather than the buoyant instability as seen in the lower branch. This can be more clearly seen in the perturbation energy budget, depicted by figure 9(a), where the kinetic energy is shown driven predominantly by shear and buoyancy production near the outer critical layer at $y\approx 8$. Unlike the buoyant instability mode shown in figure 7, for $Gr_{\infty }=1000$ and $k_\infty =0.12$, the shear production grows in strength and demonstrates a similar peak magnitude to that of the buoyancy production $\mathcal {B}$. In the region close to the base flow velocity maximum location $\delta _m$ (green dotted lines in figure 9), $\mathcal {B}$ and $\mathcal {P}$ are slightly negative – the near-wall energy is mainly sustained by the pressure transport from the outer critical layer. Near the production peaks (i.e. outer critical layer), both the base flow viscous dissipation and diffusion are relatively small; and the overall production is balanced by the pressure transport, which redistributes the excessive kinetic energy towards the near-wall and outer shear layers. A similar distribution is also seen by Versteegh & Nieuwstadt (1998) for a vertical natural convection in a differentially heated channel, where the near-wall peak of the redistributive pressure term is shown to increase the flow anisotropy near the non-slip wall. The strong pressure transport in the near-wall region (cf. figure 9a) also results in an additional peak for the temporal evolution $\partial E/\partial t$ around the near-wall critical layer, shown in figure 9(c). Figure 9(b) demonstrates the terms involving spatial fluid property variations. Similar to the buoyant instability in figure 7(b), these contributors to $\varSigma$ are at least two orders of magnitude smaller than those of the dominant terms in the kinetic energy budgets. While $\mathcal {G}_\mu$ remains negligible, an extra local maximum is developed for $\mathcal {D}_\mu$ and $\mathcal {G}_b$ close to the critical layer in the near-wall region, which is also suggestive of a change in the energy generation in the near-wall region.

Figure 8. (a) Temperature and (b) velocity perturbation contours of the marginally unstable mode $k_\infty =0.12$ at $Gr_\infty =1000$ and $\varepsilon =0.174$. Contour lines are equally spaced ($10\,\%$, $40\,\%$ and $70\,\%$ of the maximum), and the red solid, blue dashed and black dotted contour lines indicate positive, negative and zero levels. Panel (c) shows the streamwise Reynolds normal stress $\overline {\tilde {u}\tilde {u}}$ profile normalised by its instantaneous maximum $\overline {\tilde {u}\tilde {u}}_m$. Vertical dash-dotted lines indicate the maximum base flow velocity location $\delta _m$ (green) and the locations of the critical layers (violet).

Figure 9. (a) Perturbation kinetic energy budget for the marginally unstable mode $k_\infty =0.12$ at $Gr_\infty =1000$ on the upper branch of the neutral curve ($\varepsilon =0.174$); (b) the contributors involving spatial variation in fluid properties (note change in scale); (c) the rate of change $\partial E/\partial t$. All quantities are normalised in a way such that $\mathcal {D}_b=1$ at $y=0$. Vertical dash-dotted lines indicate the maximum base flow velocity location $\delta _m$ (green); and the locations of the critical layers (violet).

With the kinetic energy budgets, one is able to discern the driving mechanisms for arbitrary modes by comparing the relative contribution (defined by the $\mathcal {R}$) of shear production $\mathcal {P}$ to the buoyancy production $\mathcal {B}$ across the boundary layer,

(3.6)\begin{equation} \mathcal{R} = \frac{\int_{0}^{\infty} \mathcal{P} \,\mathrm{d} y }{\int_{0}^{\infty}\mathcal{B} \,\mathrm{d} y}. \end{equation}

The production ratio $\mathcal {R}$ conveniently distinguishes the flow regimes by the relative strength that drives the instabilities on the neutral curve: $\mathcal {R}>1$ for shear-driven instability and $\mathcal {R}<1$ for buoyancy-driven instability. As shown in figure 10(a), the shear-driven instability occurs at relatively large wavenumbers in the upper branch and it occupies increasingly greater area in the $k$–$Gr$ plane with increasing $Gr_\infty$ as the base flow gains more momentum. The buoyant instability, however, mainly happens at relatively low $Gr_{\infty }$ and small wavenumbers $k_\infty$ in the lower branch. This result is consistent with the findings in a differentially heated cavity (Suslov & Paolucci 1995b). However, we shall point out that the onset of shear and buoyant modes observed by Suslov & Paolucci (1995b) may be due to a different mechanism – Suslov & Paolucci (1995b) observed that the flow is first unstable to the shear-driven modes in the upper branch due to the breakdown of odd symmetry at the mid-span of the differentially heated cavity. For NCBLs in differentially heated cavities and channels, the shear-driven instability is typically stationary (with a phase velocity of $v_p=0$), while the buoyancy-driven modes are oscillatory in the Boussinesq limit (Suslov & Paolucci 1995b; McBain & Armfield 2003). With NOB effects, Suslov & Paolucci (1995b) showed that the resulting shear-driven mode has a negative wave speed $v_p<0$ (i.e. propagating against the base flow) for cavity flows; whereas in the present study, perturbation waves are travelling in the direction of the base flow ($v_p>0$, cf. figure 5). Notably, for our NCBL, there also exists a region at low wavenumbers where the shear production is negative (resulting in $\mathcal {R}<0$). In this region, the instability is driven entirely by buoyancy (see, e.g. the critical mode budget in figure 7a, where $\mathcal {P}$ has a peak value an order of magnitude smaller than that of $\mathcal {B}$ and acts as a sink in the outer shear layer). Figure 10(b) depicts the temperature difference effects on the $\mathcal {R}$ contours. Evidently, with increasing temperature difference $\epsilon$, the $\mathcal {R}$ contours are shifted to higher $Gr_{\infty }$ and larger $k_\infty$, effectively expanding the buoyancy-driven region in the $k$–$Gr$ plane. Such a result is not surprising as the buoyancy effect becomes increasingly strong with greater temperature differences and therefore results in a broadened buoyancy-driven instability region in figure 10. The growing strength of buoyant instability with increasing $\epsilon$ is also seen by Suslov & Paolucci (1995b), where the critical point (minimum $Gr_{\infty }$ at which the flow becomes marginally unstable to any perturbation) changes from shear driven to buoyancy driven when $\epsilon$ is large enough. However, for the temperature differences considered here up to $\epsilon =0.347$, we have not observed such a sudden change of the critical mode since the unsteady NCBL is unstable to the buoyancy-driven mode first and larger $\epsilon$ would only strengthen the buoyant instability – we conjecture that such a change in critical mode is more likely to occur in NCBLs with higher Prandtl number where the flow might become unstable to the upper shear-driven branch first.

Figure 10. (a) Ratio of integrated shear production to the buoyancy production overlaid on the neutral curve for $\epsilon =0.174$. Red dash-dotted line ($\mathcal {R}=1$) divides the shear-driven (upper) and buoyancy-driven (lower) instability regions; blue dotted line indicates $\mathcal {R}=0$, below which shear production is negative $\mathcal {P}<0$. (b) Dependence of $\mathcal {R}=1$ and $\mathcal {R}=0$ contours on the temperature differences. Colours represent the temperature differences: (green) $\epsilon =0.035$, (blue) $\epsilon =0.174$ and (red) $\epsilon =0.347$; dash-dotted lines represent $\mathcal {R}=1$ and dotted lines denote $\mathcal {R}=0$.

3.3. Neutral curves re-normalised using the film temperature

In the context of NCBL literature, flow field characterisation parameters, such as Grashof number $Gr$ and Rayleigh number $Ra$, are often evaluated at the film temperature ${T_f=(T_w+T_\infty )/2}$ for practical purposes. Earlier investigations on the NOB NCBL revealed that greater temperature differences would destabilise the gaseous boundary layer flows with a lower transition Grashof number when evaluated at $T_f$ (e.g. Chenoweth & Paolucci 1986; Sabhapathy & Cheng 1986; Jiin-Yuh & Mollendorf 1988), suggesting an earlier onset of instabilities at larger $\epsilon$ for air. As shown in figure 4(a), our linear stability results show that a higher $Gr_\infty$ is needed to trigger the instability for greater temperature differences $\epsilon$, i.e. increasing $\epsilon$ would stabilise the NCBL. However, these results do not imply that our calculations contradict earlier investigations as the neutral curves in figure 4 are normalised using the ambient temperature. Figure 4(a) can be re-normalised using fluid properties at $T_f$, as shown in figure 11. Here, $k_f$ and $Gr_f$ are the dimensionless wavenumber and the film Grashof number,

(3.7a,b)\begin{equation} k_f = k_\infty\left(\frac{L_{s,\infty}}{L_{s,f}}\right)=k_\infty \left(\frac{\kappa_\infty\rho_f}{\kappa_f\rho_\infty}\right)^{2/3},\quad Gr_f=Gr_\infty\left(\frac{\mu_\infty\rho_f}{\mu_f\rho_\infty}\right)^2. \end{equation}

As depicted by figure 11, while the neutral curves for the $\varepsilon$-range reported still share a similar shape, it is evident that with increasing $\epsilon$, the critical Grashof number decreases – indicating the air NCBL flow, when evaluated by the film temperature $T_f$, is destabilised when a larger temperature difference is applied to the base flow. Unlike figure 4(a), when normalised by $T_f$, the shear-driven instability demonstrates a greater dependence on the temperature difference $\epsilon$. For instance, the marginally unstable Grashof number for $k_f\approx 0.12$ is decreased from $Gr_f=892$ at $\epsilon =0.035$ to $Gr_{f}=564$ at $\epsilon =0.347$; however, the critical Grashof number for the buoyancy-driven instability in the lower branch decreases slightly to $Gr_{f}=355$ at $\epsilon =0.347$ from $Gr_{f}=445$ at $\epsilon =0.035$. Figure 11 indicates when evaluating at $T_f$, the NOB effects due to heating are to destabilise the air flow, prominent for the shear-driven modes but to a lesser extent for the buoyancy-driven modes. These results are consistent with the earlier linear stability investigations for NCBLs along vertical plates, e.g. Sabhapathy & Cheng (1986) and Jiin-Yuh & Mollendorf (1988) reported heating (decreasing viscosity and conductivity) would stabilise the liquid flows ($7\leq {{Pr}}\leq 100$). Note the viscosity and conductivity for liquid flows decrease with increasing temperature, which is different from the gaseous flow described by Sutherland's law.

Figure 11. Neutral curves and $\mathcal {R}=1$, $\mathcal {R}=0$ contours re-normalised using film temperature $T_f= (T_w+T_\infty )/2$. Colours represent the temperature differences: (green) $\epsilon =0.035$, (blue) $\epsilon =0.174$ and (red) $\epsilon =0.347$; dash-dotted lines represent $\mathcal {R}=1$; dotted lines denote $\mathcal {R}=0$.

3.4. Direct stability analysis

With the linear stability results obtained in § 3.1, the extent of the linear regime where a perturbation of choice, $k_\infty$, interacts linearly with the base flow can be identified by comparing the direct numerical solution to the perturbed system (2.1) in three-dimensional form. By imposing the boundary conditions as impermeable isothermal wall at $y=0$, and the quiescent ambient at $y=L_y$,

(3.8a)\begin{gather} u=v=w=0,\quad \theta = \theta_w, \quad \text{at}\ y=0, \end{gather}

(3.8b)\begin{gather} \frac{\partial{u}}{\partial y} =v= \frac{\partial{w}}{\partial y} =\theta = 0, \quad \text{at}\ y=L_y. \end{gather}

With periodic boundary condition applied in the spanwise $z$ direction, the three-dimensional direct numerical simulation (DNS) is carried out in an $L_x\times L_y \times L_z = 1035L_{s,\infty }\times 207 L_{s,\infty } \times 129 L_{s,\infty }$ computational domain with an $N_x \times N_y \times N_z = 1024\times 256 \times 128$ Cartesian grid (stretched in the $y$ direction but uniform in the homogeneous directions $x$ and $z$, with the same stretching strategy as used in the linear stability runs in § 2). The temperature perturbation, given by (2.7), is fed into the flow at $Gr_\infty =250$ with $n=8$ (so that $k_\infty =0.049$ is obtained) and $A_n=10^{-3}$. The ambient temperature is set to $T_\infty =288\,{\rm K}$ and the bulk temperature difference is set to $\Delta T=100\,{\rm K}$ to enable a direct comparison with the linear stability results for $\varepsilon =0.347$. For the DNS, the perturbation field of a flow variable $\tilde {\phi }$ is obtained by subtracting the homogeneous spatial average of a flow variable $\bar {\phi }$ from its instantaneous value $\phi$, i.e. $\tilde {\phi } = \phi -\bar {\phi }$. Figure 12(a) shows the instantaneous and the fluctuation signal of the temperature at $(x_p,y_p,z_p)=(781.6,0.25,64.5)$. Although the temperature fluctuation demonstrates a similar oscillatory behaviour as shown in the linearised solutions (e.g. figure 3), nonlinear developments become visible at $t\approx 120$ (which corresponds to $Gr_\infty \approx 2.2\times 10^4$). This can be clearly seen in figure 12(b), where the instantaneous amplitude of the temperature fluctuation signal $A_t$ is plotted against the Grashof number $Gr_\infty$. From figure 12(b), the temporal signal follows an adjustment immediately after the temperature perturbation is fed into the flow. Similar adjustments are also seen in the linearised solutions (cf. figure 3) and typically end at $t\approx 50$. During this process, the temperature fluctuation is tuned to an appropriate magnitude and phase condition in a way that it does not follow the linear amplification or decay as the boundary layer receives external artificial perturbation (also known as the receptivity process, see Reed, Saric & Arnal 1996; Saric, Reed & Kerschen 2002; Ke et al. 2019). For the unsteady DNS, a development time of $t=50$ for the base flow corresponds to $Gr_{\infty }\approx 7\times 10^3$, as shown in figure 12(b). Beyond this inevitable initial receptivity, the temperature fluctuation is amplified almost linearly in the log–log plot as the mean flow continues to develop – up until the onset of nonlinearity at $Gr_\infty \approx 1.7\times 10^4$.

Figure 12. (a) Temporal signal of the instantaneous temperature $\theta$ (left axis) and the temperature fluctuation $\tilde {\theta }$ (right axis) which is obtained by subtracting the spatial average $\bar {\theta }$ from the instantaneous value $\theta$; (b) instantaneous magnitude of the temperature fluctuation against the ambient Grashof number $Gr_{\infty }$.

This limit of the linear range can be seen in figure 13, where the amplification rate, $\sigma _k$, of the LNS is compared with the DNS results. With increasing $\epsilon$, the NOB LNS solutions are shown to have a slightly lower amplification rate than that of the OB case for $Gr_{\infty }<5\times 10^3$; however, at greater $Gr_{\infty }$, the perturbation is amplified at a much faster rate. This indicates that a larger bulk characteristic temperature difference would result in the infinitesimal perturbations being amplified more quickly as the flow continues to evolve in time, and potentially would lead to an earlier onset of nonlinearities and laminar–turbulent transition. The DNS data, however, show good agreement with the linear stability results after its initial receptivity process at $Gr_{\infty }\approx 7\times 10^3$. Notably, there are instantaneous fluctuations showing up in the DNS amplification data as there have been weakly nonlinear effects due to the interactions between the perturbation field and property/density fluctuations, all of which, being the higher order nonlinear terms, have been dropped in the linear disturbance modelling (2.6). However, this nonlinear effect is considered negligibly weak since it only has minimal impact on the amplification rates obtained by DNS. By comparing the amplification rates obtained from the DNS and the linear stability results, the linear regime for $\varepsilon =0.347$ is identified to be up until $Gr_\infty \approx 1.7\times 10^4$ (or, equivalently $Gr_{f}\approx 9.7\times 10^3$ when evaluated using $T_f$) – beyond this point, the linear stability result diverges from the DNS data since nonlinearity takes over the perturbation amplification mechanism. The onset of nonlinearity, when evaluated at film temperature, is found to be earlier than that of the OB case, whose linear range extends up to $Gr=1.3\times 10^4$ for an initial perturbation amplitude of $A_n=10^{-3}$ (Ke et al. 2019). Such a result is indeed consistent with the observation of Siebers et al. (1985), where they experimentally found that for a fixed $T_\infty$, the film Grashof number $Gr_f$ at which transition occurs decreases with increasing wall temperature $T_w$ (and thus the temperature difference $\Delta T$).

Figure 13. Amplification rate $\sigma _k$ comparison for wavenumber $k_\infty \approx 0.049$. Lines represent the OB case based on the Orr–Sommerfeld eigenvalue analysis of Ke et al. (2019) (black dotted), and linear stability (LNS) results obtained in § 3.1 for $\epsilon =0.035$, $\epsilon =0.174$ and $\epsilon =0.347$. Markers denote the direct stability results for $\varepsilon =0.347$. Vertical lines mark the end of receptivity tuning (left) and the start of nonlinear regime (right).

4. Conclusions

In this paper, the linear stability of a vertical natural convection boundary layer was investigated using both direct numerical simulation and the linearised disturbance equations. The non-Oberbeck–Boussinesq effects are accounted for by using a weakly compressible formulation (Batchelor 1953; Paolucci 1982) with the temperature-dependent fluid viscosity and conductivity modelled by Sutherland's law. The linear stability results show that with increasing temperature difference $\epsilon$, the onset of flow instability occurs at a higher ambient Grashof number $Gr_{\infty }$. However, such a neutral curve is sensitive to the reference temperature used to normalise the flow: when evaluating at the film temperature $T_f = (T_w+T_\infty )/2$, a lower critical film Grashof number $Gr_f$ is found for greater $\epsilon$, suggesting a greater temperature difference destabilises the air flow. This observation at $T_f$ is fully consistent with earlier investigations for liquid NCBL flows (Sabhapathy & Cheng 1986; Jiin-Yuh & Mollendorf 1988).

The instability dynamics of the perturbations are further explored by investigating the energy budgets for the marginally unstable modes at $Gr_\infty =600$, $k_\infty =0.049$ and $Gr_\infty =1000$, $k_\infty =0.12$, representing the lower and upper branches of the neutral curve, respectively. By comparing the relative contributions of the shear and buoyancy production to the kinetic energy, we show that the marginally unstable mode near the critical point in the lower branch ($Gr_\infty =600$, $k_\infty =0.049$) is driven by buoyancy; whereas the instability at the upper branch ($Gr_\infty =1000$, $k_\infty =0.12$) is shear driven. With increasing temperature difference, the buoyancy-driven instability region occupies an increasingly greater region in the stability diagram due to stronger buoyancy; while the shear-driven instability region expands with increasing Grashof number as the base flow gains more momentum. It is also shown that the buoyancy-driven modes are more sensitive to the temperature difference when normalised by $T_\infty$; while the shear-driven modes are more sensitive to the temperature difference when normalised by $T_f$.

Finally, the extent of the linear regime is identified by comparing the linear stability results with a three-dimensional DNS for $\epsilon =0.347$ and $k_\infty =0.049$. Although the amplification rate obtained from DNS shows slight instantaneous fluctuations due to higher order nonlinear interactions between the fluid properties and density variations, it still demonstrates good agreement with the linear stability results after its initial receptivity, suggesting the linear regime extends up to $Gr_\infty \approx 1.7\times 10^4$. When evaluating at the film temperature, this corresponds to $Gr_{f}=9.7\times 10^3$ and the onset of nonlinearities is found to be earlier than that of the OB case (Ke et al. 2019), which is consistent with earlier observations (e.g. Siebers et al. 1985, the film Grashof number $Gr_f$ at which transition occurs decreases with increasing wall temperature $T_w$).