Hostname: page-component-586b7cd67f-dsjbd Total loading time: 0 Render date: 2024-11-29T18:24:47.499Z Has data issue: false hasContentIssue false

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

Published online by Cambridge University Press:  24 July 2024

Junhao Ke*
Affiliation:
School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, Sydney,New South Wales 2006, Australia
S.W. Armfield
Affiliation:
School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, Sydney,New South Wales 2006, Australia
N. Williamson
Affiliation:
School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, Sydney,New South Wales 2006, Australia
*
Email address for correspondence: junhao.ke@sydney.edu.au

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

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2024. Published by Cambridge University Press.

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 Reference Versteegh and Nieuwstadt1999; Abedin, Tsuji & Hattori Reference Abedin, Tsuji and Hattori2009; Nakao, Hattori & Suto Reference Nakao, Hattori and Suto2017). 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 Reference Ho and Iverson2014), in which cases, the local fluid density and properties would vary significantly. In an early study, Siebers, Moffatt & Schwind (Reference Siebers, Moffatt and Schwind1985) 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é Reference Xin and Le Quéré2001; Tumin Reference Tumin2003; Xin & Le Quéré Reference Xin and Le Quéré2012), the non-Oberbeck–Boussinesq (NOB) effects on the stability for NCBL flows have received less attention. Carey & Mollendorf (Reference Carey and Mollendorf1978, Reference Carey and Mollendorf1980) 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 (Reference Carey and Mollendorf1978) 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 (Reference Sabhapathy and Cheng1986) 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 (Reference Piau1974); 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 (Reference Jiin-Yuh and Mollendorf1988) 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 (Reference Sabhapathy and Cheng1986), and showed that the glycol flow is more stable when heated than the cooled case with same temperature difference and film temperature. Chenoweth & Paolucci (Reference Chenoweth and Paolucci1985) 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 (Reference Chenoweth and Paolucci1986) who numerically investigated the laminar vertical NCBL in cavities. Using a weakly compressible formulation, Chenoweth & Paolucci (Reference Chenoweth and Paolucci1986) showed that flow instability may follow different mechanisms with increasing temperature difference. Suslov & Paolucci (Reference Suslov and Paolucci1995b) 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 (Reference Suslov and Paolucci1995a) for mixed convection flows in tall channel. More recently, Rajamanickam, Coenen & Sánchez (Reference Rajamanickam, Coenen and Sánchez2017) 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. (Reference Liu, Xia, Yan, Wan and Sun2018) 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 Reference Ahlers1980; Suslov & Paolucci Reference Suslov and Paolucci1995a,Reference Suslov and Paoluccib; Liu et al. Reference Liu, Xia, Yan, Wan and Sun2018); 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 Reference Siegel1958; Miyamoto Reference Miyamoto1977; Gebhart Reference Gebhart1988; Sammakia, Gebhart & Qureshi Reference Sammakia, Gebhart and Qureshi1982) with analytical solutions made available for the OB case (Illingworth Reference Illingworth1950; Schetz & Eichhorn Reference Schetz and Eichhorn1962). In particular, Goldstein & Briggs (Reference Goldstein and Briggs1964) and Joshi & Gebhart (Reference Joshi and Gebhart1987) 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. Reference Brooker, Patterson, Graham and Schöpf2000) and along a vertical flat plate (Joshi & Gebhart Reference Joshi and Gebhart1987; Krane & Gebhart Reference Krane and Gebhart1993; Ke et al. Reference Ke, Williamson, Armfield, McBain and Norris2019). 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 Reference Batchelor1953; Paolucci Reference Paolucci1982), 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 Reference Squire1933), 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 Reference Batchelor1953; Paolucci Reference Paolucci1982). 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 Reference Sutherland1893), 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 Reference Hilsenrath1955; White Reference White1988). 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 Reference Chenoweth and Paolucci1985). 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. Reference Bergman, Lavine, Incropera and DeWitt2011, 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 Reference Chenoweth and Paolucci1985, Reference Chenoweth and Paolucci1986; Emery & Lee Reference Emery and Lee1999, 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 Reference Suslov and Paolucci1995a,Reference Suslov and Paoluccib)

