2. Numerical method and simulation details

We simulated a fluid enclosed between two walls, kept at fixed temperatures $T_{H}$ at $y=0$ and $T_{C}=T_{H}-{\it\Delta}$ at $y=2H$ , and driven by a constant body force mimicking an imposed pressure gradient along the streamwise direction $\hat{x}$ . The equations of motion are the incompressible Navier–Stokes equation for the fluid velocity field $\boldsymbol{u}(\boldsymbol{x},t)$ in the Boussinesq approximation (namely, the density is assumed to be constant ${\it\rho}={\it\rho}_{0}$ , but for the linearized buoyancy term)

( $T_{m}=(T_{H}+T_{C})/2$ being the mean temperature) coupled with the advection–diffusion equation for the temperature field $T(\boldsymbol{x},t)$ ,

In the above equations, $P$ is the pressure field (rescaled by ${\it\rho}_{0}$ ), $\boldsymbol{g}=-g{\hat{y}}$ is the acceleration of gravity and $\boldsymbol{F}=(u_{{\it\tau}}^{2}/H)\hat{x}$ is the acceleration due to the constant body force ( $u_{{\it\tau}}$ being the friction velocity); ${\it\nu}$ , ${\it\alpha}$ and ${\it\beta}$ are the kinematic viscosity, the thermal diffusivity and thermal expansion coefficient, respectively. As a numerical scheme, we adopted a 3d LB algorithm (Benzi, Succi & Vergassola Reference Benzi, Succi and Vergassola1992; Chen & Doolen Reference Chen and Doolen1998; Aidun & Clausen Reference Aidun and Clausen2010) with two probability densities (for density/momentum and for temperature, respectively) (He, Chen & Doolen Reference He, Chen and Doolen1998). The method has been extensively used to study both thermal convection (Benzi, Toschi & Tripiccione Reference Benzi, Toschi and Tripiccione1998; Calzavarini, Toschi & Tripiccione Reference Calzavarini, Toschi and Tripiccione2002) and turbulent channel flow (Toschi et al. Reference Toschi, Amati, Succi, Benzi and Piva1999; Toschi, Lévêque & Ruiz-Chavarría Reference Toschi, Lévêque and Ruiz-Chavarría2000; Biferale et al. Reference Biferale, Lohse, Mazzitelli and Toschi2002); in particular, Lavezzo, Clercx & Toschi (Reference Lavezzo, Clercx and Toschi2011) validated the code and tested grid resolutions against previous studies with different numerical methods.

We used a computational grid $L\times 2H\times W$ of $256\times 128\times 128$ lattice points, with each run longer than $3\times 10^{6}$ time steps (in LB units), in such a way as to achieve statistically steady states of ${\sim}200T_{L}$ ( $T_{L}$ being the large-scale eddy turnover time). From (2.1) and (2.2) two dimensionless groups can be identified, giving rise to two parameters that control the dynamics, namely the shear (or friction) Reynolds number

and the Rayleigh number

quantifying, respectively, the strength of the pressure-induced shear and the strength of the buoyancy with respect to viscous dissipation.

We performed turbulent channel flow simulations with $\mathit{Re}_{{\it\tau}}\in [46,205]$ ; for each $\mathit{Re}_{{\it\tau}}$ we tuned the gravity (at fixed temperature jump), and hence the buoyancy, spanning the range $\mathit{Ra}\in [0,1.3\times 10^{7}]$ . The Prandtl number $\mathit{Pr}={\it\nu}/{\it\alpha}$ is equal to one in all cases presented.

Figure 1.

Mean profiles of the streamwise velocity $U(y)$ , normalized with the friction velocity $u_{{\it\tau}}$ , for various Rayleigh numbers at fixed shear Reynolds number $\mathit{Re}_{{\it\tau}}$ . (a)  $\mathit{Re}_{{\it\tau}}=130$ , (b)  $\mathit{Re}_{{\it\tau}}=205$ . The data have been plotted every two points (i.e. lattice spacings) for the sake of clarity of visualization.

