Turbulent mixed convection in vertical and horizontal channels

Abstract

Turbulent shear flows driven by a combination of a pressure gradient and buoyancy forcing are investigated using direct numerical simulations. Specifically, we consider the set-up of a differentially heated vertical channel subject to a Poiseuille-like horizontal pressure gradient. We explore the response of the system to its three control parameters: the Grashof number $Gr$, the Prandtl number $Pr$, and the Reynolds number $Re$ of the pressure-driven flow. From these input parameters, the relative strength of buoyancy driving to the pressure gradient can be quantified by the Richardson number $Ri=Gr/Re^2$. We compare the response of the mixed vertical convection configuration to that of mixed Rayleigh–Bénard convection, and find a nearly identical behaviour, including an increase in wall friction at higher $Gr$, and a drop in the heat flux relative to natural convection for $Ri=O(1)$. This closely matched response is despite vastly different flow structures in the systems. No large-scale organisation is visible in visualisations of mixed vertical convection – an observation that is confirmed quantitatively by spectral analysis. This analysis, combined with a statistical description of the wall heat flux, highlights how moderate shear suppresses the growth of small-scale plumes and reduces the likelihood of extreme events in the local wall heat flux. Vice versa, starting from a pure shear flow, the addition of thermal driving enhances the drag due to the emission of thermal plumes.

1. Introduction

The transport of heat by turbulent convection is integral to a wide range of natural and engineering applications, from building ventilation to the atmospheric boundary layer and the near-surface ocean. All of these examples can, under the right conditions, be classified as mixed convection. Mixed convection describes the scenario where both buoyancy and shear forces are relevant to the dynamics. This is in contrast to natural convection, where a flow is driven solely by density differences within the fluid, and forced convection, where buoyancy is negligible and the transport of heat is identical to that of a passive scalar. The relative importance of buoyancy compared to the imposed shear is quantified by the Richardson number $Ri$, with the extreme cases $Ri =\infty$ for purely thermal driving and $Ri=0$ for purely shear or pressure driving.

The foundational work on mixed convection by Obukhov (1946) was motivated by understanding the dynamics of the surface layer of the atmosphere. Obukhov supposed that the dynamics was determined solely by the surface friction (quantified by the friction velocity $u_\tau$), the surface heat flux $q$, and gravity $g$, such that dimensional analysis revealed a single possible length scale $L_O = u_\tau ^3 / g\alpha q$ that could describe the system. Using this length to rescale the problem, Monin & Obukhov (1954) derived what is now known as the Monin–Obukhov similarity theory (MOST), where ‘universal’ functions of $z/L_O$ are used to describe the mean velocity and temperature profiles in stably or unstably stratified shear flows. These universal functions are obtained by interpolating between the extreme cases of natural convection and forced convection, which were updated by Kader & Yaglom (1990) for unstable (i.e. convecting) boundary layers. A historical overview of MOST is provided in Foken (2006), and the theory has been extremely popular in atmospheric and oceanic applications. However, some of the assumptions underlying MOST have recently been coming under further scrutiny, particularly the power-law dependence of the mean profiles in the convective regime (Cheng et al. 2021).

Mixed convection has also often been studied in simple, canonical flow configurations where the system response is dependent only on a small number of dimensionless input parameters. A popular approach has been to introduce horizontal forcing into the classical Rayleigh–Bénard (RB) set-up, either through a horizontal pressure gradient (Poiseuille–Rayleigh–Bénard, P-RB) or by setting one of the boundary plates in motion (Couette–Rayleigh–Bénard, C-RB). The linear stability of the P-RB system was studied by Gage & Reid (1968), who found that streamwise perturbations are suppressed by the introduction of shear, so the fastest growing mode takes the form of convective rolls in the plane orthogonal to the mean flow. Note that the critical Rayleigh number $Ra_c=1708$ and the fastest growing wavelength $\lambda =2\sqrt {2}\,H$ do not change compared to natural RB since the linear problem is unchanged in the orthogonal plane. Such streamwise-aligned rolls were observed experimentally by Richter & Parsons (1975) in C-RB, although since their set-up was motivated by mantle convection, their working fluid had a very high Prandtl number $Pr=8600$ and low Reynolds numbers. Domaradzki & Metcalfe (1988) performed simulations of C-RB, also finding organisation into streamwise-aligned rolls but at the largest wavelength of the domain. More recent direct numerical simulations in larger domains by Pirozzoli et al. (2017) for P-RB and Blass et al. (2020,  2021) for C-RB highlight how these large-scale structures contribute a large proportion of the heat and momentum flux in mixed RB, and how their wavelength depends on the Richardson number $Ri$ of the system. Madhusudanan et al. (2022) recently reproduced the wide rolls using a linear model coupled to eddy diffusivities, showing that they are generated primarily through a classical lift-up mechanism.

The response of canonical mixed convection systems can be quantified using the friction coefficient $C_f$ and the Nusselt number $Nu$, which measure the dimensionless skin friction and heat flux. In forced convection, where buoyancy is negligible, both Poiseuille and Couette flows exhibit an identical response in $C_f$ when appropriately scaled using the centreline velocity (Orlandi, Bernardini & Pirozzoli 2015). Scagliarini et al. (2015) observed an increase in the streamwise friction coefficient in P-RB relative to pure Poiseuille flow, for which they proposed a modified formulation of the log law for the mean velocity in the presence of convection. An intriguing phenomenon of mixed RB is found in the response of the heat flux, which varies non-monotonically with Reynolds number $Re$ for fixed Rayleigh number and Prandtl number. Nusselt number $Nu$ first decreases relative to the natural convection case, before increasing at high $Re$ as the flow enters the forced convection regime (Blass et al. 2020,  2021). This behaviour, not predicted by MOST, has been attributed to the sweeping away of thermal plumes by the imposed horizontal flow. The plume sweeping concept has since been applied to form phenomenological models of the system (Scagliarini, Gylfason & Toschi 2014). Similar to the response of the friction coefficient, an identical response is found in P-RB and C-RB when appropriately scaled, and the decrease in $Nu$ has recently been shown to collapse onto a single curve when the Reynolds number of the shear flow is considered relative to the Reynolds number of the natural convection (Yerragolam et al. 2024). Yerragolam et al. (2024) also provide a theoretical estimate for this decrease in heat flux based on an extension of the Grossmann & Lohse (2000,  2001) theory for RB convection to mixed RB.

The interplay of shear and convection has an important role in another canonical natural convection problem: the differentially heated vertical channel, often simply referred to as vertical convection (VC). In this configuration, convection drives flow parallel to the boundary plates, generating a mean shear at the walls and in the bulk of the channel. An analogy can be drawn between the large-scale circulation in RB and the vertical mean flow in VC, but since the vertical buoyancy flux is not equivalent to the heat flux of interest in VC, the Grossmann & Lohse (2000) approach of linking heat flux and kinetic energy dissipation cannot be applied directly. Nevertheless, Ng et al. (2015) found similar scaling relations to RB for heat flux and dissipation rates in VC when conditionally sampling either the boundary layers or the bulk. Recent simulations at varying Prandtl number (Howland et al. 2022) have prompted renewed efforts to understand the boundary layer theory limiting the global response of the system (Ching 2023) and the dynamics setting the mean profiles in the channel centre (Li et al. 2023).

Mixed convection in a vertical channel, where an additional pressure-driven forcing is applied to the VC configuration, has been less well studied than mixed RB. The majority of studies into these flows (e.g. Kasagi & Nishimura 1997; Wetzel & Wagner 2019; Guo & Prasser 2022) impose a mean pressure gradient in the vertical direction that directly opposes or enhances the mean flow due to convection. Although this configuration is relevant to some industrial applications, from a physical perspective it breaks the symmetry of the channel, with the boundary layers at each wall subject to different shear stresses. In this study, we instead impose a horizontal pressure gradient in the channel, which leads to symmetric profiles of horizontal velocity and higher-order statistics while retaining the anti-symmetric profiles of mean vertical velocity and temperature from VC. Although we approach this configuration from a fundamental point of view, the crossflow set-up can be relevant to industrial heat exchangers in a wide range of applications. We are aware of only one other paper discussing such a system (El-Samni, Yoon & Chun 2005), which highlights tilted structures at the wall and a modification of the near-wall Reynolds stresses. However, the results of El-Samni et al. (2005) are mainly descriptive and cover a limited parameter range.