(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 Reference Otto1993; Brooker et al. Reference Brooker, Patterson, Graham and Schöpf2000; Ke et al. Reference Ke, Williamson, Armfield, McBain and Norris2019). 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 Reference Daniels and Patterson1997, Reference Daniels and Patterson2001), 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 Reference Tromans1978; Dwoyer & Hussaini Reference Dwoyer and Hussaini1987). The temporal evolution of the perturbations is then examined in a similar sense to Ke et al. (Reference Ke, Williamson, Armfield, McBain and Norris2019) 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 Reference Carey and Mollendorf1978). 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. Reference Ke, Williamson, Armfield, McBain and Norris2019, Reference Ke, Williamson, Armfield, Norris and Komiya2020, 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. (Reference Juárez, Hinojosa, Xamán and Tello2011) 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 (Reference Zhong, Yang and Lloyd1985) 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. Reference Juárez, Hinojosa, Xamán and Tello2011); 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. (Reference Zhong, Yang and Lloyd1985) 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. Reference Bergman, Lavine, Incropera and DeWitt2011) 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. Reference Xia, Wan, Liu, Wang and Sun2016; Demou, Frantzis & Grigoriadis Reference Demou, Frantzis and Grigoriadis2018) 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 Reference Huerre and Monkewitz1990). 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. Reference Ke, Williamson, Armfield, McBain and Norris2019). 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. Reference Ke, Williamson, Armfield, McBain and Norris2019). 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 (Reference Suslov and Paolucci1995b) 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 (Reference Suslov and Paolucci1995b) 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. (Reference Ke, Williamson, Armfield, McBain and Norris2019); 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 (Reference Daniels and Patterson1997, Reference Daniels and Patterson2001) 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 (Reference Armfield and Patterson1992) 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. Reference Ke, Williamson, Armfield, McBain and Norris2019) 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 Reference Versteegh and Nieuwstadt1998; Aberra et al. Reference Aberra, Armfield, Behnia and McBain2012) and water (Brooker, Patterson & Armfield Reference Brooker, Patterson and Armfield1997; Brooker et al. Reference Brooker, Patterson, Graham and Schöpf2000).

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.4ac)\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 Reference Janssen and Armfield1996; Versteegh & Nieuwstadt Reference Versteegh and Nieuwstadt1998; Ke et al. Reference Ke, Williamson, Armfield and Komiya2023). 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. Reference Leclerc, Beissner, Shaw, Den Hartog and Neumann1990, 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 Reference Mack1984; Wall & Wilson Reference Wall and Wilson1996). 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 Reference Janssen and Armfield1996), 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 (Reference Versteegh and Nieuwstadt1998) 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 Reference Suslov and Paolucci1995b). However, we shall point out that the onset of shear and buoyant modes observed by Suslov & Paolucci (Reference Suslov and Paolucci1995b) may be due to a different mechanism – Suslov & Paolucci (Reference Suslov and Paolucci1995b) 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 Reference Suslov and Paolucci1995b; McBain & Armfield Reference McBain and Armfield2003). With NOB effects, Suslov & Paolucci (Reference Suslov and Paolucci1995b) 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 (Reference Suslov and Paolucci1995b), 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 Reference Chenoweth and Paolucci1986; Sabhapathy & Cheng Reference Sabhapathy and Cheng1986; Jiin-Yuh & Mollendorf Reference Jiin-Yuh and Mollendorf1988), 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 (Reference Sabhapathy and Cheng1986) and Jiin-Yuh & Mollendorf (Reference Jiin-Yuh and Mollendorf1988) 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 Reference Reed, Saric and Arnal1996; Saric, Reed & Kerschen Reference Saric, Reed and Kerschen2002; Ke et al. Reference Ke, Williamson, Armfield, McBain and Norris2019). 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. Reference Ke, Williamson, Armfield, McBain and Norris2019). Such a result is indeed consistent with the observation of Siebers et al. (Reference Siebers, Moffatt and Schwind1985), 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. (Reference Ke, Williamson, Armfield, McBain and Norris2019) (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 Reference Batchelor1953; Paolucci Reference Paolucci1982) 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 Reference Sabhapathy and Cheng1986; Jiin-Yuh & Mollendorf Reference Jiin-Yuh and Mollendorf1988).

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. Reference Ke, Williamson, Armfield, McBain and Norris2019), which is consistent with earlier observations (e.g. Siebers et al. Reference Siebers, Moffatt and Schwind1985, the film Grashof number $Gr_f$ at which transition occurs decreases with increasing wall temperature $T_w$).

Funding

The authors gratefully acknowledge the support provided by the Australian Research Council under award DP220103209, and the computational resources provided by the Australian National Computational Infrastructure (NCI) and the University of Sydney Informatics Hub.

Declaration of interests

The authors report no conflict of interest.

References