Figure 2.

Mean profiles of the normalized temperature ${\it\Theta}=2(T-T_{m})/{\it\Delta}$ as a function of the height $y$ for various Rayleigh numbers at fixed shear Reynolds number  $\mathit{Re}_{{\it\tau}}$ . (a)  $\mathit{Re}_{{\it\tau}}=130$ , (b)  $\mathit{Re}_{{\it\tau}}=205$ ; the solid line represents the profile for the purely convective Rayleigh–Bénard case (i.e.  $\mathit{Re}_{{\it\tau}}=0$ ) with $\mathit{Ra}=8.12\times 10^{5}$ . The data have been plotted every two points (i.e. lattice spacings) for the sake of clarity of visualization.

4. Modification of the law of the wall by buoyancy

Inspection of the profiles of the mean quantities and the fluctuations showed that buoyancy alters the channel flow structure. We now focus on the mean velocity profiles. In figure 5 we report the lin–log plot of profiles, in wall units, for $\mathit{Re}_{{\it\tau}}=185$ and for various $\mathit{Ra}$ . Convective motion, realized in the form of thermal plumes rising from the walls to the bulk, disturbs the coherence of the channel flow, resulting in a source of drag that decreases the mean velocity in the channel and hence the mass throughput, as well as flattening the velocity profile. In order to quantify these observations, in the following we attempt to generalize von Kármán’s law of the wall accounting for buoyancy effects. To this aim, we propose a simple model constructed from conservation laws and phenomenological arguments, which predicts that the viscous buffer is insensitive to buoyancy (i.e.  $U^{+}=y^{+}$ ), while in the log-law region the following relation

holds (quantities have been expressed in wall units $U^{+}\equiv U/u_{{\it\tau}}$ and $y^{+}\equiv y/{\it\delta}$ , with ${\it\delta}={\it\nu}/u_{{\it\tau}}$ ). In (4.1), ${\it\kappa}$ is the von Kármán constant, $B$ is an integration constant and ${\it\kappa}_{C}$ is given by

$C_{S}$ is a phenomenological parameter. All profiles in figure 5 agree reasonably, within error bars, with our predictions: the data collapse on the function $U^{+}=y^{+}$ for $y^{+}<10$ (dashed lines) while in the logarithmic layer they can be fitted with (4.1) (dash-dotted lines), whose derivation is detailed in what follows.

Figure 4.

The skin-friction coefficient $c_{f}={\it\tau}_{w}/(1/2U_{0}^{2})$ as a function of the shear Reynolds number $\mathit{Re}_{{\it\tau}}$ for various Rayleigh numbers $\mathit{Ra}$ . The dashed lines are predictions obtained by approximating the mean velocity by a log-law over the whole channel; for the stratified cases we made use of (4.1)–(4.2) with the same fitting parameter $C_{s}=2.5$ for all $\mathit{Ra}$ (see § 4).

Figure 5.

Lin–log plot of the velocity profiles in wall units for shear Reynolds number $\mathit{Re}_{{\it\tau}}=185$ and three Rayleigh numbers; the data are shifted vertically for clarity of visualization. The dashed lines indicate the linear law, valid in the viscous layer, $U^{+}=y^{+}$ , while the dash-dotted ones correspond to the prediction (4.1)–(4.2) with ${\it\kappa}=0.42$ , $C_{S}=2.5$ and $B\in [5.8,6.7]$ .

Applying the Reynolds decomposition $\boldsymbol{u}(\boldsymbol{x},t)=U(y)\hat{x}+\tilde{\boldsymbol{u}}(\boldsymbol{x},t)$ for the velocity field and averaging the Navier–Stokes equation for the $\hat{x}$ velocity component, we get, upon integration over $y$ , an exact relation (Pope Reference Pope2001) between the mean shear $S(y)=\text{d}U(y)/\text{d}y$ and the Reynolds stress ${\it\tau}_{xy}(y)=-\overline{\tilde{u} _{x}\tilde{u} _{y}}$ , namely