In the current paper, we use direct numerical simulations to explore the dynamics of turbulent mixed convection in a vertical channel for a wide range of parameters, focusing on the transition between natural convection and forced convection. The paper is organised as follows. First, in § 2, we describe the problem set-up and details of the numerical simulations, before presenting visualisations of the resulting flow in § 3. The global response of the system is investigated in terms of the friction coefficients and the Nusselt number, and compared with the mixed RB flow in § 4. We then turn to wall-normal profiles in mixed VC in § 5, focusing primarily on the balances in the mean momentum budgets. Detailed analysis of the heat transport is then performed by spectral analysis in § 6, and through statistics of the boundary layers in § 7. Finally, our conclusions are presented in § 8, where we discuss the implications of our results and provide an outlook on future research in mixed convection.

2. Simulation set-up and numerical methods

We perform numerical simulations of the flow arising in a fluid confined between two parallel no-slip vertical walls. The walls are separated by a horizontal distance $H$ and are held at fixed temperatures, with a temperature difference $\varDelta$ between the plates. As in the schematic in figure 1(b), we take the $x$-coordinate to be horizontal and parallel to the plates, the $y$-coordinate to be normal to the boundaries, and $z$ to be in the vertical direction. We consider a fluid satisfying the Oberbeck–Boussinesq approximation, such that changes in density are significant only in the buoyancy term, and are linear with respect to changes in temperature. We therefore treat the velocity field $\boldsymbol {u}=(u,v,w)$ as divergence-free ($\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = 0$) and satisfying the Navier–Stokes momentum equation

(2.1)\begin{equation} \partial_t \boldsymbol{u} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u} ={-}\rho^{{-}1}\,\boldsymbol{\nabla} p + \nu\,\nabla^2 \boldsymbol{u} + g \alpha \theta \boldsymbol{\hat{z}} + G \boldsymbol{\hat{x}}, \end{equation}

where $\rho$ is the mean fluid density (assumed constant), $p$ is the pressure, $\nu$ is the kinematic viscosity, $g$ is the acceleration due to gravity, and $\alpha$ is the thermal expansion coefficient. A time-dependent, spatially uniform forcing $G(t)$ is applied in the streamwise ($x$) direction to maintain a constant mean flow $\langle u \rangle = U$. The magnitude of this forcing is computed at every time step to exactly cancel out any variation in the mean flow. Previous work has shown that such a forcing produces results near-identical to those of a constant pressure gradient (Quadrio, Frohnapfel & Hasegawa 2016), but allows us to use the mean flow strength as an input parameter. The temperature field $\theta$ satisfies the advection–diffusion equation

(2.2)\begin{equation} \partial_t \theta + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \theta = \kappa\,\nabla^2 \theta, \end{equation}

where $\kappa$ is the thermal diffusivity. Periodic boundary conditions are applied to both $\boldsymbol {u}$ and $\theta$ in the $x$ and $z$ directions. Unless otherwise stated, the aspect ratio of the domain in these periodic directions is taken as $\varGamma =L/H=8$, such that the length of the domain $L$ in $x$ and $z$ is equal, and eight times the plate separation distance $H$.

Figure 1. Schematics of canonical mixed convection systems. (a) Mixed RB convection with Poiseuille-type forcing as studied by e.g. Domaradzki & Metcalfe (1988), Pirozzoli et al. (2017) and Yerragolam et al. (2024). (b) Mixed VC with Poiseuille-type forcing, as used for the new simulations presented in this paper. Large arrows highlight the imposed horizontal pressure gradient, and smaller arrows highlight the large-scale flows driven by buoyancy in each configuration. The acceleration due to gravity is highlighted by a double arrow.

We perform direct numerical simulations of (2.1) and (2.2) using a highly parallelised flow solver that computes spatial derivatives using second-order central finite differences on a staggered grid configuration. The wall-normal diffusive terms are evolved in time using a Crank–Nicolson scheme, and all other terms are treated explicitly using a three-stage Runge–Kutta method. An adaptive time step is chosen using the constraint of a maximum Courant number of 1. The velocity is kept divergence-free to machine precision using a pressure correction step at each time step that is implemented with fast Fourier transforms in the periodic directions and a tridiagonal matrix solver in the wall-normal direction. A multiple-resolution technique is applied to evolve the velocity and temperature fields on independent grids, with cubic Hermite interpolation used for the buoyancy forcing and the advection of temperature. Detailed overviews of the numerical discretisation, the domain decomposition strategy and the multiple-resolution technique can be found in Verzicco & Orlandi (1996), van der Poel et al. (2015) and Ostilla-Monico et al. (2015), as well as in our software documentation.

The physical input parameters of the system are the Rayleigh number, the Prandtl number and the Reynolds number:

(2.3a–c)\begin{equation} \textit{Ra}= \frac{g\alpha\,\Delta H^3}{\nu \kappa},\quad \textit{Pr}= \frac{\nu}{\kappa},\quad \textit{Re}=\frac{U H}{\nu}. \end{equation}

When considering the strength of the flow driven by convection, it is often useful to consider the Grashof number $Gr$ as the relevant input parameter, and when comparing the relative strengths of buoyancy to pressure driving, we can construct a Richardson number as

(2.4a,b)\begin{equation} \textit{Gr} = \frac{\textit{Ra}}{\textit{Pr}} = \frac{g \alpha\,\Delta H^3}{\nu^2} = \left(\frac{U_f H}{\nu}\right)^2,\quad \textit{Ri} = \frac{\textit{Gr}}{\textit{Re}^2} = \frac{g \alpha\,\Delta H}{U^2} = \left(\frac{U_f}{U}\right)^2 . \end{equation}

These can both be considered as input parameters. Above, we have written $U_f = \sqrt {g\alpha \, \Delta H}$ as the free-fall velocity scale to give insight into the interpretation of these parameters.

In this study, we perform two sets of simulations to compare the relative impacts of the various input parameters. First, we fix $Gr=10^6$ and vary $1\leq Pr \leq 10$ along with $10^{2.5} \leq Re \leq 10^4$, which correspond to Richardson numbers $10^{-2} \leq Ri \leq 10$. For the second set, we fix $Pr=1$ and increase $Gr$ up to $10^8$ while again varying $Re$ up to $10^4$. A detailed overview of the parameters used for each simulation is given in table 1 of Appendix B. Each simulation is run to a statistically steady state, in which the flow statistics are computed and averaged for a minimum of 200 advective time units. For high $Ri\geq 1$, the relevant time unit is $H/U_f$, whereas for low $Ri\leq 1$, the relevant time unit is $H/U$. The wall-normal grid spacing is stretched following Pirozzoli & Orlandi (2021) to ensure sufficient resolution close to the wall, such that $\Delta y^+ < 0.1$ at the wall for the base velocity grid. In the periodic directions, the uniform grid spacing satisfies $\Delta x^+ \leq 5.4$ in every simulation. At the centre of the domain, the spacing of the refined grid satisfies $\Delta y_r \leq 1.05 l_B$, $\Delta x_r < 1.4 l_B$, where $l_B$ is the Batchelor scale computed using the plane-averaged turbulent kinetic energy (TKE) dissipation rate over the midplane. Full details of the grid sizes are provided in Appendix B.

3. Flow visualisation

We begin with a qualitative comparison of the simulations through visualisations of the temperature field and the local heat flux. Figure 2 displays the instantaneous local wall-normal heat flux at the boundary plate $y=0$ for cases with fixed $Pr=1$, and ranges $10^6\leq Gr\leq 10^8$, $10^3\leq Re\leq 10^4$. The dimensionless heat flux plotted here is defined as