Abedin, M.Z., Tsuji, T. & Hattori, Y. 2009 Direct numerical simulation for a time-developing natural-convection boundary layer along a vertical flat plate. Intl J. Heat Mass Transfer 52 (19–20), 45254534.CrossRefGoogle Scholar
Aberra, T., Armfield, S.W., Behnia, M. & McBain, G.D. 2012 Boundary layer instability of the natural convection flow on a uniformly heated vertical plate. Intl J. Heat Mass Transfer 55 (21–22), 60976108.CrossRefGoogle Scholar
Ahlers, G. 1980 Effect of departures from the Oberbeck–Boussinesq approximation on the heat transport of horizontal convecting fluid layers. J. Fluid Mech. 98 (1), 137148.CrossRefGoogle Scholar
Armfield, S.W. & Patterson, J.C. 1992 Wave properties of natural-convection boundary layers. J. Fluid Mech. 239, 195211.CrossRefGoogle Scholar
Batchelor, G.K. 1953 The conditions for dynamical similarity of motions of a frictionless perfect-gas atmosphere. Q. J. R. Meteorol. Soc. 79 (340), 224235.CrossRefGoogle Scholar
Bergman, T.L., Lavine, A.S., Incropera, F.P. & DeWitt, D.P. 2011 Fundamentals of Heat and Mass Transfer. Wiley.Google Scholar
Brooker, A.M.H., Patterson, J.C. & Armfield, S.W. 1997 Non-parallel linear stability analysis of the vertical boundary layer in a differentially heated cavity. J. Fluid Mech. 352, 265281.CrossRefGoogle Scholar
Brooker, A.M.H., Patterson, J.C., Graham, T. & Schöpf, W. 2000 Convective instability in a time-dependent buoyancy driven boundary layer. Intl J. Heat Mass Transfer 43 (2), 297310.CrossRefGoogle Scholar
Carey, V.P. & Mollendorf, J.C. 1978 Natural convection in liquids with temperature dependent viscosity. In International Heat Transfer Conference Digital Library. Begel House.CrossRefGoogle Scholar
Carey, V.P. & Mollendorf, J.C. 1980 Variable viscosity effects in several natural convection flows. Intl J. Heat Mass Transfer 23 (1), 95109.CrossRefGoogle Scholar
Chenoweth, D.R. & Paolucci, S. 1985 Gas flow in vertical slots with large horizontal temperature differences. Phys. Fluids 28 (8), 23652374.CrossRefGoogle Scholar
Chenoweth, D.R. & Paolucci, S. 1986 Natural convection in an enclosed vertical air layer with large horizontal temperature differences. J. Fluid Mech. 169, 173210.CrossRefGoogle Scholar
Daniels, P.G. & Patterson, J.C. 1997 On the long-wave instability of natural-convection boundary layers. J. Fluid Mech. 335, 5773.CrossRefGoogle Scholar
Daniels, P.G. & Patterson, J.C. 2001 On the short-wave instability of natural convection boundary layers. Proc. R. Soc. Lond. A 457 (2007), 519538.CrossRefGoogle Scholar
Demou, A.D., Frantzis, C. & Grigoriadis, D.G.E. 2018 A numerical methodology for efficient simulations of non-Oberbeck–Boussinesq flows. Intl J. Heat Mass Transfer 125, 11561168.CrossRefGoogle Scholar
Dwoyer, D.L. & Hussaini, M.Y. 1987 Stability of Time Dependent and Spatially Varying Flows. Springer.CrossRefGoogle Scholar
Emery, A.F. & Lee, J.W. 1999 The effects of property variations on natural convection in a square enclosure. Trans. ASME J. Heat Transfer 121, 5762.CrossRefGoogle Scholar
Gebhart, B. 1988 Buoyancy-Induced Flows and Transport. Hemisphere Publishing.Google Scholar
Goldstein, R.J. & Briggs, D.G. 1964 Transient free convection about vertical plates and circular cylinders. Trans. ASME J. Heat Transfer 86 (4), 490500.CrossRefGoogle Scholar
Hilsenrath, J. 1955 Tables of Thermal Properties of Gases: Comprising Tables of Thermodynamic and Transport Properties of Air, Argon, Carbon Dioxide, Carbon Monoxide, Hydrogen, Nitrogen, Oxygen, and Steam, vol. 564. US Department of Commerce, National Bureau of Standards.Google Scholar
Ho, C.K. & Iverson, B.D. 2014 Review of high-temperature central receiver designs for concentrating solar power. Renew. Sustain. Energy Rev. 29, 835846.CrossRefGoogle Scholar
Huerre, P. & Monkewitz, P.A. 1990 Local and global instabilities in spatially developing flows. Annu. Rev. Fluid Mech. 22 (1), 473537.CrossRefGoogle Scholar
Illingworth, C.R. 1950 Unsteady laminar flow of gas near an infinite flat plate. In Math. Proc. Camb. Philos. Soc., vol. 46, pp. 603–613. Cambridge University Press.CrossRefGoogle Scholar
Janssen, R. & Armfield, S.W. 1996 Stability properties of the vertical boundary layers in differentially heated cavities. Intl J. Heat Fluid Flow 17 (6), 547556.CrossRefGoogle Scholar
Jiin-Yuh, J. & Mollendorf, J.C. 1988 The stability of a vertical natural convection boundary layer with temperature dependent viscosity. Intl J. Engng Sci. 26 (1), 112.CrossRefGoogle Scholar
Joshi, Y. & Gebhart, B. 1987 Transition of transient vertical natural-convection flows in water. J. Fluid Mech. 179, 407438.CrossRefGoogle Scholar
Juárez, J.O., Hinojosa, J.F., Xamán, J.P. & Tello, M.P. 2011 Numerical study of natural convection in an open cavity considering temperature-dependent fluid properties. Intl J. Therm. Sci. 50 (11), 21842197.CrossRefGoogle Scholar
Ke, J., Williamson, N., Armfield, S.W. & Komiya, A. 2023 The turbulence development of a vertical natural convection boundary layer. J. Fluid Mech. 964, A24.CrossRefGoogle Scholar
Ke, J., Williamson, N., Armfield, S.W., McBain, G.D. & Norris, S.E. 2019 Stability of a temporally evolving natural convection boundary layer on an isothermal wall. J. Fluid Mech. 877, 11631185.CrossRefGoogle Scholar
Ke, J., Williamson, N., Armfield, S.W., Norris, S.E. & Komiya, A. 2020 Law of the wall for a temporally evolving vertical natural convection boundary layer. J. Fluid Mech. 902, A31.CrossRefGoogle Scholar
Krane, M.J.M. & Gebhart, B. 1993 The hydrodynamic stability of a one-dimensional transient buoyancy-induced flow. Intl J. Heat Mass Transfer 36 (4), 977988.CrossRefGoogle Scholar
Leclerc, M.Y., Beissner, K.C., Shaw, R.H., Den Hartog, G. & Neumann, H.H. 1990 The influence of atmospheric stability on the budgets of the Reynolds stress and turbulent kinetic energy within and above a deciduous forest. J. Appl. Meteorol. 30, 916933.2.0.CO;2>CrossRefGoogle Scholar
Liu, S., Xia, S.-N., Yan, R., Wan, Z.-H. & Sun, D.-J. 2018 Linear and weakly nonlinear analysis of Rayleigh–Bénard convection of perfect gas with non-Oberbeck–Boussinesq effects. J. Fluid Mech. 845, 141169.CrossRefGoogle Scholar
Mack, L.M. 1984 Boundary-layer linear stability theory. Agard Rep. 709 (3), 13.Google Scholar
McBain, G.D. & Armfield, S.W. 2003 Natural convection in a vertical slot: accurate solution of the linear stability equations. ANZIAM J. 45, C92C105.CrossRefGoogle Scholar
Miyamoto, M. 1977 Influence of variable properties upon transient and steady-state free convection. Intl J. Heat Mass Transfer 20, 12581261.CrossRefGoogle Scholar
Nakao, K., Hattori, Y. & Suto, H. 2017 Numerical investigation of a spatially developing turbulent natural convection boundary layer along a vertical heated plate. Intl J. Heat Fluid Flow 63, 128138.CrossRefGoogle Scholar
Otto, S.R. 1993 On the stability of a time dependent boundary layer. NASA Tech. Note, pp. CR–191542. Institute for Computer Applications in Science and Engineering, NASA Langley Research Center.Google Scholar
Paolucci, S. 1982 On the filtering of sound from the Navier–Stokes equations. Tech. Rep., SAND82-8257. Analytical Thermal/Fluid Mechanics Division, Sandia National Laboratories, Livermore.Google Scholar
Piau, J.-M. 1974 Influence des variations des propriétés physiques et de la stratification en convection naturelle. Intl J. Heat Mass Transfer 17 (4), 465476.CrossRefGoogle Scholar
Rajamanickam, P., Coenen, W. & Sánchez, A.L. 2017 Non-Boussinesq stability analysis of natural-convection gaseous flow on inclined hot plates. Intl J. Heat Mass Transfer 109, 949957.CrossRefGoogle Scholar
Reed, H.L., Saric, W.S. & Arnal, D. 1996 Linear stability theory applied to boundary layers. Annu. Rev. Fluid Mech. 28 (1), 389428.CrossRefGoogle Scholar
Sabhapathy, P. & Cheng, K.C. 1986 The effects of temperature-dependent viscosity and coefficient of thermal expansion on the stability of laminar, natural convective flow along an isothermal, vertical surface. Intl J. Heat Mass Transfer 29 (10), 15211529.CrossRefGoogle Scholar
Sammakia, B., Gebhart, B. & Qureshi, Z.H. 1982 Measurements and calculations of transient natural convection in water. Trans. ASME J. Heat Transfer 104 (4), 644648.CrossRefGoogle Scholar
Saric, W.S., Reed, H.L. & Kerschen, E.J. 2002 Boundary-layer receptivity to freestream disturbances. Annu. Rev. Fluid Mech. 34 (1), 291319.CrossRefGoogle Scholar
Schetz, J.A. & Eichhorn, R. 1962 Unsteady natural convection in the vicinity of a doubly infinite vertical plate. Trans. ASME J. Heat Transfer 84 (4), 334338.CrossRefGoogle Scholar
Siebers, D.L., Moffatt, R.F. & Schwind, R.G. 1985 Experimental, variable properties natural convection from a large, vertical, flat surface. Trans. ASME J. Heat Transfer 107 (1), 124132.CrossRefGoogle Scholar
Siegel, R. 1958 Transient free convection from a vertical flat plate. Trans. Am. Soc. Mech. Engrs 80 (2), 347357.CrossRefGoogle Scholar
Squire, H.B. 1933 On the stability for three-dimensional disturbances of viscous fluid flow between parallel walls. Proc. R. Soc. Lond. A 142 (847), 621628.Google Scholar
Suslov, S.A. & Paolucci, S. 1995 a Stability of mixed-convection flow in a tall vertical channel under non-Boussinesq conditions. J. Fluid Mech. 302, 91115.CrossRefGoogle Scholar
Suslov, S.A. & Paolucci, S. 1995 b Stability of natural convection flow in a tall vertical enclosure under non-Boussinesq conditions. Intl J. Heat Mass Transfer 38 (12), 21432157.CrossRefGoogle Scholar
Sutherland, W. 1893 LII. The viscosity of gases and molecular force. Lond. Edinb. Dubl. Phil. Mag. J. Sci. 36 (223), 507531.CrossRefGoogle Scholar
Tromans, P.S. 1978 Stability and transition of periodic pipe flows. PhD thesis, University of Cambridge.Google Scholar
Tumin, A. 2003 The spatial stability of natural convection flow on inclined plates. Trans. ASME J. Fluids Engng 125 (3), 428437.CrossRefGoogle Scholar
Versteegh, T.A.M. & Nieuwstadt, F.T.M. 1998 Turbulent budgets of natural convection in an infinite, differentially heated, vertical channel. Intl J. Heat Fluid Flow 19 (2), 135149.CrossRefGoogle Scholar
Versteegh, T.A.M. & Nieuwstadt, F.T.M. 1999 A direct numerical simulation of natural convection between two infinite vertical differentially heated walls scaling laws and wall functions. Intl J. Heat Mass Transfer 42 (19), 36733693.CrossRefGoogle Scholar
Wall, D.P. & Wilson, S.K. 1996 The linear stability of channel flow of fluid with temperature-dependent viscosity. J. Fluid Mech. 323, 107132.CrossRefGoogle Scholar
White, F.M. 1988 Heat and Mass Transfer. Addison-Wesley.Google Scholar
Xia, S.-N., Wan, Z.-H., Liu, S., Wang, Q. & Sun, D.-J. 2016 Flow reversals in Rayleigh–Bénard convection with non-Oberbeck–Boussinesq effects. J. Fluid Mech. 798, 628642.CrossRefGoogle Scholar
Xin, S. & Le Quéré, P. 2001 Linear stability analyses of natural convection flows in a differentially heated square cavity with conducting horizontal walls. Phys. Fluids 13 (9), 25292542.CrossRefGoogle Scholar
Xin, S. & Le Quéré, P. 2012 Stability of two-dimensional (2D) natural convection flows in air-filled differentially heated cavities: 2D/3D disturbances. Fluid Dyn. Res. 44 (3), 031419.CrossRefGoogle Scholar
Zhong, Z.Y., Yang, K.T. & Lloyd, J.R. 1985 Variable property effects in laminar natural convection in a square enclosure. Trans. ASME J. Heat Transfer 107 (1), 133138.CrossRefGoogle Scholar
Figure 0

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).

Figure 1

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

Figure 2

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 3

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

Figure 4

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

Figure 5

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).

Figure 6

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).

Figure 7

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 8

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).

Figure 9

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

Figure 10

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

Figure 11

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

Figure 12

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).