or, in wall units,

Notice that (4.3)–(4.4) remain valid also in the presence of buoyancy, since the forcing term ${\it\beta}g(T-T_{m}){\hat{y}}$ does not act along the streamwise direction. For very large $\mathit{Re}_{{\it\tau}}$ (in principle $\mathit{Re}_{{\it\tau}}\rightarrow \infty$ ), (4.4) gives $S^{+}+{\it\tau}_{xy}^{+}=1$ ; close to the wall (for $y^{+}<10$ ), where viscous terms dominate, $S^{+}\gg {\it\tau}_{xy}^{+}$ whence $S^{+}\approx 1$ , which implies $U^{+}\simeq y^{+}$ , as anticipated above. Conversely, away from the wall (typically for $y^{+}\gtrsim 30$ ), turbulent fluctuations dominate over viscous processes and ${\it\tau}_{xy}^{+}\gg S^{+}$ , i.e. (L’vov et al. Reference L’vov, Pomyalov, Procaccia and Tiberkevich2004)

Such a result alone, however, does not give any further insight into the behaviour of the velocity profile in the corresponding region of the channel. To this aim, we need to consider the budget equation of turbulent kinetic energy $E_{K}(y)=\overline{|\tilde{\boldsymbol{u}}|^{2}}/2$ (Pope Reference Pope2001). In the log-law layer the energy production $\mathscr{P}={\it\tau}_{xy}(y)S(y)+{\it\beta}gQ(y)$ (where $Q(y)=\overline{u_{y}{\it\theta}}$ is the turbulent heat flux and ${\it\theta}$ are the temperature fluctuations) is balanced by dissipation ${\it\varepsilon}(y)=2{\it\nu}\overline{(\partial _{i}\tilde{u} _{j})^{2}}$ (Pope Reference Pope2001); the latter cannot be calculated exactly but only inferred by means of phenomenological arguments: for large Reynolds number and outside the viscous boundary layer, the dissipation equals the turbulent energy flux, which can be estimated as $E_{K}(y)/{\it\tau}(y)$ , where ${\it\tau}(y)\propto y/\sqrt{E_{K}(y)}$ is the typical eddy turnover time at $y$ (L’vov et al. Reference L’vov, Pomyalov, Procaccia and Tiberkevich2004). Thus, the energy balance equation reads

where $a$ is a non-dimensional parameter of order unity.

Equation (4.6) is not yet closed: a relation between $E_{K}(y)$ and ${\it\tau}_{xy}(y)$ is required. Dimensional analysis suggests that the two quantities should be proportional to each other, i.e.  $E_{K}(y)\propto {\it\tau}_{xy}(y)$ (L’vov et al. Reference L’vov, Pomyalov, Procaccia and Tiberkevich2004); such an assumption can be readily tested in the numerical simulations: in the main panel of figure 6 we show that the ratio $E_{K}(y)/{\it\tau}_{xy}(y)$ indeed approaches a constant value for $y^{+}\gtrsim 40$ . We can then rewrite (4.6) as (in view also of (4.5))

with $c$ yet another dimensionless number. The term for energy production by buoyancy, viz. the heat flux, needs also to be modelled; although more refined closures can be employed (Johansson & Wikström Reference Johansson and Wikström1999; Hattori et al. Reference Hattori, Morita and Nagano2006) involving tensorial eddy diffusivities and coupling with gradients of the temperature field in all directions, we adopt a simple standard mixing length ansatz, i.e.