(3.1)\begin{equation} q_\theta(x,z,t) ={-}\frac{H}{\varDelta} \left.\frac{\partial \theta}{\partial y}\right|_{y=0}, \end{equation}

such that its time- and plane-averaged value is equivalent to the Nusselt number, $Nu=\overline {q_\theta }$.

Figure 2. Vertical planes of the instantaneous local wall-normal heat flux at the boundary $y=0$. All simulations shown are with fixed $Pr=1$. The Richardson number $Ri=Gr/Re^2$ is comparable along diagonals from the upper left to the lower right, with the largest value ($Ri=100$) in the top right and the lowest value ($Ri=0.01$) in the bottom left. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-2.ipynb.

As mentioned above, the relative strength of convection to the horizontal flow can be characterised by the Richardson number $Ri=Gr/Re^2$, which is constant along diagonals in figure 2. At high $Ri$, as in figure 2(c), the horizontal flow has little impact on the local distribution of the wall heat flux. The near-wall temperature structures are the same as in the case of no crossflow, with regions of large local heat flux (dark spots) separated by longer, streaky structures aligned in the vertical (Howland et al. 2022). As $Ri$ decreases, such as in figures 2(e,i) where $Ri=1$, the prominence of the large heat flux regions diminishes, and the streaks become aligned in the diagonal. This visualises the local direction of the flow along the wall, which at $Ri=1$ is due to a combination of the VC and the horizontal pressure gradient. At lower $Ri$, these structures become more aligned with the horizontal, eventually spanning the domain as in figure 2(g), which is reminiscent of classical low-speed streaks in turbulent channel flow (Kline et al. 1967; Antonia, Abe & Kawamura 2009). A more quantitative analysis of the change in the near-wall heat flux distribution will be provided in § 7.

In figure 3, we present visualisations of the same simulations but now at the midplane of the simulation domain $y=H/2$. As would be expected, the fields at higher $Gr$ and $Re$ exhibit structures with a wider range of spatial scales. Aside from this dynamical range, the effect of increasing $Re$ is less noticeable at the midplane than at the wall in figure 2. At the centre of a turbulent channel, the mean profile of horizontal velocity is relatively flat, with zero mean shear and a local minimum in TKE. By contrast in VC, a mean shear in the vertical velocity drives the generation of turbulence in the bulk, argued by Li et al. (2023) to follow a mixing-layer-like behaviour. In figure 3, greater mixing at higher $Gr$ leads to smaller values of the temperature perturbations in the midplane, although the same effect is not evident as $Re$ increases. Compared to the mixed RB system, where gravity is orthogonal to the plates and the temperature field organises into large, streamwise-aligned coherent rolls, the fields in mixed VC appear rather featureless. A hint of such large-scale rolls is noticeable only for the cases dominated by strong pressure driving with high $Re$ and low $Ri$ (figures 3g,h), where the temperature structures appear more aligned along the streamwise ($x$) axis. This contrast in the organisation of mixed convection systems will be investigated more quantitatively in § 6 through spectral analysis, but we first turn to the global responses of the two mixed convection systems in the next section.

Figure 3. Vertical planes of the temperature field at the channel centre $y=H/2$. See the caption of figure 2 for details of the simulations presented. Here, $\theta _0$ is the arbitrary reference temperature that is the midpoint of the two boundary values. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-3.ipynb.

4. Global response quantities in mixed convection systems

In this section, we compare the responses of the mixed VC set-up with those of mixed RB using the simulation data of Yerragolam et al. (2024). Those simulations cover comparable ranges of parameters $10^6\leq Gr \leq 10^8$, $0.5\leq Pr \leq 5$ and $Re\leq 10^4$ in a domain of streamwise aspect ratio $\varGamma _x=8$ and spanwise aspect ratio $\varGamma _y=4$. The flow solver shares an identical code base except for the multiple-resolution technique, which is used only in the newly reported mixed VC simulations.

4.1. Friction coefficients

A key global response parameter in shear-driven flows is the friction coefficient

(4.1)\begin{equation} C_f = \frac{2 \tau}{\rho U^2}, \end{equation}

where $\tau$ is the wall shear stress and $U$ is the velocity magnitude in the bulk. Since the friction coefficient is determined solely by the velocity field, we expect the dependence on the Prandtl number to be relatively weak, and consider a relationship $C_f(Re)$. Power-law scalings for $C_f(Re)$ can be derived for laminar boundary layer flows, with e.g. $C_f\sim Re^{-1}$ applicable to Couette or Poiseuille flows, and $C_f \sim Re^{-1/2}$ arising from the classical Blasius boundary layer solution (Schlichting & Gersten 2016). For a turbulent boundary layer in the sense of Prandtl and von Kármán, the friction coefficient satisfies the relation

(4.2)\begin{equation} \sqrt{\frac{2}{C_f}} = \frac{1}{\kappa_u} \log{\left( \textit{Re}\,\sqrt{\frac{C_f}{8}} \right) } + B, \end{equation}

known as the Prandtl friction law after Prandtl (1932). The von Kármán constant $\kappa _u$ typically takes a value of approximately 0.4, and the intercept $B$ is close to 4, but the exact values, their universality and the way in which they are fit to data remain an active topic of research (Monkewitz & Nagib 2023). Due to our similar set-up and numerical methods, we take the values suggested by Pirozzoli, Bernardini & Orlandi (2014), namely $\kappa _u = 0.41$ and $B=5$.

Since the mixed VC flow is driven in orthogonal directions by the pressure gradient and by buoyancy, we can construct separate friction coefficients for each component of the wall shear stress. Understanding the response of the friction coefficients in this context relies on choosing an appropriate Reynolds number for each component of the flow. For the streamwise ($x$) direction in which the mean flow is imposed by a pressure gradient, this Reynolds number is simply the input parameter defined in (2.3a–c).

We consider the response of the streamwise friction coefficient in figure 4, where only the streamwise component of the shear stress ${\tau =\rho \nu \,\partial _y \bar {u}}$ is applied to the definition (4.1). The global response of mixed VC is near identical to that of mixed RB, with a transition from a laminar power-law scaling to the Prandtl friction law of (4.2). For comparable parameter values, the largest difference in $C_f$ between the mixed RB and mixed VC cases is 16 %. In the laminar scaling regime, stronger buoyancy driving (characterised by higher $Gr$) leads to an increase in the streamwise skin friction for a given $Re$. As suggested earlier, the dependence of $C_f$ on $Pr$ is very weak compared to the other control parameters. Unlike in standard Poiseuille or Couette flow, where a subcritical transition arises due to instability of the laminar base flow and leads to a jump in $C_f$, the streamwise boundary layer transition in mixed convection systems appears smooth. We anticipate that the laminar scaling regime remains relevant until its intersection with the friction law (4.2). The increase of $C_f$ with $Gr$ in the laminar regime would therefore delay the transition to this ‘fully turbulent’ Prandtl friction law (4.2) to higher $Re$. At low $Re$, although the relationship exhibits a laminar-like scaling, one should recall that the convection flow in the interior remains turbulent. In figure 4(b), we focus on clarifying this low $Re$ regime and the increase in $C_f$ with stronger convection. Across both mixed convection systems and for a range of $Pr$, we find a collapse of the data upon rescaling by $Gr^{1/4}$. The $Re^{-1}$ scaling that arises from laminar profiles in Couette/Poiseuille flow appears somewhat too steep to describe the data accurately. Blass et al. (2020) reported a scaling $C_f\sim Re^{-0.90}$ in C-RB, but at this time there is no theoretical basis for such a scaling. Note that one could equivalently express the simplified $Re^{-1}$ collapse as $C_f \sim Ri^{1/4}\,Re^{-1/2}$ using the definitions of (2.3a–c). In the case of mixed convection, the buoyancy-driven flow generates non-zero Reynolds stresses in the equation for the mean profile, which can be expected to modify the mean momentum budget and lead to such an increase in $C_f$. This will be analysed in further detail for mixed VC in § 5.

Figure 4. (a) Friction coefficients calculated for the streamwise component of velocity as a function of Reynolds number. (b) Friction coefficients normalised by $Gr^{1/4}$ to collapse the data at low $Re$. Large markers are used for the mixed VC simulations, whereas small dots show the data from the mixed RB cases of Yerragolam et al. (2024). Grashof numbers are highlighted by the colour of the markers, and for the mixed VC cases, different markers signify different $Pr$. Black dashed lines show a scaling law of $Re^{-1}$, and the black solid line marks the Prandtl friction law of (4.2). This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-4.ipynb.

We now turn to the friction coefficient associated with the buoyancy-driven component of the flow. For the mixed VC system, we can simply take the peak velocity $W_{max}$ of the mean vertical velocity $\bar {w}(y)$ as the relevant velocity scale, and directly measure the mean vertical shear stress at the wall, $\tau = \rho \nu \,\partial _y \bar {w}$. In the RB configuration, the convection has no preferential direction along the walls, resulting in zero mean shear stress. However, we can still construct a friction coefficient associated with the persistent large-scale circulation by using the root mean square (r.m.s.) horizontal velocity profile $u_H(y) = (\overline {u^2 + v^2})^{1/2}$. When a horizontal crossflow is added to the RB system, the large-scale circulation aligns itself perpendicular to the imposed flow (Pirozzoli et al. 2017; Yerragolam et al. 2022), so for the mixed RB cases, we construct a friction coefficient using only the spanwise r.m.s. velocity. Defining the friction coefficient in this way will be appropriate only for cases where the convectively driven flow is stronger than the spanwise velocity fluctuations induced by the turbulent shear flow, i.e. for $Ri\geq O(1)$.

In terms of the Reynolds number, the plate separation $H$ is no longer the appropriate length scale for describing the boundary layer dynamics of the convectively driven flow. As shown in the inset of figure 5(a), the mean profile $\bar {w}(y)$ of the vertical velocity in (mixed) VC reaches its peak value at a certain wall-normal distance $\delta$. From this, we can define a boundary layer Reynolds number

(4.3)\begin{equation} \textit{Re}_\delta = \frac{W_{max} \delta}{\nu} \end{equation}

that drives the behaviour at the wall. We construct an analogous Reynolds number for the mixed RB system using the spanwise r.m.s. velocity profile. The friction coefficients of the convective flow component are plotted against $Re_\delta$ in figure 5. Note that $Re_\delta$ is not known a priori, but is itself a response parameter of the system that varies with $Gr$, $Pr$ and $Re$. Similar to the low-$Re$ regime for the streamwise friction, we observe a power-law scaling close to $C_f\sim Re_\delta ^{-1}$.

Figure 5. (a) Friction coefficients associated with the convectively driven flow as a function of the boundary layer Reynolds number $Re_\delta =W_{max} \delta /\nu$, where $W_{max}$ and $\delta$ are computed as in the mean velocity profile shown in the inset. (b) The same data normalised against the values from the corresponding natural convection system. For the mixed VC cases, different symbols denote different $Pr$ values as outlined in the legend of figure 4. The colour of the markers represents $Gr$, and simulations with $Ri<1/4$ are plotted as semi-transparent. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-5.ipynb.

This is made clearer in figure 5(b), where we collapse the data using the friction coefficient $C_0$ and Reynolds number $Re_{\delta,0}$ obtained from the corresponding natural convection flows, matching $Ra$ and $Pr$. Deviations from the scaling relation are observed for cases where $Ri<1/4$, highlighted by translucent symbols in the figure. As mentioned above, for the mixed RB system this is likely an artefact of using the spanwise r.m.s. velocity to construct the friction coefficient. However, we also observe the discrepancy for low $Ri$ in mixed VC, suggesting that at low Richardson numbers, the turbulence generated by the imposed horizontal flow disrupts the near-wall vertical velocity. Within the range of parameters explored here, the pressure-driven horizontal flow does not modify the vertical $Re_\delta$ by more than a factor 4, suggesting that even in the case of mixed convection, the primary control parameters determining $Re_\delta$ are $Gr$ and $Pr$. As discussed in the Appendix of Howland et al. (2022), a ‘fully turbulent’ transition of this boundary layer may be possible, analogous to the so-called ‘ultimate regime’ in RB convection (Lohse & Shishkina 2023), but only at very high $Gr$. A more in-depth analysis of the mean vertical momentum budget in mixed VC will be presented in § 5.

4.2. Nusselt number

The dimensionless heat flux through the system is characterised by the Nusselt number, defined as

(4.4a,b)\begin{equation} \textit{Nu} = \frac{F_\theta}{\kappa \varDelta / H},\quad F_\theta ={-}\kappa\,\frac{\partial \bar{\theta}}{\partial y} + \overline{v'\theta'} . \end{equation}

Here, $F_\theta$ is the horizontal heat flux through the system (normalised by the specific heat $\rho c_p$), and the overbar denotes averaging in time and in directions parallel to the plates. Integration of the mean temperature equation (3.1) shows that $F_\theta$ is constant across the domain in a statistically steady state.

In figure 6(a), we plot the Nusselt number compensated by $Pr^{1/2}$. This prefactor seems appropriate for the $Pr$ dependence of passive scalar transport in turbulent boundary layers when $Pr\lesssim O(1)$ (Kays et al. 2005), although at higher $Pr$, one expects a transition towards a $Pr^{1/3}$ dependence (Kader & Yaglom 1972; Alcántara-Ávila & Hoyas 2021). The data of both systems converge towards the expression

(4.5)\begin{equation} \textit{Nu}\approx \frac{C_f(\textit{Re})}{4}\,\textit{Re}\,\textit{Pr}^{1/2}, \end{equation}

where $C_f(\textit {Re})$ follows the Prandtl law of (4.2). This expression draws a parallel between the transport of heat and momentum at the wall, known as the Reynolds analogy, and describes the heat transport in ‘forced convection’ when buoyancy no longer affects the flow. From these data, we anticipate that the forced convection expression (4.5) applies for Reynolds numbers greater than that where the friction coefficient begins to follow the turbulent friction law shown in figure 4. At low $Re$ (or more precisely high $Ri$), the Nusselt number responses of the two systems (VC and RB) do not match as precisely as the friction coefficients. Indeed, in the absence of an external flow, VC and RB do not exhibit the same $Nu(Ra,Pr)$ response due to the lack of coupling between the kinetic energy budget and the heat flux in VC (Ng et al. 2015).

Figure 6. Nusselt numbers plotted as functions of Reynolds number. (a) The Nusselt number is normalised by $Pr^{1/2}$ to focus on the response in the shear-dominated regime at high $Re$. The black dashed line marks the Reynolds analogy $Nu\approx \frac {1}{4}C_f\,Re\,Pr^{1/2}$, where $C_f$ satisfies the Prandtl friction law (4.2). (b) The data are normalised by the values associated with natural convection. Here, the black dashed line marks the recently proposed scaling relation $Nu/Nu_0 \sim (\sqrt {1 + (Re/Re_0)^2})^{-1/5}$ from Yerragolam et al. (2024). This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-6.ipynb.

However, we observe a more universal behaviour when the Nusselt numbers in the mixed convection systems are normalised by the values for the equivalent natural convection systems $Nu_0(Gr, Pr, Re) = Nu(Gr, Pr)|_{Re=0}$. In figure 6(b), the normalised Nusselt numbers are plotted as a function of the input Reynolds number normalised by the Reynolds number of the natural convection case $Re_0=W_0H/\nu$. For the mixed VC cases, we define $W_0(Gr, Pr, Re) = W_{max}(Gr, Pr)|_{Re=0}$ as the peak velocity of the natural VC flow, as highlighted in the inset of figure 5(a). For mixed RB, we follow Yerragolam et al. (2024) in using the volume-averaged r.m.s. velocity for $W_0$, and note that the results are insensitive to this choice of velocity scale in describing the ‘wind’ of the large-scale circulation. For $Re/Re_0=O(1)$, all the data from both configurations collapse onto a single curve, showing a drop in the heat flux of up to 25 %. Given this collapse, the $Re_0$ of the natural convection appears to be a critical $Re$ above which the Nusselt number is significantly affected by the horizontal crossflow. Yerragolam et al. (2024) provide an estimate for the drop in $Nu$, derived from the kinetic energy balance in mixed RB, but this balance cannot be related to the horizontal heat flux that is relevant for mixed VC.