A first-order closure like (4.8) must not be expected to work well for second-order quantities like temperature fluctuations, Reynolds stresses, etc., but, as we will show, mean streamwise velocity profiles are satisfactorily reproduced through such a model; according to Prandtl’s hypothesis the mixing length $\ell _{m}$ is proportional to the distance from the wall, i.e.  $\ell _{m}\propto y$ , whence

where $b$ is a numerical constant. The mean temperature gradient itself will, in principle, depend on the shear profile; however, in a perturbative spirit, we postulate here, for simplicity, a logarithmic form, such that

where $T_{\ast }={\it\alpha}{\it\Delta}/(2Hu_{{\it\tau}})$ is a scale for the dynamic temperature. The functional form (4.10) somehow interpolates between the passive scalar case (Johansson & Wikström Reference Johansson and Wikström1999) and the natural convection (Ahlers et al. Reference Ahlers, Bodenschatz, Funfschilling, Grossmann, He, Lohse, Stevens and Verzicco2012). Inserted in the expression for the heat flux (4.9), (4.10) gives

Figure 6.

Main panel: ratio of the turbulent kinetic energy $E_{K}$ to the Reynolds stress ${\it\tau}_{xy}$ (for the case $\mathit{Re}_{{\it\tau}}=205$ , $\mathit{Ra}=8.12\times 10^{5}$ ): the plateau attained for $y^{+}\gtrsim 40$ suggests that the two quantities are proportional to each other in the log-law region. Inset: ratio of the turbulent heat flux $Q(y)$ to the mean shear $S(y)$ ; the solid line indicates a linear dependence on $y^{+}$ and represents a numerical check of the validity of the assumptions leading to (4.11).

For the sake of validation of our arguments, we check (4.11) (which predicts $Q(y)\propto yS(y)$ ) against the numerics in the inset of figure 6, finding a reasonably good agreement. Inserting (4.11) into (4.7) provides

which can be recast, introducing the wall units and the definitions of $\mathit{Re}_{{\it\tau}}$ and $\mathit{Ra}$ , in the following form

It must be noticed that the two parameters ${\it\kappa}$ and $C_{S}$ are combinations of dimensionless quantities emerging in the derivation, which cannot, however, be derived from first principles; fits of the numerical data indicate that good estimates are the values ${\it\kappa}=0.42$ and $C_{S}=2.5$ . Making explicit the shear term in (4.13), we get

whose integration finally yields expressions (4.1) and (4.2) for the velocity profile. The robustness of the model is confirmed in figure 7 where we plot the fitted values of the parameter ${\it\kappa}_{C}$ as a function of $\mathit{Ra}$ for fixed $\mathit{Re}_{{\it\tau}}$ (main panel) and as a function of $\mathit{Re}_{{\it\tau}}$ for fixed $\mathit{Ra}$ (inset), together with the predictions of (4.2) (dashed lines).

Figure 7.

Plots of the parameter ${\it\kappa}_{C}(\mathit{Ra},\mathit{Re}_{{\it\tau}})$ of the model as function of $\mathit{Ra}$ for fixed $\mathit{Re}_{{\it\tau}}$ ’s (main panel) and as function of $\mathit{Re}_{{\it\tau}}$ for fixed $\mathit{Ra}$ ’s (inset), whose values are $\mathit{Ra}=8.12\times 10^{5}$ (●), $\mathit{Ra}=6.5\times 10^{6}$ (▪) and $\mathit{Ra}=1.3\times 10^{7}$ (▾). The dashed lines represent (4.2) with constant $C_{S}$ fixed at $C_{S}=2.5$ for all cases.