In summary, in this section we have demonstrated the universality in the global response parameters of mixed RB and mixed VC, namely in the friction coefficient $C_f$ and the Nusselt number, and the limitations of this universality. In the following sections, we will compare more local quantities, starting with the wall-normal profiles.

5. Wall-normal profiles in mixed VC

We now turn to the first- and second-order statistics, averaged parallel to the plates, to further investigate the dynamics behind the observed global responses. For clarity, we focus solely on the new simulations of mixed VC, and study the variation across the three-parameter space of $Gr$, $Pr$ and $Re$.

We begin with the response of the mean streamwise velocity $\bar {u}(y)$ in figure 7. For a fixed $Gr=10^6$, as in figures 7(a,b), the effect of increasing $Re$ can be seen most clearly when the mean velocity profiles are scaled by viscous wall units in figure 7(b). As $Re$ increases, the velocity profile tends towards the classical log-law profile, with the case $Re=10^4$ closely matching the profile of turbulent Poiseuille flow, as in e.g. Lee & Moser (2015). The effect of stronger thermal convection on the mean profile is also similar to that observed in other mixed convection systems in the literature (Scagliarini et al. 2015; Blass et al. 2020). In figure 7(c), where $Re=10^3$ is fixed and $Gr$ varies between $10^6$ and $10^8$, higher $Gr$ leads to a flatter mean profile in the bulk of the channel. This is illustrated further in wall units in figure 7(d), where a significant drop in $\bar {u}^+$ is observed for $y^+ =O(10)$. Plus symbols denote scaling in viscous wall units, with velocity $u_\tau = \sqrt {\tau _u/\rho }$ and length $\nu /u_\tau$. Such a drop is consistent with the previous findings of Scagliarini et al. (2015) for mixed RB, who proposed a modified log law based on mixing length arguments coupled to the temperature field.

Figure 7. Mean profiles of the streamwise velocity $\bar {u}(y)$. (a,b) Data from simulations at a fixed $Gr=10^6$, varying $Re$ and $Pr$. The $Re$ variation is denoted by the line colour, and the $Pr$ variation is denoted by the line style. (c,d) Present data from simulations at a fixed $Re=10^3$, varying $Gr$ and $Pr$. (a,c) The profiles are normalised by the imposed bulk velocity $\langle u\rangle = U$ and the plate separation $H$. (b,d) The profiles are presented in viscous wall units, where $\bar {u}^+ = \bar {u}/{u_\tau }$ and $y^+=y u_\tau /\nu$. Dashed black lines represent the linear relation $u^+=y^+$ and the logarithmic region $u^+ = \kappa _u^{-1} \log {y^+} + B$, where the von Kármán constant $\kappa _u = 0.41$ and $B=5$ are taken from Pirozzoli et al. (2014). This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-7.ipynb.

Further insight for the streamwise velocity can be gained from the appropriate component of the Reynolds stress. Considering a statistically steady state, we average the streamwise component of the momentum equation (2.1) to obtain

(5.1)\begin{equation} \partial_y \overline{uv} = \nu\,\partial_{yy}\bar{u} + \bar{G}, \end{equation}

where an overbar denotes an average in the periodic ($x, z$) directions and in time. Note that wall-normal advective fluxes in incompressible channel flows are purely turbulent, i.e. $\overline {uv}\equiv \overline {u'v'}$ since $\bar {v}\equiv 0$. From volume averaging, we can also relate the mean pressure gradient forcing to the mean wall shear stress through $\bar {G}=2\tau _u/\rho H$. The first integral of (5.1) can therefore be written as

(5.2)\begin{equation} \partial_{y^+} \bar{u}^+- \overline{uv}^+= 1 - \frac{2y}{H} . \end{equation}

From (5.1) and (5.2), the close coupling of the mean streamwise velocity and the Reynolds stress $\overline {uv}$ is evident. We therefore present the profiles of $\overline {uv}(y)$ in figure 8. As expected, the Reynolds stress dominates the viscous contribution to (5.2) away from the walls, leading to a balance of $-\overline {uv}^+\approx 1 - 2y/H$ as shown in figures 8(a,c). Relative to $H$, the boundary layer in which the viscous term is relevant becomes thinner as both $Re$ and $Gr$ increase. The near-wall behaviour of $\overline {uv}$ exhibits a remarkable collapse when scaled by $Ri^{1/4}$ as in figures 8(b,d), except for the highest $Re$ case with $Ri=0.01$.

Figure 8. Wall-normal profiles of the streamwise Reynolds stress component $\overline {uv}(y)$ scaled by the streamwise friction velocity $\overline {uv}^+ = \overline {uv}/u_\tau ^2$. As in figure 7, (a,b) are at fixed $Gr=10^6$ and (c,d) are at fixed $Re=10^3$. Colours and line styles are as detailed in the caption of figure 7. (a,c) Profiles relative to the plate separation $H$; (b,d) profiles in viscous wall units scaled by $Ri^{1/4}$. The inset of (d) presents the data of (b,d) with a logarithmic scale on both axes to highlight the near-wall collapse of the data. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-8.ipynb.

The additional factor $Ri^{1/4}$ suggests that the appropriate near-wall length scale for the Reynolds stress is modified from the standard viscous wall unit as

(5.3a,b)\begin{equation} \textit{Ri}^{1/4} y^+= \frac{y}{l},\quad l =\frac{\nu}{\sqrt{\dfrac{U_f}{U}\,\dfrac{\tau_u}{\rho}}}. \end{equation}

The additional prefactor $U_f/U$ in front of the shear stress suggests that the vertical, convectively driven component of shear cannot be neglected when considering the streamwise Reynolds stress. An improved collapse to that seen in figures 8(c,d) can be found by computing a viscous length scale using the total shear stress at the wall, $\tau = \sqrt {\tau _u^2 + \tau _w^2}$. As shown explicitly in Appendix A, with this scaling, the Reynolds stress for the highest $Re$ case also matches the other curves. The presence of the convectively driven flow increases the near-wall Reynolds stress, which in turn leads to the flattened mean velocity profiles observed in figure 7(d). This result explains qualitatively the origin of the change in $C_f$ with $Gr$ seen in figure 4, where enhanced buoyancy driving led to a larger skin friction. The increase in near-wall Reynolds stress that arises due to convection thins the boundary layer of the mean horizontal velocity, which in turn produces a larger mean gradient at the wall and a larger friction coefficient $C_f$.

The effect of convection on the shear is by no means a one-way interaction. This is evident from the modification of the mean vertical velocity profile by the imposed horizontal flow. In figures 9(a,b), vertical velocity profiles are shown for fixed $Gr=10^6$ and varying $Re$, with the reference natural convection case ($Re=0$) highlighted in blue for comparison. Compared to the natural VC case, the introduction of horizontal driving at moderate $Re$ leads to an increase in the peak vertical velocity, and hence an increase in the mean shear both in the bulk and at the walls. At the highest $Re=10^4$ (the darkest red line in figure 9), a subsequent decrease is observed in the peak vertical velocity, as well as a nonlinear profile in the bulk. None of the cases studied here exhibits a log layer in the vertical velocity, with the largest $\bar {w}^+$ being approximately 6.5 for the most strongly convective case, $Gr=10^8$. For fixed $Ri=1$ (shown in figures 9c,d), all cases show a similar increase in the peak velocity, and the mean gradient in the bulk appears largely independent of $Gr$ and $Pr$. The distance of the velocity peak from the wall (compared to the channel width $H$) decreases for larger $Gr$, but the value of the peak velocity in free-fall units does not depend strongly on $Gr$. This similarity at constant $Ri$ is suggestive that the vertical velocity modification is determined primarily by the relative strength of the horizontal flow to convection.

Figure 9. Mean profiles of the vertical velocity $\bar {w}(y)$. As in figure 7, (a,b) are at fixed $Gr=10^6$, but here (c,d) are at fixed $Ri=1$. (a,c) Profiles are normalised by the free-fall velocity $U_f = \sqrt {g\alpha \,\Delta H}$ and the plate separation $H$. (b,d) Profiles are presented in viscous wall units, where $\bar {w}^+ = \bar {w}/{w_\tau }$ and $y^+=y w_\tau /\nu$. Reference data for natural VC (with $Gr=10^6$, $Re=0$) are shown in blue. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-9.ipynb.

To further investigate the behaviour of the vertical velocity profiles, we now turn to the mean vertical momentum budget. Unlike the streamwise velocity in (5.1)–(5.2), the vertical velocity is not only tied to the Reynolds stress profile. Rather, the mean vertical momentum equation reads

(5.4)\begin{equation} \partial_y \overline{vw} = \nu\,\partial_{yy} \bar{w} + g\alpha \bar{\theta} . \end{equation}

Close to the wall, we expect the Reynolds stress to be negligible and a balance to arise between buoyancy and viscosity. Due to the symmetry of the boundary conditions, all three terms must be zero at the channel centre. In the bulk of the flow, the mean velocity is approximately linear, so we expect a balance between buoyancy and Reynolds stress. These features are present in each of the simulations highlighted in figure 10. As $Re$ increases, the key modification to the budget arises in the Reynolds stress term $-\partial _y \overline {vw}$. In natural VC, there is a minuscule positive peak in the Reynolds stress term close to the wall, but its amplitude is so small that it is indistinguishable in figures 10(a,f). This peak grows with $Re$, becoming visible at $Re=10^{3.5}$ in figures 10(d,i), coinciding with a flattened profile of the viscous term (shown in green). By $Re=10^4$, at which the Reynolds stresses are significantly energised by the horizontal forcing, the nonlinear term peak exceeds the contribution from the buoyancy term. This leads to a significant drop in the viscous term at $y\approx 10^{-2}H$, which determines the modified mean velocity observed in figures 9(a,b).

Figure 10. Profiles of the mean vertical momentum budget terms for $Pr=1$ and (a–e) $Gr=10^6$ or ( f–j) $Gr=10^7$. Reynolds numbers vary from (a,f) $Re=0$ to (e,j) $Re=10^4$. Colours denote the budget term being plotted as detailed in the legend, and each budget term is plotted normalised against the buoyancy scale $g\alpha \varDelta$. Within each column, semi-transparent lines are added to show the profiles for the Reynolds numbers highlighted by the other rows as comparison. The wall-normal coordinate $y$ is plotted on a logarithmic axis to highlight the variation in the near-wall region. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-10.ipynb.

Changes in the mean temperature profile $\bar {\theta }$ in figure 10 are more subtle, with the most obvious feature being a slight drop in figure 10(e), coinciding with the Reynolds stress peak. A clearer picture of the mean temperature response can be found in figure 11, where profiles are plotted both on linear axes and in wall units. Since the dimensionless wall temperature gradient is equivalent to the Nusselt number, which is most strongly dependent on $Ra$, we compare the temperature profiles in figures 11(a,b) at varying $Re$ and $Pr$ but fixed $Ra=10^7$. As shown in these plots, the temperature gradient in the bulk increases with $Re$, and any possible log-law profile does not collapse to a universal slope or coefficient. For fixed $Ri=1$, however, shown in figures 11(c,d), the bulk gradient appears independent of $Gr$, suggesting that $Ri$ and $Pr$ are the key control parameters determining the flow properties away from the walls.

Figure 11. Mean temperature profiles for (a,b) fixed $Ra=10^7$ and (c,d) fixed $Ri=1$. As in figures 7–9, colours denote changes in $Re$ and $Gr$, whereas line styles show changes in $Pr$. (a,c) Temperature is normalised by the difference across the plates $\varDelta$. (b,d) Temperature is shown in terms of wall units $\theta ^+=(\theta _w - \theta ) U_\tau / F_\theta$, where $\theta _w$ is the temperature at the wall $y=0$. The $U_\tau$ value used here and in $y^+=yU_\tau /\nu$ is calculated using both components of the wall shear stress, as defined in (A1). This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-11.ipynb.

The mean profiles presented here do not further clarify the origin of the non-monotonic behaviour of $Nu$ with $Re$ seen in figure 6. When normalised by the wall heat flux, the mean temperature profiles in figure 11(b) show a monotonic increasing trend with $Re$ away from the wall. The non-monotonic behaviour observed in the vertical velocity profiles of figures 9(a,b) suggests an intrinsic connection between the vertical flow and the wall-normal heat transport, although the very subtle changes in the vertical momentum budget remain difficult to interpret in the context of the non-monotonic Nusselt number variation. The later sections on spectral analysis and statistics of the boundary layer aim to shed more light on this global heat transport.

To conclude our analysis of the mean profiles, we compare the new simulations of mixed VC to the universal functions of the MOST (Monin & Obukhov 1954; Foken 2006). As outlined in the Introduction, the theory considers a region of the flow where viscosity and molecular diffusion are negligible, and the appropriate length scale is the Obukhov length

(5.5)\begin{equation} L_O=\frac{u_\tau^3}{g\alpha F_\theta}. \end{equation}

On dimensional grounds, universal functions that depend only on $y/L_O$ for the first- and second-order statistics can then be constructed for the temperature and velocity fields. One key assumption used to derive the universal functions for velocity is that the turbulent momentum flux $\rho \,\overline {uv}$ is constant over the region of interest. However, as shown in figure 8, this is not true in a turbulent channel flow driven by a pressure gradient, where the Reynolds stress follows a linear relation in the bulk. We can therefore expect the universal functions of Monin & Obukhov (1954) to hold for only a small range of $y$ in our simulations.