Before concluding, we assimilate our results, mutatis mutandis, to those obtained in the literature for PBL. Our set-up differs from PBL in that we simulate a channel flow with fixed temperature (and zero velocity) at the two walls; unlike the classical MO54 approach (Monin & Obukhov Reference Monin and Obukhov1954; Dyer Reference Dyer1974), as written in (4.8) and (4.11), we do not assume the heat flux to be constant as a sort of boundary condition (Monin & Obukhov Reference Monin and Obukhov1954) (which is an approximation, as pointed out in Mellor & Yamada (Reference Mellor and Yamada1982)). However, some analogies might be considered. With respect to other models (Mellor & Yamada Reference Mellor and Yamada1982; Dubrulle et al. Reference Dubrulle, Laval, Sullivan and Werne2002b ), our approach enjoys some peculiar features: (i) it relies on one single empirical input (the parameter $C_{S}$ ); (ii) it provides an explicit analytic expression for the velocity profile and for its dependence on $\mathit{Ra}$ , $\mathit{Re}_{{\it\tau}}$ and  $\mathit{Pr}$ , which makes it particularly suitable for comparison against data from channel flow simulations where one has direct control of these parameters. On the other hand, our model must not be expected to be fully trustworthy for high values of the stability parameter ${\it\zeta}=y/L_{O}$ , where $L_{O}=u_{{\it\tau}}^{3}/({\it\kappa}{\it\beta}gQ_{w})$ is the Obukhov length (Monin & Obukhov Reference Monin and Obukhov1954) and $Q_{w}$ the heat flux at the wall. For a further comparison, it is worth rewriting (4.14) in terms of the dimensionless mean wind gradient ${\it\phi}_{U}$ as

where the characteristic length $L_{S}$ is defined as

Figure 8.

Temperature fluctuations ${\it\Theta}_{\ast }\equiv \overline{(T-\overline{T})^{2}}^{1/2}/(Q_{w}/u_{{\it\tau}})$ from our LB simulations (symbols) and RDT prediction (4.19), with the expression (4.17) for $L_{O}$ (dashed lines), as functions of the wall-normal direction coordinate in wall units (the range shown is restricted to the logarithmic layer). The data have been fitted with the theoretical curves, using for the parameters $({\it\alpha}_{1},{\it\alpha}_{2})$ the values: $(0.26,0.6)$ for $\mathit{Ra}=6.5\times 10^{6}$ and $(0.22,0.9)$ for $\mathit{Ra}=1.3\times 10^{7}$ .

If we introduce the Nusselt number $\mathit{Nu}=Q_{w}/({\it\alpha}{\it\Delta}/2H)$ , the Obukhov length $L_{O}$ , in turn, can be expressed in terms of $(\mathit{Ra},\mathit{Re}_{{\it\tau}},\mathit{Pr})$ as

where ${\it\lambda}_{{\it\theta}}$ is the boundary layer width (the last equality is only approximate since, as mentioned before, in certain cases there is no clear thermal shortcut (Grossmann & Lohse Reference Grossmann and Lohse2000)). A direct look at (4.16) and (4.17) suggests that

but ${\it\lambda}_{{\it\theta}}$ is fixed by the heat flux from the wall and, hence, it depends in a non-trivial way on the control parameters. In the form (4.15)–(4.16), our model appears, then, as an equivalent similarity theory for unstably stratified turbulent channel flows. Furthermore, the expression for the Obukhov length $L_{O}$ (4.17) in terms of the controlling parameters allows us to check the consistency of our data, in the logarithmic layer, with previous theoretical studies (Dubrulle et al. Reference Dubrulle, Laval, Sullivan and Werne2002b ) and, as a consequence, with experimental data (Businger et al. Reference Businger, Wyngaard, Izumi and Bradley1971). In figure 8 we compare the temperature fluctuations ${\it\Theta}_{\ast }\equiv \overline{(T-\overline{T})^{2}}^{1/2}/(Q_{w}/u_{{\it\tau}})$ from our simulations with the RDT prediction (Dubrulle et al. Reference Dubrulle, Laval, Sullivan and Werne2002b ):

The numerical data and the theoretical result (dashed line) appear to be in fair agreement, within error bars; in particular, the dependence on $(\mathit{Ra},\mathit{Re}_{{\it\tau}})$ as in (4.17) is fulfilled.