Considering these caveats, we present a collection of mean profiles from the mixed VC simulations as a function of $y/L_O$ in figure 12. Following Pirozzoli et al. (2017), who present a similar figure (their figure 16) for mixed RB convection, we emphasise only the region of the domain where the turbulent fluxes of heat and momentum are greater than 80 % of their maximum value, with the rest of the profiles made semi-transparent. We compare our results with the updated scaling theories of Kader & Yaglom (1990), and for the mean temperature and velocity profiles, the commonly used Businger–Dyer relations outlined in Foken (2006). For the highlighted region where the fluxes are close to their peak, the majority of the results collapse onto a single function of $y/L_O$ close to these profiles. In figure 12(a), there are two notable discrepancies from the theoretical estimates for the mean temperature profile. First, the values at transitional $y/L_O\approx 10^{-1}$ are approximately a factor 2 larger than the suggested profiles. Second, the effective power-law scaling observed for strong buoyancy driving $y\geq L_O$ appears steeper than the $-1/3$ proposed by Kader & Yaglom (1990) (dashed black lines), and closer to the $-1/2$ scaling associated with the asymptotic limit of the Businger–Dyer relation (solid blue line). A discrepancy in the two theoretical predictions for the mean velocity profile at $y\gg L_O$, highlighted in figure 12(b), cannot be resolved from the results presented here, although we note that the data here agree well with the mixed RB data of Pirozzoli et al. (2017). The variation in $Pr$ across the simulations collapses for most quantities in figure 12, except for the r.m.s. temperature shown in figure 12(e). The origin of this discrepancy with $Pr$ is currently unclear, but understanding how the $Pr$ dependence may affect the Monin–Obukhov profiles is important in the context of their applications outside the atmosphere, for example in the surface layer of the upper ocean (Zheng, Harcourt & D'Asaro 2021).

Figure 12. Comparison of the mixed VC simulation results with the ‘universal functions’ of the MOST for first- and second-order statistics. Profiles are shown as semi-transparent unless the advective wall-normal fluxes of heat $\overline {v \theta }$ and momentum $\overline {vu}$ are within 80 % of their maximum value. Consistent with the Monin–Obukhov formulation, all quantities here are normalised by the horizontal component of the friction velocity $u_\tau$ and the wall-normal heat flux $F_\theta$. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-12.ipynb.

6. Spectral analysis

We now investigate the scale dependence of the thermal structures and heat flux in mixed convection through analysis of the power spectra. To ensure that we capture the full range of dynamical scales, we perform further simulations in extended aspect ratio domains, with $L_x=L_z=24 H$. The details of these simulations are outlined in table 2 of Appendix B.

In unsheared RB convection, Krug, Lohse & Stevens (2020) show that large-scale patterns, also known as superstructures, can be identified from a low-wavenumber peak in the power spectrum of the temperature field and in the co-spectrum of the wall-normal heat flux. This peak denotes the scale of the large-scale circulation or ‘wind’ of convection that is key to theories describing RB convection. Since there is no preferential horizontal flow direction in RB convection, one-dimensional spectra can be analysed, but in mixed VC, the two directions parallel to the plates must be considered separately.

We therefore compute time-averaged spectra using two-dimensional Fourier transforms in the periodic directions. For visualisation purposes, we present separate one-dimensional spectra for the streamwise ($x$) and spanwise ($z$) wavenumbers, which are computed by integrating over the other wavenumber. Precisely, the two- and one-dimensional co-spectra of any two variables $f$ and $g$ are defined as

(6.1a,b)\begin{equation} \varPhi_{fg}(k_x, y, k_z) = {\rm Re}[\,\hat{f}^* \hat{g}],\quad \varPhi_{fg}(k_x,y) = \int {\rm Re}[\,\hat{f}^* \hat{g}]\, \mathrm{d}k_z, \end{equation}

where $\hat {f}$ denotes the two-dimensional Fourier transform of $f$ in $x$ and $z$, and ${\rm Re}$ denotes the real part. The Fourier transforms are normalised such that integrating the spectrum over wavenumber space recovers the corresponding volume-averaged quantity $\langle\, f g \rangle = \iint \varPhi _{fg} \,\mathrm {d}k_x \,\mathrm {d}k_z$.

The power spectrum of temperature $\varPhi _{\theta \theta }$ is presented in figure 13 for the four extended simulations. In figure 13(a), for the standard RB configuration, we see behaviour similar to that in Krug et al. (2020), with a distinct, sharp low-wavenumber peak at $k H\approx 1$, and a more broad peak at smaller scales. As mentioned above, the lack of a preferential direction means that the two directional spectra are virtually identical. When shear is added to the RB system in figure 13(b), the horizontal isotropy in the system is destroyed. As expected from previous work on mixed RB (Pirozzoli et al. 2017; Blass et al. 2020,  2021; Yerragolam et al. 2022), coherent rolls aligned with the streamwise axis dominate the signal. Note the difference in $y$ axis between figures 13(a) and 13(b). In the streamwise spectrum $\varPhi _{\theta \theta }(k_z,y)$, a sharp peak at $kH\approx 2$ is visible at all wall-normal locations, but is particularly prominent in the near-wall region (light green) where there is a maximum in $\overline {\theta ^2}(y)$. The corresponding wavelength of these structures is $\lambda = 2{\rm \pi} H/k_z \approx 3 H$, which is consistent with that observed previously in the mixed convection literature at $Ri=1$ (e.g. Pirozzoli et al. 2017). The alignment of the flow structures in the streamwise direction leads to a broad contribution to the streamwise spectrum at low $k_x$. This behaviour is purely a result of the regular alignment, and is not related to any domain size effects.

Figure 13. One-dimensional power spectra $\varPhi _{\theta \theta }(k,y)$ of temperature for the extended domain simulations described in table 2. Solid lines denote spectra as functions of $k=k_x$, and dashed lines are functions of $k=k_z$. The colour of each line is determined by its wall normal location $y$, with the specific locations highlighted on the colour bar. Each spectrum is multiplied by the wavenumber $k$ such that area under each curve is representative of the relative contribution in wavenumber space. (a,c) Simulations of natural convection ($Re=0$). (b,d) Simulations of mixed convection at $Ri=1$. (a,b) The RB cases with gravity in the wall-normal ($y$) direction. (c,d) The VC cases with gravity parallel to the wall in z. Note that P-VC here denotes the Poiseuille-VC configuration to compare with P-RB. All simulations have fixed $Gr=10^7$, $Pr=1$. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-13.ipynb.

In figures 13(c,d), we present the same analysis but for natural and mixed vertical convection. The contrast to the RB cases is immediately apparent, with no sharp peaks at any wavenumber or wall-normal position. This confirms the earlier visual observation of figure 3, where the instantaneous snapshots showed no clear coherent length scale. In natural VC, shown in figure 13(c), the two one-dimensional spectra do not overlap exactly as in the RB case due to the buoyancy-driven mean flow in the vertical. The $k_x$ spectrum shows a small peak at $kH\approx 3$ that is absent from the $k_z$ spectrum and becomes more prominent towards the channel centre. Comparing figures 13(c) and 13(d), we see that this peak is suppressed by the addition of external shear. Furthermore, all of the broad mid-range peaks in the spectra exhibit a decrease in amplitude and a shift to higher wavenumbers.

The heat flux co-spectra $\varPhi _{v\theta }$, shown in figure 14, exhibit similar features to the power spectra. In figures 14(a,b), a sharp low-wavenumber peak is once again visible for RB and mixed RB systems, and for the heat flux, this peak increases towards the channel centre. A broad high-wavenumber peak is also observed in the spectrum close to the wall. This peak flattens and shifts to lower wavenumbers (i.e. larger scales) as distance from the wall increases, which can be interpreted as the emergence and coalescence of small-scale plume structures. A similar behaviour can be found in the VC configurations of figures 14(c,d), with the broad peak in both the $k_x$ and $k_z$ spectra shifting to lower wavenumbers for larger $y$. However, in the $k_x$ spectrum for natural VC (where $x$ is perpendicular to gravity), this peak also increases in amplitude towards the channel centre, highlighting that the majority of the wall-normal heat transport at the channel centre occurs due to structures of size $\lambda \approx H$. Comparing this to the mixed VC case in figure 14(d), the mid-scale peak appears significantly suppressed when the external shear is imposed, suggesting that the coalescence of plumes in the bulk is disrupted by the shear. This is consistent with the interpretation of Domaradzki & Metcalfe (1988) and Scagliarini et al. (2014), who argued that the drop in heat flux in mixed RB can be related phenomenologically to disruption of the organisation of small-scale convective plumes.

Figure 14. One-dimensional co-spectra $\varPhi _{v\theta }(k,y)$ of the wall-normal heat flux for the extended domain simulations described in table 2. See the caption of figure 13 for more details on the meanings of line styles and colours. The spectra are normalised by $\sqrt {Ra\,Pr}$ such that integration over $k$ recovers the dimensionless advective heat flux $Nu-1 = \langle v\theta \rangle /(\kappa \varDelta /H)$. As in figure 13, the configurations presented are (a) RB, (b) mixed RB, (c) VC, (d) mixed VC. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-14.ipynb.

7. Heat flux statistics in the boundary layer

Whereas the previous section focused on the advective heat flux in the bulk and its modification due to shear, we now turn to the conductive heat flux that dominates the transport in the boundary layer. Specifically, we investigate the statistical distribution of the local conductive heat flux at the boundary plates. As defined earlier, in (3.1), we consider the local dimensionless heat flux $q_\theta$ whose time and plane average is equivalent to the Nusselt number. The pointwise data of $q_\theta (x,z,t)$, also shown in figure 2, are collected over time and space to construct histograms that are normalised to produce probability density functions (p.d.f.s) in figure 15. The data shown here are only for $Pr=1$ and $Gr=10^6, 10^7$ since these cases cover the widest range of Richardson numbers, highlighting the transition from natural convection to forced convection.

Figure 15. Probability density functions (p.d.f.s) of the dimensionless local heat flux $q_\theta$ as defined in (3.1) for (a,c) $Gr=10^6$ and (b,d) $Gr=10^7$. Line colours highlight variation in the Reynolds number, as shown by the colour bar, with the blue dashed line plotting the data for natural convection with $Re=0$. Plots (a,b) are presented with linear $y$ axes, and (c,d) present the same data on logarithmic $y$ axes to highlight the exponential tails. This figure is also available as an interactive JFM notebook: https://www.cambridge.org/S0022112024005986/JFM-Notebooks/files/Figure-15.ipynb.

In figures 15(a,b), where the data are presented on linear axes, we find the same effect of increasing shear at both Grashof numbers. As $Re$ increases from zero (shown by progressively darker lines, with $Re=0$ marked in blue for comparison), the peak of the p.d.f. increases in amplitude. In statistical terms, this means that the most common values of local heat flux become more common as shear is increased, up to $Re\approx 10^{3.5}$. Due to the skewed nature of the distributions, these increasingly common values of heat flux are below the mean, and although the peak shifts slightly to the right as $Re$ increases, these parameters are associated with the drop in heat flux observed in figure 6. This corresponds to the visualisation of figure 2, where the streaks of low heat flux span a larger proportion of the boundary at $Ri=1$ (figures 2e,i) than at high $Ri$ (figures 2b,c). Once $Re$ is sufficiently high, which in these cases is for $Re> 10^{3.5}$, the whole distribution shifts more significantly to higher $q_\theta$ as the transport becomes dominated by the pressure-driven shear flow.

At moderate $Re$, another key modification is to the tails of the distribution. These are visualised most clearly in figures 15(c,d), where the heat flux distributions are plotted on a logarithmic axis. In this representation, it is clear that the right tails of large heat flux decay exponentially for $Gr=10^7$, and close to exponentially for $Gr=10^6$. At both Grashof numbers, the tails are reduced as $Re$ increases relative to the natural convection case, showing that the probability of extreme local heat flux values is reduced due to the introduction of a mean horizontal flow. Again, this is consistent with what was observed in the visualisation of figure 2, where dark patches associated with large local heat flux became less prevalent as $Re$ increases. In natural VC, Pallares et al. (2010) used conditional sampling to show that these patches are most often associated with instantaneous flow reversals at the walls, which are likely a result of impacting plumes originating from the opposing wall.

8. Conclusion and outlook

In this paper, we have investigated mixed convection in a differentially heated vertical channel subject to a horizontal pressure gradient, referred to as mixed vertical convection (VC), through direct numerical simulations. By simulating the system across the parameter ranges $10^6 \leq Gr \leq 10^8$, ${1 \leq Pr \leq 10}$ and ${10^{2.5} \leq Re \leq 10^4}$, we have explored the transition from natural convection to forced convection, characterised by the Richardson number $10^{-2} \leq Ri \leq 10^2$. Across this parameter range, the response of the streamwise skin friction is identical to the response seen in mixed Rayleigh–Bénard (RB) convection, with a power-law $Re$ dependence for $C_f$ giving way to the Prandtl friction law (4.2) at sufficiently high $Re$. The presence of convection acts to increase the skin friction for a given $Re$, due to the thermal plumes emitted from the plates, with a collapse observed for all the data in the power-law regime $C_f \sim Gr^{1/4} Re^{-\gamma }$ for both configurations across the entire range of $Pr$. For mixed VC, the streamwise momentum budgets show an enhanced impact of Reynolds stresses close to the wall at higher $Ri$, which modify the shape of the mean velocity profile. The Reynolds stresses in the boundary layer are driven by the combined shear stress of the horizontal pressure-driven flow and the VC flow.

Friction coefficients could also be obtained for the flow component associated with buoyancy driving. For cases with significant buoyancy effects at $Ri > 1/4$, a reasonable collapse was found for the laminar-like scaling $C_f\sim Re_\delta ^{-1}$ where $Re_\delta$ is the boundary layer Reynolds number of the convection flow. The introduction of the horizontal pressure gradient leads to significant modification of the mean vertical velocity, with $Re_\delta$ increasing up to three times its value for natural convection, and a corresponding increase in the mean shear in the channel bulk for moderate $Re$. For $Ri=O(1)$, the response of the Nusselt number $Nu$ to the introduction of shear in mixed VC is also identical to that in mixed RB. Before the flow undergoes a transition to a forced convection regime at high $Re$, the response can be expressed as $Nu/Nu_0 = f(Re/Re_0)$, where $Nu_0$ and $Re_0$ are the Nusselt number and Reynolds number associated with natural convection. The data of both mixed convection systems collapse onto this single curve, which describes the drop in $Nu$ as $Re$ increases.

The near identical quantitative response of mixed RB and mixed VC is observed despite striking qualitative differences between the two configurations in terms of large-scale flow organisation. Whereas mixed RB features large convective rolls oriented along the streamwise axis, no such structures form in mixed VC except at low $Ri$ when the dynamics are dominated solely by the pressure gradient forcing. The absence of a low-wavenumber peak in the heat flux co-spectrum for mixed VC confirms that the advective heat flux across the domain is transported by eddies or plumes with a wide range of scales rather than by large coherent rolls. Comparing the spectra from natural VC with mixed VC reveals that the organisation of plume structures with a horizontal scale $\lambda \approx H$ is suppressed by the horizontal mean flow. This is reflected also in the distribution of local heat flux at the boundaries, which shows that instantaneous events of extreme heat flux become less likely as $Re$ increases. As the boundaries become dominated by streaky structures, the formation of localised plume structures is disrupted, consistent with the earlier interpretations of reduced heat flux in mixed RB (Domaradzki & Metcalfe 1988; Scagliarini et al. 2014).

The striking agreement between the two channel configurations compared here, regardless of the gravity direction, opens up the question of how universal such a response in skin friction and heat flux is for other mixed convection systems. The independence to the gravity direction suggests that our results may be applicable more generally to inclined layer convection, as studied by Daniels & Bodenschatz (2002) at low $Gr$, subject to a horizontal pressure gradient at comparable $Gr$ and $Re$. However, it is less clear how directly applicable the findings here are to a wall plume subject to a crossflow. Such a scenario is relevant to environmental applications such as at a melting ice face subject to ambient ocean currents (Jackson et al. 2020). Whereas both boundaries impact the transport in mixed RB and mixed VC, a wall plume constantly entrains fluid from the ambient at its outer edge. Recent work has provided an innovative way to simulate wall plumes at high $Gr$, through studying temporally growing boundary layers (Ke et al. 2023; Wells 2023), and it would be fruitful to understand how the growth of these boundary layers is affected by the presence of a turbulent crossflow. For ice–ocean interactions, the picture is further complicated through multicomponent transport (Howland, Verzicco & Lohse 2023) and the development of rough boundaries that are closely coupled to the flow structures (Couston et al. 2021; Ravichandran, Toppaladoddi & Wettlaufer 2022).

From a more theoretical standpoint, our current work has also emphasised how certain flow properties are noticeably modified under the transition from natural VC to forced convection. These include the sign of the near-wall Reynolds stress in figure 10 and the distribution of the local wall heat flux in figure 15. In both RB and VC, predictions have been made for a transition to the ‘ultimate’ regime at sufficiently high buoyancy driving, where the boundary layers behave as turbulent boundary layers (Lohse & Shishkina 2023), although this regime has thus far been inaccessible to three-dimensional numerical simulation. Analysis of the aforementioned statistical quantities in natural convection at high $Gr$ may help in identifying key markers of such a transition in natural convection systems.